#############################################################
############### Resource allocation - Reproductive effort ###
#############################################################


setwd("//hest.nas.ethz.ch/green_users_all_homes/schoebch/Documents/ETHZ/Collaborations/Jianguo Chen/resource allocation crops")
#setwd("H:\\Work\\Work in 2019\\Resource allocation")

Sys.setlocale(locale = "English")

# libraries
.libPaths("C:/Users/schoebch/Documents/Rlibraries")

library(nlme)
library(lme4)
library(lmerTest)
library(effects)
library(lattice)
library(ggplot2)
library(sjPlot)
library(ggeffects)
library(plyr)
library(MuMIn)
library(see)
library(car)
library(Hmisc)
library(ggpubr)
library(rockchalk)
library(lavaan)
library(lavaanPlot)
library(dplyr)
library(pdiv)
library(tidyr)
library(performance)
library(forecast)

### read individual level data and define variables
Crop<-read.table("traitdata_CHSP2018_20200404.txt",header=T)
str(Crop)
names(Crop)
head(Crop)
Crop$RV <- Crop$yield/Crop$biomass
Crop$RDV <- Crop$yield/(Crop$biomass+Crop$yield)
Crop$SPNUMB <- as.factor(Crop$spnumber)
Crop$BED <- as.factor(Crop$bed)
Crop$bedId <- as.factor(paste(Crop$country,Crop$BED))
Crop$DIV <- as.factor(paste(Crop$diversity,Crop$SPNUMB))
Crop <- Crop %>% 
  rename(
    biomass_ind = biomass,
    seedmass_ind = yield
  )

Crop$totmass_ind <- Crop$seedmass_ind+Crop$biomass_ind

single <- subset(Crop, SPNUMB=="0")
means_single <- aggregate(cbind(totmass_ind, biomass_ind, seedmass_ind)~country+ecotype+fertilised+species, single,mean, na.action=na.omit)

Crop1 <- merge(Crop, means_single, by=c("country","ecotype","fertilised","species"), all.x=T, suffixes = c("",".single"))

Crop1$RII_tot <- (Crop1$totmass_ind-Crop1$totmass_ind.single)/(Crop1$totmass_ind+Crop1$totmass_ind.single)
Crop1$RII_bio <- (Crop1$biomass_ind-Crop1$biomass_ind.single)/(Crop1$biomass_ind+Crop1$biomass_ind.single)
Crop1$RII_seed <- (Crop1$seedmass_ind-Crop1$seedmass_ind.single)/(Crop1$seedmass_ind+Crop1$seedmass_ind.single)
Crop1$NIntC_tot <- 2*((Crop1$totmass_ind-Crop1$totmass_ind.single)/(Crop1$totmass_ind.single+Crop1$totmass_ind+(abs(Crop1$totmass_ind-Crop1$totmass_ind.single))))
Crop1$NIntC_bio <- 2*((Crop1$biomass_ind-Crop1$biomass_ind.single)/(Crop1$biomass_ind.single+Crop1$biomass_ind+(abs(Crop1$biomass_ind-Crop1$biomass_ind.single))))
Crop1$NIntC_seed <- 2*((Crop1$seedmass_ind-Crop1$seedmass_ind.single)/(Crop1$seedmass_ind.single+Crop1$seedmass_ind+(abs(Crop1$seedmass_ind-Crop1$seedmass_ind.single))))



#### subset of complete data (i.e. individuals with yield and biomass data)
datRDV <- subset(Crop1, RDV>0)
datRDV <- subset(datRDV, biomass_ind>0)

#### define further variables
datRDV$Home <- factor(as.numeric(datRDV$country=="Spain" & datRDV$ecotype =="ES" | datRDV$country=="Switzerland" & datRDV$ecotype =="CH"))
datRDV$factor <- as.factor(paste(datRDV$species,datRDV$country,datRDV$ecotype,datRDV$fertilised,datRDV$DIV))

datRDV$fertilised <- revalue(datRDV$fertilised, c("no"="without fertiliser", "yes"="with fertiliser"))

datRDV$targetFG <- vector(mode="character",length=nrow(datRDV))
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Lentil") {
    datRDV$targetFG[i]= "legume"
  }
}
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Lupinus") {
    datRDV$targetFG[i]= "legume"
  }
}
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Avena") {
    datRDV$targetFG[i]= "cereal"
  }
}
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Triticum") {
    datRDV$targetFG[i]= "cereal"
  }
}
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Camelina") {
    datRDV$targetFG[i]= "oldHerb"
  }
}
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Linum") {
    datRDV$targetFG[i]= "oldHerb"
  }
}
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Coriander") {
    datRDV$targetFG[i]= "newHerb"
  }
}
for(i in 1:nrow(datRDV)){
  if (datRDV$species[i] == "Quinoa") {
    datRDV$targetFG[i]= "newHerb"
  }
}

datRDV$targetFG <- as.factor(datRDV$targetFG)

datRDV$cereal <- vector(mode="character",length=nrow(datRDV))
for(i in 1:nrow(datRDV)){
  if (datRDV$targetFG[i] == "cereal") {
    datRDV$cereal[i]= "1"
  } else {
    datRDV$cereal[i]= "0"
  }
}

datRDV$oldHerb <- vector(mode="character",length=nrow(datRDV))
for(i in 1:nrow(datRDV)){
  if (datRDV$targetFG[i] == "oldHerb") {
    datRDV$oldHerb[i]= "1"
  } else {
    datRDV$oldHerb[i]= "0"
  }
}

datRDV$legume <- vector(mode="character",length=nrow(datRDV))
for(i in 1:nrow(datRDV)){
  if (datRDV$targetFG[i] == "legume") {
    datRDV$legume[i]= "1"
  } else {
    datRDV$legume[i]= "0"
  }
}

datRDV$newHerb <- vector(mode="character",length=nrow(datRDV))
for(i in 1:nrow(datRDV)){
  if (datRDV$targetFG[i] == "newHerb") {
    datRDV$newHerb[i]= "1"
  } else {
    datRDV$newHerb[i]= "0"
  }
}

datRDV$targetFG <- revalue(datRDV$targetFG, c("cereal"="Cereals", "legume"="Legumes","newHerb"="Superasterids","oldHerb"="Superrosids"))
datRDV$species <- revalue(datRDV$species, c("Avena"="oat", "Camelina"="camelina","Coriander"="coriander","Lentil"="lentil","Linum"="linseed", "Lupinus"="lupin","Quinoa"="quinoa","Triticum"="wheat"))
datRDV$SPNUMB <- revalue(datRDV$SPNUMB, c("0"="single plant", "1"="monoculture","2"="2-species mixture","4"="4-species mixture"))

comm <- subset(datRDV, !community=="single")
comm <- droplevels(comm)


######### read plot-level data
plotlevel <- read.table("plot_level_data1.txt",header=T)
str(plotlevel)
plotlevel$plotId <- as.factor(paste(plotlevel$Country,plotlevel$plot, sep="_"))
plotlevel$species <- revalue(plotlevel$species, c("Avena"="oat", "Camelina"="camelina","Coriander"="coriander","Lentil"="lentil","Linum"="linseed", "Lupinus"="lupin","Quinoa"="quinoa","Triticum"="wheat"))
plotlevel <- plotlevel %>% 
  rename(
    biomass_sp = biomass,
    seedmass_sp = seedmass
  )
plotlevel$totmass_sp <- plotlevel$biomass_sp+plotlevel$seedmass_sp


### figure: species yield per plot
plotlevel$multiplier <- vector(mode="numeric", length = nrow(plotlevel))

for(i in 1:nrow(plotlevel)){
  if (plotlevel$sp_numb[i] =="one") {plotlevel$multiplier[i] = 4}
  else if (plotlevel$sp_numb[i] =="two") {plotlevel$multiplier[i] = 8}
  else {plotlevel$multiplier[i] = 16}
}

plotlevel$biomass_sp_m <- plotlevel$biomass_sp*plotlevel$multiplier
plotlevel$seedmass_sp_m <- plotlevel$seedmass_sp*plotlevel$multiplier
plotlevel$totmass_sp_m <- plotlevel$totmass_sp*plotlevel$multiplier

