--- title: "Increasing crop rotational diversity can enhance cereal yields" author: "Alessio Costa and Monique Smith" date: '2023-02-27' output: html_document: df_print: paged pdf_document: default word_document: default --- > sessionInfo() R version 4.2.2 (2022-10-31) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Ventura 13.2.1 Matrix products: default LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] parameters_0.20.2 ggpubr_0.6.0 emmeans_1.8.4-1 merTools_0.5.2 arm_1.13-1 [6] MASS_7.3-58.1 readr_2.1.4 readxl_1.4.2 ggrepel_0.9.3 gridExtra_2.3 [11] tidyr_1.3.0 plyr_1.8.8 effects_4.2-2 glmmTMB_1.1.5 DHARMa_0.4.6 [16] MuMIn_1.47.1 sjstats_0.18.2 dplyr_1.1.0 reshape2_1.4.4 plotly_4.10.1 [21] ggeffects_1.1.5 car_3.1-1 carData_3.0-5 xlsx_0.6.5 lattice_0.20-45 [26] ggplot2_3.4.1 listviewer_3.0.0 lmerTest_3.1-3 lme4_1.1-31 Matrix_1.5-1 loaded via a namespace (and not attached): [1] minqa_1.2.5 colorspace_2.1-0 ggsignif_0.6.4 ellipsis_0.3.2 [5] sjlabelled_1.2.0 estimability_1.4.1 rstudioapi_0.14 listenv_0.9.0 [9] furrr_0.3.1 fansi_1.0.4 mvtnorm_1.1-3 codetools_0.2-18 [13] splines_4.2.2 knitr_1.42 sjmisc_2.8.9 jsonlite_1.8.4 [17] nloptr_2.0.3 rJava_1.0-6 broom_1.0.3 broom.mixed_0.2.9.4 [21] shiny_1.7.4 compiler_4.2.2 httr_1.4.4 backports_1.4.1 [25] fastmap_1.1.0 lazyeval_0.2.2 survey_4.1-1 cli_3.6.0 [29] later_1.3.0 htmltools_0.5.4 tools_4.2.2 coda_0.19-4 [33] gtable_0.3.1 glue_1.6.2 Rcpp_1.0.10 cellranger_1.1.0 [37] vctrs_0.5.2 nlme_3.1-160 iterators_1.0.14 insight_0.19.0 [41] xfun_0.37 stringr_1.5.0 globals_0.16.2 xlsxjars_0.6.1 [45] mime_0.12 lifecycle_1.0.3 rstatix_0.7.2 future_1.31.0 [49] scales_1.2.1 promises_1.2.0.1 hms_1.1.2 parallel_4.2.2 [53] TMB_1.9.2 yaml_2.3.7 stringi_1.7.12 bayestestR_0.13.0 [57] foreach_1.5.2 blme_1.0-5 boot_1.3-28 rlang_1.0.6 [61] pkgconfig_2.0.3 evaluate_0.20 purrr_1.0.1 htmlwidgets_1.6.1 [65] cowplot_1.1.1 tidyselect_1.2.0 parallelly_1.34.0 magrittr_2.0.3 [69] R6_2.5.1 generics_0.1.3 DBI_1.1.3 pillar_1.8.1 [73] withr_2.5.0 survival_3.4-0 datawizard_0.6.5 abind_1.4-5 [77] nnet_7.3-18 tibble_3.1.8 performance_0.10.2 modelr_0.1.10 [81] utf8_1.2.3 tzdb_0.3.0 rmarkdown_2.20 grid_4.2.2 [85] data.table_1.14.8 forcats_1.0.0 digest_0.6.31 xtable_1.8-4 [89] httpuv_1.6.9 numDeriv_2016.8-1.1 stats4_4.2.2 munsell_0.5.0 [93] viridisLite_0.4.1 mitools_2.4 egg_0.4.5 ```{r libraries} library(listviewer) library(ggplot2) library(lattice) library(xlsx) library(car) library(lme4) #citation(package = "lme4") library(lmerTest) library(ggeffects) #citation(package = "ggeffects") library(plotly) library(reshape2) library(dplyr) library(sjstats) library(MuMIn) library(DHARMa) #citation(package = "DHARMa") library(glmmTMB) library(effects) library(plyr) library(tidyr) library(gridExtra) library(ggrepel) library(readxl) library(readr) library(merTools) library(emmeans) library(ggpubr) library(parameters) ``` #Data preparation ```{r data, message=FALSE, warning=FALSE, include=FALSE} # load data # setwd("~/Library...) yield.df <- as.data.frame(read.delim("yield_anomalies_publication.tsv", header = TRUE, sep = "\t", dec = ".")) yield.df[sapply(yield.df, is.character)] <- lapply(yield.df[sapply(yield.df, is.character)], as.factor) yield.df$year.f=factor(yield.df$year.f) yield.df$F_types=factor(yield.df$F_types) yield.df$year.f <- as.factor(yield.df$year) #### yield standardisation #### # we have not included the raw yields but you can see below how the standardised yields were calculated. # yield.df$site.crop <- paste(yield.df$site, yield.df$crop,sep=".") # yield.df$site.crop <- as.factor(yield.df$site.crop) # # sym=aggregate(yield.df$yield, list(yield.df$site.crop), mean) # # yield.df$s.yield=0 # # # for (i in 1:length(sym[,1])) # { # yield.df[which(yield.df$site.crop==sym[i,1]),"s.yield"]=yield.df[which(yield.df$site.crop==sym[i,1]),"yield"]-sym[i,2] # } # # # names(yield.df)[names(yield.df) == 'yield'] <- 'r.yield' # # nrow(yield.df) # yield observations 27460 # #site years # yield.df$site.time <- paste(yield.df$site, yield.df$time,sep=".") # yield.df$site.time <- as.factor(yield.df$site.time) # str(yield.df$site.time) #957 #### Split the data into Spring and Winter crops #### Spring.df <- subset(yield.df, phenology == "spring") Winter.df <- subset(yield.df, phenology == "winter") # Subset Spring data into US and Europe #### Spring.df.US <- subset(Spring.df, country == "North America") Spring.df.EU <- subset(Spring.df, country != "North America") ``` #Model selection species diversity ```{r} #detach("package:lmerTest", unload = TRUE) #necessary to look at AICs #Maize model Spring.df.US$fert=relevel(Spring.df.US$fert, ref="low") #Setting fertiliser reference level as low mX<- lmer(s.yield ~ D*time+I(D^2)*time+D*I(time^2)+I(D^2)*I(time^2)+fert+D:fert+I(D^2)*fert + #Full model (1|site:group) + (1|site:year.f), data=Spring.df.US) drop1(mX, test="Chisq") #all model terms significantly improved the model, no removal necessary library(lmerTest) #Saving the final model mXf= lmer(s.yield ~ D*time+I(D^2)*time+D*I(time^2)+I(D^2)*I(time^2)+fert+D:fert+I(D^2)*fert + (1|site:group) + (1|site:year.f), data=Spring.df.US) #Spring small grain cereal model Spring.df.EU$fert=relevel(Spring.df.EU$fert, ref="low") detach("package:lmerTest", unload = TRUE) sX<- lmer(s.yield ~ D*time+I(D^2)*time+D*I(time^2)+I(D^2)*I(time^2)+D*fert+I(D^2)*fert + (1|site:group) + (1|site:year.f), data=Spring.df.EU) drop1(sX, test="Chisq") #Removing several terms improved the model, the interaction of the highest order will be removed sX1<-update (sX, .~. -I(D^2):I(time^2)) #updating the model by removing unnecessary terms drop1(sX1, test="Chisq") #removing D:I(time^2) improved the model sX2<-update (sX1, .~. -D:I(time^2) ) drop1(sX2, test="Chisq") #removing I(time^2) improved the model sX3<-update (sX2, .~. -I(time^2) ) drop1(sX3, test="Chisq") #no further removal necessary library(lmerTest) #Saving the final model sX= lmer(s.yield ~ D*time+I(D^2)*time+D*I(time^2)+I(D^2)*I(time^2)+fert+D:fert+I(D^2)*fert + (1|site:group) + (1|site:year.f), data=Spring.df.EU) sXf<-update (sX, .~. -I(D^2):I(time^2) -D:I(time^2) -I(time^2) ) #Winter small grain cereal model Winter.df$fert=relevel(Winter.df$fert, ref="low") detach("package:lmerTest", unload = TRUE) wX<- lmer(s.yield ~ D*time+I(D^2)*time+D*I(time^2)+I(D^2)*I(time^2)+D*fert+I(D^2)*fert + (1|site:group) + (1|site:year.f), data=Winter.df) drop1(wX, test="Chisq") #Removing several terms improved the model, the interaction of the highest order will be removed wX1<-update (wX, .~. #updating the model by removing unnecessary terms -I(D^2):I(time^2) ) drop1(wX1, test="Chisq") #no further removal necessary wXf=wX1 #final model D + time + I(D^2) + I(time^2) + fert + D:time + time:I(D^2) + D:I(time^2) + D:fert + I(D^2):fert + (1 | site:group) + (1 | site:year.f) library(lmerTest) #Saving the final model wX<- lmer(s.yield ~ D*time+I(D^2)*time+D*I(time^2)+I(D^2)*I(time^2)+D*fert+I(D^2)*fert + (1|site:group) + (1|site:year.f), data=Winter.df) wXf<-update (wX, .~. #updating the model by removing unnecessary terms -I(D^2):I(time^2) ) ``` #Model diagnostic species diversity Checking assumptions ```{r} #Maize model simulationOutput <- simulateResiduals(fittedModel = mXf, n=1000) plot(simulationOutput, asFactor =F) ``` ```{r} #Spring small grain cereals model simulationOutput <- simulateResiduals(fittedModel = sXf, n=1000) plot(simulationOutput, asFactor =F) ``` Significant deviation from KS test, however the plot suggests that there is a good fit between observed and expected quantiles. Since sample size is large, even slight deviations can be detected as significant. No action taken. ```{r} #Winter small grain cereals model simulationOutput <- simulateResiduals(fittedModel = wXf, n=1000) plot(simulationOutput, asFactor =F) ``` Significant deviation from KS test, however the plot suggests that there is a good fit between observed and expected quantiles. No action taken. #Model selection functional richness ```{r} #Maize model detach("package:lmerTest", unload = TRUE) Spring.df.US$fert=relevel(Spring.df.US$fert, ref="low") mX<- lmer(s.yield ~ F_types*time+F_types*I(time^2)+fert+F_types:fert+ (1|site:group) + (1|site:year.f), data=Spring.df.US) drop1(mX, test="Chisq") #no removal necessary library(lmerTest) #saving the final model mXf_f=lmer(s.yield ~ F_types*time+F_types*I(time^2)+fert+F_types:fert+ (1|site:group) + (1|site:year.f), data=Spring.df.US) #Spring small grain cereal model detach("package:lmerTest", unload = TRUE) Spring.df.EU$fert=relevel(Spring.df.EU$fert, ref="low") sX<- lmer(s.yield ~ F_types*time+F_types*I(time^2)+fert+F_types:fert+ (1|site:group) + (1|site:year.f), data=Spring.df.EU) drop1(sX, test="Chisq") #removing F_types:I(time^2) improved the model sX1<-update (sX, .~. -F_types:I(time^2) ) drop1(sX1, test="Chisq") sX2<-update (sX1, .~. #removing I(time^2) improved the model -I(time^2) ) drop1(sX2, test="Chisq") #no further removal necessary library(lmerTest) #saving the final model sX<- lmer(s.yield ~ F_types*time+F_types*I(time^2)+fert+F_types:fert+ (1|site:group) + (1|site:year.f), data=Spring.df.EU) sXf_f=update (sX, .~. -F_types:I(time^2) -I(time^2) ) #Winter small grain cereals detach("package:lmerTest", unload = TRUE) Winter.df$fert=relevel(Winter.df$fert, ref="low") wX<- lmer(s.yield ~F_types*time+F_types*I(time^2)+fert+F_types:fert+ (1|site:group) + (1|site:year.f), data=Winter.df) drop1(wX, test="Chisq") #no term removal necessary #saving the final model library(lmerTest) wXf_f=lmer(s.yield ~F_types*time+F_types*I(time^2)+fert+F_types:fert+ (1|site:group) + (1|site:year.f), data=Winter.df) ``` #Model diagnostic functional richness Checking assumptions ```{r} #Maize model simulationOutput <- simulateResiduals(fittedModel = mXf_f, n=1000) plot(simulationOutput, asFactor =F) ``` No issue detected. ```{r} #Spring small grain cereals model simulationOutput <- simulateResiduals(fittedModel = sXf_f, n=1000) plot(simulationOutput, asFactor =F) ``` Significant deviation from KS test, however the plot suggests that there is a good fit between observed and expected quantiles. Since sample size is large, even slight deviations can be detected as significant. No action taken. ```{r} #Winter small grain cereals model simulationOutput <- simulateResiduals(fittedModel = wXf_f, n=1000) plot(simulationOutput, asFactor =F) ``` Significant deviation from KS test, however the plot suggests that there is a good fit between observed and expected quantiles. No action taken. There was a significant deviation from the outlier test. Identifying observations that produced a standard deviation higher than 2.5. ```{r} res1=resid(wXf_f, type = "pearson") Winter.df[which(abs(res1) > 2.5),] #which(abs(res1) > 2.5) returns the row number of the observations that produced sd>2.5 ``` This resulted in 29 outliers from English, Italian, Polish and Swedish sites. Below you can see the models summary with and without outliers. ```{r} W.NO=Winter.df[-which(abs(res1) > 2.5),] WIIf=lmer(s.yield ~F_types*time+F_types*I(time^2)+fert+F_types:fert+ (1|site:group) + (1|site:year.f), data=W.NO) model_parameters(wXf_f, effects="fixed", summary=T) model_parameters(WIIf, effects="fixed", summary=T) ``` Removing the outliers does not lead to dramatic changes in the estimates. The outliers will be kept in the analyses. # Main figures ## Figure 1 - maps- not included ##Figure 2 ```{r} #Species diversity Spring.US.effect1 <- ggeffect(mXf, terms = c("D[1:4.6, by=0.1]", "time [5,20,35]","fert")) Spring.US.ref1 <- ggeffect(mXf, terms = c("D[1,2]", "time [0,1]","fert")) Spring.US.effect1$conf.low <-Spring.US.effect1$conf.low-Spring.US.ref1$predicted[1] Spring.US.effect1$conf.high <-Spring.US.effect1$conf.high-Spring.US.ref1$predicted[1] Spring.US.effect1$predicted <-Spring.US.effect1$predicted-Spring.US.ref1$predicted[1] Winter.effect1 <- ggeffect(wXf, terms = c("D[1:6, by=0.1]", "time [5,20,35]","fert")) Winter.ref1 <- ggeffect(wXf, terms = c("D[1,2]", "time [0,1]","fert")) Winter.effect1$conf.low <-Winter.effect1$conf.low-Winter.ref1$predicted[1] Winter.effect1$conf.high <-Winter.effect1$conf.high-Winter.ref1$predicted[1] Winter.effect1$predicted <-Winter.effect1$predicted-Winter.ref1$predicted[1] Spring.EU.effect1 <- ggeffect(sXf, terms = c("D[1:6, by=0.1]", "time [5,20,35]","fert")) Spring.EU.ref1 <- ggeffect(sXf, terms = c("D[1,2]", "time [0,1]","fert")) Spring.EU.effect1$conf.low <-Spring.EU.effect1$conf.low-Spring.EU.ref1$predicted[1] Spring.EU.effect1$conf.high <-Spring.EU.effect1$conf.high-Spring.EU.ref1$predicted[1] Spring.EU.effect1$predicted <-Spring.EU.effect1$predicted-Spring.EU.ref1$predicted[1] #Functional richness Spring.US.effect2 <- ggeffect(mXf_f, terms = c("F_types", "time [5,20,35]","fert")) Spring.US.ref2 <- ggeffect(mXf_f, terms = c("F_types[1,2]", "time [0,1]","fert")) Spring.US.effect2$conf.low <-Spring.US.effect2$conf.low-Spring.US.ref2$predicted[1] Spring.US.effect2$conf.high <-Spring.US.effect2$conf.high-Spring.US.ref2$predicted[1] Spring.US.effect2$predicted <-Spring.US.effect2$predicted-Spring.US.ref2$predicted[1] Spring.EU.effect2 <- ggeffect(sXf_f, terms = c("F_types", "time [5,20,35]","fert")) Spring.EU.ref2 <- ggeffect(sXf_f, terms = c("F_types[1,2]", "time [0,1]","fert")) Spring.EU.effect2$conf.low <-Spring.EU.effect2$conf.low-Spring.EU.ref2$predicted[1] Spring.EU.effect2$conf.high <-Spring.EU.effect2$conf.high-Spring.EU.ref2$predicted[1] Spring.EU.effect2$predicted <-Spring.EU.effect2$predicted-Spring.EU.ref2$predicted[1] Winter.effect2 <- ggeffect(wXf_f, terms = c("F_types", "time [5,20,35]","fert")) Winter.ref2 <- ggeffect(wXf_f, terms = c("F_types[1,2]", "time [0,1]","fert")) Winter.effect2$conf.low <-Winter.effect2$conf.low-Winter.ref2$predicted[1] Winter.effect2$conf.high <-Winter.effect2$conf.high-Winter.ref2$predicted[1] Winter.effect2$predicted <-Winter.effect2$predicted-Winter.ref2$predicted[1] fert_labeller <- as_labeller(c("low" = "Low fertilisation", #Needed for fertilisation facet label "high" = "High fertilisation" )) Spring.EU.a=ggplot(Spring.EU.effect1, aes(x = x, y = predicted, colour = group)) + geom_line(aes(linetype = group), size = 0.8, show.legend = FALSE) + geom_ribbon( aes(ymin = conf.low, ymax = conf.high, fill = group, color = NULL), alpha = .15, show.legend = FALSE) + theme_bw() + scale_color_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'))+ scale_fill_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'), name="Time (yr)") + scale_linetype_manual(values=c("twodash", "dotted","solid"), name="Time (yr)")+ labs(colour = "Time (yr)", x = "Species Diversity", y = "Yield benefit (t/ha)", title = "Spring small grain cereals" )+ ylim(-0.5,3)+ facet_grid(~factor(facet,levels=c("low","high")), labeller = fert_labeller) Spring.US.a=ggplot(Spring.US.effect1, aes(x = x, y = predicted, colour = group)) + geom_line(aes(linetype = group), size = 0.8, show.legend = FALSE) + geom_ribbon( aes(ymin = conf.low, ymax = conf.high, fill = group, color = NULL), alpha = .15, show.legend = FALSE) + theme_bw() + scale_color_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'))+ scale_fill_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'), name="Time (yr)") + scale_linetype_manual(values=c("twodash", "dotted","solid"), name="Time (yr)")+ labs(colour = "Time (yr)", x = "Species Diversity", y = "Yield benefit (t/ha)", title = "Maize" )+ xlim(1,6)+ # ylim(-0.5,4.5)+ facet_grid(~factor(facet,levels=c("low","high")), labeller = fert_labeller) Winter.a=ggplot(Winter.effect1, aes(x = x, y = predicted, colour = group)) + geom_line(aes(linetype = group), size = 0.8) + geom_ribbon( aes(ymin = conf.low, ymax = conf.high, fill = group, color = NULL), alpha = .15) + theme_bw() + scale_color_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'))+ scale_fill_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'), name="Time (yr)") + scale_linetype_manual(values=c("twodash", "dotted","solid"), name="Time (yr)")+ labs(colour = "Time (yr)", x = "Species Diversity", y = "Yield benefit (t/ha)", title = "Winter small grain cereals" )+ ylim(-0.5,3)+ facet_grid(~factor(facet,levels=c("low","high")), labeller = fert_labeller) pd <- position_dodge(width=0.2) Spring.EU.b=ggplot(Spring.EU.effect2, aes(x = x, y = predicted, colour = group)) + geom_point(position=pd, show.legend = FALSE) + geom_errorbar(aes(ymin = conf.low, ymax=conf.high, color = group), width=.15,size=.5,position=pd, show.legend = FALSE) + theme_bw() + scale_color_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'))+ scale_fill_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'), name="Time (yr)") + labs( colour = "Time (yr)", x = "Functional Richness", y = "Yield benefit (t/ha)", title = "" ) + theme( plot.title = element_text(vjust=-52,hjust=1.3)# move title to bottom right corner )+ ylim(-0.4,3)+ facet_grid(~factor(facet,levels=c("low","high")), labeller = fert_labeller) blank=ggplot(Spring.EU.effect2, aes(x = x, y = predicted, colour = group)) + geom_blank() Spring.US.b=blank+ geom_point(Spring.US.effect2, inherit.aes = TRUE, mapping=aes(x = x, y = predicted, colour = group),position=pd, show.legend = FALSE) + geom_errorbar(Spring.US.effect2, inherit.aes = TRUE,mapping=aes(ymin = conf.low, ymax=conf.high, color = group), width=.15,size=.5,position=pd, show.legend = FALSE) + theme_bw() + scale_color_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'))+ scale_fill_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'), name="Time (yr)") + labs( colour = "Time (yr)", x = "Functional Richness", y = "Yield benefit (t/ha)", title = "" ) + theme( plot.title = element_text(vjust=-52,hjust=1.3)# move title to bottom right corner )+ ylim(-0.4,5.5)+ facet_grid(~factor(facet,levels=c("low","high")), labeller = fert_labeller) blank=ggplot(Spring.EU.effect2, aes(x = x, y = predicted, colour = group))+ facet_grid(~factor(facet,levels=c("low","high")), labeller = fert_labeller)+ geom_blank() Winter.b=blank+ geom_point(Winter.effect2, inherit.aes = TRUE, mapping=aes(x = x, y = predicted, colour = group),position=pd) + geom_errorbar(Winter.effect2, inherit.aes = TRUE,mapping=aes(ymin = conf.low, ymax=conf.high, color = group), width=.15,size=.5,position=pd) + theme_bw() + scale_color_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'))+ scale_fill_manual(values=c('#66c2a5','#fc8d62', '#8da0cb'), name="Time (yr)") + labs( colour = "Time (yr)", x = "Functional Richness", y = "Yield benefit (t/ha)", title = "" ) + ylim(-0.4,3)+ theme( plot.title = element_text(vjust=-52,hjust=1.3)# move title to bottom right corner ) egg::ggarrange(Spring.EU.a,Spring.US.a,Winter.a,Spring.EU.b,Spring.US.b,Winter.b,nrow=2, ncol=3) ``` ##Fig 3 ###Preparing required data subset - Species Diversity ```{r} for (t in c(5,20,35)) #loop code to generate lsmeans estimates at Time snapshot of 5, 20, 35 years { #generating lsmean estimates for every 0.01 point of D - needed to extract yield-maximizing diversity at selected time snapshot and HIGH fertilisation print("Spring small grain cereals") hsls=as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(3,3.5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans) hsls=rbind(hsls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(3.5,4,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hsls=rbind(hsls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(4.,4.5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hsls=rbind(hsls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(4.5,5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hsls=rbind(hsls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(5,5.5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hsls=rbind(hsls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(5.5,6,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hsmax=hsls[which(hsls$lsmean==max(hsls$lsmean)),] #extracting the lsmean estimate at the yield-maximizing diversity print("Maize") #repeating the same procedure for maize hmls=as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=seq(3,3.5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans) hmls=rbind(hmls,as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=seq(3.5,4,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hmls=rbind(hmls,as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=seq(4.,4.57,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hmmax=hmls[which(hmls$lsmean==max(hmls$lsmean)),] print("Winter small grain cereals") #repeating the same procedure for winter cereals hwls=as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(3,3.5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans) hwls=rbind(hwls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(3.5,4,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hwls=rbind(hwls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(4.,4.5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hwls=rbind(hwls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(4.5,5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hwls=rbind(hwls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(5,5.5,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hwls=rbind(hwls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(5.5,6,by=0.01), time=t,fert=c("high")), adjust="TUKEY")$lsmeans)) hwmax=hwls[which(hwls$lsmean==max(hwls$lsmean)),] #Repeating the procedures (all 3 indicator crops) under LOW fertilisation print("Spring small grain cereals") sls=as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(3,3.5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans) sls=rbind(sls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(3.5,4,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) sls=rbind(sls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(4.,4.5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) sls=rbind(sls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(4.5,5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) sls=rbind(sls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(5,5.5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) sls=rbind(sls,as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=seq(5.5,6,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) smax=sls[which(sls$lsmean==max(sls$lsmean)),] print("Maize") mls=as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=seq(3,3.5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans) mls=rbind(mls,as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=seq(3.5,4,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) mls=rbind(mls,as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=seq(4.,4.57,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) mmax=mls[which(mls$lsmean==max(mls$lsmean)),] print("Winter small grain cereals") wls=as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(3,3.5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans) wls=rbind(wls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(3.5,4,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) wls=rbind(wls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(4.,4.5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) wls=rbind(wls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(4.5,5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) wls=rbind(wls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(5,5.5,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) wls=rbind(wls,as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=seq(5.5,6,by=0.01), time=t,fert=c("low")), adjust="TUKEY")$lsmeans)) wmax=wls[which(wls$lsmean==max(wls$lsmean)),] #extracting lsmeans estimate for the baseline monoculture (sref,mref,wref), and for high-input monocultures (shigh,mhigh,whigh) sref=as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=1, time=0,fert=c("low")), adjust="TUKEY")$lsmeans) shigh=as.data.frame(lsmeans(sXf, pairwise~D*time+D*fert, at=list(D=1, time=t,fert=c("high")), adjust="TUKEY")$lsmeans) mref=as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=1, time=0,fert=c("low")), adjust="TUKEY")$lsmeans) mhigh=as.data.frame(lsmeans(mXf, pairwise~D*time+D*fert, at=list(D=1, time=t,fert=c("high")), adjust="TUKEY")$lsmeans) wref=as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=1, time=0,fert=c("low")), adjust="TUKEY")$lsmeans) whigh=as.data.frame(lsmeans(wXf, pairwise~D*time+D*fert, at=list(D=1, time=t,fert=c("high")), adjust="TUKEY")$lsmeans) #Dataset assemblage#### stest=rbind( cbind( "SSGC", smax[,c("D","time","fert")], smax[,c("lsmean","asymp.LCL","asymp.UCL")]-sref[,"lsmean"] #Each reported estimate is compared with the baseline (sref,mref,wref) ), cbind( "SSGC", hsmax[,c("D","time","fert")], hsmax[,c("lsmean","asymp.LCL","asymp.UCL")]-sref[,"lsmean"] ), cbind("SSGC", shigh[,c("D","time","fert")], shigh[,c("lsmean","asymp.LCL","asymp.UCL")]-sref[,"lsmean"] ) ) colnames(stest)[1]="ind" mtest=rbind( cbind( "Maize", mmax[,c("D","time","fert")], mmax[,c("lsmean","asymp.LCL","asymp.UCL")]-mref[,"lsmean"] ), cbind( "Maize", hmmax[,c("D","time","fert")], hmmax[,c("lsmean","asymp.LCL","asymp.UCL")]-mref[,"lsmean"] ), cbind( "Maize", mhigh[,c("D","time","fert")], mhigh[,c("lsmean","asymp.LCL","asymp.UCL")]-mref[,"lsmean"] ) ) colnames(mtest)[1]="ind" wtest=rbind( cbind( "WSGC", wmax[,c("D","time","fert")], wmax[,c("lsmean","asymp.LCL","asymp.UCL")]-wref[,"lsmean"] ), cbind( "WSGC", hwmax[,c("D","time","fert")], hwmax[,c("lsmean","asymp.LCL","asymp.UCL")]-wref[,"lsmean"] ), cbind( "WSGC", whigh[,c("D","time","fert")], whigh[,c("lsmean","asymp.LCL","asymp.UCL")]-wref[,"lsmean"] ) ) colnames(wtest)[1]="ind" if (t==5) { finaltest=rbind(stest,mtest,wtest) } else {finaltest=rbind(finaltest,stest,mtest,wtest)} } colnames(finaltest)[2]="CRD_val" fert_est=cbind(finaltest,"SD",c("CRD", "CRD_Fert","Fertilisation")) colnames(fert_est)[8]="CRD" colnames(fert_est)[9]="effect.of" ``` ###Preparing required data subset - Functional Richness ```{r} #Repeating the previous procedure for functional richness for (t in c(5,20,35)) { print("Spring small grain cereals") sls=as.data.frame(lsmeans(sXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types=c("1","2","3","4"), time=t,fert=c("low")), adjust="TUKEY")$lsmeans) hsls=as.data.frame(lsmeans(sXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types=c("1","2","3","4"), time=t,fert=c("high")), adjust="TUKEY")$lsmeans) smax=sls[which(sls$lsmean==max(sls$lsmean)),] hsmax=hsls[which(hsls$lsmean==max(hsls$lsmean)),] print("Maize") mls=as.data.frame(lsmeans(mXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types=c("1","2","3"), time=t,fert=c("low")), adjust="TUKEY")$lsmeans) hmls=as.data.frame(lsmeans(mXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types=c("1","2","3"), time=t,fert=c("high")), adjust="TUKEY")$lsmeans) mmax=mls[which(mls$lsmean==max(mls$lsmean)),] hmmax=hmls[which(hmls$lsmean==max(hmls$lsmean)),] print("Winter small grain cereals") wls=as.data.frame(lsmeans(wXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types=c("1","2","3"), time=t,fert=c("low")), adjust="TUKEY")$lsmeans) hwls=as.data.frame(lsmeans(wXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types=c("1","2","3"), time=t,fert=c("high")), adjust="TUKEY")$lsmeans) wmax=wls[which(wls$lsmean==max(wls$lsmean)),] hwmax=hwls[which(hwls$lsmean==max(hwls$lsmean)),] sref=as.data.frame(lsmeans(sXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types="1", time=0,fert=c("low")), adjust="TUKEY")$lsmeans) shigh=as.data.frame(lsmeans(sXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types="1", time=t,fert=c("high")), adjust="TUKEY")$lsmeans) mref=as.data.frame(lsmeans(mXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types="1", time=0,fert=c("low")), adjust="TUKEY")$lsmeans) mhigh=as.data.frame(lsmeans(mXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types="1", time=t,fert=c("high")), adjust="TUKEY")$lsmeans) wref=as.data.frame(lsmeans(wXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types="1", time=0,fert=c("low")), adjust="TUKEY")$lsmeans) whigh=as.data.frame(lsmeans(wXf_f, pairwise~F_types*time+F_types*fert, at=list(F_types="1", time=t,fert=c("high")), adjust="TUKEY")$lsmeans) stest=rbind( cbind( "SSGC", smax[,c("F_types","time","fert")], smax[,c("lsmean","asymp.LCL","asymp.UCL")]-sref[,"lsmean"] ), cbind( "SSGC", hsmax[,c("F_types","time","fert")], hsmax[,c("lsmean","asymp.LCL","asymp.UCL")]-sref[,"lsmean"] ), cbind("SSGC", shigh[,c("F_types","time","fert")], shigh[,c("lsmean","asymp.LCL","asymp.UCL")]-sref[,"lsmean"] ) ) colnames(stest)[1]="ind" mtest=rbind( cbind( "Maize", mmax[,c("F_types","time","fert")], mmax[,c("lsmean","asymp.LCL","asymp.UCL")]-mref[,"lsmean"] ), cbind( "Maize", hmmax[,c("F_types","time","fert")], hmmax[,c("lsmean","asymp.LCL","asymp.UCL")]-mref[,"lsmean"] ), cbind( "Maize", mhigh[,c("F_types","time","fert")], mhigh[,c("lsmean","asymp.LCL","asymp.UCL")]-mref[,"lsmean"] ) ) colnames(mtest)[1]="ind" wtest=rbind( cbind( "WSGC", wmax[,c("F_types","time","fert")], wmax[,c("lsmean","asymp.LCL","asymp.UCL")]-wref[,"lsmean"] ), cbind( "WSGC", hwmax[,c("F_types","time","fert")], hwmax[,c("lsmean","asymp.LCL","asymp.UCL")]-wref[,"lsmean"] ), cbind( "WSGC", whigh[,c("F_types","time","fert")], whigh[,c("lsmean","asymp.LCL","asymp.UCL")]-wref[,"lsmean"] ) ) colnames(wtest)[1]="ind" if (t==5) { finaltest=rbind(stest,mtest,wtest) } else {finaltest=rbind(finaltest,stest,mtest,wtest)} } colnames(finaltest)[2]="CRD_val" FR_Fert_est=cbind(finaltest,"FR",c("CRD", "CRD_Fert","Fertilisation")) colnames(FR_Fert_est)[8]="CRD" colnames(FR_Fert_est)[9]="effect.of" ``` ###Producing the graph ```{r} est.df.all=rbind(fert_est,FR_Fert_est) #combining SD and FR subsets est.df.all$ind <- factor(est.df.all$ind, levels=c("SSGC", "Maize", "WSGC")) est.df.all$effect.of <- factor(est.df.all$effect.of, levels=c("Fertilisation", "CRD", "CRD_Fert")) pd <- position_dodge(0.6) plot.est.all <- ggplot(est.df.all, aes(x=ind, y=lsmean, colour=effect.of, shape = CRD)) + geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL), width=.1, position=pd, size = 1) + geom_point(position=pd, size = 4) + theme_bw() + scale_color_manual(name = "With reference to low\nN fertilisation monoculture,\nyield gain by increasing\n ", values=c('#a6cee3','#1f78b4','#b2df8a'), labels=c("N fertilisation (low to high)\n ", "CRD (mono to\nyield-max CRD)\n ", "Both N fertilisation & CRD\n(low to high;\nmono to yield-max CRD" ))+ scale_shape_manual(name = "CRD metric", values=c(15, 17))+ labs(x = "Indicator crop", y = "Yield benefit (t/ha)") plot.est.all ``` # Extended data figures ##Fig S1 ```{r} yield.df$F_comb <- factor(yield.df$F_comb, levels = c("C", "C,B", "C,LEG", "C,LEY", "C,B,LEG","C,B,LEY", "C,LEG,LEY", "C,B,LEG,LEY")) ggplot(data=yield.df, aes(x = F_types, fill=F_comb))+ geom_bar(col=I("black"), alpha=0.7) + labs(x="Functional Richness", y="Number of observations")+ theme_bw() ``` ##Fig S2 ```{r} #### Figure S2 - FR vs SD #### diversity_summary_plot1=ggplot(yield.df, aes(x = D, y = F_types)) + #geom_point()+ geom_count() + theme_bw() + labs(x = "Species diversity", y = "Functional richness"#, #title = "Species diversity vs functional richness" ) diversity_summary_plot1 ``` ##Fig S3 ```{r} # create the data needed # winter small grain cereals Winter.df=Winter.df %>% mutate(t_s=case_when( time<=quantile(time,.25) ~ "<=9 yr", time>=quantile(time,.75) ~ ">=28 yr", timequantile(time,.25) ~ "9 yr=28 yr")) Yield_means_fit_winter <- Winter.df %>% group_by(lte, fert, D, t_s) %>% dplyr::summarize(meanY=mean(s.yield), sdY = sd(s.yield)) Winter.effect1 <- ggeffect(wXf, terms = c("D[1:6, by=0.1]", "time [5,20,35]","fert")) Winter.effect1$t_s <- ifelse(Winter.effect1$group== '5', "<=9 yr", ifelse(Winter.effect1$group== '20', "9 yr=28 yr")) Winter.effect1$t_s <- ordered(Winter.effect1$t_s, levels = c("<=9 yr", "9 yr=28 yr")) Winter.effect1$s.yield <- Winter.effect1$predicted Winter.effect1$D <- Winter.effect1$x Winter.effect1$fert <- Winter.effect1$facet #spring small grain cereals Spring.df.EU=Spring.df.EU %>% mutate(t_s=case_when( time<=quantile(time,.25) ~ "<=9 yr", time>=quantile(time,.75) ~ ">=28 yr", timequantile(time,.25) ~ "9 yr=28 yr")) Yield_means_fit_spring <- Spring.df.EU %>% group_by(lte, fert, D, t_s) %>% dplyr::summarize(meanY=mean(s.yield), sdY = sd(s.yield)) Spring.EU.effect1 <- ggeffect(sXf, terms = c("D[1:6, by=0.1]", "time [5,20,35]","fert")) Spring.EU.effect1$t_s <- ifelse(Spring.EU.effect1$group== '5', "<=9 yr", ifelse(Spring.EU.effect1$group== '20', "9 yr=28 yr")) Spring.EU.effect1$t_s <- ordered(Spring.EU.effect1$t_s, levels = c("<=9 yr", "9 yr=28 yr")) Spring.EU.effect1$s.yield <- Spring.EU.effect1$predicted Spring.EU.effect1$D <- Spring.EU.effect1$x Spring.EU.effect1$fert <- Spring.EU.effect1$facet # maize Spring.df.US=Spring.df.US %>% mutate(t_s=case_when( time<=quantile(time,.25) ~ "<=9 yr", time>=quantile(time,.75) ~ ">=28 yr", timequantile(time,.25) ~ "9 yr=28 yr")) Yield_means_fit_maize <- Spring.df.US %>% group_by(lte, fert, D, t_s) %>% dplyr::summarize(meanY=mean(s.yield), sdY = sd(s.yield)) Spring.US.effect1 <- ggeffect(mXf, terms = c("D[1:4.6, by=0.1]", "time [5,20,35]","fert")) Spring.US.effect1$t_s <- ifelse(Spring.US.effect1$group== '5', "<=9 yr", ifelse(Spring.US.effect1$group== '20', "9 yr=28 yr")) Spring.US.effect1$t_s <- ordered(Spring.US.effect1$t_s, levels = c("<=9 yr", "9 yr=28 yr")) Spring.US.effect1$s.yield <- Spring.US.effect1$predicted Spring.US.effect1$D <- Spring.US.effect1$x Spring.US.effect1$fert <- Spring.US.effect1$facet # set up plots fert_time_labeller <- as_labeller(c("low" = "Low fertilisation", #Needed for fertilisation facet label "high" = "High fertilisation", "<=9 yr"="<=9 yr", "9 yr=28 yr"=">=28 yr" )) pd <- position_dodge(0.2) # make plots winter.sites3 <- ggplot(Yield_means_fit_winter, aes(x = D, y=meanY)) + geom_jitter(data=Winter.df,aes(D, s.yield,col=lte,group=lte),alpha = 0.3, size=1,width = 0.05, show.legend = F)+ geom_errorbar(aes(ymin=meanY-sdY, ymax=meanY+sdY,group=lte), width=.1, position=pd) + geom_point(aes(fill=lte), position=pd, size=2, shape=21) + geom_line(data = Winter.effect1, aes(x=D, y=s.yield), colour="black", show.legend = F) + # geom_ribbon(data = Spring.US.effect1, aes(ymin = conf.low, ymax = conf.high), # show.legend=FALSE, alpha = .5) + facet_grid(t_s~factor(fert,levels=c("low","high")), labeller = fert_time_labeller) + theme_bw() + xlim(1,6)+ labs(fill = "", x = "Species Diversity", y = "Standardised yield (t/ha)", title = "Winter small grain cereals" )+guides(fill = guide_legend(override.aes = list(shape = 21, size=3))) spring.sites3 <- ggplot(Yield_means_fit_spring, aes(x = D, y=meanY)) + geom_jitter(data=Spring.df.EU,aes(D, s.yield,col=lte,group=lte),alpha = 0.3, size=1,width = 0.05, show.legend = F)+ geom_errorbar(aes(ymin=meanY-sdY, ymax=meanY+sdY,group=lte), width=.1, position=pd) + geom_point(aes(fill=lte), position=pd, size=2, shape=21) + geom_line(data = Spring.EU.effect1, aes(x=D, y=s.yield), colour="black", show.legend = F) + # geom_ribbon(data = Spring.US.effect1, aes(ymin = conf.low, ymax = conf.high), # show.legend=FALSE, alpha = .5) + facet_grid(t_s~factor(fert,levels=c("low","high")), labeller = fert_time_labeller) + theme_bw() + xlim(1,6)+ labs(fill = "", x = "Species Diversity", y = "Standardised yield (t/ha)", title = "Spring small grain cereals" )+guides(fill = guide_legend(override.aes = list(shape = 21, size=3))) maize.sites4 <- ggplot(Yield_means_fit_maize, aes(x = D, y=meanY)) + geom_jitter(data=Spring.df.US,aes(D, s.yield,col=lte,group=lte),alpha = 0.3, size=1,width = 0.05, show.legend = F)+ geom_errorbar(aes(ymin=meanY-sdY, ymax=meanY+sdY,group=lte), width=.1, position=pd) + geom_point(aes(fill=lte), position=pd, size=2, shape=21) + geom_line(data = Spring.US.effect1, aes(x=D, y=s.yield), colour="black", show.legend = F) + # geom_ribbon(data = Spring.US.effect1, aes(ymin = conf.low, ymax = conf.high), # show.legend=FALSE, alpha = .5) + facet_grid(t_s~factor(fert,levels=c("low","high")), labeller = fert_time_labeller) + theme_bw() + xlim(1,6)+ labs(fill = "", x = "Species Diversity", y = "Standardised yield (t/ha)", title = "Maize" )+guides(fill = guide_legend(override.aes = list(shape = 21, size=3))) # combine plots # ggarrange(spring.sites3, maize.sites4, winter.sites3, # ncol = 3, legend = "top") combined_fig <- ggarrange(spring.sites3 + rremove("ylab") + rremove("xlab"), maize.sites4 + rremove("ylab") + rremove("xlab"), winter.sites3 + rremove("ylab") + rremove("xlab"), labels = NULL, ncol = 3, legend = "top") require(grid) # for the textGrob() function annotate_figure(combined_fig, left = textGrob("Mean-centred yield (t/ha)", rot = 90, vjust = 1, gp = gpar(cex = 1)), bottom = textGrob("Species diversity", gp = gpar(cex = 1))) ``` ## Supplementary Tables ##Tab S4 and S5 ```{r} #Spring small grain cereals, species diversity model_parameters(sXf, effects="fixed", summary=T) #Spring small grain cereals, functional richness model_parameters(sXf_f, effects="fixed", summary=T) #Maize, species diversity model_parameters(mXf, effects="fixed", summary=T) #Maize, functional richness model_parameters(mXf_f, effects="fixed", summary=T) #Winter small grain cereals, species diversity model_parameters(wXf, effects="fixed", summary=T) #Winter small grain cereals, functional richness model_parameters(wXf_f, effects="fixed", summary=T) ``` #Tab S5 ```{r} mod.winter<- lmer(s.yield ~ ley+leg+bl+ (1|site:group) + (1|year.f), data=Winter.df) mod.spring<- lmer(s.yield ~ ley+leg+bl+ (1|site:group) + (1|year.f), data=Spring.df.EU) mod.maize<- lmer(s.yield ~ ley+leg+bl+ (1|site:group) + (1|year.f), data=Spring.df.US) model_parameters(mod.spring, effects="fixed", summary=T) model_parameters(mod.maize, effects="fixed", summary=T) model_parameters(mod.winter, effects="fixed", summary=T) ```