PaNDiv experiment, competition network analysis
Here will be presented the data analysis attached with “Fast-slow
traits predict competition network structure and its response to
resources and enemies.”
In order to run this file, the
attached R script “functions_competition_network.R” must be present in
the “code” subfolder of the working directory, and all data must be
present in a subfolder called “data”. An R environment
“all_competition_matrices.Rdata” is available in the “results”
subfolder, for easy loading, but each matrix is also individually
available if needed.
To avoid long running time, several Rdata
environments were also saved in the “result” folder, they can be loaded
in the corresponding parts of the script.
WARNING: before knitting, make sure the knit directory
is the project directory and not the document directory, it can be
changed in the knit options.
0. Loading packages, working directory
library("rmdformats") # this package needs to be loaded to knit
library("tidyr") # data wrangling
library("dplyr") # data wrangling
library("readxl") # useful to open data for supp fig. 8
library("Rmisc") # data wrangling for supp figure 8
library("stringr") # data wrangling
library("ggplot2") # plotting package
library("ggthemes") # plotting package
library("ggpubr") # plotting package, ggarrange and getting intercepts from models
library("reshape") # data wrangling, for melt.array function
library("reshape2") # data wrangling
library("igraph") # package for plotting networks in a igraph object
library("doParallel") # for parallel computing
library("parallel") # for parallel computing
library("moments") # skewness package
library("bipartite") # useful for network handling
library("GiniWegNeg") # best definition of Gini with negative values
library("lme4") # linear model package
library("jtools") # good presentation of data
library("ciTools") # used for prediction of lmer with random effects
library("sjPlot") # package for tab_model function
library("modelsummary") # another package for summarising model outputs because sjPlot decided one day to stop working1. Example figures of network metrics
First, we computed examples of matrix and igraph corresponding to a high value of each network metric (see Fig. 1).
a. Diagonal dominance (D)
Strong diagonal dominance was simulated with high values of intraspecific competition (in the range of our observed competition coefficients), then we built a matrix and an igraph in order to illustrate the example.
rich<-6 # species richness level used for example
diagexample<-matrix(-runif(324),
nrow = 18) # sampling random numbers for interaction coefficients
diag(diagexample)<-sample(-12:-3, 18, replace = TRUE) # sampling random numbers for
# strong intraspecific competition
diagexample<-diagexample[1:rich,1:rich] # building the matrix
rownames(diagexample)<-c("A","B","C","D","E","F") # naming the example species
colnames(diagexample)<-c("A","B","C","D","E","F")
plot_diagonal<-melt.array(diagexample) # formatting in the long format for ggplot
plot_diagonal$value<-c(diagexample)
diagplot<-ggplot(plot_diagonal,
aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",
mid = "white",
high="red",
limits=c(-12,4)) +
labs(x="species",
y="species",
title="Diagonal dominance example") +
theme_bw() +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank())+
scale_y_discrete(limits = rev)+
scale_x_discrete(position = "top") # plotting the diagonal dominance matrix example
# with an explicit color code : purple = negative, orange = positive interactions
#example igraph
species<-colnames(diagexample) # species in a variable
weight<-c(t(diagexample)) # define the weight of the paths
species2<-rep(species,lenght.out = length(weight))
species1<-rep(species,each=rich)
network<-data.frame(species1,species2) # explicitly building all paths names
g<-graph_from_data_frame(network)
E(g)$weight<-weight
E(g)$color<-"mediumorchid4"
E(g)[weight>=0]$color<-"sienna1"plot.igraph(g,
edge.width=abs(weight),
edge.arrow.size=0.1,
edge.arrow.width=1.5,
edge.curved=0.2,
edge.color=E(g)$color,
vertex.size=35,
vertex.label.cex=1,
vertex.label.color="black",
vertex.color="white",
layout=layout.circle,
mark.border = "white")
b. Asymmetry (A)
We then computed examples following the same method for all other metrics.
asymexample<-matrix(runif(324,-12,-5),nrow = 18)
diag(asymexample)<--0.00001
asymexample[lower.tri(asymexample, diag = F)]<-runif(153,-1,0)
asymexample<-asymexample[1:rich,1:rich]
rownames(asymexample)<-c("A","B","C","D","E","F")
colnames(asymexample)<-c("A","B","C","D","E","F")
plot_asymmetry<-melt.array(asymexample)
plot_asymmetry$value<-c(asymexample)
asymplot<-ggplot(plot_asymmetry,
aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",
mid = "white",
high="red",
limits=c(-12,4)) +
labs(x="species",
y="species",
title="Asymmetry example") +
theme_bw() +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank())+
scale_y_discrete(limits = rev)+
scale_x_discrete(position = "top")
#example igraph
species<-colnames(asymexample)
weight<-c(t(asymexample))
sp2<-rep(species,lenght.out = length(weight))
sp1<-rep(species,each=rich)
network<-data.frame(sp1,sp2)
g<-graph_from_data_frame(network)
E(g)$weight<-weight
E(g)$color<-"mediumorchid4"
E(g)[weight>=0]$color<-"sienna1"
plot.igraph(g,
edge.width=abs(weight),
edge.arrow.size=0.1,
edge.arrow.width=1.5,
edge.curved=0.2,
edge.color=E(g)$color,
vertex.size=35,
vertex.label.cex=1,
vertex.label.color="black",
vertex.color="white",
layout=layout.circle,
mark.border = "white")
c. Skewness (γ)
skewexample<-matrix(runif(324, -12,4),nrow = 18)
skewexample<-skewexample[1:6,1:6]
rownames(skewexample)<-c("A","B","C","D","E","F")
colnames(skewexample)<-c("A","B","C","D","E","F")
plot_skewness<-melt.array(skewexample)
plot_skewness$value<-c(skewexample)
skewplot<-ggplot(plot_skewness,
aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",
mid = "white",
high="red",
limits=c(-12,4)) +
labs(x="species",
y="species",
title="Skewness example") +
theme_bw() + theme(axis.text.x=element_blank(),
axis.text.y=element_blank())+
scale_y_discrete(limits = rev)+
scale_x_discrete(position = "top")
#example igraph
species<-colnames(skewexample)
weight<-c(t(skewexample))
sp2<-rep(species,
lenght.out = length(weight))
sp1<-rep(species,
each=rich)
network<-data.frame(sp1,sp2)
g<-graph_from_data_frame(network)
E(g)$weight<-weight
E(g)$color<-"mediumorchid4"
E(g)[weight>=0]$color<-"sienna1"
plot.igraph(g,
edge.width=abs(weight),
edge.arrow.size=0.1,
edge.arrow.width=1.5,
edge.curved=0.2,
edge.color=E(g)$color,
vertex.size=35,
vertex.label.cex=1,
vertex.label.color="black",
vertex.color="white",
layout=layout.circle,
mark.border = "white")
d. Evenness (1-G)
giniexample<-matrix(-runif(324),nrow = 18)
giniexample<-giniexample[1:6,1:6]
rownames(giniexample)<-c("A","B","C","D","E","F")
colnames(giniexample)<-c("A","B","C","D","E","F")
plot_gini<-melt.array(giniexample)
plot_gini$value<-c(giniexample)
giniplot<-ggplot(plot_gini, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",
mid = "white",
high="red",
limits=c(-12,4)) +
labs(x="species",
y="species",
title="Gini (evenness) example") +
theme_bw() + theme(axis.text.x=element_blank(),
axis.text.y=element_blank())+
scale_y_discrete(limits = rev)+
scale_x_discrete(position = "top")
#example igraph
species<-colnames(giniexample)
weight<-c(t(giniexample))
sp2<-rep(species,lenght.out = length(weight))
sp1<-rep(species,each=rich)
network<-data.frame(sp1,sp2)
g<-graph_from_data_frame(network)
E(g)$weight<-weight
E(g)$color<-"mediumorchid4"
E(g)[weight>=0]$color<-"sienna1"
plot.igraph(g,
edge.width=abs(weight),
edge.arrow.size=0.1,
edge.arrow.width=1.5,
edge.curved=0.2,
edge.color=E(g)$color,
vertex.size=35,
vertex.label.cex=1,
vertex.label.color="black",
vertex.color="white",
layout=layout.circle,
mark.border = "white")
e. Modularity (M)
moduexample<-matrix(runif(324,-12,-2),nrow = 18)
moduexample[1,]<-runif(18, -0.5,0)
moduexample[,1]<-runif(18, -0.5,0)
moduexample[6,]<-runif(18, -0.5,0)
moduexample[,6]<-runif(18, -0.5,0)
diag(moduexample)<--0.0001
moduexample<-moduexample[1:6,1:6]
rownames(moduexample)<-c("A","B","C","D","E","F")
colnames(moduexample)<-c("A","B","C","D","E","F")
plot_modularity<-melt.array(moduexample)
plot_modularity$value<-c(moduexample)
moduplot<-ggplot(plot_modularity, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",
mid = "white",
high="red",
limits=c(-12,4)) +
labs(x="species",
y="species",
title="Modularity example") +
theme_bw() +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank())+
scale_y_discrete(limits = rev)+
scale_x_discrete(position = "top")
#example igraph
species<-colnames(moduexample)
weight<-c(t(moduexample))
sp2<-rep(species,lenght.out = length(weight))
sp1<-rep(species,each=rich)
network<-data.frame(sp1,sp2)
g<-graph_from_data_frame(network)
E(g)$weight<-weight
# define best clusters using cluster_otpima function
clu<-cluster_optimal(g,weights = abs(E(g)$weight))
indclass<-data.frame(values=species,ind=as.character(clu$membership)) # save cluster number for each species
cols <- c("1"="salmon2", "2"="seashell2") # assign colors for the two clusters
V(g)$color <- cols[indclass$ind[match(V(g)$name,indclass$values)]] # cluster colors goes in a variable
E(g)$color<-"mediumorchid4"
E(g)[weight>=0]$color<-"sienna1"
plot.igraph(g,
edge.width=abs(weight),
edge.arrow.size=0.1,
edge.arrow.width=1.5,
edge.curved=0.2,
edge.color=E(g)$color,
vertex.size=35,
vertex.label.cex=1,
vertex.label.color="black",
mark.groups=clu,
mark.col=cols,
layout=layout.circle,
mark.border = "white")
2. Network metrics for communities of 5 species
First, we computed our 5 network metrics for every community of 5 PaNDiv plant species. We then ran multimembership models and extracted effects of treatment, SLA and SLA variance on the structure of the networks.
a. Loading competition matrices
Here is the list of the competition matrices and their transformation (all coefficients were divided by the standard errors obtained from the glmmTMB models), used for the rest of the analysis.
# open matrices and set dataframe
alpha_june_control<-read.table("data/biommatrix_june_2020_control.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_june_control<-read.table("data/SE_biommatrix_june_2020_control.txt",header = T, sep = "\t", row.names = 1)
alpha_june_control<-alpha_june_control/SE_alpha_june_control
alpha_june_control<-as.data.frame(alpha_june_control)
alpha_june_nitrogen<-read.table("data/biommatrix_june_2020_nitrogen.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_june_nitrogen<-read.table("data/SE_biommatrix_june_2020_nitrogen.txt",header = T, sep = "\t", row.names = 1)
alpha_june_nitrogen<-alpha_june_nitrogen/SE_alpha_june_nitrogen
alpha_june_nitrogen<-as.data.frame(alpha_june_nitrogen)
alpha_june_fungicide<-read.table("data/biommatrix_june_2020_fungicide.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_june_fungicide<-read.table("data/SE_biommatrix_june_2020_fungicide.txt",header = T, sep = "\t", row.names = 1)
alpha_june_fungicide<-alpha_june_fungicide/SE_alpha_june_fungicide
alpha_june_fungicide<-as.data.frame(alpha_june_fungicide)
alpha_june_combined<-read.table("data/biommatrix_june_2020_combined.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_june_combined<-read.table("data/SE_biommatrix_june_2020_combined.txt",header = T, sep = "\t", row.names = 1)
alpha_june_combined<-alpha_june_combined/SE_alpha_june_combined
alpha_june_combined<-as.data.frame(alpha_june_combined)
alpha_june_control<-as.matrix(alpha_june_control)
alpha_june_nitrogen<-as.matrix(alpha_june_nitrogen)
alpha_june_fungicide<-as.matrix(alpha_june_fungicide)
alpha_june_combined<-as.matrix(alpha_june_combined)
alpha_august_control<-read.table("data/biommatrix_august_2020_control.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_august_control<-read.table("data/SE_biommatrix_august_2020_control.txt",header = T, sep = "\t", row.names = 1)
alpha_august_control<-alpha_august_control/SE_alpha_august_control
alpha_august_control<-as.data.frame(alpha_august_control)
alpha_august_nitrogen<-read.table("data/biommatrix_august_2020_nitrogen.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_august_nitrogen<-read.table("data/SE_biommatrix_august_2020_nitrogen.txt",header = T, sep = "\t", row.names = 1)
alpha_august_nitrogen<-alpha_august_nitrogen/SE_alpha_august_nitrogen
alpha_august_nitrogen<-as.data.frame(alpha_august_nitrogen)
alpha_august_fungicide<-read.table("data/biommatrix_august_2020_fungicide.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_august_fungicide<-read.table("data/SE_biommatrix_august_2020_fungicide.txt",header = T, sep = "\t", row.names = 1)
alpha_august_fungicide<-alpha_august_fungicide/SE_alpha_august_fungicide
alpha_august_fungicide<-as.data.frame(alpha_august_fungicide)
alpha_august_combined<-read.table("data/biommatrix_august_2020_combined.txt",header = T, sep = "\t", row.names = 1)
SE_alpha_august_combined<-read.table("data/SE_biommatrix_august_2020_combined.txt",header = T, sep = "\t", row.names = 1)
alpha_august_combined<-alpha_august_combined/SE_alpha_august_combined
alpha_august_combined<-as.data.frame(alpha_august_combined)
alpha_august_control<-as.matrix(alpha_august_control)
alpha_august_nitrogen<-as.matrix(alpha_august_nitrogen)
alpha_august_fungicide<-as.matrix(alpha_august_fungicide)
alpha_august_combined<-as.matrix(alpha_august_combined)
save.image(file = "results/all_competition_matrices.Rdata")For an easy load, you can just load “all_competition_matrices.Rdata”
We then plotted the interaction matrices in order to obtain Fig. 1 (here we present August matrices as well).
# simple re ordering of the species
alpha_june_control<-alpha_june_control[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_june_control<-alpha_june_control[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
alpha_june_nitrogen<-alpha_june_nitrogen[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_june_nitrogen<-alpha_june_nitrogen[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
alpha_june_fungicide<-alpha_june_fungicide[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_june_fungicide<-alpha_june_fungicide[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
alpha_june_combined<-alpha_june_combined[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_june_combined<-alpha_june_combined[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
alpha_august_control<-alpha_august_control[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_august_control<-alpha_august_control[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
alpha_august_nitrogen<-alpha_august_nitrogen[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_august_nitrogen<-alpha_august_nitrogen[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
alpha_august_fungicide<-alpha_august_fungicide[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_august_fungicide<-alpha_august_fungicide[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
alpha_august_combined<-alpha_august_combined[c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt"),]
alpha_august_combined<-alpha_august_combined[,c("Am","Dc","Fr","Pm","Be","Hp","Sp","Pg","Cj","Ao",
"To","Cb","Lp","Ra","Dg","Hl","Ga","Pt")]
#define the color of each species to differentiate fast from slow growing species
colorsp<-c("turquoise4","turquoise4","turquoise4","turquoise4","turquoise4","turquoise4","turquoise4",
"turquoise4","turquoise4","turquoise4","chartreuse4","chartreuse4","chartreuse4","chartreuse4","chartreuse4","chartreuse4",
"chartreuse4","chartreuse4")
# first, format the matrix in the long format for ggplot2
plot_alpha_june_control<-melt.array(alpha_june_control)
plot_alpha_june_control$value<-c(alpha_june_control)
plot_alpha_june_control<-plot_alpha_june_control[plot_alpha_june_control$value!=0,]
# then we plotted the interaction matrix for June Control
contju<-ggplot(plot_alpha_june_control,
aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",
mid = "white",
high="red",
limits=c(-12,4)) +
labs(x="species",
y="species",
title="June Control") +
theme_bw() +
theme(axis.text.x=element_text(size=9,
angle=0,
vjust=0.3,
color = colorsp),
axis.text.y=element_text(size=9,
color = rev(colorsp)),
plot.title=element_text(size=11))+
scale_y_discrete(limits = rev)+
scale_x_discrete(position = "top")
# then, similarly for the other treatments and season...
plot_alpha_june_nitrogen<-melt.array(alpha_june_nitrogen)
plot_alpha_june_nitrogen<-plot_alpha_june_nitrogen[plot_alpha_june_nitrogen$value!=0,]
niju<-ggplot(plot_alpha_june_nitrogen, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",mid = "white", high="red", limits=c(-12,4)) +
labs(x="species", y="species", title="June Nitrogen") +
theme_bw() + theme(axis.text.x=element_text(size=9,
angle=0,
vjust=0.3,
color = colorsp),
axis.text.y=element_text(size=9,
color = rev(colorsp)),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.title=element_text(size=11),)+
scale_y_discrete(limits = rev)+scale_x_discrete(position = "top")
plot_alpha_june_fungicide<-melt.array(alpha_june_fungicide)
plot_alpha_june_fungicide<-plot_alpha_june_fungicide[plot_alpha_june_fungicide$value!=0,]
funju<-ggplot(plot_alpha_june_fungicide, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",mid = "white", high="red", limits=c(-12,4)) +
labs(x="species", y="species", title="June Fungicide") +
theme_bw() + theme(axis.text.x=element_text(size=9,
angle=0,
vjust=0.3,
color = colorsp),
axis.ticks.x = element_blank(),
axis.text.y=element_text(size=9, color = rev(colorsp)),
axis.title.x = element_blank(),
plot.title=element_text(size=11))+
scale_y_discrete(limits = rev)+scale_x_discrete(position = "top")
plot_alpha_june_combined<-melt.array(alpha_june_combined)
plot_alpha_june_combined<-plot_alpha_june_combined[plot_alpha_june_combined$value!=0,]
combju<-ggplot(plot_alpha_june_combined, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",mid = "white", high="red", limits=c(-12,4)) +
labs(x="species", y="species", title="June Combined (Nitrogen + Fungicide)") +
theme_bw() + theme(axis.text.x=element_text(size=9,
angle=0,
vjust=0.3,
color = colorsp),
axis.text.y=element_text(size=9,
color = rev(colorsp)),
axis.ticks.y = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title=element_text(size=11))+
scale_y_discrete(limits = rev)+scale_x_discrete(position = "top")
plot_alpha_august_control<-melt.array(alpha_august_control)
plot_alpha_august_control<-plot_alpha_august_control[plot_alpha_august_control$value!=0,]
contau<-ggplot(plot_alpha_august_control, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",mid = "white", high="red", limits=c(-12,4)) +
labs(x="species", y="species", title="August control") +
theme_bw() + theme(axis.text.x=element_text(size=9, angle=0, vjust=0.3, color = colorsp),
axis.text.y=element_text(size=9, color = rev(colorsp)),
plot.title=element_text(size=11))+
scale_y_discrete(limits = rev)+scale_x_discrete(position = "top")
plot_alpha_august_nitrogen<-melt.array(alpha_august_nitrogen)
plot_alpha_august_nitrogen<-plot_alpha_august_nitrogen[plot_alpha_august_nitrogen$value!=0,]
niau<-ggplot(plot_alpha_august_nitrogen, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",mid = "white", high="red", limits=c(-12,4)) +
labs(x="species", y="species", title="August nitrogen") +
theme_bw() + theme(axis.text.x=element_text(size=9, angle=0, vjust=0.3, color = colorsp),
axis.text.y=element_text(size=9, color = rev(colorsp)),
axis.title.y = element_blank(),
plot.title=element_text(size=11))+
scale_y_discrete(limits = rev)+scale_x_discrete(position = "top")
plot_alpha_august_fungicide<-melt.array(alpha_august_fungicide)
plot_alpha_august_fungicide<-plot_alpha_august_fungicide[plot_alpha_august_fungicide$value!=0,]
funau<-ggplot(plot_alpha_august_fungicide, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",mid = "white", high="red", limits=c(-12,4)) +
labs(x="species", y="species", title="August fungicide") +
theme_bw() + theme(axis.text.x=element_text(size=9, angle=0, vjust=0.3, color = colorsp),
axis.text.y=element_text(size=9, color = rev(colorsp)),
axis.title.x = element_blank(),
plot.title=element_text(size=11))+
scale_y_discrete(limits = rev)+scale_x_discrete(position = "top")
plot_alpha_august_combined<-melt.array(alpha_august_combined)
plot_alpha_august_combined<-plot_alpha_august_combined[plot_alpha_august_combined$value!=0,]
combau<-ggplot(plot_alpha_august_combined, aes(x = X1, y = X2)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient2(low="darkblue",mid = "white", high="red", limits=c(-12,4)) +
labs(x="species", y="species", title="August combined (nitrogen + fungicide)") +
theme_bw() + theme(axis.text.x=element_text(size=9, angle=0, vjust=0.3, color = colorsp),
axis.text.y=element_text(size=9, color = rev(colorsp)),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title=element_text(size=11))+
scale_y_discrete(limits = rev)+scale_x_discrete(position = "top") 
Using the same data on the raw plant interactions, we can plot supplementary figure 2 and 3, the histograms of plant interactions, respectively competitive effects and responses (all treatments combined).
#save the effects together in a dataframe
all_matrices_effect<-as.data.frame(rbind(t(alpha_june_control),t(alpha_june_nitrogen),t(alpha_june_fungicide),t(alpha_june_combined),
t(alpha_august_control),t(alpha_august_nitrogen),t(alpha_august_fungicide),t(alpha_august_combined)))
rownames(all_matrices_effect)<-NULL
all_matrices_effect<-all_matrices_effect %>%
pivot_longer(colnames(all_matrices_effect))
#highlight the extremes in effect (high)
all_matrices_effect$extreme<-with(all_matrices_effect, ifelse(name=="Hl", 0,1))
all_matrices_effect$extreme<-as.factor(all_matrices_effect$extreme)
hist_comp_effect<-ggplot(all_matrices_effect, aes(x=value, fill=extreme))+
geom_histogram(color="white")+
facet_wrap(~ name)+
geom_vline(xintercept = mean(all_matrices_effect$value), linetype="dashed", linewidth= 1, color="blue")+
theme_classic()+
theme(legend.position = "none")
hist_comp_effect
#save the responses together in a dataframe
all_matrices_response<-as.data.frame(rbind(alpha_june_control,alpha_june_nitrogen,alpha_june_fungicide,alpha_june_combined,
alpha_august_control,alpha_august_nitrogen,alpha_august_fungicide,alpha_august_combined))
rownames(all_matrices_response)<-NULL
all_matrices_response<-all_matrices_response %>%
pivot_longer(colnames(all_matrices_response))
#highlight the extremes in response (low)
all_matrices_response$extreme<-with(all_matrices_response, ifelse(name=="Ra"|name=="Dc", 0,1))
all_matrices_response$extreme<-as.factor(all_matrices_response$extreme)
hist_comp_response<-ggplot(all_matrices_response, aes(x=value, fill=extreme))+
geom_histogram(color="white")+
facet_wrap(~ name)+
geom_vline(xintercept = mean(all_matrices_response$value), linetype="dashed", linewidth= 1, color="blue")+
theme_classic()+
theme(legend.position = "none")
hist_comp_response
b. Loading SLA data
In order to distinguish fast growing species from slow growing species, we used data on the specific leaf area gathered on PaNDiv experiment monocultures in 2019. We log transformed the value in order to remove the correlation between mean SLA and variance SLA.
#SLA for continuous value of fast versus slow growing species
sla<- read.table("data/SLA_PaNDiv_2016-2019.txt",header = T)
sla<-cbind(sla,str_split_fixed(sla$sampling.period, "/", 3)) # explicit date of sampling period
colnames(sla)[12:14]<-c("month","day","year")
sla$day<-NULL # remove useless data
sla2019<-na.omit(sla)# average sla across the 4 years #subset(sla,sla$year=="2019") # only using data collected in 2019
meansla2019nonlog<-sla2019 %>% group_by(Species) %>% summarise_at(vars(mean_SLA), list(value=mean)) #averaging both June and August SLA
meansla2019<-meansla2019nonlog # saving non log values in a variable, just in case
meansla2019$value<-log(meansla2019$value) # log transform SLA valuesc. All possible networks of 5 species
Here is a list of all PaNDiv species abbreviations. From there, we calculated all possible combo for a given network size (5 species). The network size is flexible, but the closer to 9 species, the longer the code takes to run due to the high number of possible combos.
sp<-c("Hp","Be","Ao","Dg","Fr","Hl","Lp","Pt","Am","Cj","Cb","To","Dc","Sp","Pg","Ra","Pm","Ga") # list of all species
slowsp<-c("Hp","Be", "Ao", "Fr", "Am", "Cj", "Dc", "Sp", "Pg", "Pm") # list of all slow growing species
fastsp<-c("Cb","Dg","Ga","Hl","Lp", "Pt", "Ra", "To") # list of all fast growing species
rich<-5 #assign a number between 1 and 18 for species richness of the sub-matrix
slowdat<-as.data.frame(combn(slowsp,rich)) # all slow combination of 5
fastdat<-as.data.frame(combn(fastsp,rich)) # all fast combination of 5
mixdat<-as.data.frame(combn(sp,rich)) # all combinations of 5 speciesThe following function file can be found in the “code” subfolder, it needs to be loaded for the rest of the analysis.
In order to optimize the time of code running, we removed the unnecessary combinations (only fast or only slow in mixed data), and we set parallel computing (packages “doparallel” and “parallel”).
#remove all fast and all slow network from mixed communities
mixdat<-removeidentical(mixdat,slowdat) # personal function from "function_competition_network.R"
mixdat<-removeidentical(mixdat,fastdat)
#set parallel computation
numCores <- detectCores() #number of cores in your computer
registerDoParallel(numCores/2-1) #important: define parallel computation to numCoresd. Calculate network metrics for each community
The following code calculates all network metrics for each combo of 5 species (or any other input of richness value). This code takes more than 30 min to run, if you are only interested in the output, you can load it directly from below.
# list of competition matrices for all treatments and seasons
all<-list(alpha_june_control,alpha_june_nitrogen,alpha_june_fungicide,alpha_june_combined,alpha_august_control,alpha_august_nitrogen,alpha_august_fungicide,alpha_august_combined)
names(all)<-c("June C", "June N", "June F", "June NF", "August C", "August N", "August F", "August NF")
#compute all network metrics for networks of mixed SLA
mixnet<-networkmetrics(mixdat,rich) # /!\ WARNING : takes 30 minutes to run (output directly available below)
speciesmixed<-mixnet[[2]]
mixnet<-mixnet[[1]]
mixnet$Communities<-"Mixed"
#compute all network metrics for networks of low SLA
slownet<-networkmetrics(slowdat,rich) # networkmetrics() is a personal function defined in "function_competition_network.R"
speciesslow<-slownet[[2]]
slownet<-slownet[[1]]
slownet$Communities<-"Slow"
#compute all network metrics for networks of high SLA
fastnet<-networkmetrics(fastdat,rich)
speciesfast<-fastnet[[2]]
fastnet<-fastnet[[1]]
fastnet$Communities<-"Fast"
save.image(file = "results/network_metrics_5_species.Rdata")Here is the environment where the output is saved for networks of 5 species:
e. Scaling variables, formatting treatments, and matrix of species presence/absence
Once we obtained the data frame with all network metrics for every community, we scaled the variable in order to compare their effect sizes in the multimembership models.
savecomp<-comp # save the competition data frame before scaling just in case
#scaling response variables
comp$SLAv<-scale(comp$SLAv, center = T)
comp$SLA<-scale(comp$SLA, center = T)We also transformed treatment in a two level factor (for both Nitrogen and Fungicide).
comp$Nitrogen<-NA
comp$Fungicide<-NA
for (i in 1:nrow(comp)) {
if (comp$Treatment[i] == "N" | comp$Treatment[i] =="NF") {
comp$Nitrogen[i]<- 1
}
else {
comp$Nitrogen[i]<- 0
}
if (comp$Treatment[i] == "F" | comp$Treatment[i] =="NF") {
comp$Fungicide[i]<- 1
}
else {
comp$Fungicide[i]<- 0
}
}
comp$Nitrogen<-as.factor(comp$Nitrogen)
comp$Fungicide<-as.factor(comp$Fungicide)And finally, we built a presence / absence matrix in order to run multimemberships models with presence / absence data as random factors.
speciesmat<-speciesmixed[rep(rownames(speciesmixed),8),] # in speciesmixed was saved the community species names, it is repeated 8 times because we have 4 treatment x 2 seasons = 8 matrices
fastmat<-speciesfast[rep(rownames(speciesfast),8),]
slowmat<-speciesslow[rep(rownames(speciesslow),8),]
speciesmat<-rbind(slowmat,speciesmat,fastmat)
rownames(speciesmat)<-NULL
# bind the presence/absence matrix to the overall data frame
comp<-cbind(comp,speciesmat)f. Multi-membership modelling applied to the communities
We then applied multimembership models to our “comp” dataset. To do so, we need to create a variable that will be a placeholder random factor for our presence / abscence matrix
#fake random factor with 18 levels, for 18 species
Species <- rep(LETTERS[1:18], length.out= nrow(comp))
comp$Species<-Species
mycol<-c("#000000" , "#6A9400", "#26CDFF", "#962EDC")Running the models and their prediction takes about 45 min, if you are just interested in the output, you can load the image file below.
This code presents the average effect across June and August 2020, but you can also look at the results on the two sampling periods separately and obtain the Supplementary figures 6 and 7.
- Diagonal dominance (D)
#Diagonal dominance
# WARNING : running multi-membership model + prediction takes ~ 10 min for each variable
# input from model formula, here my response is Diagonal dominance, SLA and SLA variance + treatment as explanatory variables, Species is random
lmod <- lFormula(Diagstrength~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
# I modified the random factor by introducing "speciesmat" which is a presence/absence matrix of PaNDiv species in my communities of 5 plant species.
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
# preparing to run the optimisation of the lmer
devfun <- do.call(mkLmerDevfun, lmod)
# optimising per se
opt <- optimizeLmer(devfun)
# run the optimisation and save the model into a variable called m1, m2,... to m5 -> important !
m2 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
season_datapred<-slapred(m2,comp) # takes two-three minutes to run, personal function defined in "function_competition_network.R"
#average across seasons for more clarity, remove to look at the two sampling periods separately
datapred<-aggregate(cbind(fit,lwr,upr)~SLA+Treatment,FUN=mean,data= season_datapred)
season_datavpred<-slavpred(m2,comp) # takes two-three minutes to run
#average across seasons for more clarity, remove to look at the two sampling periods separately
datavpred<-aggregate(cbind(fit,lwr,upr)~SLAv+Treatment,FUN=mean,data= season_datavpred)
season_diagpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=season_datapred)+
facet_grid(~Season)+ #look at the two sampling periods separately
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Diagonal dominance")+
ylim(c(-9,-4.5))+ #might need to change it if scaled
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
diagpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=datapred)+
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Diagonal dominance")+
ylim(c(-9,-4.5))+ #might need to change it if scaled
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
season_diagvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=season_datavpred)+
facet_grid(~Season)+ #look at the two sampling periods separately
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(-9,-4.5))+ #might need to change it if scaled
xlab(" ")+
ylab("Diagonal dominance")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
diagvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=datavpred)+
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(-9,-4.5))+ #might need to change it if scaled
xlab(" ")+
ylab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())- Asymmetry (A)
#Asymmetry
lmod <- lFormula(Asymmetry~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m3 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
season_datapred<-slapred(m3,comp)
#average across seasons for more clarity, remove to look at the two sampling periods separately
datapred<-aggregate(cbind(fit,lwr,upr)~SLA+Treatment,FUN=mean,data= season_datapred)
season_datavpred<-slavpred(m3,comp) #takes two-three minutes each to run
#average across seasons for more clarity, remove to look at the two sampling periods separately
datavpred<-aggregate(cbind(fit,lwr,upr)~SLAv+Treatment,FUN=mean,data= season_datavpred)
datapredb<-subset(datapred, datapred$Treatment=="Control"|datapred$Treatment=="Nitrogen"|datapred$Treatment=="Fungicide")
season_asympred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=season_datapred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(1.3,2.0))+ #might need to change it if scaled
ylab("Asymmetry")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
asympred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=datapred)+
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(1.3,2.0))+ #might need to change it if scaled
ylab("Asymmetry")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
season_asymvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=season_datavpred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(1.3,2.0))+ #might need to change it if scaled
ylab("Asymmetry")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
asymvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=datavpred)+
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(1.3,2.0))+ #might need to change it if scaled
ylab(" ")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())- Skewness (γ)
#Skewness
lmod <- lFormula(Skewness~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m4 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
season_datapred<-slapred(m4,comp)
#average across seasons for more clarity, remove to look at the two sampling periods separately
datapred<-aggregate(cbind(fit,lwr,upr)~SLA+Treatment,FUN=mean,data= season_datapred)
season_datavpred<-slavpred(m4,comp) #takes two-three minutes each to run
#average across seasons for more clarity, remove to look at the two sampling periods separately
datavpred<-aggregate(cbind(fit,lwr,upr)~SLAv+Treatment,FUN=mean,data= season_datavpred)
season_skewpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=season_datapred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(-1.5,0.5))+ #might need to change it if scaled
ylab("Skewness")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
skewpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=datapred)+
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(-1.5,0.5))+ #might need to change it if scaled
ylab("Skewness")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
season_skewvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=season_datavpred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(-2,0.5))+ #might need to change it if scaled
xlab(" ")+
ylab("Skewness")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
skewvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=datavpred)+
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylim(c(-2,0.5))+ #might need to change it if scaled
xlab(" ")+
ylab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())- Evenness (1-G)
#Evenness
lmod <- lFormula(1-Gini~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m5 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
season_datapred<-slapred(m5,comp)
#average across seasons for more clarity, remove to look at the two sampling periods separately
datapred<-aggregate(cbind(fit,lwr,upr)~SLA+Treatment,FUN=mean,data= season_datapred)
season_datavpred<-slavpred(m5,comp) #takes two-three minutes each to run
#average across seasons for more clarity, remove to look at the two sampling periods separately
datavpred<-aggregate(cbind(fit,lwr,upr)~SLAv+Treatment,FUN=mean,data= season_datavpred)
season_ginipred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=season_datapred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
ylim(c(0.4,0.65))+ #might need to change it if scaled
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Evenness")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
ginipred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=datapred)+
ylim(c(0.4,0.65))+ #might need to change it if scaled
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Evenness")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
season_ginivpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=season_datavpred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
ylim(c(0.4,0.65))+ #might need to change it if scaled
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Evenness")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
ginivpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=datavpred)+
ylim(c(0.4,0.65))+ #might need to change it if scaled
geom_line(size=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab(" ")+
xlab(" ")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())- Modularity (M)
#Modularity
lmod <- lFormula(Modularity~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
season_datapred<-slapred(m1,comp)
#average across seasons for more clarity, remove to look at the two sampling periods separately
datapred<-aggregate(cbind(fit,lwr,upr)~SLA+Treatment,FUN=mean,data=season_datapred)
season_datavpred<-slavpred(m1,comp) #takes two-three minutes each to run
#average across seasons for more clarity, remove to look at the two sampling periods separately
datavpred<-aggregate(cbind(fit,lwr,upr)~SLAv+Treatment,FUN=mean,data= season_datavpred)
season_modupred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=season_datapred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
#ylim(c(-2,2))+
geom_line(linewidth=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Modularity")+
ylim(c(0.05,0.2))+ #might need to change it if scaled
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
modupred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=datapred)+
#ylim(c(-2,2))+
geom_line(linewidth=0.8)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Modularity")+
ylim(c(0.05,0.2))+ #might need to change it if scaled
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
season_moduvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=season_datavpred)+
facet_grid(~Season)+ #add to look at the two sampling periods separately
geom_line(linewidth=0.8)+
#ylim(c(-2,2))+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Modularity")+
xlab("SLA variance")+
ylim(c(0.05,0.2))+ #might need to change it if scaled
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
moduvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=datavpred)+
geom_line(linewidth=0.8)+
#ylim(c(-2,2))+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab(" ")+
xlab("SLA variance")+
ylim(c(0.05,0.2))+ #might need to change it if scaled
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
#saving the models and predictions in an R environment
save.image(file = "results/multimembership_results_5_species.Rdata")Here is the direct output of the multimembership models and predictions:
With the model predictions we can plot Fig.2, the overall yearly effect of treatment and SLA on our network metrics.
We can also plot the two seasons separately to obtain figs. S6 and
S7:
Then, with the model outputs we can also plot the effects sizes (Fig. 3)
#scaling all variables to compare effect sizes
comp$Modularity<-scale(comp$Modularity, center = T)
comp$Diagstrength<-scale(comp$Diagstrength, center = T)
comp$Asymmetry<-scale(comp$Asymmetry, center = T)
comp$Skewness<-scale(comp$Skewness, center = T)
comp$Gini<-scale(comp$Gini, center = T)
#re run models with scaled variables
#Diagonal dominance
lmod <- lFormula(Diagstrength~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m2 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Asymmetry
lmod <- lFormula(Asymmetry~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m3 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Skewness
lmod <- lFormula(Skewness~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m4 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Evenness
lmod <- lFormula(1-Gini~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m5 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Modularity
lmod <- lFormula(Modularity~Nitrogen*Fungicide*(SLA+SLAv)+Season*(SLA+SLAv)+(1|Species), data=comp)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
# summarise the model outputs
models<- list(m2,m3,m4,m5,m1)
models<-dvnames(models)
# rename models with explicit network metric names
names(models)<-c("Diagonal dominance", "Asymmetry", "Skewness", "Evenness", "Modularity")
# re-ordering the explanatory variables and intercept and giving them better names
cm<-c("SLAv:SeasonAugust" = "SLA variance x August",
"SLA:SeasonAugust" = "SLA x August",
"SeasonAugust" ="August",
"Nitrogen1:Fungicide1:SLAv" ="Nitrogen x Fungicide x SLA variance",
"Fungicide1:SLAv" = "Fungicide x SLA variance",
"Nitrogen1:SLAv" = "Nitrogen x SLA variance",
"SLAv"= "SLA variance",
"Nitrogen1:Fungicide1:SLA" = "Nitrogen x Fungicide x SLA",
"Fungicide1:SLA" = "Fungicide x SLA",
"Nitrogen1:SLA" = "Nitrogen x SLA",
"SLA"="SLA",
"Nitrogen1:Fungicide1" = "Nitrogen x Fungicide",
"Fungicide1"="Fungicide",
"Nitrogen1"="Nitrogen")
# define color vector for plotting
cbp2 <- c( "black", "black", "black",
"black","black")
b <- list(geom_vline(xintercept = 0, color = "white"),
annotate("rect", alpha = .05,
xmin = -1, xmax = 1,
ymin = -Inf, ymax = Inf))
# saving plot in a g1 ggplot
g1<-modelplot(models,
coef_omit = "Interc",
coef_map = cm,
size=1,
background = b)+
aes(alpha = ifelse(estimate > 0.1 | estimate <(-0.1),
"Strong effect ( | x | > 0.1)",
"Weak effect ( | x | < 0.1)"),
shape = ifelse(p.value > 0.05,
"Not significant",
"Significant")) +
scale_colour_manual(values=cbp2)+
facet_grid(~model)+
geom_vline(xintercept=0)+
scale_alpha_manual(values = c(1,0.25))+
scale_shape_manual(values = c(4,19))+
xlim(c(-1,1))+
guides(color = "none")+
theme(text=element_text(size=20),
axis.text.y=element_text(size=15),
#panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.text.x=element_text(size=9),
axis.title.x.bottom = element_text(size=14, margin= margin(8,0,15,0)),
axis.title=element_text(size=18),
plot.title=element_text(size=20),
#legend.text=element_text(size=20),
#legend.title=element_text(size=20)
)
g1
modelist<- list(m2,m3,m4,m5,m1)
names(modelist) <- c("Diagonal dominance", "Asymmetry", "Skewness", "Evenness", "Modularity")
modelsummary(modelist, statistic = NULL,
estimate = "{estimate}{stars} ({p.value})",
fmt = fmt_statistic(estimate = 2, p.value = 3),
coef_rename = c("SLAv" = "SLA variance", "Nitrogen1" = "Nitrogen", "Fungicide1" = "Fungicide", "SeasonAugust"="Season"),
coef_omit = "(Intercept)",
stars = T)| Diagonal dominance | Asymmetry | Skewness | Evenness | Modularity | |
|---|---|---|---|---|---|
| Nitrogen | −0.25*** (0.000) | 0.19*** (0.000) | 0.23*** (0.000) | 0.31*** (0.000) | 0.68*** (0.000) |
| Fungicide | −0.28*** (0.000) | 0.00 (0.938) | 0.65*** (0.000) | 0.96*** (0.000) | −0.35*** (0.000) |
| SLA | 0.05 (0.792) | 0.21+ (0.072) | 0.09 (0.584) | −0.17 (0.235) | −0.23* (0.021) |
| SLA variance | 0.03*** (0.000) | −0.05*** (0.000) | −0.16*** (0.000) | 0.10*** (0.000) | 0.34*** (0.000) |
| Season | −0.01** (0.008) | −0.01* (0.049) | −0.04*** (0.000) | 0.07*** (0.000) | 0.78*** (0.000) |
| Nitrogen:Fungicide | 0.24*** (0.000) | −0.15*** (0.000) | −0.45*** (0.000) | −0.82*** (0.000) | −0.35*** (0.000) |
| Nitrogen:SLA | −0.02* (0.018) | −0.25*** (0.000) | −0.26*** (0.000) | 0.09*** (0.000) | 0.04*** (0.000) |
| Nitrogen:SLA variance | −0.16*** (0.000) | −0.11*** (0.000) | 0.33*** (0.000) | 0.11*** (0.000) | −0.32*** (0.000) |
| Fungicide:SLA | 0.06*** (0.000) | −0.06*** (0.000) | −0.20*** (0.000) | 0.00 (0.981) | 0.29*** (0.000) |
| Fungicide:SLA variance | −0.12*** (0.000) | 0.03** (0.001) | 0.43*** (0.000) | 0.11*** (0.000) | −0.37*** (0.000) |
| SLA:Season | 0.10*** (0.000) | −0.09*** (0.000) | 0.08*** (0.000) | −0.06*** (0.000) | 0.10*** (0.000) |
| SLA variance:Season | 0.14*** (0.000) | 0.03*** (0.000) | −0.11*** (0.000) | −0.29*** (0.000) | −0.05*** (0.000) |
| Nitrogen:Fungicide:SLA | 0.06*** (0.000) | 0.29*** (0.000) | 0.13*** (0.000) | 0.08*** (0.000) | 0.02 (0.175) |
| Nitrogen:Fungicide:SLA variance | −0.04*** (0.000) | 0.15*** (0.000) | −0.36*** (0.000) | −0.09*** (0.000) | 0.17*** (0.000) |
| SD (Observations) | 0.62 | 0.87 | 0.70 | 0.70 | 0.71 |
| Num.Obs. | 68976 | 68976 | 68976 | 68976 | 68976 |
| R2 Marg. | 0.089 | 0.038 | 0.143 | 0.225 | 0.381 |
| R2 Cond. | 0.360 | 0.112 | 0.316 | 0.355 | 0.434 |
| AIC | 130209.3 | 176524.5 | 146637.3 | 147047.4 | 148657.7 |
| BIC | 130364.7 | 176679.9 | 146792.7 | 147202.9 | 148813.1 |
| ICC | 0.3 | 0.1 | 0.2 | 0.2 | 0.1 |
| RMSE | 0.62 | 0.87 | 0.70 | 0.70 | 0.71 |
3. Other supplementary figures
a. Network metrics for communities of 5, 7, 11 and 15 species
In order to check whether the size of the network will impact the effect of treatment, season or SLA and SLA variance on our network metrics, we looked at the effect of several network size (5, 7, 11, or 15 species). This code is heavy to run and will take at least 4 hours to be fully computed. Once again, you can simply load the output in the R environment “results/richness_data_5-7-11-15-18_species.Rdata”.
n<-list(5,7,11,15) #assign numbers between 1 and 18 for sp richness of sub matrix
allrich<-NA
allrich<-lapply(n,networkrich) #takes one day to run
#saving the models and predictions in an R environment
save.image(file = "results/richness_data_5-7-11-15-18_species.Rdata")Now we have computed all possible combination of species for these 4 levels of species richness, we will once again run multimembership model, with network size being computed as an ordered factor.
This code is also long to run, around 5 hours, the output can be loaded directly below in the “multimembership_results_richness_5-7-11-15_species.Rdata”.
#format the table
allrich_data<-melt(allrich, id=colnames(allrich[[1]]))
allrich_data$L1<-NULL
#scale the response variables
allrich_data$SLA<-scale(allrich_data$SLA, center = T)
allrich_data$SLAv<-scale(allrich_data$SLAv, center = T)
#adding a two level factor for treatment
allrich_data$Nitrogen<-NA
allrich_data$Fungicide<-NA
for (i in 1:nrow(allrich_data)) {
if (allrich_data$Treatment[i] == "N" | allrich_data$Treatment[i] =="NF") {
allrich_data$Nitrogen[i]<- 1
}
else {
allrich_data$Nitrogen[i]<- 0
}
if (allrich_data$Treatment[i] == "F" | allrich_data$Treatment[i] =="NF") {
allrich_data$Fungicide[i]<- 1
}
else {
allrich_data$Fungicide[i]<- 0
}
}
allrich_data$Nitrogen<-as.factor(allrich_data$Nitrogen)
allrich_data$Fungicide<-as.factor(allrich_data$Fungicide)
allrich_data<-subset(allrich_data, allrich_data$Richness != 18)
speciesmat<-allrich_data[,12:29]
#make an ordered factor out of the richness
allrich_data$network.size<-as.factor(as.numeric(allrich_data$Richness))
levels(allrich_data$network.size)<-c(5,7,11,15)
allrich_data$network.size<-ordered(allrich_data$network.size)
#multi membership initialisation
Species <- rep(LETTERS[1:18], length.out= nrow(allrich_data))
allrich_data$Species<-Species
#Diagonal dominance
lmod <- lFormula(Diagstrength~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m2 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
datadiagpred<-slarichpred(m2,allrich_data,n)
datadiagvpred <-slarichvpred(m2, allrich_data,n)
diagrichpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=datadiagpred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Diagstrength")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
diagrichvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=datadiagvpred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Diagstrength")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
#Asymmetry
lmod <- lFormula(Asymmetry~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m3 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
dataasympred<-slarichpred(m3,allrich_data,n)
dataasymvpred <-slarichvpred(m3, allrich_data,n)
asymrichpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=dataasympred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Asymmetry")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
asymrichvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=dataasymvpred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Asymmetry")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
#Skewness
lmod <- lFormula(Skewness~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m4 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
dataskewpred<-slarichpred(m4,allrich_data,n)
dataskewvpred <-slarichvpred(m4, allrich_data,n)
skewrichpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=dataskewpred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Skewness")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
skewrichvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=dataskewvpred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Skewness")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
#Evenness
lmod <- lFormula(1-Gini~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m5 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
dataginipred<-slarichpred(m5,allrich_data,n)
dataginivpred <-slarichvpred(m5, allrich_data,n)
ginirichpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=dataginipred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Gini")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
ginirichvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=dataginivpred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Gini")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
#Modularity
lmod <- lFormula(Modularity~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
datamodupred<-slarichpred(m1,allrich_data,n)
datamoduvpred <-slarichvpred(m1, allrich_data,n)
modurichpred<-ggplot(aes(x=SLA, y=fit, colour=Treatment), data=datamodupred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Modularity")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
modurichvpred<-ggplot(aes(x=SLAv, y=fit, colour=Treatment), data=datamoduvpred)+
facet_grid(~network.size)+
geom_smooth(method = lm)+
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1)+
ylab("Modularity")+
theme_bw()+
scale_color_manual(values = mycol)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
save.image(file = "results/multimembership_results_richness_5-7-11-15-18_species.Rdata")
From the results of these multimembership model, we can plot the
Supplementary figure 8, which present the impact of treatment and SLA on
mixed communities of different size across June and August 2020.
And the Supplementary figure 9, which present the impact of
treatment and SLA variance on mixed communities of different size across
June and August 2020.
Then, we can also look at the model output directly and plot the Supplementary table 3,
#scaling all variables
allrich_data$Modularity<-scale(allrich_data$Modularity, center = T)
allrich_data$Diagstrength<-scale(allrich_data$Diagstrength, center = T)
allrich_data$Asymmetry<-scale(allrich_data$Asymmetry, center = T)
allrich_data$Skewness<-scale(allrich_data$Skewness, center = T)
allrich_data$Gini<-scale(allrich_data$Gini, center = T)
#re running models with scaled variables to compare effect sizes
#Diagonal dominance
lmod <- lFormula(Diagstrength~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m2 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Asymmetry
lmod <- lFormula(Asymmetry~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m3 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Skewness
lmod <- lFormula(Skewness~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m4 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Evenness
lmod <- lFormula(1-Gini~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m5 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
#Modularity
lmod <- lFormula(Modularity~Nitrogen*Fungicide*network.size+Nitrogen*Fungicide*(SLA+SLAv)+(1|Species), data=allrich_data)
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(speciesmat))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
modelist<- list(m2,m3,m4,m5,m1)
names(modelist) <- c("Diagonal dominance", "Asymmetry", "Skewness", "Evenness", "Modularity")
modelsummary(modelist, statistic = NULL,
estimate = "{estimate}{stars} ({p.value})",
fmt = fmt_statistic(estimate = 2, p.value = 3),
coef_rename = c("SLAv" = "SLA variance", "network.size.L"="Network size (linear)", "Nitrogen1" = "Nitrogen", "Fungicide1" = "Fungicide"),
coef_omit = "(Intercept)|network.size.C|network.size.Q",
stars = T)| Diagonal dominance | Asymmetry | Skewness | Evenness | Modularity | |
|---|---|---|---|---|---|
| Nitrogen | −0.24*** (0.000) | 0.31*** (0.000) | 0.64*** (0.000) | 0.43*** (0.000) | 0.64*** (0.000) |
| Fungicide | −0.22*** (0.000) | 0.00 (0.999) | 1.00*** (0.000) | 1.31*** (0.000) | −0.65*** (0.000) |
| Network size (linear) | −2.69*** (0.000) | 0.00 (1.000) | −0.38 (0.514) | −0.25 (0.605) | −0.63* (0.024) |
| SLA | −0.02*** (0.000) | 0.17*** (0.000) | 0.12*** (0.000) | −0.22*** (0.000) | −0.15*** (0.000) |
| SLA variance | 0.06*** (0.000) | −0.04*** (0.000) | −0.22*** (0.000) | −0.02*** (0.000) | 0.19*** (0.000) |
| Nitrogen:Fungicide | 0.17*** (0.000) | −0.25*** (0.000) | −1.05*** (0.000) | −1.20*** (0.000) | −0.25*** (0.000) |
| Nitrogen:Network size (linear) | −0.25*** (0.000) | 0.00 (1.000) | 0.55*** (0.000) | 0.02 (0.230) | −0.53*** (0.000) |
| Fungicide:Network size (linear) | −0.20*** (0.000) | 0.00 (1.000) | 0.27*** (0.000) | 0.00 (0.870) | −0.17*** (0.000) |
| Nitrogen:SLA | −0.01*** (0.000) | −0.25*** (0.000) | −0.14*** (0.000) | 0.09*** (0.000) | 0.01*** (0.001) |
| Nitrogen:SLA variance | −0.05*** (0.000) | −0.13*** (0.000) | 0.25*** (0.000) | 0.10*** (0.000) | −0.23*** (0.000) |
| Fungicide:SLA | 0.02*** (0.000) | −0.07*** (0.000) | −0.08*** (0.000) | 0.00+ (0.060) | 0.25*** (0.000) |
| Fungicide:SLA variance | −0.04*** (0.000) | 0.03*** (0.000) | 0.34*** (0.000) | 0.08*** (0.000) | −0.30*** (0.000) |
| Nitrogen:Fungicide:Network size (linear) | 0.13*** (0.000) | 0.00 (1.000) | −0.70*** (0.000) | −0.12*** (0.000) | 0.37*** (0.000) |
| Nitrogen:Fungicide:SLA | 0.03*** (0.000) | 0.29*** (0.000) | −0.03*** (0.000) | 0.09*** (0.000) | 0.03*** (0.000) |
| Nitrogen:Fungicide:SLA variance | −0.03*** (0.000) | 0.19*** (0.000) | −0.26*** (0.000) | −0.08*** (0.000) | 0.11*** (0.000) |
| SD (Observations) | 0.24 | 0.85 | 0.63 | 0.66 | 0.73 |
| Num.Obs. | 584256 | 584256 | 584256 | 584256 | 584256 |
| R2 Marg. | 0.905 | 0.053 | 0.238 | 0.340 | 0.391 |
| R2 Cond. | 0.936 | 0.122 | 0.399 | 0.437 | 0.417 |
| AIC | −981.3 | 1465266.6 | 1120415.0 | 1167962.9 | 1295528.0 |
| BIC | −688.0 | 1465559.8 | 1120708.3 | 1168256.1 | 1295821.2 |
| ICC | 0.3 | 0.1 | 0.2 | 0.1 | 0.0 |
| RMSE | 0.24 | 0.85 | 0.63 | 0.66 | 0.73 |
b. Individual species effect
In order to interpret in details the results of our multimembership models, we ran some extra analysis to detect any species effects.
- insight on the treatment effect at the species level
We wanted to test whether treatment impacted competition coefficients at the species level. For this, we plotted Supplementary figure 4, the change in competitive effect (response is also presented here) with treatment for every species. The black line represents the values that would be obtained if there was no impact of the treatments.The dotted colored lines represent the obtained values by treatment (the difference was not significant).
#difference of competitive effect by species with treatment
bmatrix<-rowSums(alpha_june_control)
bmatrixn<-rowSums(alpha_june_nitrogen)
bmatrixf<-rowSums(alpha_june_fungicide)
bmatrixnf<-rowSums(alpha_june_combined)
bmatrixall<-data.frame(bmatrix,bmatrixn,bmatrixf, bmatrixnf)
matsp<-rownames(bmatrixall)
bmatrixall<-data.frame(bmatrix,bmatrixn,bmatrixf,bmatrixnf,matsp)
#bmatrixall[6,]<-NA #check if Holcus lanatus is driving the pattern
effect<-ggplot(data=bmatrixall,aes(x=bmatrix,y=bmatrixn))+
geom_smooth(aes(x=bmatrix,y=bmatrixn),method="lm",color="#7CAE00", alpha=0.3, fill="lightgray", linetype="dashed")+
geom_smooth(aes(x=bmatrix,y=bmatrixf),method="lm",color="#00BFC4", fill="lightgray", alpha=0.3, linetype="dashed")+
geom_smooth(aes(x=bmatrix,y=bmatrixnf),method="lm",color="#C77CFF", fill="lightgray", alpha=0.3, linetype="dashed")+
geom_point(color="#7CAE00")+
geom_text(aes(label=matsp),color="#7CAE00", hjust=1.5)+
geom_point(aes(x=bmatrix,y=bmatrixf), color="#00BFC4")+
geom_text(aes(x=bmatrix,y=bmatrixf,label=matsp), color="#00BFC4", hjust=1.5)+
geom_point(aes(x=bmatrix,y=bmatrixnf), color="#C77CFF")+
geom_text(aes(x=bmatrix,y=bmatrixnf,label=matsp), color="#C77CFF", hjust=1.5)+
geom_abline(intercept = 0, slope=1)+
xlim(c(-90,10))+ylim(c(-90,10))+
labs(title = "Competitive effect control vs. treatment June 2020",
y="Competitive coefficient with treatment", x="Competitive effect in control conditions")+
scale_fill_manual(values = c("#7CAE00","#00BFC4"), labels=c("Nitrogen","Fungicide"),drop=FALSE)+
theme_clean()+
theme(plot.background=element_blank())
amatrix<-colSums(alpha_june_control)
amatrixn<-colSums(alpha_june_nitrogen)
amatrixf<-colSums(alpha_june_fungicide)
amatrixnf<-colSums(alpha_june_combined)
amatrixall<-data.frame(amatrix,amatrixn,amatrixf,amatrixnf)
matsp<-rownames(amatrixall)
amatrixall<-data.frame(amatrix,amatrixn,amatrixf, amatrixnf, matsp)
response<-ggplot(data=amatrixall,aes(x=amatrix,y=amatrixn))+
geom_smooth(method="lm",aes(color="Nitrogen"), fill="lightgray", alpha=0.3, linetype="dashed")+
geom_smooth(aes(x=amatrix,y=amatrixf,color="Fungicide"),method="lm", fill="lightgray", alpha=0.3, linetype="dashed")+
geom_smooth(aes(x=amatrix,y=amatrixnf,color="Combined"),method="lm", fill="lightgray", alpha=0.3, linetype="dashed")+
geom_point(color="#7CAE00")+
geom_text(aes(label=matsp),color="#7CAE00", hjust=1.5)+
geom_point(aes(x=amatrix,y=amatrixf), color="#00BFC4")+
geom_text(aes(x=amatrix,y=amatrixf,label=matsp), color="#00BFC4", hjust=1.5)+
geom_point(aes(x=amatrix,y=amatrixnf), color="#C77CFF")+
geom_text(aes(x=amatrix,y=amatrixnf,label=matsp), color="#C77CFF", hjust=1.5)+
geom_abline(intercept = 0, slope=1)+
xlim(c(-90,10))+ylim(c(-90,10))+
labs(title = "Competitive response control vs. treatment June 2020",
y=" ", x="Competitive response in control conditions")+
scale_color_manual(name="Treatments",values = c(Nitrogen="#7CAE00",Fungicide="#00BFC4",Combined="#C77CFF"), breaks=c("Nitrogen","Fungicide", "Combined"))+
theme_clean()+
theme(plot.background=element_blank())
- comparison to the pairwise interactions level, effect of nitrogen vs. fungicide
Thirdly, we wanted to investigate whether the pairwise competitive coefficients change similarly across nitrogen and fungicide. For this, we calculated the absolute differences between our interaction coefficients in June control vs. nitrogen, and June control vs. fungicide. The 1.1 line represent what we would have expected if changes were correlated. We highlighted 4 pairwise interaction to illustrate our result in supplementary figure 5:
#absolute value between nitrogen and control in June
delta_contni_matrix<-abs(alpha_june_control-alpha_june_nitrogen)
plot_alpha_delta_contni<-melt.array(delta_contni_matrix)
#absolute value between fungicide and control in June
delta_contfun_matrix<-abs(alpha_june_control-alpha_june_fungicide)
plot_alpha_delta_contfun<-melt.array(delta_contfun_matrix)
#gathering the two in one data.frame
plot_alpha_delta_contnifun<-plot_alpha_delta_contni
plot_alpha_delta_contnifun$value2<-plot_alpha_delta_contfun$value
colnames(plot_alpha_delta_contnifun)<-c("sp1","sp2","control.vs.nitrogen","control.vs.fungicide")
#selecting 4 pairwise interaction to illustrate the result
selection<-subset(plot_alpha_delta_contnifun, plot_alpha_delta_contnifun$sp1=="Am"| plot_alpha_delta_contnifun$sp2=="Am")
selection<-selection[selection$sp2=="Hl"|selection$sp2=="Ga"|selection$sp2 == "Ao" | selection$sp1 == "Ga",]
rownames(selection)<-NULL
figs9a<-ggplot(plot_alpha_delta_contnifun, aes(x=control.vs.nitrogen,y=control.vs.fungicide))+
geom_point(alpha=0.2, shape=16)+
geom_point(data=selection, aes(x=control.vs.nitrogen,y=control.vs.fungicide),
color="red", shape=16, size=2)+
annotate("text", x=selection$control.vs.nitrogen[1], y=selection$control.vs.fungicide[1]+0.15,
label=c("Ga_Am"), color="red")+
annotate("text", x=selection$control.vs.nitrogen[2], y=selection$control.vs.fungicide[2]+0.15,
label=c("Am_Ao"), color="red")+
annotate("text", x=selection$control.vs.nitrogen[3], y=selection$control.vs.fungicide[3]+0.15,
label=c("Am_Hl"), color="red")+
annotate("text", x=selection$control.vs.nitrogen[4], y=selection$control.vs.fungicide[4]+0.15,
label=c("Am_Ga"), color="red")+
geom_abline(intercept = 0, slope = 1, color = "grey")+
theme_bw()
c. PaNDiv experiment whole field cover data
- insight on the fungicide effect
We used PaNDiv experiment cover data from 2020, to fit a linear model between the species average competitive response and the change in their cover with fungicide addition in June 2020.
#first, we need to do a little bit of data formatting from long to large cover data format
dat2020 <- read.table("data/202005_perco.txt", sep = "\t", header = T)
#changing names to same format and deleting useless columns
names(dat2020)[names(dat2020)=="Datum"]<-"datum"
names(dat2020)[names(dat2020)=="Plot_Nr"]<-"plot"
names(dat2020)[names(dat2020)=="Weed_species"]<-"species"
dat2020$Person<-NULL
dat2020$Block<-NULL
dat2020$original_species_in_non_VIP<-NULL
dat2020$Comments<-NULL
#simplifying dates (all are June cover)
dat2020$datum<-2020
#2020 table
datot20<-matrix(data=0, nrow=336, ncol=26, byrow = F, dimnames = NULL)
colnames(datot20)<-c("year","plot","Am","Ao","As","Be","Cb","Cj","Dc","Dg","Fr","Ga",
"Hl","Hp","Hs","Lp","Pg","Pm","Pt","Ra","Sp","To","bg",
"legume","herb","grass")
datot20<-as.data.frame(datot20)
dat2020<-as.data.frame(dat2020)
dat2020<-dat2020[order(dat2020$plot),]
#loop for filling 2020 cover table (large format)
i<-1
j<-1
k<-1
datot20$plot[i]<-i
datot20$year[i]<-2020
datot20<-tryCatch({
for (i in 1:nrow(datot20)){
for (k in 1:nrow(dat2020)) {
while (datot20$plot[i]==dat2020$plot[k]) {
for (j in 1:ncol(datot20)) {
if (colnames(datot20)[j]==dat2020$species[k]){
datot20[i,j]<-dat2020$X.weeded[k]
j<-j+1
}
else{
j<-j+1
}
}
k<-k+1
j<-1
}
}
i<-i+1
datot20$plot[i]<-i
datot20$year[i]<-2020
}
},
error = function (e) {
return(datot20)
})
datot20$sum<-rowSums(datot20[,3:26])
rel.cov20<-datot20[,3:26]
rel.cov20<-rel.cov20/datot20$sum
rel.cov20$plot<-datot20$plot
#adding treatment info to the plot
info<-read.table("data/plot_all.txt", sep = "\t", header = T)
rel.cov20$Nitrogen<-info$N
rel.cov20$Fungicide<-info$F
#creating dataset for testing the fungicide effect
test<-rel.cov20[1:24]
test[test==0]<-NA
test$Fungicide<-rel.cov20$Fungicide
test$Nitrogen<-rel.cov20$Nitrogen
# separating competitive effect from response
amatrix<-colSums(alpha_june_control)
bmatrix<-rowSums(alpha_june_control)
ctest<-subset(test, test$Fungicide==0)
cm<-colMeans(ctest, na.rm=T)
ftest<-subset(test, test$Fungicide==1)
fm<-colMeans(ftest, na.rm=T)
#creating the variable of the change in percentage cover with fungicide
diff<-fm-cm
#removing species that are not in the phytometers data set
diff$Hs<-NA
diff$As<-NA
diff<-diff[!is.na(diff)]
diff<-unlist(diff,use.names=T)
diff<-diff[1:18]
#order them in the same order than competition matrix
diff<-diff[order(names(diff))]
amatrix<-amatrix[order(names(amatrix))]
bmatrix<-bmatrix[order(names(bmatrix))]
#creating final dataframe for testing the linear model
x<-data.frame(diff=diff,amatrix=amatrix,bmatrix=bmatrix)
#x[10,]<-NA #remove the 2 extremes? still significant
#x[1,]<-NA
x<-na.omit(x)
model<-lm(data=x,diff~bmatrix)
summary(model)##
## Call:
## lm(formula = diff ~ bmatrix, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03337 -0.02122 -0.01010 0.01316 0.10188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0393661 0.0141631 2.779 0.0134 *
## bmatrix 0.0010785 0.0004811 2.242 0.0395 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03483 on 16 degrees of freedom
## Multiple R-squared: 0.239, Adjusted R-squared: 0.1915
## F-statistic: 5.025 on 1 and 16 DF, p-value: 0.03951We can then plot Supplementary figure 11,
Lastly, we used PaNDiv SLA data in order to check the plots shift in SLA with nitrogen treatment (monocultures were excluded). A negative Delta_SLA means that communities are shifting toward lower values of SLA, and that they become more slow growing.
- insight on the nitrogen effect and SLA shift in PaNDiv plots
SLA_data <- read.table("data/SLA shift.txt", header=T, sep="\t")
SLA_data$Delta_SLA<-as.numeric(SLA_data$Delta_SLA)
SLA_data <- na.omit(SLA_data)
compo<-read.table("data/info_plot.txt",sep="\t", header=T)
compo$Plot<-compo$Plot_Nr
compo$Plot_Nr<-NULL
compo$Block<-NULL
SLA_data<-merge(SLA_data,compo, by = "Plot")
#add the treatment in the right format (4-level factor)
SLA_data$Treatment<-NA
for (i in 1:nrow(SLA_data)) {
if (SLA_data$Nitrogen[i] == 0 & SLA_data$Fungicide[i] == 0) {
SLA_data$Treatment[i]<- "C"
}
if (SLA_data$Nitrogen[i] == 0 & SLA_data$Fungicide[i] == 1) {
SLA_data$Treatment[i]<- "F"
}
if (SLA_data$Nitrogen[i] == 1 & SLA_data$Fungicide[i] == 0) {
SLA_data$Treatment[i]<- "N"
}
if (SLA_data$Nitrogen[i] == 1 & SLA_data$Fungicide[i] == 1) {
SLA_data$Treatment[i]<- "NF"
}
}
SLA_data$Treatment<-as.factor(SLA_data$Treatment)
SLA_data$Treatment<-factor(SLA_data$Treatment, levels = c("C","N","F","NF"))
SLA_datac <- Rmisc::summarySE(SLA_data, measurevar="Delta_SLA", groupvars=c("Treatment","Year")) ##sd,seֵ
#some mixed modelling
#modify year for model
SLA_data$Year[SLA_data$Year==2016]<--3
SLA_data$Year[SLA_data$Year==2017]<--2
SLA_data$Year[SLA_data$Year==2018]<--1
SLA_data$Year[SLA_data$Year==2019]<-0
SLA_data$Year[SLA_data$Year==2020]<-1
SLA_data$Year[SLA_data$Year==2021]<-2
SLA_data$Year[SLA_data$Year==2022]<-3
SLA_data$composition<-as.factor(SLA_data$composition)
#mixed modelling
sla_shift_model<-lmer(Delta_SLA~Nitrogen*Fungicide*Year+(1|Plot/Block)+(1|Year)+(1|composition), data=SLA_data)
m<-modelsummary(list("(b)"=sla_shift_model),statistic = NULL,
estimate = "{estimate}{stars} ({p.value})",
fmt = fmt_statistic(estimate = 2, p.value = 3),
coef_omit = NULL,
stars = T)
m|
|
|
|---|---|
| (Intercept) | −0.23 (0.435) |
| Nitrogen | −0.52** (0.003) |
| Fungicide | −0.19 (0.281) |
| Year | 0.29* (0.011) |
| Nitrogen × Fungicide | 0.71** (0.005) |
| Nitrogen × Year | −0.10 (0.260) |
| Fungicide × Year | 0.13 (0.141) |
| Nitrogen × Fungicide × Year | −0.09 (0.462) |
| SD (Intercept Plot) | 0.00 |
| SD (Intercept BlockPlot) | 0.00 |
| SD (Intercept composition) | 1.32 |
| SD (Intercept Year) | 0.51 |
| SD (Observations) | 2.42 |
| Num.Obs. | 1494 |
| R2 Marg. | 0.064 |
| R2 Cond. | |
| AIC | 7051.0 |
| BIC | 7120.0 |
| RMSE | 2.37 |
#removing useless treatments, remove the following lines if you also want to look at fungicide and combined effect
#SLA_datac[SLA_datac$Treatment=="F",]<-NA
#SLA_datac<-na.omit(SLA_datac)
#SLA_datac[SLA_datac$Treatment=="NF",]<-NA
#SLA_datac<-na.omit(SLA_datac)
pd <- position_dodge(0.3)
a1 <- ggplot(SLA_datac, aes(x=Year, y=Delta_SLA, colour=Treatment, group=Treatment)) +
geom_errorbar(aes(ymin=Delta_SLA-se, ymax=Delta_SLA+se), colour="black", width=.1, position=pd)+
geom_point(position=pd, size=3)+
theme_bw()+
#scale_colour_manual(values = c("#F8766D", "#7CAE00"))+
ggtitle("(a)")+
scale_color_manual(values = mycol)+
theme(panel.grid = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(colour = "black") )