sp.mean <- aggregate(cbind(biomass_sp_m,seedmass_sp_m,totmass_sp_m)~Country+ecotype+fertilizer+species+sp_numb, plotlevel, mean, na.action=na.omit)
sp.mean1 <- aggregate(cbind(biomass_sp_m,seedmass_sp_m,totmass_sp_m)~Country+species+sp_numb, plotlevel, mean, na.action=na.omit)
sp.sd1 <- aggregate(cbind(biomass_sp_m,seedmass_sp_m,totmass_sp_m)~Country+species+sp_numb, plotlevel, sd, na.action=na.omit)
sp.length1 <- aggregate(cbind(biomass_sp_m,seedmass_sp_m,totmass_sp_m)~Country+species+sp_numb, plotlevel, length, na.action=na.omit)
sp.mean1$biomass_sp_m_se <- sp.sd1$biomass_sp_m/sqrt(sp.length1$biomass_sp_m)
sp.mean1$seedmass_sp_m_se <- sp.sd1$seedmass_sp_m/sqrt(sp.length1$seedmass_sp_m)
sp.mean1$totmass_sp_m_se <- sp.sd1$totmass_sp_m/sqrt(sp.length1$totmass_sp_m)
sp.mean2 <- sp.mean1[!sp.mean1$sp_numb=="single",]
sp.mean2 <- droplevels(sp.mean2)
sp.mean2$sp_numb <- factor(sp.mean2$sp_numb,levels = c("one", "two", "four"))  
sp.mean2$SPNUMB <- revalue(sp.mean2$sp_numb, c("one"="monoculture","two"="2-species mixture","four"="4-species mixture"))


pd <- position_dodge(0.1) # move them .05 to the left and right
png('ploY_CoDiSp.png', width = 5, height = 3, units = 'in', res = 300)
ploY1 <- ggplot(sp.mean2, aes(x = SPNUMB, y = seedmass_sp_m, colour = species),expand.grid=T)+
  labs(x="",y=expression("Seed yield (g"~m^-2*")")) + 
  geom_errorbar(aes(ymin=seedmass_sp_m-seedmass_sp_m_se, ymax=seedmass_sp_m+seedmass_sp_m_se), width=.1, position = pd) +
  geom_point(position = pd)+
  scale_color_brewer(palette="Set1")+
  theme(legend.position = "right", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))+
  facet_wrap(~Country)
ploY1
dev.off()

png('ploB_CoDiSp.png', width = 5, height = 3, units = 'in', res = 300)
ploB1 <- ggplot(sp.mean2, aes(x = SPNUMB, y = biomass_sp_m, colour = species),expand.grid=T)+
  labs(x="",y=expression("Vegetative biomass (g"~m^-2*")")) + 
  geom_errorbar(aes(ymin=biomass_sp_m-biomass_sp_m_se, ymax=biomass_sp_m+biomass_sp_m_se), width=.1, position = pd) +
  geom_point(position = pd)+
  scale_color_brewer(palette="Set1")+
  theme(legend.position = "right", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))+
  facet_wrap(~Country)
ploB1
dev.off()

## combine plots
jpeg('ED_Fig2.jpg', width = 8, height = 4, units = 'in', res = 300)
ggarrange(ploY1+theme(plot.margin = margin(l=10)), ploB1+theme(plot.margin = margin(l=10)), 
          labels = c("a", "b"),
          hjust = 0,
          ncol = 2, nrow = 1,
          common.legend = T,
          widths = c(2,2)
)
dev.off()


### calculate plot-level productivity and yield
plotlevel2 <- aggregate(cbind(biomass_sp,seedmass_sp,totmass_sp)~plotId, plotlevel, sum, na.action=na.omit)
plotlevel2 <- aggregate(cbind(biomass_sp,seedmass_sp,totmass_sp)~plotId, plotlevel, sum, na.action=na.omit)


### add plot-level productivity and yield data to individual-level database
comm2 <- merge(comm,plotlevel2, by="plotId", all.x = T)
names(comm2)[64] <- "totmass_plot"
names(comm2)[63] <- "seedmass_plot"
names(comm2)[62] <- "biomass_plot"
comm2$indId <- as.factor(paste(comm2$country,comm2$label2,sep = "_"))
str(comm2)

#########################
#### plotlevel yield ####
#########################

plotyield <- aggregate(cbind(seedmass_plot,biomass_plot,totmass_plot)~country+ecotype+Home+fertilised+diversity+SPNUMB+spcomp+bedId+plotId, comm2, mean)
str(plotyield)
plotyield$RE <- plotyield$seedmass_plot/plotyield$totmass_plot

ploY <- lmer(sqrt(seedmass_plot)~
               (country+ecotype+Home)*
               fertilised*
               (diversity+SPNUMB)+
               (1|spcomp)+
               (1|bedId),
             data=plotyield)
plot(ploY)
hist(resid(ploY))
check_model(ploY)
anova(ploY, type=1)
summary(ploY)
r.squaredGLMM(ploY)

ploB <- lmer(sqrt(biomass_plot)~
               (country+ecotype+Home)*
               fertilised*
               (diversity+SPNUMB)+
               (1|spcomp)+
               (1|bedId),
             data=plotyield)
plot(ploB)
hist(resid(ploB))
check_model(ploB)
anova(ploB, type=1)
summary(ploB)
r.squaredGLMM(ploB)

# figure
# Country, DIVERSITY

jpeg("ploY_CoDi.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("ploY_CoDi.pdf", width = 4, height = 5)
png('ploY_CoDi.png', width = 2.5, height = 3, units = 'in', res = 300)
plot_model(ploY, type = "emm", terms = c("SPNUMB","country"), title="", axis.title = c("",expression("Seed yield"~plot^-1~"(g)")), expand.grid=T)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
dev.off()

ploYvals <- ggemmeans(ploY, terms=c("SPNUMB","country"))

pd <- position_dodge(0.1) # move them .05 to the left and right
png('ploY_CoDi1_m2.png', width = 2.5, height = 3, units = 'in', res = 300)
pdf('Fig.1.pdf', width = 3.5, height = 4)
ggplot(ploYvals, aes(x = x, y = predicted^2*4, colour = group),expand.grid=T)+
  labs(x="",y=expression("Seed yield (g"~m^-2*")")) + 
  geom_errorbar(aes(ymin=conf.low^2*4, ymax=conf.high^2*4), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
dev.off()

ploBvals <- ggemmeans(ploB, terms=c("SPNUMB","country"))

pd <- position_dodge(0.1) # move them .05 to the left and right
png('ploB_CoDi1_m2.png', width = 2.5, height = 3, units = 'in', res = 300)
jpeg('ED_Fig1.jpg', width = 2.5, height = 3, units = 'in', res = 300)
ggplot(ploBvals, aes(x = x, y = predicted^2*4, colour = group),expand.grid=T)+
  labs(x="",y=expression("Vegetative biomass (g"~m^-2*")")) + 
  geom_errorbar(aes(ymin=conf.low^2*4, ymax=conf.high^2*4), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
dev.off()

########## plotlevel RE
ploRE <- lmer((RE)~
               (country+ecotype+Home)*
               fertilised*
               (diversity+SPNUMB)+
               (1|spcomp)+
               (1|bedId),
             data=plotyield)
plot(ploRE)
hist(resid(ploRE))
check_model(ploRE)
anova(ploRE, type=1)
summary(ploRE)
r.squaredGLMM(ploRE)

png('ploRE_CoDi_HI.png', width = 2.5, height = 3, units = 'in', res = 300)
jpeg('ED_Fig7.jpg', width = 2.5, height = 3, units = 'in', res = 300)
plot_model(ploRE, type = "emm", terms = c("SPNUMB","country"), title="", axis.title = c("",expression("Harvest index"~plot^-1~"(g)")), expand.grid=T)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
dev.off()

ploREvals <- ggemmeans(ploRE, terms=c("SPNUMB","country"))



#####################################################
########### Plant interaction - partitioning ########
#####################################################
comm3 <- merge(comm2,plotlevel[,cbind("biomass_sp","seedmass_sp","totmass_sp","species","plotId")], by=c("species","plotId"))
str(comm3)

###################### SPECIES LEVEL
##### NIntC net
net <- aggregate(cbind(NIntC_seed,NIntC_bio,NIntC_tot, height, leafarea, SLA, LDMC, mass_per_seed, RDV)~plotId+species+country+ecotype+fertilised+community+diversity+SPNUMB+spcomp+BED+bedId+DIV+Home,comm3,mean, na.action=na.pass) # 1796 species per plot
str(net)
net2 <- merge(net,plotlevel[,cbind("biomass_sp","seedmass_sp","totmass_sp","species","plotId")], by=c("species","plotId"))

#### NIntC fac
pos_seed <- subset(net, NIntC_seed>0, select = c(plotId,species,NIntC_seed)) # 82 species per plot (4.6% of the species per plot)
pos_seed <- droplevels(pos_seed)
pos_bio <- subset(net, NIntC_bio>0, select = c(plotId,species,NIntC_bio)) # 28 species per plot (1.6% of the species per plot)
pos_bio <- droplevels(pos_bio)
pos_tot <- subset(net, NIntC_tot>0, select = c(plotId,species,NIntC_tot)) # 26 species per plot (1.4% of the species per plot)
pos_tot <- droplevels(pos_tot)

RII2 <- merge(net2,pos_seed, by=c("plotId","species"), all.x=T, suffixes = c("","_pos"))
RII2 <- merge(RII2,pos_bio, by=c("plotId","species"), all.x=T, suffixes = c("","_pos"))
RII2 <- merge(RII2,pos_tot, by=c("plotId","species"), all.x=T, suffixes = c("","_pos"))
str(RII2)

##### NIntC comp
neg_seed <- subset(net, NIntC_seed<0, select = c(plotId,species,NIntC_seed)) #1714 species per plot (95.4% of the species per plot)
neg_seed <- droplevels(neg_seed)
neg_bio <- subset(net, NIntC_bio<0, select = c(plotId,species,NIntC_bio)) #1768 species per plot (98.4% of the species per plot)
neg_bio <- droplevels(neg_bio)
neg_tot <- subset(net, NIntC_tot<0, select = c(plotId,species,NIntC_tot)) #1770 species per plot (98.6% of the species per plot)
neg_tot <- droplevels(neg_tot)

RII3 <- merge(RII2,neg_seed, by=c("plotId","species"), all.x=T, suffixes = c("","_neg"))
RII3 <- merge(RII3,neg_bio, by=c("plotId","species"), all.x=T, suffixes = c("","_neg"))
RII3 <- merge(RII3,neg_tot, by=c("plotId","species"), all.x=T, suffixes = c("","_neg"))
str(RII3)

# replace NAs in neg and pos with 0
RII3$NIntC_seed_neg[is.na(RII3$NIntC_seed_neg)] <- 0
RII3$NIntC_seed_pos[is.na(RII3$NIntC_seed_pos)] <- 0
RII3$NIntC_bio_neg[is.na(RII3$NIntC_bio_neg)] <- 0
RII3$NIntC_bio_pos[is.na(RII3$NIntC_bio_pos)] <- 0
RII3$NIntC_tot_neg[is.na(RII3$NIntC_tot_neg)] <- 0
RII3$NIntC_tot_pos[is.na(RII3$NIntC_tot_pos)] <- 0

#summary(RII3$NIntC_seed_neg+RII3$NIntC_seed_pos-RII3$NIntC_seed)
#summary(RII3$NIntC_bio_neg+RII3$NIntC_bio_pos-RII3$NIntC_bio)
#summary(RII3$NIntC_tot_neg+RII3$NIntC_tot_pos-RII3$NIntC_tot)

########### plots
bwplot(NIntC_seed~SPNUMB|species, data=RII3)
bwplot(NIntC_bio~SPNUMB|species, data=RII3)
bwplot(NIntC_tot~SPNUMB|species, data=RII3)


##################################################
######### mixtures only ##########################
# subset mixtures
mix <- subset(RII3, diversity=="mix") # 1671 species per plot
mix <- droplevels(mix)
mono <- subset(RII3, diversity=="mono") # 125 species per plot
mono <- droplevels(mono)

monoMean <- aggregate(cbind(NIntC_seed,NIntC_bio,NIntC_tot,seedmass_sp,biomass_sp,totmass_sp, RDV, height, leafarea,SLA,LDMC)~species+country+ecotype+fertilised, mono, mean, na.action=na.omit)
mix2 <- merge(mix, monoMean, by=c("species","country","ecotype","fertilised"), all.x = T, suffixes = c("","_mono"))
str(mix2)

mix2$mixNIntC_seed = mix2$NIntC_seed - mix2$NIntC_seed_mono
mix2$mixNIntC_bio = mix2$NIntC_bio - mix2$NIntC_bio_mono
mix2$mixNIntC_tot = mix2$NIntC_tot - mix2$NIntC_tot_mono

mix2$mixNIntC_seed_pos = mix2$mixNIntC_seed
mix2$mixNIntC_seed_pos[mix2$mixNIntC_seed_pos<0] <- 0
mix2$mixNIntC_seed_neg = mix2$mixNIntC_seed
mix2$mixNIntC_seed_neg[mix2$mixNIntC_seed_neg>0] <- 0

mix2$mixNIntC_bio_pos = mix2$mixNIntC_bio
mix2$mixNIntC_bio_pos[mix2$mixNIntC_bio_pos<0] <- 0
mix2$mixNIntC_bio_neg = mix2$mixNIntC_bio
mix2$mixNIntC_bio_neg[mix2$mixNIntC_bio_neg>0] <- 0

mix2$mixNIntC_tot_pos = mix2$mixNIntC_tot
mix2$mixNIntC_tot_pos[mix2$mixNIntC_tot_pos<0] <- 0
mix2$mixNIntC_tot_neg = mix2$mixNIntC_tot
mix2$mixNIntC_tot_neg[mix2$mixNIntC_tot_neg>0] <- 0

mix2$sampling_sign_seed = mix2$mixNIntC_seed
mix2$sampling_sign_seed[mix2$sampling_sign_seed < 0] <- -1
mix2$sampling_sign_seed[mix2$sampling_sign_seed > 0] <- 1
mix2$sampling_sign_seed <- mix2$sampling_sign_seed*mix2$seedmass_sp_mono

mix2$sampling_sign_bio = mix2$mixNIntC_bio
mix2$sampling_sign_bio[mix2$sampling_sign_bio < 0] <- -1
mix2$sampling_sign_bio[mix2$sampling_sign_bio > 0] <- 1
mix2$sampling_sign_bio <- mix2$sampling_sign_bio*mix2$biomass_sp_mono

mix2$sampling_sign_tot = mix2$mixNIntC_tot
mix2$sampling_sign_tot[mix2$sampling_sign_tot < 0] <- -1
mix2$sampling_sign_tot[mix2$sampling_sign_tot > 0] <- 1
mix2$sampling_sign_tot <- mix2$sampling_sign_tot*mix2$totmass_sp_mono

mix2$mixRDV = mix2$RDV-mix2$RDV_mono
mix2$mixleafarea = mix2$leafarea-mix2$leafarea_mono
mix2$mixSLA = mix2$SLA-mix2$SLA_mono
mix2$mixLDMC = mix2$LDMC-mix2$LDMC_mono
mix2$mixmass_per_seed = mix2$mass_per_seed-mix2$seedmass_sp_mono
mix2$mixheight = mix2$height-mix2$height_mono

###
modmixRDV <- lmer(mixRDV~species*
                    (mixNIntC_tot+
                    mixheight+mixleafarea+mixSLA+mixLDMC)+ #
                    (1|bedId)+
                    (1|plotId),
                  data=mix2)
plot(modmixRDV)
hist(resid(modmixRDV))
check_model(modmixRDV)
anova(modmixRDV)
summary(modmixRDV)
r.squaredGLMM(modmixRDV)

png('mixRDV_NIntCtot_HI.png', width = 3, height = 3, units = 'in', res = 300)
RDVplot <- plot_model(modmixRDV, type = "emm", terms = c("mixNIntC_tot","species"), title="", axis.title = c(expression(NIntC[mixture]~-~NIntC[monoculture]),expression("Harvest"~index[mixture]~"- Harvest"~index[monoculture])), expand.grid=T, colors = "Set1")+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
RDVplot
dev.off()

png('mixRDV_height_HI.png', width = 3, height = 3, units = 'in', res = 300)
RDVplot0 <- plot_model(modmixRDV, type = "emm", terms = c("mixheight","species"), title="", axis.title = c(expression(Height[mixture]~-~Height[monoculture]),expression("Harvest"~index[mixture]~"- Harvest"~index[monoculture])), expand.grid=T, colors = "Set1")+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
RDVplot0
dev.off()

png('mixRDV_LA_HI.png', width = 3, height = 3, units = 'in', res = 300)
RDVplot1a <- plot_model(modmixRDV, type = "emm", terms = c("mixleafarea","species"), title="", axis.title = c(expression(Leaf~area[mixture]~-~Leaf~area[monoculture]),expression("Harvest"~index[mixture]~"- Harvest"~index[monoculture])), expand.grid=T, colors = "Set1")+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
RDVplot1a
dev.off()

png('mixRDV_SLA_HI.png', width = 3, height = 3, units = 'in', res = 300)
RDVplot1 <- plot_model(modmixRDV, type = "emm", terms = c("mixSLA","species"), title="", axis.title = c(expression(SLA[mixture]~-~SLA[monoculture]),expression("Harvest"~index[mixture]~"- Harvest"~index[monoculture])), expand.grid=T, colors = "Set1")+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
RDVplot1
dev.off()

png('mixRDV_LDMC_HI.png', width = 3, height = 3, units = 'in', res = 300)
RDVplot2 <- plot_model(modmixRDV, type = "emm", terms = c("mixLDMC","species"), title="", axis.title = c(expression(LDMC[mixture]~-~LDMC[monoculture]),expression("Harvest"~index[mixture]~"- Harvest"~index[monoculture])), expand.grid=T, colors = "Set1")+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
RDVplot2
dev.off()

png('mixRDV_traits1_HI.png', width = 9, height = 9, units = 'in', res = 300)
pdf('Fig.4.pdf', width = 9, height = 9)
ggarrange(RDVplot+theme(plot.margin = margin(l=10)), RDVplot0+theme(axis.title.y = element_blank()),RDVplot1a+theme(axis.title.y = element_blank()),RDVplot1+theme(plot.margin = margin(l=10)),RDVplot2+theme(axis.title.y = element_blank()), 
          labels = c("a", "b", "c", "d", "e"),
          hjust = 0,
          ncol = 3, nrow = 2,
          common.legend = T)
dev.off()



##### plot level
# complementarity effect = mixNIntC = NIntCmix - NIntCmono
# sampling effect = NIntC_sampling = min(mixNIntC_pos,abs_mixNIntC_neg)
plot <- aggregate(cbind(mixNIntC_seed,NIntC_seed,NIntC_seed_pos,NIntC_seed_neg, NIntC_seed_mono,mixNIntC_bio,NIntC_bio,NIntC_bio_pos,NIntC_bio_neg, NIntC_bio_mono,mixNIntC_tot,NIntC_tot,NIntC_tot_pos,NIntC_tot_neg, NIntC_tot_mono)~plotId+country+ecotype+fertilised+SPNUMB+spcomp+BED+bedId+DIV+Home, mix2, mean)
str(plot) #637 plots
plot1 <- aggregate(cbind(mixNIntC_seed_pos,mixNIntC_seed_neg, sampling_sign_seed,mixNIntC_bio_pos,mixNIntC_bio_neg, sampling_sign_bio,mixNIntC_tot_pos,mixNIntC_tot_neg, sampling_sign_tot)~plotId+country+ecotype+fertilised+SPNUMB+spcomp+BED+bedId+DIV+Home, mix2, sum)
str(plot1)

plot2 <- merge(plot, plot1, by=c("plotId","country","ecotype","fertilised","SPNUMB","spcomp","BED","bedId","DIV","Home"), all.x=T)
plot2$abs_mixNIntC_seed_neg <- abs(plot2$mixNIntC_seed_neg)
plot2$abs_mixNIntC_bio_neg <- abs(plot2$mixNIntC_bio_neg)
plot2$abs_mixNIntC_tot_neg <- abs(plot2$mixNIntC_tot_neg)
str(plot2)
plot2$sampling_sign_seed[plot2$sampling_sign_seed > 0] <- 1
plot2$sampling_sign_seed[plot2$sampling_sign_seed < 0] <- -1
plot2$sampling_sign_bio[plot2$sampling_sign_bio > 0] <- 1
plot2$sampling_sign_bio[plot2$sampling_sign_bio < 0] <- -1
plot2$sampling_sign_tot[plot2$sampling_sign_tot > 0] <- 1
plot2$sampling_sign_tot[plot2$sampling_sign_tot < 0] <- -1

plot2$NIntC_seed_sampling <- apply(cbind(plot2$mixNIntC_seed_pos, plot2$abs_mixNIntC_seed_neg), 1, min)*plot2$sampling_sign_seed
plot2$NIntC_bio_sampling <- apply(cbind(plot2$mixNIntC_bio_pos, plot2$abs_mixNIntC_bio_neg), 1, min)*plot2$sampling_sign_bio
plot2$NIntC_tot_sampling <- apply(cbind(plot2$mixNIntC_tot_pos, plot2$abs_mixNIntC_tot_neg), 1, min)*plot2$sampling_sign_tot
str(plot2)

plot(plot2$NIntC_bio_sampling~plot2$NIntC_seed_sampling)
plot(plot2$NIntC_tot_sampling~plot2$NIntC_seed_sampling)
plot(plot2$NIntC_bio_sampling~plot2$NIntC_tot_sampling)

plot(plot2$mixNIntC_bio~plot2$mixNIntC_seed)
plot(plot2$mixNIntC_tot~plot2$mixNIntC_seed)
plot(plot2$mixNIntC_bio~plot2$mixNIntC_tot)

###
bwplot(NIntC_seed_sampling~SPNUMB|country, data=plot2)
bwplot(mixNIntC_seed~SPNUMB|country, data=plot2)

#########################
### biodiversity effects
#########################
####*** diversity effects manually ***#### CHRISTIAN
mix2$spnumb <- vector("numeric", length = nrow(mix2))
mix2$spnumb <- ifelse(mix2$SPNUMB=="2-species mixture",2,4)

# species-level
mix2$RYoi_seed <- mix2$seedmass_sp/mix2$seedmass_sp_mono
mix2$RYei_seed <- 1/mix2$spnumb
mix2$Yei_seed <- mix2$RYei_seed*mix2$seedmass_sp_mono
mix2$deltaRYi_seed <- mix2$RYoi_seed-mix2$RYei_seed

mix2$RYoi_bio <- mix2$biomass_sp/mix2$biomass_sp_mono
mix2$RYei_bio <- 1/mix2$spnumb
mix2$Yei_bio <- mix2$RYei_bio*mix2$biomass_sp_mono
mix2$deltaRYi_bio <- mix2$RYoi_bio-mix2$RYei_bio

mix2$RYoi_tot <- mix2$totmass_sp/mix2$totmass_sp_mono
mix2$RYei_tot <- 1/mix2$spnumb
mix2$Yei_tot <- mix2$RYei_tot*mix2$totmass_sp_mono
mix2$deltaRYi_tot <- mix2$RYoi_tot-mix2$RYei_tot

# sampling effect
SE <- as.data.frame(as.factor(levels(mix2$plotId)))
SE$SE_seed <- vector(mode="numeric", length=nlevels(mix2$plotId))
SE$SE_bio <- vector(mode="numeric", length=nlevels(mix2$plotId))
SE$SE_tot <- vector(mode="numeric", length=nlevels(mix2$plotId))
names(SE)[1] <- "plotId"
plotID <- as.factor(levels(mix2$plotId))

for (i in 1:nlevels(mix2$plotId)){
  dat <- subset(mix2, plotId==plotID[i])
  SE$SE_seed[i] <- mean(dat$spnumb)*(cov(dat$deltaRYi_seed,dat$seedmass_sp_mono))*((mean(dat$spnumb)-1)/mean(dat$spnumb))
  SE$SE_bio[i] <- mean(dat$spnumb)*(cov(dat$deltaRYi_bio,dat$biomass_sp_mono))*((mean(dat$spnumb)-1)/mean(dat$spnumb))
  SE$SE_tot[i] <- mean(dat$spnumb)*(cov(dat$deltaRYi_tot,dat$totmass_sp_mono))*((mean(dat$spnumb)-1)/mean(dat$spnumb))
}

# plot-level
plot10 <- aggregate(cbind(RYoi_seed, Yei_seed, seedmass_sp,RYoi_bio, Yei_bio, biomass_sp,RYoi_tot, Yei_tot, totmass_sp)~plotId+country+ecotype+fertilised+diversity+SPNUMB+spnumb+spcomp+BED+bedId+DIV+Home, mix2, sum)
names(plot10)[13] <- "RYo_seed"
names(plot10)[14] <- "Ye_seed"
names(plot10)[15] <- "Yo_seed"
names(plot10)[16] <- "RYo_bio"
names(plot10)[17] <- "Ye_bio"
names(plot10)[18] <- "Yo_bio"
names(plot10)[19] <- "RYo_tot"
names(plot10)[20] <- "Ye_tot"
names(plot10)[21] <- "Yo_tot"
str(plot10)

plot11 <- aggregate(cbind(deltaRYi_seed, seedmass_sp_mono,deltaRYi_bio, biomass_sp_mono,deltaRYi_tot, totmass_sp_mono)~plotId+country+ecotype+fertilised+diversity+SPNUMB+spnumb+spcomp+BED+bedId+DIV+Home, mix2, mean) 
names(plot11)[13] <- "meandeltaRY_seed"
names(plot11)[14] <- "meanM_seed"
names(plot11)[15] <- "meandeltaRY_bio"
names(plot11)[16] <- "meanM_bio"
names(plot11)[17] <- "meandeltaRY_tot"
names(plot11)[18] <- "meanM_tot"
str(plot11)

plot10$deltaY_seed <- plot10$Yo_seed-plot10$Ye_seed # net effect
plot10$deltaY_bio <- plot10$Yo_bio-plot10$Ye_bio # net effect
plot10$deltaY_tot <- plot10$Yo_tot-plot10$Ye_tot # net effect

plot12 <- merge(plot10, plot11, by=c("plotId","country","ecotype","fertilised","diversity","SPNUMB","spnumb","spcomp","BED","bedId","DIV","Home"))
plot13 <- merge(plot12,SE, by=("plotId"))

# Relative Yield Total (RYT)
plot13$RYT_seed <- 1+(plot13$spnumb*plot13$meandeltaRY_seed)
plot13$RYT_bio <- 1+(plot13$spnumb*plot13$meandeltaRY_bio)
plot13$RYT_tot <- 1+(plot13$spnumb*plot13$meandeltaRY_tot)

# complementarity effect
plot13$CE_seed <- plot13$spnumb*plot13$meandeltaRY_seed*plot13$meanM_seed
plot13$CE_bio <- plot13$spnumb*plot13$meandeltaRY_bio*plot13$meanM_bio
plot13$CE_tot <- plot13$spnumb*plot13$meandeltaRY_tot*plot13$meanM_tot

plot(plot13$SE_bio~plot13$SE_seed)
plot(plot13$SE_bio~plot13$SE_tot)
plot(plot13$SE_seed~plot13$SE_tot)

plot(plot13$CE_bio~plot13$CE_seed)
plot(plot13$CE_bio~plot13$CE_tot)
plot(plot13$CE_seed~plot13$CE_tot)

plot20 <- merge(plot13, plot2, by=c("plotId","country","ecotype","fertilised","SPNUMB","spcomp","BED","bedId","DIV","Home"))
str(plot20)

plot(RYT_bio~RYT_seed, data=plot20)
abline(a=0, b=1)


plot(CE_seed~mixNIntC_seed, data=plot20)
plot(SE_seed~NIntC_seed_sampling, data=plot20)
cor.test(plot20$SE_seed,plot20$NIntC_seed_sampling, na.rm=T)
cor.test(plot20$CE_seed,plot20$mixNIntC_seed, na.rm=T)

plot(CE_bio~mixNIntC_bio, data=plot20)
plot(SE_bio~NIntC_bio_sampling, data=plot20)
cor.test(plot20$SE_bio,plot20$NIntC_bio_sampling, na.rm=T)
cor.test(plot20$CE_bio,plot20$mixNIntC_bio, na.rm=T)

plot(CE_tot~mixNIntC_tot, data=plot20)
plot(SE_tot~NIntC_tot_sampling, data=plot20)
cor.test(plot20$SE_tot,plot20$NIntC_tot_sampling, na.rm=T)
cor.test(plot20$CE_tot,plot20$mixNIntC_tot, na.rm=T)

#### wide to long
plot20long <- rbind(plot20[,1:12],plot20[,1:12],plot20[,1:12])
plot20long$level <- factor(c(rep("vegetative biomass",nrow(plot20)),rep("seed yield",nrow(plot20)),rep("totmass",nrow(plot20))))
plot20long$Yo_tot <- c(plot20$Yo_tot,plot20$Yo_tot,plot20$Yo_tot)
plot20long$RYo <- c(plot20$RYo_bio,plot20$RYo_seed,plot20$RYo_tot)
plot20long$Yo <- c(plot20$Yo_bio,plot20$Yo_seed,plot20$Yo_tot)
plot20long$Ye <- c(plot20$Ye_bio,plot20$Ye_seed,plot20$Ye_tot)
plot20long$deltaY <- c(plot20$deltaY_bio,plot20$deltaY_seed,plot20$deltaY_tot)
plot20long$SE <- c(plot20$SE_bio,plot20$SE_seed,plot20$SE_tot)
plot20long$CE <- c(plot20$CE_bio,plot20$CE_seed,plot20$CE_tot)
plot20long$RYT <- c(plot20$RYT_bio,plot20$RYT_seed,plot20$RYT_tot)
plot20long$mixNIntC <- c(plot20$mixNIntC_bio,plot20$mixNIntC_seed,plot20$mixNIntC_tot)
plot20long$mixNIntC_pos <- c(plot20$mixNIntC_bio_pos,plot20$mixNIntC_seed_pos,plot20$mixNIntC_tot_pos)
plot20long$mixNIntC_neg <- c(plot20$mixNIntC_bio_neg,plot20$mixNIntC_seed_neg,plot20$mixNIntC_tot_neg)
plot20long$NIntC <- c(plot20$NIntC_bio,plot20$NIntC_seed,plot20$NIntC_tot)
plot20long$NIntC_pos <- c(plot20$NIntC_bio_pos,plot20$NIntC_seed_pos,plot20$NIntC_tot_pos)
plot20long$NIntC_neg <- c(plot20$NIntC_bio_neg,plot20$NIntC_seed_neg,plot20$NIntC_tot_neg)
plot20long$NIntC_sampling <- c(plot20$NIntC_bio_sampling,plot20$NIntC_seed_sampling,plot20$NIntC_tot_sampling)
str(plot20long)

plot21long <- subset(plot20long, !level=="totmass")
plot21long <- droplevels(plot21long)
plot21long$sign_deltaY <- plot21long$deltaY
plot21long$sign_deltaY[plot21long$sign_deltaY<0] <- -1
plot21long$sign_deltaY[plot21long$sign_deltaY>0] <- 1
plot21long$sqrt_deltaY <- sqrt(abs(plot21long$deltaY))*plot21long$sign_deltaY
plot21long$sign_CE <- plot21long$CE
plot21long$sign_CE[plot21long$sign_CE<0] <- -1
plot21long$sign_CE[plot21long$sign_CE>0] <- 1
plot21long$sqrt_CE <- sqrt(abs(plot21long$CE))*plot21long$sign_CE
plot21long$sign_SE <- plot21long$SE
plot21long$sign_SE[plot21long$sign_SE<0] <- -1
plot21long$sign_SE[plot21long$sign_SE>0] <- 1
plot21long$sqrt_SE <- sqrt(abs(plot21long$SE))*plot21long$sign_SE


plot(plot21long$mixNIntC~plot21long$NIntC_pos)
plot(plot21long$mixNIntC~plot21long$NIntC_neg)


## stats
# net
modNet <- lmer(sqrt_deltaY~level*
                ((country+ecotype+Home)*
                   fertilised*
                   (SPNUMB))+
                #(1|spcomp)+
                (1|bedId)+
                (1|plotId),
              data=plot21long)
plot(modNet)
hist(resid(modNet))
anova(modNet, type=1)
summary(modNet)
r.squaredGLMM(modNet)

jpeg("net_CoDi.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("net_CoDi.pdf", width = 4, height = 5)
png('net_CoDi.png', width = 3, height = 3, units = 'in', res = 300)
netplot <- plot_model(modNet, type = "emm", terms = c("SPNUMB","level","country"), title="", axis.title = c("",expression("Net effect"~plot^-1~"(g)")), expand.grid=T)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
netplot
dev.off()

netplotvals <- ggemmeans(modNet, c("SPNUMB","level","country")) 

png('net_CoDi1_m2.png', width = 3, height = 3, units = 'in', res = 300)
netplot1 <- ggplot(netplotvals, aes(x = x, y = predicted^2*4, colour = group),expand.grid=T)+
  labs(x="",y=expression("Net effect (g"~m^-2*")")) + 
  geom_errorbar(aes(ymin=conf.low^2*4, ymax=conf.high^2*4), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))+
  facet_wrap(~facet)
netplot1
dev.off()

ggemmeans(modNet, c("level","country"))
ggemmeans(modNet, c("level","SPNUMB"))
ggemmeans(modNet, c("level","SPNUMB","country")) 

ggeffect(modNet, c("level","SPNUMB","country"), transform="exp") #


# complementarity
modCE <- lmer(sqrt_CE~level*
                ((country+ecotype+Home)*
                   fertilised*
                   (SPNUMB))+
                #(1|spcomp)+
                (1|bedId)+
                (1|plotId),
              data=plot21long)
plot(modCE)
hist(resid(modCE))
anova(modCE, type=1)
summary(modCE)
r.squaredGLMM(modCE)

jpeg("CE_CoDi.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("CE_CoDi.pdf", width = 4, height = 5)
png('CE_CoDi.png', width = 3, height = 3, units = 'in', res = 300)
CEplot <- plot_model(modCE, type = "emm", terms = c("SPNUMB","level","country"), title="", axis.title = c("",expression("Complementarity effect"~plot^-1~"(g)")), expand.grid=T)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
CEplot
dev.off()

CEplotvals <- ggemmeans(modCE, c("SPNUMB","level","country"))

png('CE_CoDi1_m2.png', width = 3, height = 3, units = 'in', res = 300)
CEplot1 <- ggplot(CEplotvals, aes(x = x, y = predicted^2*4, colour = group),expand.grid=T)+
  labs(x="",y=expression("Complementarity effect (g"~m^-2*")")) + 
  geom_errorbar(aes(ymin=conf.low^2*4, ymax=conf.high^2*4), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))+
  facet_wrap(~facet)
CEplot1
dev.off()

ggemmeans(modCE, c("level","country"))
ggemmeans(modCE, c("SPNUMB","country"))

# sampling
modSE <- lmer(sqrt_SE~level*
                ((country+ecotype+Home)*
                   fertilised*
                   (SPNUMB))+
                #(1|spcomp)+
                (1|bedId)+
                (1|plotId),
              data=plot21long)
plot(modSE)
hist(resid(modSE))
check_model(modSE)
r2(modSE)
anova(modSE, type=1)
summary(modSE)
r.squaredGLMM(modSE)

jpeg("SE_CoDi.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("SE_CoDi.pdf", width = 4, height = 5)
png('SE_CoDi.png', width = 3, height = 3, units = 'in', res = 300)
SEplot <- plot_model(modSE, type = "emm", terms = c("SPNUMB","level","country"), title="", axis.title = c("",expression("Sampling effect"~plot^-1~"(g)")), expand.grid=T)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
SEplot
dev.off

SEplotvals <- ggemmeans(modSE, c("SPNUMB","level","country"))

png('SE_CoDi1_m2.png', width = 3, height = 3, units = 'in', res = 300)
SEplot1 <- ggplot(SEplotvals, aes(x = x, y = predicted^2*4, colour = group),expand.grid=T)+
  labs(x="",y=expression("Sampling effect (g"~m^-2*")")) + 
  geom_errorbar(aes(ymin=conf.low^2*4, ymax=conf.high^2*4), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))+
  facet_wrap(~facet)
SEplot1
dev.off()

ggemmeans(modSE, c("level","country"))
ggemmeans(modSE, c("level","SPNUMB","country"))
ggemmeans(modSE, c("SPNUMB","country"))
ggemmeans(modSE, c("SPNUMB"))
ggemmeans(modSE, c("SPNUMB","level"))


# RYT
modRYT <- lmer(RYT~level*
                ((country+ecotype+Home)*
                   fertilised*
                   (SPNUMB))+
                #(1|spcomp)+
                (1|bedId)+
                (1|plotId),
              data=plot21long)
plot(modRYT)
hist(resid(modRYT))
anova(modRYT, type=1)
summary(modRYT)
r.squaredGLMM(modRYT)

jpeg("RYT_CoDi.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("RYT_CoDi.pdf", width = 4, height = 5)
png('RYT_CoDi.png', width = 3, height = 3, units = 'in', res = 300)
RYTplot <- plot_model(modRYT, type = "emm", terms = c("SPNUMB","level","country"), title="", axis.title = c("",expression("Relative Yield Total"~plot^-1)), expand.grid=T)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
RYTplot
dev.off()

RYTplotvals <- ggemmeans(modRYT, c("SPNUMB","level","country"))

png('RYT_CoDi1_new.png', width = 3, height = 3, units = 'in', res = 300)
RYTplot1 <- ggplot(RYTplotvals, aes(x = x, y = predicted, colour = group),expand.grid=T)+
  labs(x="",y=expression("Relative Yield Total")) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "top", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))+
  facet_wrap(~facet)
RYTplot1
dev.off()

ggemmeans(modCE, c("level","country"))
ggemmeans(modCE, c("SPNUMB","country"))


# combine plots: net & RYT
png('netRYT_CoDi_m2.png', width = 3, height = 6, units = 'in', res = 300)
pdf('Fig.2.pdf', width = 3.5, height = 7)
ggarrange(netplot1+theme(axis.text.x = element_blank(),plot.margin = margin(l=10)), RYTplot1+theme(plot.margin = margin(l=10)), 
          labels = c("a", "b"),
          hjust = 0,
          ncol = 1, nrow = 2,
          common.legend = T,
          heights = c(2,2.9)
          )
dev.off()

# combine plots: CE & SE
png('CESE_CoDi_m2.png', width = 3, height = 6, units = 'in', res = 300)
jpeg('ED_Fig.3.jpg', width = 3, height = 6, units = 'in', res = 300)
ggarrange(CEplot1+theme(axis.text.x = element_blank(),plot.margin = margin(l=10)), SEplot1+theme(plot.margin = margin(l=10)), 
          labels = c("a", "b"),
          hjust = 0,
          ncol = 1, nrow = 2,
          common.legend = T,
          heights = c(2,2.9)
)
dev.off()

# combine plots: net, CE & SE
png('addpart_CoDi1.png', width = 3, height = 9, units = 'in', res = 300)
ggarrange(netplot1+theme(axis.text.x = element_blank(),plot.margin = margin(l=10)), CEplot1+theme(axis.text.x = element_blank(),plot.margin = margin(l=10)), SEplot1+theme(plot.margin = margin(l=10)), 
          labels = c("a", "b", "c"),
          hjust = 0,
          ncol = 1, nrow = 3,
          common.legend = T,
          heights = c(2,2,2.9)
)
dev.off()

#########################################################
##################################### RE (RDV)
########################################################

#### full model
modK00 <- lmer(sqrt(RDV)~(species)+
                (country+ecotype+Home)+
                fertilised+
                (diversity+SPNUMB)+
                (species):(country+ecotype+Home)+
                (species):fertilised+
                (species):(diversity+SPNUMB)+
                (country+ecotype+Home):fertilised+
                (country+ecotype+Home):(diversity+SPNUMB)+
                fertilised:(diversity+SPNUMB)+
                (species):(country+ecotype+Home):fertilised+
                (species):(country+ecotype+Home):(diversity+SPNUMB)+
                (country+ecotype+Home):fertilised:(diversity+SPNUMB)+
                (species):(country+ecotype+Home):fertilised:(diversity+SPNUMB)+
                (1|bedId)+(1|plotId)+
                (1|spcomp)
              , comm2)
summary(modK00)
anova(modK00, type=1)
r.squaredGLMM(modK00)
plot(modK00)
hist(resid(modK00))



########## for figures: get rid of NA estimates, i.e. nested factors: diversity and SPNUMB
modL1a <- lm(sqrt(RDV)~(species)+
             (country+ecotype+Home)+
             fertilised+
             (SPNUMB)+
             (species):(country+ecotype+Home)+
             (species):fertilised+
             (species):(SPNUMB)+
             (country+ecotype+Home):fertilised+
             (country+ecotype+Home):(SPNUMB)+
             fertilised:(SPNUMB)+
             (species):(country+ecotype+Home):fertilised+
             (species):(country+ecotype+Home):(SPNUMB)+
             (country+ecotype+Home):fertilised:(SPNUMB)+
             (species):(country+ecotype+Home):fertilised:(SPNUMB)
           , 
           data=comm2, na.action=na.omit
)
anova(modL1a)
summary(modL1a)
coef(modL1a)


############# figures
# Country, DIVERSITY
jpeg("RE_CoDi.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("RE_CoDi.pdf", width = 4, height = 5)
png('RE_CoDi_HI.png', width = 3, height = 5, units = 'in', res = 300)
plot1 <- plot_model(modL1a, type = "emm", terms = c("SPNUMB","country"), title="", axis.title = c("","Harvest index"), expand.grid=T)+
  theme(legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1),
        legend.position="top")
plot1
dev.off()


# Country, DIVERSITY, species
bwplot(pred~community|targetFG:country, data=datRDV)

jpeg("RE_CoDiS_comm.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("RE_CoDiS_comm.pdf", width = 4, height = 5)
png('RE_CoDiS_comm_HI.png', width = 5, height = 5, units = 'in', res = 300)
plot2 <- plot_model(modL1a, type = "emm", terms = c("SPNUMB","species","country"), title="", axis.title = c("","Harvest index"), expand.grid=T, colors = "Set1")+
  geom_line()+
  theme(legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
plot2
dev.off()

### combine country x diversity and country x diversity x species
pdf('Fig.3.pdf', width = 7, height = 4)
ggarrange(plot1+theme(plot.margin = margin(l=10)), plot2+theme(plot.margin = margin(l=10)), 
          labels = c("a", "b"),
          hjust = 0,
          ncol = 2, nrow = 1,
          common.legend = F,
          widths = c(2,3.5)
)
dev.off()


# Country, DIVERSITY, fertiliser, ecotype
bwplot(pred~community|targetFG:country, data=datRDV)

jpeg("RE_CoDiFEco_comm.jpeg", width = 500, height = 500, units = "px", pointsize = 12, quality = 300)
pdf("RE_CoDiFEco_comm.pdf", width = 4, height = 5)
png('RE_CoDiFEco_comm_HI.png', width = 5, height = 5, units = 'in', res = 300)
plot_model(modL1a, type = "emm", terms = c("SPNUMB","country","Home","fertilised"), title="", axis.title = c("","Harvest index"), expand.grid=T, colors = "Set1")+
  geom_line()+
  theme(legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust=1))
dev.off()


ggemmeans(modL1a, terms = c("Home","species"))
ggemmeans(modL1a, terms = c("SPNUMB","country"))
ggemmeans(modL1a, terms = c("SPNUMB","country","species"))
ggemmeans(modL1a, terms = c("SPNUMB","Home","fertilised"))

REplotvals1 <- ggemmeans(modL1a, terms = c("Home","SPNUMB","country","fertilised"))
REplotvals2 <- REplotvals1[!REplotvals1$group=="monoculture",]
REplotvals2$x <- revalue(REplotvals2$x, c("0"="Away", "1"="Home"))

jpeg('ED_Fig4.jpg', width = 5, height = 4, units = 'in', res = 300)
REplot1 <- ggplot(REplotvals2, aes(x = x, y = predicted, colour = group),expand.grid=T)+
  labs(x="",y=expression("Harvest index")) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "top", legend.title = element_blank())+
  facet_grid(facet~panel)
REplot1
dev.off()

REplotvals5 <- ggemmeans(modL1a, terms = c("Home","species"))
REplotvals5$x <- revalue(REplotvals5$x, c("0"="Away", "1"="Home"))

jpeg('ED_Fig6.jpg', width = 5, height = 3, units = 'in', res = 300)
REplot2 <- ggplot(REplotvals5, aes(x = x, y = predicted, colour = group),expand.grid=T)+
  scale_color_brewer(palette="Set1")+
    labs(x="",y=expression("Harvest index")) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.1, position = pd) +
  geom_point(position = pd)+
  theme(legend.position = "right", legend.title = element_blank())
REplot2
dev.off()


#RE for each species composition
#Avena
avena <- subset(comm3, species=="oat")
avena <- droplevels(avena)
head(avena)

plotav<-ggplot(avena,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                       width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("")+coord_flip()+theme(legend.position = "none")+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.2,0.5), breaks=seq(0.2,0.5,0.1))

plotav

#Triticum
triticum <- subset(comm3, species=="wheat")
triticum <- droplevels(triticum)

plottr<-ggplot(triticum,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                          width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("")+coord_flip()+theme(legend.position = "none")+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.2,0.4), breaks=seq(0.2,0.4,0.1))
plottr


#Camelina
camelina <- subset(comm3,species=="camelina")
camelina <- droplevels(camelina)

plotca<-ggplot(camelina,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                          width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("")+coord_flip()+theme(legend.position = "none")+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.2,0.5), breaks=seq(0.2,0.5,0.1))

plotca


#Coriander
coriander <- subset(comm3, species=="coriander")
coriander <- droplevels(coriander)

plotco<-ggplot(coriander,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                            width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("Harvest index")+coord_flip()+theme(legend.position = "none")+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.3,0.5), breaks=seq(0.3,0.5,0.1))
plotco

#Lentil
lentil <- subset(comm3, species=="lentil")
lentil <- droplevels(lentil)

plotle<-ggplot(lentil,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                        width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("")+coord_flip()+theme(legend.position = "none")+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.3,0.5), breaks=seq(0.3,0.5,0.1))

plotle

#Lupin
lupin <- subset(comm3, species=="lupin")
lupin <- droplevels(lupin)

plotlu<-ggplot(lupin,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                        width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("")+coord_flip()+theme(legend.position = "none")+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.3,0.5), breaks=seq(0.3,0.5,0.1))

plotlu

#Linum
linum <- subset(comm3, species=="linseed")
linum <- droplevels(linum)

plotli<-ggplot(linum,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                        width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("Harvest index")+coord_flip()+theme(legend.position = "none")+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.2,0.4), breaks=seq(0.2,0.4,0.1))

plotli

#Quinoa
quinoa <- subset(comm3, species=="quinoa")
quinoa <- droplevels(quinoa)

plotqu<-ggplot(quinoa,aes(x=spcomp,y=RDV))+stat_summary(geom = "errorbar", 
                                                        width = 0.1,fun.data = "mean_se", fun.args = list(mult = 1),position = position_dodge(0), aes(color = country)) + stat_summary(geom = "line", fun = mean, aes(color = country)) + 
  stat_summary(geom = "point", size = 3, fun = mean, aes(color = country))+xlab("Species composition")+ylab("Harvest index")+coord_flip()+
  theme(panel.background=element_rect(fill='transparent',color='black'))+
  theme(axis.text=element_text(size=8,color='black'))+theme(axis.title=element_text(size=8,face="bold"))+
  scale_y_continuous(limits=c(0.2,0.5), breaks=seq(0.2,0.5,0.1))
plotqu

### combine plots
jpeg('ED_Fig5.jpg', width = 10, height = 7, units = 'in', res = 300)
ggarrange(plotav+theme(plot.margin = margin(l=11)), plottr+theme(plot.margin = margin(l=11)), 
          plotle+theme(plot.margin = margin(l=11)), plotlu+theme(plot.margin = margin(l=11)),
          plotca+theme(plot.margin = margin(l=11)), plotli+theme(plot.margin = margin(l=11)),
          plotco+theme(plot.margin = margin(l=11)), plotqu+theme(plot.margin = margin(l=11)),
          labels = c("oat", "wheat","lentil","lupin","camelina","linseed","coriander","quinoa"),
          font.label = list(size =7),
          label.x = 0.25,
        #  hjust=-2,
          align = "hv",
          ncol = 3, nrow = 3,
          common.legend = T
          )
dev.off()



#####################################################################
########## plant development ########################################
#####################################################################
pwheat<-read.table("wheat.txt",header=T)
poat<-read.table("oat.txt",header=T)
plentil<-read.table("lentil.txt",header=T)
plupin<-read.table("lupin.txt",header=T)
pquinoa<-read.table("quinoa.txt",header=T)
pcoriander<-read.table("coriander.txt",header=T)
plinseed<-read.table("linseed.txt",header=T)
pcamelina<-read.table("camelina.txt",header=T)

# wheat
pwheat2 <- merge(plotlevel,pwheat, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
pwheat2$flower <- as.Date(pwheat2$flower, "%d.%m.%Y" )
pwheat2$startTillering <- as.Date(pwheat2$startTillering, "%d.%m.%Y" )
pwheat2$milk <- as.Date(pwheat2$milk, "%d.%m.%Y" )

str(pwheat2)
summary(pwheat2)
pwheat2$sp_numb <- factor(pwheat2$sp_numb , levels=c("single", "one", "two", "four"))
pwheat2$SP_NUMB <- revalue(pwheat2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

till <- bwplot(SP_NUMB~startTillering|Country, data=pwheat2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Wheat")
flow <- bwplot(SP_NUMB~flower|Country, data=pwheat2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Wheat")
frui <- bwplot(SP_NUMB~milk|Country, data=pwheat2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Wheat")
png('Phen_wheat.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_wheat.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow, Tillering=till, layout =c(2,3))
dev.off()

# oat
poat2 <- merge(plotlevel,poat, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
summary(poat2)

poat2$inflo <- as.Date(poat2$startInflorescence, "%d.%m.%Y" )
poat2$startTillering <- as.Date(poat2$startTillering, "%d.%m.%Y" )
poat2$milk <- as.Date(poat2$milk, "%d.%m.%Y" )

str(poat2)

poat2$sp_numb <- factor(poat2$sp_numb , levels=c("single", "one", "two", "four"))
poat2$SP_NUMB <- revalue(poat2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

till <- bwplot(SP_NUMB~startTillering|Country, data=poat2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Oat")
flow <- bwplot(SP_NUMB~inflo|Country, data=poat2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Oat")
frui <- bwplot(SP_NUMB~milk|Country, data=poat2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Oat")
png('Phen_oat.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_oat.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow, Tillering=till, layout =c(2,3))
dev.off()

# lupin
plupin2 <- merge(plotlevel,plupin, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
summary(plupin2)

plupin2$stemElongation <- as.Date(plupin2$stemElongation, "%d.%m.%Y" )
plupin2$flowering <- as.Date(plupin2$flowering, "%d.%m.%Y" )
plupin2$podGreen <- as.Date(plupin2$podGreen, "%d.%m.%Y" )

str(plupin2)

plupin2$sp_numb <- factor(plupin2$sp_numb , levels=c("single", "one", "two", "four"))
plupin2$SP_NUMB <- revalue(plupin2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

elon <- bwplot(SP_NUMB~stemElongation|Country, data=plupin2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Lupin")
flow <- bwplot(SP_NUMB~flowering|Country, data=plupin2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Lupin")
frui <- bwplot(SP_NUMB~podGreen|Country, data=plupin2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Lupin")
png('Phen_lupin.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_lupin.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow, Elongation=elon, layout =c(2,3))
dev.off()

# lentil  
plentil2 <- merge(plotlevel,plentil, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
summary(plentil2)

plentil2$firstLeaf <- as.Date(plentil2$firstLeaf, "%d.%m.%Y" )
plentil2$flowering <- as.Date(plentil2$flowering, "%d.%m.%Y" )
plentil2$podsFull <- as.Date(plentil2$podsFull, "%d.%m.%Y" )

str(plentil2)

plentil2$sp_numb <- factor(plentil2$sp_numb , levels=c("single", "one", "two", "four"))
plentil2$SP_NUMB <- revalue(plentil2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

leaf <- bwplot(SP_NUMB~firstLeaf|Country, data=plentil2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Lentil")
flow <- bwplot(SP_NUMB~flowering|Country, data=plentil2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Lentil")
frui <- bwplot(SP_NUMB~podsFull|Country, data=plentil2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Lentil")
png('Phen_lentil.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_lentil.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow, Leaf=leaf, layout =c(2,3))
dev.off()

# camelina  
pcamelina2 <- merge(plotlevel,pcamelina, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
summary(pcamelina2)

pcamelina2$moreLeaves <- as.Date(pcamelina2$moreLeaves, "%d.%m.%Y" )
pcamelina2$startInflorescence <- as.Date(pcamelina2$startInflorescence, "%d.%m.%Y" )
pcamelina2$startGrain <- as.Date(pcamelina2$startGrain, "%d.%m.%Y" )

str(pcamelina2)

pcamelina2$sp_numb <- factor(pcamelina2$sp_numb , levels=c("single", "one", "two", "four"))
pcamelina2$SP_NUMB <- revalue(pcamelina2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

leaf <- bwplot(SP_NUMB~moreLeaves|Country, data=pcamelina2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Camelina")
flow <- bwplot(SP_NUMB~startInflorescence|Country, data=pcamelina2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Camelina")
frui <- bwplot(SP_NUMB~startGrain|Country, data=pcamelina2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Camelina")
png('Phen_camelina.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_camelina.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow,Leaf=leaf,  layout =c(2,3))
dev.off()

# linseed  
plinseed2 <- merge(plotlevel,plinseed, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
summary(plinseed2)

plinseed2$sideShoot <- as.Date(plinseed2$sideShoot, "%d.%m.%Y" )
plinseed2$flowering <- as.Date(plinseed2$flowering, "%d.%m.%Y" )
plinseed2$endFruiting <- as.Date(plinseed2$endFruiting, "%d.%m.%Y" )

str(plinseed2)

plinseed2$sp_numb <- factor(plinseed2$sp_numb , levels=c("single", "one", "two", "four"))
plinseed2$SP_NUMB <- revalue(plinseed2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

shoot <- bwplot(SP_NUMB~sideShoot|Country, data=plinseed2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Linseed")
flow <- bwplot(SP_NUMB~flowering|Country, data=plinseed2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Linseed")
frui <- bwplot(SP_NUMB~endFruiting|Country, data=plinseed2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Linseed")
png('Phen_linseed.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_linseed.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow,Shoot=shoot,  layout =c(2,3))
dev.off()

# quinoa  
pquinoa2 <- merge(plotlevel,pquinoa, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
summary(pquinoa2)

pquinoa2$firstPairLeaves <- as.Date(pquinoa2$firstPairLeaves, "%d.%m.%Y" )
pquinoa2$flowerPyramids <- as.Date(pquinoa2$flowerPyramids, "%d.%m.%Y" )
pquinoa2$startGrain <- as.Date(pquinoa2$startGrain, "%d.%m.%Y" )

str(pquinoa2)

pquinoa2$sp_numb <- factor(pquinoa2$sp_numb , levels=c("single", "one", "two", "four"))
pquinoa2$SP_NUMB <- revalue(pquinoa2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

leaf <- bwplot(SP_NUMB~firstPairLeaves|Country, data=pquinoa2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Quinoa")
flow <- bwplot(SP_NUMB~flowerPyramids|Country, data=pquinoa2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Quinoa")
frui <- bwplot(SP_NUMB~startGrain|Country, data=pquinoa2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Quinoa")
png('Phen_quinoa.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_quinoa.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow, Leaf=leaf, layout =c(2,3))
dev.off()

# coriander  
pcoriander2 <- merge(plotlevel,pcoriander, by.x=c("plot","Country"), by.y=c("plot","country"), all=F)
summary(pcoriander2)

pcoriander2$firstTrueLeaf <- as.Date(pcoriander2$firstTrueLeaf, "%d.%m.%Y" )
pcoriander2$flowering <- as.Date(pcoriander2$flowering, "%d.%m.%Y" )
pcoriander2$fruiting <- as.Date(pcoriander2$fruiting, "%d.%m.%Y" )

str(pcoriander2)

pcoriander2$sp_numb <- factor(pcoriander2$sp_numb , levels=c("single", "one", "two", "four"))
pcoriander2$SP_NUMB <- revalue(pcoriander2$sp_numb, c("single"="single plant", "one"="monoculture", "two"="2-species mixture", "four"="4-species mixture"))

leaf <- bwplot(SP_NUMB~firstTrueLeaf|Country, data=pcoriander2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Coriander")
flow <- bwplot(SP_NUMB~flowering|Country, data=pcoriander2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Coriander")
frui <- bwplot(SP_NUMB~fruiting|Country, data=pcoriander2, xlim=as.Date(c("2018-02-01", "2018-08-30")), xlab = "", main="Coriander")
png('Phen_coriander.png', width = 7, height = 5, units = 'in', res = 300)
jpeg('Phen_coriander.jpg', width = 7, height = 5, units = 'in', res = 300)
c(Fruiting=frui, Flowering=flow, Leaf=leaf, layout =c(2,3))
dev.off()