####################################################################################################################################### ################################################################# Header ############################################################## ####################################################################################################################################### # Purpose of this code: Re-run statistical analyses and produce main results (tables, figures) for Jochum et al. 2021, Oikos # Author of script: Malte Jochum (malte.jochum@idiv.de) # # Paper title: Earthworm invasion causes declines across soil fauna size classes and biodiversity facets in northern North American forests # First author: Malte Jochum # Journal: Oikos # Year of publication: 2021 # DOI: 10.1111/oik.07867 # # Last time performed: by Malte Jochum, 09.02.2021 # R version 3.6.3 (2020-02-29) -- "Holding the Windsock" # Copyright (C) 2020 The R Foundation for Statistical Computing # Platform: x86_64-pc-linux-gnu (64-bit) ### # What has been done in preparation of this analysis and is not covered in this script: # Loading and aggregating community-level data for earthworms, macro-, meso-, and microfauna (nematodes) ### # What this script does: # This script takes aggregated datasets and produces main figures and statistical output tables for the main analyses in the paper # If you would like to see details of other steps described in the manuscript (e.g. aggregation of data), please let me know and I will # happily provide further information (contact me at malte.jochum@idiv.de). ### # This script has been streamlined and annotated to the best of my knowledge. Please direct questions and reports of bugs or missing information to me via malte.jochum@idiv.de. ####################################################################################################################################### ########################################################### Table of contents ######################################################### ####################################################################################################################################### # This script contains the following sections: # I. Setup and head # II. Model the effect of earthworm invasion status on soil fauna communities of different size groups (macrofauna, mesofauna, nematodes) (Figure 1, Table 1) # III. Calculate and plot LNRR effect sizes with bootstrapped confidence intervals for relationships in section II (Figure 2) # IV. Model the effect of earthworm invasion intensity (earthworm fresh biomass) on soil fauna communities of different size groups (macrofauna, mesofauna, nematodes) (Figure 3, Table 2) # V. Run a model-selection procedure to test if biodiversity facet (community property) identity and size group affect the response of soil fauna to earthworm invasion (Supporting Information Table S6) # To navigate from one section to the next, best search this script for "I.", "II.", "III.", "IV.", or "V.". ####################################################################################################################################### ############### I. Setup and head ################### ####################################################################################################################################### rm(list=ls()) library(nlme) setwd(...) # set local working directory ####################################################################################################################################### ########## Load data sets plotdata_macro <- read.csv("./plotdata_macro_constrained_20210105.csv", sep = "," ) # Macrofauna community data and design variables (including earthworm invasion status and biomass) plotdata_meso <- read.csv("./plotdata_meso_constrained_20210105.csv", sep = "," ) # Mesofauna community data and design variables (including earthworm invasion status and biomass) plotdata_nema <- read.csv("./plotdata_nema_constrained_20210105.csv", sep = "," ) # Nematode (microfauna) community data and design variables (including earthworm invasion status and biomass) str(plotdata_macro) str(plotdata_meso) str(plotdata_nema) # For units, please refer to the README.txt file. For details on sampling, please refer to the methods section in Jochum et al. 2021 (Oikos). ####################################################################################################################################### ############### II. Model the effect of earthworm invasion status on soil fauna communities of different size groups (macrofauna, mesofauna, nematodes) (Figure 1, Table 1) ################### ####################################################################################################################################### ################################################################################################################################## # Figure 1: 12-panel figure: abundance, biomass, richness, shannon all in one; facets in columns, macro, meso, micro (nematodes) in rows, from top par(mar=c(1,3.75,1,1), oma=c(2,1.5,1,0.5)) layout(matrix(c(1,4,7,10, 2,5,8,11, 3,6,9,12), 3, 4, byrow = TRUE)) #################### abundance - all size classes ################################################## ###### #macro mod1.5_ma <- lme(log10(abundance.m2)~worms, random=~1|forest, data=plotdata_macro) summary(mod1.5_ma) boxplot(log10(plotdata_macro$abundance.m2[plotdata_macro$worms==0]),log10(plotdata_macro$abundance.m2[plotdata_macro$worms==1]), names=NA, col=c("#8C04CC","#5A2971"), cex.axis=1.2) mtext(bquote(log[10]~" abundance (ind. " ~ m^-2~ ")"), side=2, line=2.5, cex=1) text(1.5,2.45,"***", cex=2) mtext(expression(bold("Abundance")), side=3, line=0.4, cex=1) ###### # meso mod1.5_me <- lme(log10(abundance.m2)~worms, random=~1|forest, data=plotdata_meso) summary(mod1.5_me) boxplot(log10(plotdata_meso$abundance.m2[plotdata_meso$worms==0]),log10(plotdata_meso$abundance.m2[plotdata_meso$worms==1]), names=NA, col=c("#8C04CC","#5A2971"), cex.axis=1.2) mtext(bquote(log[10]~" abundance (ind. " ~ m^-2~ ")"), side=2, line=2.5, cex=1) text(1.5,5.8,"***", cex=2) ###### # micro mod1.5_mi <- lme(log10(abundance)~worms, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod1.5_mi) boxplot(log10(plotdata_nema$abundance[plotdata_nema$worms==0]),log10(plotdata_nema$abundance[plotdata_nema$worms==1]), col=c("#8C04CC","#5A2971"), cex.axis=1.2, names=c("","")) mtext(bquote(log[10]~" abundance (ind. " ~ m^-2~ ")"), side=2, line=2.5, cex=1) mtext(c("Low", "High Invasion"), at=c(1,2), side=1, line=1.5, cex=1) text(1.5,7,"***", cex=2) #################### biomass - all size classes ################################################## ###### #macro mod2.5_ma <- lme(log10(biomass.g.m2)~worms, random=~1|forest, data=plotdata_macro) summary(mod2.5_ma) boxplot(log10(plotdata_macro$biomass.g.m2[plotdata_macro$worms==0]),log10(plotdata_macro$biomass.g.m2[plotdata_macro$worms==1]), names=NA, col=c("#FF7400","#AA6C39"), cex.axis=1.2) mtext(bquote(log[10]~" fresh mass (g " ~ m^-2~ ")"), side=2, line=2.5, cex=1) text(1.5,0.73,"n.s.", cex=1.5) mtext(expression(bold("Biomass")), side=3, line=0.4, cex=1) ###### # meso mod2.5_me <- lme(log10(biomass.m2)~worms, random=~1|forest, data=plotdata_meso) summary(mod2.5_me) boxplot(log10(plotdata_meso$biomass.m2[plotdata_meso$worms==0]),log10(plotdata_meso$biomass.m2[plotdata_meso$worms==1]), names=NA, col=c("#FF7400","#AA6C39"), cex.axis=1.2) mtext(bquote(log[10]~" fresh mass (g " ~ m^-2~ ")"), side=2, line=2.5, cex=1) text(1.5,0.6,"**", cex=2) ###### # micro mod2.5_mi <- lme(log10(biomass.g.per.m2)~worms, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod2.5_mi) boxplot(log10(plotdata_nema$biomass.g.per.m2[plotdata_nema$worms==0]),log10(plotdata_nema$biomass.g.per.m2[plotdata_nema$worms==1]), col=c("#FF7400","#AA6C39"), cex.axis=1.2, names=c("","")) mtext(bquote(log[10]~" fresh mass (g " ~ m^-2~ ")"), side=2, line=2.5, cex=1) mtext(c("Low", "High Invasion"), at=c(1,2), side=1, line=1.5, cex=1) text(1.5,0.8,"***", cex=2) #################### richness - all size classes ################################################## ###### #macro mod3.5_ma <- lme(richness~worms, random=~1|forest, data=plotdata_macro) summary(mod3.5_ma) boxplot(plotdata_macro$richness[plotdata_macro$worms==0],plotdata_macro$richness[plotdata_macro$worms==1], names=NA, col=c("#00C2C2","#226666"), cex.axis=1.2) mtext(bquote("'species' richness"), side=2, line=2.5, cex=1) text(1.5,43,"***", cex=2) mtext(expression(bold("Richness")), side=3, line=0.4, cex=1) ###### # meso mod3.5_me <- lme(richness~worms, random=~1|forest, data=plotdata_meso) summary(mod3.5_me) boxplot(plotdata_meso$richness[plotdata_meso$worms==0],plotdata_meso$richness[plotdata_meso$worms==1], names=NA, col=c("#00C2C2","#226666"), cex.axis=1.2) mtext(bquote("'taxon' richness"), side=2, line=2.5, cex=1) text(1.5,17.5,"***", cex=2) ###### # micro mod3.5_mi <- lme(log10(richness)~worms, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod3.5_mi) boxplot(log10(plotdata_nema$richness[plotdata_nema$worms==0]),log10(plotdata_nema$richness[plotdata_nema$worms==1]), col=c("#00C2C2","#226666"), cex.axis=1.2, names=c("","")) mtext(bquote(log[10]~" genus richness"), side=2, line=2.5, cex=1) mtext(c("Low", "High Invasion"), at=c(1,2), side=1, line=1.5, cex=1) text(1.5,1.56,"***", cex=2) #################### shannon - all size classes ################################################## ###### #macro mod5.5_ma <- lme(shannon~worms, random=~1|forest, data=plotdata_macro) summary(mod5.5_ma) boxplot(plotdata_macro$shannon[plotdata_macro$worms==0],plotdata_macro$shannon[plotdata_macro$worms==1], names=NA, col=c("#9CEE00","#7A9E35"), cex.axis=1.2) mtext(bquote("'species' Shannon"), side=2, line=2.5, cex=1) text(1.5,3.4,"*", cex=2) mtext(expression(bold("Shannon")), side=3, line=0.4, cex=1) ###### # meso mod5.5_me <- lme(shannon~worms, random=~1|forest, data=plotdata_meso) summary(mod5.5_me) boxplot(plotdata_meso$shannon[plotdata_meso$worms==0],plotdata_meso$shannon[plotdata_meso$worms==1], names=NA, col=c("#9CEE00","#7A9E35"), cex.axis=1.2) mtext(bquote("'taxon' Shannon"), side=2, line=2.5, cex=1) text(1.5,2.5,"n.s", cex=1.5) ###### # micro mod5.5_mi <- lme(shannon~worms, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod5.5_mi) boxplot(plotdata_nema$shannon[plotdata_nema$worms==0],plotdata_nema$shannon[plotdata_nema$worms==1], col=c("#9CEE00","#7A9E35"), cex.axis=1.2, names=c("","")) mtext(bquote("genus Shannon"), side=2, line=2.5, cex=1) mtext(c("Low", "High Invasion"), at=c(1,2), side=1, line=1.5, cex=1) text(1.5,3.1,"***", cex=2) #dev.off() ################################################################################################################################## # Extracting model parameters (n, slope estimate, confint, p-value, R2) model_names <- c("macro_abundance", "macro_biomass", "macro_richness", "macro_shannon","meso_abundance", "meso_biomass", "meso_richness", "meso_shannon","micro_abundance", "micro_biomass", "micro_richness", "micro_shannon") #### #sample size model_ns <- c(summary(mod1.5_ma)$dims$N, summary(mod2.5_ma)$dims$N, summary(mod3.5_ma)$dims$N, summary(mod5.5_ma)$dims$N, summary(mod1.5_me)$dims$N, summary(mod2.5_me)$dims$N, summary(mod3.5_me)$dims$N, summary(mod5.5_me)$dims$N, summary(mod1.5_mi)$dims$N, summary(mod2.5_mi)$dims$N, summary(mod3.5_mi)$dims$N, summary(mod5.5_mi)$dims$N) #### #slopes model_slopes <- c(coef(mod1.5_ma)[1,2], coef(mod2.5_ma)[1,2], coef(mod3.5_ma)[1,2], coef(mod5.5_ma)[1,2], coef(mod1.5_me)[1,2], coef(mod2.5_me)[1,2], coef(mod3.5_me)[1,2], coef(mod5.5_me)[1,2], coef(mod1.5_mi)[1,2], coef(mod2.5_mi)[1,2], coef(mod3.5_mi)[1,2], coef(mod5.5_mi)[1,2]) #### #confidence intervals (using intervals() approximation as it is lme's) model_intlow <- c(intervals(mod1.5_ma)$fixed[2,1], intervals(mod2.5_ma)$fixed[2,1], intervals(mod3.5_ma)$fixed[2,1], intervals(mod5.5_ma)$fixed[2,1], intervals(mod1.5_me)$fixed[2,1], intervals(mod2.5_me,which=c("fixed"))$fixed[2,1], intervals(mod3.5_me,which=c("fixed"))$fixed[2,1], intervals(mod5.5_me)$fixed[2,1], # some weird issue otherwise... intervals(mod1.5_mi)$fixed[2,1], intervals(mod2.5_mi)$fixed[2,1], intervals(mod3.5_mi)$fixed[2,1], intervals(mod5.5_mi)$fixed[2,1]) model_intupp <- c(intervals(mod1.5_ma)$fixed[2,3], intervals(mod2.5_ma)$fixed[2,3], intervals(mod3.5_ma)$fixed[2,3], intervals(mod5.5_ma)$fixed[2,3], intervals(mod1.5_me)$fixed[2,3], intervals(mod2.5_me,which=c("fixed"))$fixed[2,3], intervals(mod3.5_me,which=c("fixed"))$fixed[2,3], intervals(mod5.5_me)$fixed[2,3], # some weird issue otherwise... intervals(mod1.5_mi)$fixed[2,3], intervals(mod2.5_mi)$fixed[2,3], intervals(mod3.5_mi)$fixed[2,3], intervals(mod5.5_mi)$fixed[2,3]) #### #F-values model_Fvals <- c(anova(mod1.5_ma)[2,3] , anova(mod2.5_ma)[2,3] , anova(mod3.5_ma)[2,3] , anova(mod5.5_ma)[2,3] , anova(mod1.5_me)[2,3] , anova(mod2.5_me)[2,3] , anova(mod3.5_me)[2,3] , anova(mod5.5_me)[2,3] , anova(mod1.5_mi)[2,3] , anova(mod2.5_mi)[2,3] , anova(mod3.5_mi)[2,3] , anova(mod5.5_mi)[2,3] ) #### #p-values model_pvals <- c(summary(mod1.5_ma)$tTable[2,5], summary(mod2.5_ma)$tTable[2,5], summary(mod3.5_ma)$tTable[2,5], summary(mod5.5_ma)$tTable[2,5], summary(mod1.5_me)$tTable[2,5], summary(mod2.5_me)$tTable[2,5], summary(mod3.5_me)$tTable[2,5], summary(mod5.5_me)$tTable[2,5], summary(mod1.5_mi)$tTable[2,5], summary(mod2.5_mi)$tTable[2,5], summary(mod3.5_mi)$tTable[2,5], summary(mod5.5_mi)$tTable[2,5]) #### library(piecewiseSEM) model_mar_r2 <- c(rsquared(mod1.5_ma)[1,5], rsquared(mod2.5_ma)[1,5], rsquared(mod3.5_ma)[1,5], rsquared(mod5.5_ma)[1,5], rsquared(mod1.5_me)[1,5], rsquared(mod2.5_me)[1,5], rsquared(mod3.5_me)[1,5], rsquared(mod5.5_me)[1,5], rsquared(mod1.5_mi)[1,5], rsquared(mod2.5_mi)[1,5], rsquared(mod3.5_mi)[1,5], rsquared(mod5.5_mi)[1,5]) model_con_r2 <- c(rsquared(mod1.5_ma)[1,6], rsquared(mod2.5_ma)[1,6], rsquared(mod3.5_ma)[1,6], rsquared(mod5.5_ma)[1,6], rsquared(mod1.5_me)[1,6], rsquared(mod2.5_me)[1,6], rsquared(mod3.5_me)[1,6], rsquared(mod5.5_me)[1,6], rsquared(mod1.5_mi)[1,6], rsquared(mod2.5_mi)[1,6], rsquared(mod3.5_mi)[1,6], rsquared(mod5.5_mi)[1,6]) model_transf <- c("log10", "log10", "none", "none", "log10", "log10", "none", "none", "log10", "log10", "log10", "none") #### put all together model_params <- data.frame(model=model_names, n=model_ns, resp_transf=model_transf, slope=round(model_slopes,2), low=round(model_intlow,2), upp=round(model_intupp,2), Fval=round(model_Fvals,2), pval=round(model_pvals,2), pval_holm=round(p.adjust(model_pvals, method="holm"),2), # Holm's adjustment to p-values when running multiple tests pseudoR2mar=round(model_mar_r2,2), pseudoR2con=round(model_con_r2,2)) ############### Table 1 model_params ####################################################################################################################################### ############### III. Calculate and plot LNRR effect sizes with bootstrapped confidence intervals for relationships in section II ################### ####################################################################################################################################### ################################################################################################################################## # Extracting raw data for calculating LNRR instead of percentage change #Abundance lnrr_macro_abund <- ddply(plotdata_macro, .(forest,worms), summarise, mean = mean(abundance.m2), sd = sd(abundance.m2), n = length(abundance.m2[is.na(abundance.m2)==FALSE])) lnrr_meso_abund <- ddply(plotdata_meso, .(forest,worms), summarise, mean = mean(abundance.m2), sd = sd(abundance.m2), n = length(abundance.m2[is.na(abundance.m2)==FALSE])) lnrr_micro_abund <- ddply(plotdata_nema, .(forest,worms), summarise, mean = mean(abundance, na.rm=TRUE), sd = sd(abundance, na.rm=TRUE), n = length(abundance[is.na(abundance)==FALSE])) lnrr_macro_biom <- ddply(plotdata_macro, .(forest,worms), summarise, mean = mean(biomass.g.m2), sd = sd(biomass.g.m2), n = length(biomass.g.m2[is.na(biomass.g.m2)==FALSE])) lnrr_meso_biom <- ddply(plotdata_meso, .(forest,worms), summarise, mean = mean(biomass.m2), sd = sd(biomass.m2), n = length(biomass.m2[is.na(biomass.m2)==FALSE])) lnrr_micro_biom <- ddply(plotdata_nema, .(forest,worms), summarise, mean = mean(biomass.g.per.m2, na.rm=TRUE), sd = sd(biomass.g.per.m2, na.rm=TRUE), n = length(biomass.g.per.m2[is.na(biomass.g.per.m2)==FALSE])) lnrr_macro_rich <- ddply(plotdata_macro, .(forest,worms), summarise, mean = mean(richness), sd = sd(richness), n = length(richness[is.na(richness)==FALSE])) lnrr_meso_rich <- ddply(plotdata_meso, .(forest,worms), summarise, mean = mean(richness), sd = sd(richness), n = length(richness[is.na(richness)==FALSE])) lnrr_micro_rich <- ddply(plotdata_nema, .(forest,worms), summarise, mean = mean(richness, na.rm=TRUE), sd = sd(richness, na.rm=TRUE), n = length(richness[is.na(richness)==FALSE])) lnrr_macro_shan <- ddply(plotdata_macro, .(forest,worms), summarise, mean = mean(shannon), sd = sd(shannon), n = length(shannon[is.na(shannon)==FALSE])) lnrr_meso_shan <- ddply(plotdata_meso, .(forest,worms), summarise, mean = mean(shannon), sd = sd(shannon), n = length(shannon[is.na(shannon)==FALSE])) lnrr_micro_shan <- ddply(plotdata_nema, .(forest,worms), summarise, mean = mean(shannon, na.rm=TRUE), sd = sd(shannon, na.rm=TRUE), n = length(shannon[is.na(shannon)==FALSE])) response <- c(rep("abundance",24),rep("biomass",24),rep("richness",24),rep("shannon",24)) # add column "response" sizeclass <- c(rep(c(rep("macro",8),rep("meso",8),rep("micro",8)),4)) # add column "response" lnrr_data <- rbind(lnrr_macro_abund, lnrr_meso_abund, lnrr_micro_abund, lnrr_macro_biom, lnrr_meso_biom, lnrr_micro_biom, lnrr_macro_rich, lnrr_meso_rich, lnrr_micro_rich, lnrr_macro_shan, lnrr_meso_shan, lnrr_micro_shan) lnrr_data <- cbind(response, sizeclass, lnrr_data) ################################################################################################################################## # Calculating LNRR based on Olga's code: # first, make some changes to the data structure: # concatenate response and size group lnrr_data$res_size <- paste(lnrr_data$response,lnrr_data$sizeclass, sep="_") lnrr_data_control <- subset(lnrr_data, lnrr_data$worms==0) names(lnrr_data_control)[5:7] <- c("mean_control", "sd_control", "N_control") lnrr_data_treat <- subset(lnrr_data, lnrr_data$worms==1) names(lnrr_data_treat)[5:7] <- c("mean_treat", "sd_treat", "N_treat") lnrr_data_fin <- cbind(lnrr_data_control, lnrr_data_treat[5:7]) library(metafor) #Build subset for each measure ab.macro<-subset(lnrr_data_fin,res_size=="abundance_macro") ab.meso<-subset(lnrr_data_fin,res_size=="abundance_meso") ab.micro<-subset(lnrr_data_fin,res_size=="abundance_micro") bi.macro<-subset(lnrr_data_fin,res_size=="biomass_macro") bi.meso<-subset(lnrr_data_fin,res_size=="biomass_meso") bi.micro<-subset(lnrr_data_fin,res_size=="biomass_micro") ri.macro<-subset(lnrr_data_fin,res_size=="richness_macro") ri.meso<-subset(lnrr_data_fin,res_size=="richness_meso") ri.micro<-subset(lnrr_data_fin,res_size=="richness_micro") sh.macro<-subset(lnrr_data_fin,res_size=="shannon_macro") sh.meso<-subset(lnrr_data_fin,res_size=="shannon_meso") sh.micro<-subset(lnrr_data_fin,res_size=="shannon_micro") #Calculate effect sizes, ROM= Log response ratio es1<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=ab.macro) es2<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=ab.meso) es3<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=ab.micro) es4<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=bi.macro) es5<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=bi.meso) es6<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=bi.micro) es7<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=ri.macro) es8<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=ri.meso) es9<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=ri.micro) es10<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=sh.macro) es11<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=sh.meso) es12<-escalc(measure="ROM", m1i=mean_treat, m2i=mean_control, sd1i=sd_treat, sd2i=sd_control, n1i=N_treat, n2i=N_control, data=sh.micro) ##Conduct meta-regression with forest as random factor meta_reg.ab.macro<-rma.mv(yi, vi, random=~1|forest, data=es1) meta_reg.ab.meso<-rma.mv(yi, vi, random=~1|forest, data=es2) meta_reg.ab.micro<-rma.mv(yi, vi, random=~1|forest, data=es3) meta_reg.bi.macro<-rma.mv(yi, vi, random=~1|forest, data=es4) meta_reg.bi.meso<-rma.mv(yi, vi, random=~1|forest, data=es5) meta_reg.bi.micro<-rma.mv(yi, vi, random=~1|forest, data=es6) meta_reg.ri.macro<-rma.mv(yi, vi, random=~1|forest, data=es7) meta_reg.ri.meso<-rma.mv(yi, vi, random=~1|forest, data=es8) meta_reg.ri.micro<-rma.mv(yi, vi, random=~1|forest, data=es9) meta_reg.sh.macro<-rma.mv(yi, vi, random=~1|forest, data=es10) meta_reg.sh.meso<-rma.mv(yi, vi, random=~1|forest, data=es11) meta_reg.sh.micro<-rma.mv(yi, vi, random=~1|forest, data=es12) meta_reg.ab.macro meta_reg.ab.meso meta_reg.ab.micro meta_reg.bi.macro meta_reg.bi.meso meta_reg.bi.micro meta_reg.ri.macro meta_reg.ri.meso meta_reg.ri.micro meta_reg.sh.macro meta_reg.sh.meso meta_reg.sh.micro ########## Extracting parameters ma_a_lnrr <- meta_reg.ab.macro$b ma_a_lnrr_ci_low <- meta_reg.ab.macro$ci.lb ma_a_lnrr_ci_upp <- meta_reg.ab.macro$ci.ub ma_a_lnrr_pval <- meta_reg.ab.macro$pval me_a_lnrr <- meta_reg.ab.meso$b me_a_lnrr_ci_low <- meta_reg.ab.meso$ci.lb me_a_lnrr_ci_upp <- meta_reg.ab.meso$ci.ub me_a_lnrr_pval <- meta_reg.ab.meso$pval mi_a_lnrr <- meta_reg.ab.micro$b mi_a_lnrr_ci_low <- meta_reg.ab.micro$ci.lb mi_a_lnrr_ci_upp <- meta_reg.ab.micro$ci.ub mi_a_lnrr_pval <- meta_reg.ab.micro$pval ma_b_lnrr <- meta_reg.bi.macro$b ma_b_lnrr_ci_low <- meta_reg.bi.macro$ci.lb ma_b_lnrr_ci_upp <- meta_reg.bi.macro$ci.ub ma_b_lnrr_pval <- meta_reg.bi.macro$pval me_b_lnrr <- meta_reg.bi.meso$b me_b_lnrr_ci_low <- meta_reg.bi.meso$ci.lb me_b_lnrr_ci_upp <- meta_reg.bi.meso$ci.ub me_b_lnrr_pval <- meta_reg.bi.meso$pval mi_b_lnrr <- meta_reg.bi.micro$b mi_b_lnrr_ci_low <- meta_reg.bi.micro$ci.lb mi_b_lnrr_ci_upp <- meta_reg.bi.micro$ci.ub mi_b_lnrr_pval <- meta_reg.bi.micro$pval ma_r_lnrr <- meta_reg.ri.macro$b ma_r_lnrr_ci_low <- meta_reg.ri.macro$ci.lb ma_r_lnrr_ci_upp <- meta_reg.ri.macro$ci.ub ma_r_lnrr_pval <- meta_reg.ri.macro$pval me_r_lnrr <- meta_reg.ri.meso$b me_r_lnrr_ci_low <- meta_reg.ri.meso$ci.lb me_r_lnrr_ci_upp <- meta_reg.ri.meso$ci.ub me_r_lnrr_pval <- meta_reg.ri.meso$pval mi_r_lnrr <- meta_reg.ri.micro$b mi_r_lnrr_ci_low <- meta_reg.ri.micro$ci.lb mi_r_lnrr_ci_upp <- meta_reg.ri.micro$ci.ub mi_r_lnrr_pval <- meta_reg.ri.micro$pval ma_s_lnrr <- meta_reg.sh.macro$b ma_s_lnrr_ci_low <- meta_reg.sh.macro$ci.lb ma_s_lnrr_ci_upp <- meta_reg.sh.macro$ci.ub ma_s_lnrr_pval <- meta_reg.sh.macro$pval me_s_lnrr <- meta_reg.sh.meso$b me_s_lnrr_ci_low <- meta_reg.sh.meso$ci.lb me_s_lnrr_ci_upp <- meta_reg.sh.meso$ci.ub me_s_lnrr_pval <- meta_reg.sh.meso$pval mi_s_lnrr <- meta_reg.sh.micro$b mi_s_lnrr_ci_low <- meta_reg.sh.micro$ci.lb mi_s_lnrr_ci_upp <- meta_reg.sh.micro$ci.ub mi_s_lnrr_pval <- meta_reg.sh.micro$pval ########## Plotting LNRR results ####################################################### ###### Bootstrap the confidence intervals: #Please note that bootstrapped values differ between individual runs and may thus deviate slightly from what is shown in Fig. 2 of the paper. ##Non-parametric bootstrapping (here, only for macrofauna abundance) library(boot) #macrofauna abundance boot.func.ab.macro<-function(es1, indices) { res.ab.macro<- try(rma(yi, vi, data=es1, subset=indices), silent=TRUE) if (is.element("try-error", class(res.ab.macro))) { c(NA,NA,NA,NA) } else { c(coef(res.ab.macro), vcov(res.ab.macro), res.ab.macro$tau2, res.ab.macro$se.tau2^2) } } res.boot.ab.macro <- boot(es1, boot.func.ab.macro, R=1000) ab_macro_boot_ci <- boot.ci(res.boot.ab.macro, index=1:2) ## Use Bias Corrected Bootstrap CIs ma_a_lnrr_ci_low2 <- ab_macro_boot_ci$basic[4] ma_a_lnrr_ci_upp2 <- ab_macro_boot_ci$basic[5] #mesofauna abundance boot.func.ab.meso<-function(es2, indices) { res.ab.meso<- try(rma(yi, vi, data=es2, subset=indices), silent=TRUE) if (is.element("try-error", class(res.ab.meso))) { c(NA,NA,NA,NA) } else { c(coef(res.ab.meso), vcov(res.ab.meso), res.ab.meso$tau2, res.ab.meso$se.tau2^2) } } res.boot.ab.meso <- boot(es2, boot.func.ab.meso, R=1000) ab_meso_boot_ci <- boot.ci(res.boot.ab.meso, index=1:2) ## Use Bias Corrected Bootstrap CIs me_a_lnrr_ci_low2 <- ab_meso_boot_ci$basic[4] me_a_lnrr_ci_upp2 <- ab_meso_boot_ci$basic[5] #microfauna abundance boot.func.ab.micro<-function(es3, indices) { res.ab.micro<- try(rma(yi, vi, data=es3, subset=indices), silent=TRUE) if (is.element("try-error", class(res.ab.micro))) { c(NA,NA,NA,NA) } else { c(coef(res.ab.micro), vcov(res.ab.micro), res.ab.micro$tau2, res.ab.micro$se.tau2^2) } } res.boot.ab.micro <- boot(es3, boot.func.ab.micro, R=1000) ab_micro_boot_ci <- boot.ci(res.boot.ab.micro, index=1:2) ## Use Bias Corrected Bootstrap CIs mi_a_lnrr_ci_low2 <- ab_micro_boot_ci$basic[4] mi_a_lnrr_ci_upp2 <- ab_micro_boot_ci$basic[5] ######### #macrofauna biomass boot.func.bi.macro<-function(es4, indices) { res.bi.macro<- try(rma(yi, vi, data=es4, subset=indices), silent=TRUE) if (is.element("try-error", class(res.bi.macro))) { c(NA,NA,NA,NA) } else { c(coef(res.bi.macro), vcov(res.bi.macro), res.bi.macro$tau2, res.bi.macro$se.tau2^2) } } res.boot.bi.macro <- boot(es4, boot.func.bi.macro, R=1000) bi_macro_boot_ci <- boot.ci(res.boot.bi.macro, index=1:2) ## Use Bias Corrected Bootstrap CIs ma_b_lnrr_ci_low2 <- bi_macro_boot_ci$basic[4] ma_b_lnrr_ci_upp2 <- bi_macro_boot_ci$basic[5] #mesofauna biomass boot.func.bi.meso<-function(es5, indices) { res.bi.meso<- try(rma(yi, vi, data=es5, subset=indices), silent=TRUE) if (is.element("try-error", class(res.bi.meso))) { c(NA,NA,NA,NA) } else { c(coef(res.bi.meso), vcov(res.bi.meso), res.bi.meso$tau2, res.bi.meso$se.tau2^2) } } res.boot.bi.meso <- boot(es5, boot.func.bi.meso, R=1000) bi_meso_boot_ci <- boot.ci(res.boot.bi.meso, index=1:2) ## Use Bias Corrected Bootstrap CIs me_b_lnrr_ci_low2 <- bi_meso_boot_ci$basic[4] me_b_lnrr_ci_upp2 <- bi_meso_boot_ci$basic[5] #microfauna biomass boot.func.bi.micro<-function(es6, indices) { res.bi.micro<- try(rma(yi, vi, data=es6, subset=indices), silent=TRUE) if (is.element("try-error", class(res.bi.micro))) { c(NA,NA,NA,NA) } else { c(coef(res.bi.micro), vcov(res.bi.micro), res.bi.micro$tau2, res.bi.micro$se.tau2^2) } } res.boot.bi.micro <- boot(es6, boot.func.bi.micro, R=1000) bi_micro_boot_ci <- boot.ci(res.boot.bi.micro, index=1:2) ## Use Bias Corrected Bootstrap CIs mi_b_lnrr_ci_low2 <- bi_micro_boot_ci$basic[4] mi_b_lnrr_ci_upp2 <- bi_micro_boot_ci$basic[5] ######### #macrofauna richness boot.func.ri.macro<-function(es7, indices) { res.ri.macro<- try(rma(yi, vi, data=es7, subset=indices), silent=TRUE) if (is.element("try-error", class(res.ri.macro))) { c(NA,NA,NA,NA) } else { c(coef(res.ri.macro), vcov(res.ri.macro), res.ri.macro$tau2, res.ri.macro$se.tau2^2) } } res.boot.ri.macro <- boot(es7, boot.func.ri.macro, R=1000) ri_macro_boot_ci <- boot.ci(res.boot.ri.macro, index=1:2) ## Use Bias Corrected Bootstrap CIs ma_r_lnrr_ci_low2 <- ri_macro_boot_ci$basic[4] ma_r_lnrr_ci_upp2 <- ri_macro_boot_ci$basic[5] #mesofauna richness boot.func.ri.meso<-function(es8, indices) { res.ri.meso<- try(rma(yi, vi, data=es8, subset=indices), silent=TRUE) if (is.element("try-error", class(res.ri.meso))) { c(NA,NA,NA,NA) } else { c(coef(res.ri.meso), vcov(res.ri.meso), res.ri.meso$tau2, res.ri.meso$se.tau2^2) } } res.boot.ri.meso <- boot(es8, boot.func.ri.meso, R=1000) ri_meso_boot_ci <- boot.ci(res.boot.ri.meso, index=1:2) ## Use Bias Corrected Bootstrap CIs me_r_lnrr_ci_low2 <- ri_meso_boot_ci$basic[4] me_r_lnrr_ci_upp2 <- ri_meso_boot_ci$basic[5] #microfauna richness boot.func.ri.micro<-function(es9, indices) { res.ri.micro<- try(rma(yi, vi, data=es9, subset=indices), silent=TRUE) if (is.element("try-error", class(res.ri.micro))) { c(NA,NA,NA,NA) } else { c(coef(res.ri.micro), vcov(res.ri.micro), res.ri.micro$tau2, res.ri.micro$se.tau2^2) } } res.boot.ri.micro <- boot(es9, boot.func.ri.micro, R=1000) ri_micro_boot_ci <- boot.ci(res.boot.ri.micro, index=1:2) ## Use Bias Corrected Bootstrap CIs mi_r_lnrr_ci_low2 <- ri_micro_boot_ci$basic[4] mi_r_lnrr_ci_upp2 <- ri_micro_boot_ci$basic[5] ######### #macrofauna shannon boot.func.sh.macro<-function(es10, indices) { res.sh.macro<- try(rma(yi, vi, data=es10, subset=indices), silent=TRUE) if (is.element("try-error", class(res.sh.macro))) { c(NA,NA,NA,NA) } else { c(coef(res.sh.macro), vcov(res.sh.macro), res.sh.macro$tau2, res.sh.macro$se.tau2^2) } } res.boot.sh.macro <- boot(es10, boot.func.sh.macro, R=1000) sh_macro_boot_ci <- boot.ci(res.boot.sh.macro, index=1:2) ## Use Bias Corrected Bootstrap CIs ma_s_lnrr_ci_low2 <- sh_macro_boot_ci$basic[4] ma_s_lnrr_ci_upp2 <- sh_macro_boot_ci$basic[5] #mesofauna shannon boot.func.sh.meso<-function(es11, indices) { res.sh.meso<- try(rma(yi, vi, data=es11, subset=indices), silent=TRUE) if (is.element("try-error", class(res.sh.meso))) { c(NA,NA,NA,NA) } else { c(coef(res.sh.meso), vcov(res.sh.meso), res.sh.meso$tau2, res.sh.meso$se.tau2^2) } } res.boot.sh.meso <- boot(es11, boot.func.sh.meso, R=1000) sh_meso_boot_ci <- boot.ci(res.boot.sh.meso, index=1:2) ## Use Bias Corrected Bootstrap CIs me_s_lnrr_ci_low2 <- sh_meso_boot_ci$basic[4] me_s_lnrr_ci_upp2 <- sh_meso_boot_ci$basic[5] #microfauna shannon boot.func.sh.micro<-function(es12, indices) { res.sh.micro<- try(rma(yi, vi, data=es12, subset=indices), silent=TRUE) if (is.element("try-error", class(res.sh.micro))) { c(NA,NA,NA,NA) } else { c(coef(res.sh.micro), vcov(res.sh.micro), res.sh.micro$tau2, res.sh.micro$se.tau2^2) } } res.boot.sh.micro <- boot(es12, boot.func.sh.micro, R=1000) sh_micro_boot_ci <- boot.ci(res.boot.sh.micro, index=1:2) ## Use Bias Corrected Bootstrap CIs mi_s_lnrr_ci_low2 <- sh_micro_boot_ci$basic[4] mi_s_lnrr_ci_upp2 <- sh_micro_boot_ci$basic[5] #################### plot LNRR with bootstrapped CI!! x <- c(15,20) y <- c(1,2) plot(y,x, type="n", ylim=c(-2,1), xlim=c(0.5,12.5), xaxt="n", ylab="", cex.axis=0.8, # yaxt="n", # ylab="slope", xlab="", frame=F ) mtext("LNRR +/- bootstrapped CI", side=2, line=2, cex=0.8) # polygons for better readability (maybe) polygon(c(0.05,0.05,14,14),c(0,0.1,0.1,0), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(0.2,0.3,0.3,0.2), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(0.40,0.50,0.50,0.40), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(0.60,0.70,0.70,0.60), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(0.80,0.90,0.90,0.80), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(0,-0.10,-0.10,0), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-0.20,-0.30,-0.30,-0.20), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-0.40,-0.50,-0.50,-0.40), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-0.60,-0.70,-0.70,-0.60), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-0.80,-0.90,-0.90,-0.80), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-1,-1.10,-1.10,-1), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-1.20,-1.30,-1.30,-1.20), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-1.40,-1.50,-1.50,-1.40), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-1.60,-1.70,-1.70,-1.60), col= "gray90", border=FALSE) polygon(c(0.05,0.05,14,14),c(-1.80,-1.90,-1.90,-1.80), col= "gray90", border=FALSE) abline(h=0, lty=2) #abline(v=4.5, lty=1) # vertical lines to separate the fauna groups #abline(v=8.5, lty=1) points(1,ma_a_lnrr, col = "#6D1099", pch=19, cex=1) points(2,ma_b_lnrr, col = "#E5700E", pch=19, cex=1) points(3,ma_r_lnrr, col = "#098A8A", pch=19, cex=1) points(4,ma_s_lnrr, col = "#8BC817", pch=19, cex=1) points(5,me_a_lnrr, col = "#6D1099", pch=19, cex=1) points(6,me_b_lnrr, col = "#E5700E", pch=19, cex=1) points(7,me_r_lnrr, col = "#098A8A", pch=19, cex=1) points(8,me_s_lnrr, col = "#8BC817", pch=19, cex=1) points(9,mi_a_lnrr, col = "#6D1099", pch=19, cex=1) points(10,mi_b_lnrr, col = "#E5700E", pch=19, cex=1) points(11,mi_r_lnrr, col = "#098A8A", pch=19, cex=1) points(12,mi_s_lnrr, col = "#8BC817", pch=19, cex=1) arrows(1, ma_a_lnrr_ci_low2, 1, ma_a_lnrr_ci_upp2, length=0, col="#6D1099", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(2, ma_b_lnrr_ci_low2, 2, ma_b_lnrr_ci_upp2, length=0, col="#E5700E", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(3, ma_r_lnrr_ci_low2, 3, ma_r_lnrr_ci_upp2, length=0, col="#098A8A", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(4, ma_s_lnrr_ci_low2, 4, ma_s_lnrr_ci_upp2, length=0, col="#8BC817", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(5, me_a_lnrr_ci_low2, 5, me_a_lnrr_ci_upp2, length=0, col="#6D1099", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(6, me_b_lnrr_ci_low2, 6, me_b_lnrr_ci_upp2, length=0, col="#E5700E", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(7, me_r_lnrr_ci_low2, 7, me_r_lnrr_ci_upp2, length=0, col="#098A8A", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(8, me_s_lnrr_ci_low2, 8, me_s_lnrr_ci_upp2, length=0, col="#8BC817", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(9, mi_a_lnrr_ci_low2, 9, mi_a_lnrr_ci_upp2, length=0, col="#6D1099", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(10, mi_b_lnrr_ci_low2, 10, mi_b_lnrr_ci_upp2, length=0, col="#E5700E", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(11, mi_r_lnrr_ci_low2, 11, mi_r_lnrr_ci_upp2, length=0, col="#098A8A", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, arrows(12, mi_s_lnrr_ci_low2, 12, mi_s_lnrr_ci_upp2, length=0, col="#8BC817", lwd=2) # from x, from y, to x, to y; length 0 says no head! otherwise do: length=0.025, angle=90, code=3, text(seq(1,12,by=1), par("usr")[3]-0.25, srt = 45, adj= 1, xpd = TRUE, labels = c("Abundance","Biomass","Richness","Shannon", "Abundance","Biomass","Richness","Shannon", "Abundance","Biomass","Richness","Shannon"), cex=0.8) #dev.off() ############## # save these percentage change values for the manuscript text: name <- c(rep("macro",4),rep("meso",4),rep("micro",4)) name2 <- rep(c("abund", "biomass", "richness", "shannon"),3) lnrr_means <- c(ma_a_lnrr, ma_b_lnrr, ma_r_lnrr, ma_s_lnrr,me_a_lnrr, me_b_lnrr, me_r_lnrr, me_s_lnrr,mi_a_lnrr, mi_b_lnrr, mi_r_lnrr, mi_s_lnrr) lnrr_lower <- c(ma_a_lnrr_ci_low2, ma_b_lnrr_ci_low2, ma_r_lnrr_ci_low2, ma_s_lnrr_ci_low2,me_a_lnrr_ci_low2, me_b_lnrr_ci_low2, me_r_lnrr_ci_low2, me_s_lnrr_ci_low2,mi_a_lnrr_ci_low2, mi_b_lnrr_ci_low2, mi_r_lnrr_ci_low2, mi_s_lnrr_ci_low2) lnrr_upper <- c(ma_a_lnrr_ci_upp2, ma_b_lnrr_ci_upp2, ma_r_lnrr_ci_upp2, ma_s_lnrr_ci_upp2,me_a_lnrr_ci_upp2, me_b_lnrr_ci_upp2, me_r_lnrr_ci_upp2, me_s_lnrr_ci_upp2,mi_a_lnrr_ci_upp2, mi_b_lnrr_ci_upp2, mi_r_lnrr_ci_upp2, mi_s_lnrr_ci_upp2) lnrr_summary <- data.frame(sizegroup=name, facet=name2, lnrr=round(lnrr_means,2), lower=round(lnrr_lower,2), upper=round(lnrr_upper,2)) ####################################################################################################################################### ############### IV. Model the effect of earthworm invasion intensity (earthworm fresh biomass) on soil fauna communities of different size groups (macrofauna, mesofauna, nematodes) (Figure 3, Table 2) ################### ####################################################################################################################################### ### add worm biomass data to plotdata_nema plotdata_nema$worm_biomass.m2 <- plotdata_macro$worm_biomass.m2 ####################################### par(mar=c(1.5,3.75,1,1), oma=c(3,1.5,1,0.5)) layout(matrix(c(1,4,7,10, 2,5,8,11, 3,6,9,12), 3, 4, byrow = TRUE)) # the ### ABUNDANCE ### Macro # make my variables full-on columns in the data frame (no log transformation in model call) plotdata_macro$log_worm_bmass_plusone <- log10(plotdata_macro$worm_biomass.m2+1) log_worm_bmass_plusone <- log10(plotdata_macro$worm_biomass.m2+1) plotdata_macro$log_abundance.m2 <- log10(plotdata_macro$abundance.m2) mod1.5_ma_w <- lme(log_abundance.m2~log_worm_bmass_plusone, random=~1|forest, data=plotdata_macro) summary(mod1.5_ma_w) plot(plotdata_macro$log_worm_bmass_plusone,plotdata_macro$log_abundance.m2, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(1.2,2.6), pch=16, frame=F ) mtext(bquote(log[10]~" abundance (ind. " ~ m^-2~ ")"), side=2, line=2.5, cex=1) points(log10(plotdata_macro$worm_biomass.m2[plotdata_macro$worms==0]+1),log10(plotdata_macro$abundance.m2[plotdata_macro$worms==0]), col="#8C04CC", pch=21) points(log10(plotdata_macro$worm_biomass.m2[plotdata_macro$worms==1]+1),log10(plotdata_macro$abundance.m2[plotdata_macro$worms==1]), col="#5A2971", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod1.5_ma_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod1.5_ma_w)$tTable[2,5]<=0.05){lines(x,y, col="#8C04CC", lwd=2.5)} if(summary(mod1.5_ma_w)$tTable[2,5]>=0.05){lines(x,y, col="#8C04CC", lwd=2.5, lty=2)} legend("topright", legend = c("low invasion","high invasion", "n.s."), bty = "n", col = "grey", lty=c(NA,NA,2), pch=c(21,16,NA), lwd=2, cex=1) mtext(expression(bold("Abundance")), side=3, line=0.4, cex=1) ### Meso plotdata_meso$log_worm_bmass_plusone <- log10(plotdata_meso$worm_biomass.m2+1) log_worm_bmass_plusone <- log10(plotdata_meso$worm_biomass.m2+1) plotdata_meso$log_abundance.m2 <- log10(plotdata_meso$abundance.m2) mod1.5_me_w <- lme(log_abundance.m2~log_worm_bmass_plusone, random=~1|forest, data=plotdata_meso) summary(mod1.5_me_w) plot(plotdata_meso$log_worm_bmass_plusone,plotdata_meso$log_abundance.m2, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(3.5,6), pch=16, frame=F ) mtext(bquote(log[10]~" abundance (ind. " ~ m^-2~ ")"), side=2, line=2.5, cex=1) points(log10(plotdata_meso$worm_biomass.m2[plotdata_meso$worms==0]+1),log10(plotdata_meso$abundance.m2[plotdata_meso$worms==0]), col="#8C04CC", pch=21) points(log10(plotdata_meso$worm_biomass.m2[plotdata_meso$worms==1]+1),log10(plotdata_meso$abundance.m2[plotdata_meso$worms==1]), col="#5A2971", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod1.5_me_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod1.5_me_w)$tTable[2,5]<=0.05){lines(x,y, col="#8C04CC", lwd=2.5)} if(summary(mod1.5_me_w)$tTable[2,5]>=0.05){lines(x,y, col="#8C04CC", lwd=2.5, lty=2)} ### Nema ## plotdata_nema$log_worm_bmass_plusone <- log10(plotdata_nema$worm_biomass.m2+1) log_worm_bmass_plusone <- log10(plotdata_nema$worm_biomass.m2+1) plotdata_nema$log_abundance.m2 <- log10(plotdata_nema$abundance) mod1.5_mi_w <- lme(log_abundance.m2~log_worm_bmass_plusone, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod1.5_mi_w) plot(plotdata_nema$log_worm_bmass_plusone,plotdata_nema$log_abundance.m2, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(5.6,7.2), pch=16, frame=F ) mtext(bquote(log[10]~" abundance (ind. " ~ m^-2~ ")"), side=2, line=2.5, cex=1) mtext(bquote(log[10]~"earthworm biomass (g"~ m^-2~ "+1)"), side=1, line=3.5, cex=1) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==0],plotdata_nema$log_abundance.m2[plotdata_nema$worms==0], col="#8C04CC", pch=21) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==1],plotdata_nema$log_abundance.m2[plotdata_nema$worms==1], col="#5A2971", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod1.5_mi_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod1.5_mi_w)$tTable[2,5]<=0.05){lines(x,y, col="#8C04CC", lwd=2.5)} if(summary(mod1.5_mi_w)$tTable[2,5]>=0.05){lines(x,y, col="#8C04CC", lwd=2.5, lty=2)} ### BIOMASS ### Macro plotdata_macro$log_biomass.g.m2 <- log10(plotdata_macro$biomass.g.m2) mod2.5_ma_w <- lme(log_biomass.g.m2~log_worm_bmass_plusone, random=~1|forest, data=plotdata_macro) summary(mod2.5_ma_w) plot(plotdata_macro$log_worm_bmass_plusone,plotdata_macro$log_biomass.g.m2, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(-1,1), pch=16, frame=F ) mtext(bquote(log[10]~" fresh mass (g " ~ m^-2~ ")"), side=2, line=2.5, cex=1) points(plotdata_macro$log_worm_bmass_plusone[plotdata_macro$worms==0],plotdata_macro$log_biomass.g.m2[plotdata_macro$worms==0], col="#FF7400", pch=21) points(plotdata_macro$log_worm_bmass_plusone[plotdata_macro$worms==1],plotdata_macro$log_biomass.g.m2[plotdata_macro$worms==1], col="#AA6C39", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod2.5_ma_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod2.5_ma_w)$tTable[2,5]<=0.05){lines(x,y, col="#FF7400", lwd=2.5)} if(summary(mod2.5_ma_w)$tTable[2,5]>=0.05){lines(x,y, col="#FF7400", lwd=2.5, lty=2)} mtext(expression(bold("Biomass")), side=3, line=0.4, cex=1) ### Meso plotdata_meso$log_biomass.m2 <- log10(plotdata_meso$biomass.m2) mod2.5_me_w <- lme(log_biomass.m2~log_worm_bmass_plusone, random=~1|forest, data=plotdata_meso) summary(mod2.5_me_w) plot(plotdata_meso$log_worm_bmass_plusone,plotdata_meso$log_biomass.m2, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(-1.5,1), pch=16, frame=F ) mtext(bquote(log[10]~" fresh mass (g " ~ m^-2~ ")"), side=2, line=2.5, cex=1) points(plotdata_meso$log_worm_bmass_plusone[plotdata_meso$worms==0],plotdata_meso$log_biomass.m2[plotdata_meso$worms==0], col="#FF7400", pch=21) points(plotdata_meso$log_worm_bmass_plusone[plotdata_meso$worms==1],plotdata_meso$log_biomass.m2[plotdata_meso$worms==1], col="#AA6C39", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod2.5_me_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod2.5_me_w)$tTable[2,5]<=0.05){lines(x,y, col="#FF7400", lwd=2.5)} if(summary(mod2.5_me_w)$tTable[2,5]>=0.05){lines(x,y, col="#FF7400", lwd=2.5, lty=2)} ### Nema plotdata_nema$log_biomass.g.per.m2 <- log10(plotdata_nema$biomass.g.per.m2) mod2.5_mi_w <- lme(log_biomass.g.per.m2~log_worm_bmass_plusone, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod2.5_mi_w) plot(plotdata_nema$log_worm_bmass_plusone,plotdata_nema$log_biomass.g.per.m2, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(-2.5,1), pch=16, frame=F ) mtext(bquote(log[10]~" fresh mass (g " ~ m^-2~ ")"), side=2, line=2.5, cex=1) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==0],plotdata_nema$log_biomass.g.per.m2[plotdata_nema$worms==0], col="#FF7400", pch=21) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==1],plotdata_nema$log_biomass.g.per.m2[plotdata_nema$worms==1], col="#AA6C39", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod2.5_mi_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod2.5_mi_w)$tTable[2,5]<=0.05){lines(x,y, col="#FF7400", lwd=2.5)} if(summary(mod2.5_mi_w)$tTable[2,5]>=0.05){lines(x,y, col="#FF7400", lwd=2.5, lty=2)} ### RICHNESS ### Macro plotdata_macro$log_richness <- log10(plotdata_macro$richness) mod3.5_ma_w <- lme(log_richness~log_worm_bmass_plusone, random=~1|forest, data=plotdata_macro) summary(mod3.5_ma_w) plot(plotdata_macro$log_worm_bmass_plusone,plotdata_macro$log_richness, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(0.8,1.8), pch=16, frame=F ) mtext(bquote(log[10]~"'species' richness"), side=2, line=2.5, cex=1) points(plotdata_macro$log_worm_bmass_plusone[plotdata_macro$worms==0],plotdata_macro$log_richness[plotdata_macro$worms==0], col="#00C2C2", pch=21) points(plotdata_macro$log_worm_bmass_plusone[plotdata_macro$worms==1],plotdata_macro$log_richness[plotdata_macro$worms==1], col="#226666", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod3.5_ma_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod3.5_ma_w)$tTable[2,5]<=0.05){lines(x,y, col="#00C2C2", lwd=2.5)} if(summary(mod3.5_ma_w)$tTable[2,5]>=0.05){lines(x,y, col="#00C2C2", lwd=2.5, lty=2)} mtext(expression(bold("Richness")), side=3, line=0.4, cex=1) ### Meso mod3.5_me_w <- lme(richness~log_worm_bmass_plusone, random=~1|forest, data=plotdata_meso) summary(mod3.5_me_w) plot(plotdata_meso$log_worm_bmass_plusone,plotdata_meso$richness, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(2,20), pch=16, frame=F ) mtext(bquote("'taxon' richness"), side=2, line=2.5, cex=1) points(plotdata_meso$log_worm_bmass_plusone[plotdata_meso$worms==0],plotdata_meso$richness[plotdata_meso$worms==0], col="#00C2C2", pch=21) points(plotdata_meso$log_worm_bmass_plusone[plotdata_meso$worms==1],plotdata_meso$richness[plotdata_meso$worms==1], col="#226666", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod3.5_me_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod3.5_me_w)$tTable[2,5]<=0.05){lines(x,y, col="#00C2C2", lwd=2.5)} if(summary(mod3.5_me_w)$tTable[2,5]>=0.05){lines(x,y, col="#00C2C2", lwd=2.5, lty=2)} ### Nema mod3.5_mi <- lme(log10(richness)~worms, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod3.5_mi) mod3.5_mi_w <- lme(richness~log_worm_bmass_plusone, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod3.5_mi_w) plot(plotdata_nema$log_worm_bmass_plusone,plotdata_nema$richness, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(10,40), pch=16, frame=F ) mtext(bquote("genus richness"), side=2, line=2.5, cex=1) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==0],plotdata_nema$richness[plotdata_nema$worms==0], col="#00C2C2", pch=21) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==1],plotdata_nema$richness[plotdata_nema$worms==1], col="#226666", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod3.5_mi_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod3.5_mi_w)$tTable[2,5]<=0.05){lines(x,y, col="#00C2C2", lwd=2.5)} if(summary(mod3.5_mi_w)$tTable[2,5]>=0.05){lines(x,y, col="#00C2C2", lwd=2.5, lty=2)} ### SHANNON ### Macro plotdata_macro$log_shannon <- log10(plotdata_macro$shannon) mod5.5_ma_w <- lme(shannon~log_worm_bmass_plusone, random=~1|forest, data=plotdata_macro) summary(mod5.5_ma_w) plot(plotdata_macro$log_worm_bmass_plusone,plotdata_macro$shannon, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(1,3.5), pch=16, frame=F ) mtext(bquote("'species' Shannon"), side=2, line=2.5, cex=1) points(plotdata_macro$log_worm_bmass_plusone[plotdata_macro$worms==0],plotdata_macro$shannon[plotdata_macro$worms==0], col="#9CEE00", pch=21) points(plotdata_macro$log_worm_bmass_plusone[plotdata_macro$worms==1],plotdata_macro$shannon[plotdata_macro$worms==1], col="#7A9E35", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod5.5_ma_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod5.5_ma_w)$tTable[2,5]<=0.05){lines(x,y, col="#9CEE00", lwd=2.5)} if(summary(mod5.5_ma_w)$tTable[2,5]>=0.05){lines(x,y, col="#9CEE00", lwd=2.5, lty=2)} mtext(expression(bold("Shannon")), side=3, line=0.4, cex=1) ### Meso mod5.5_me_w <- lme(shannon~log_worm_bmass_plusone, random=~1|forest, data=plotdata_meso) summary(mod5.5_me_w) plot(plotdata_meso$log_worm_bmass_plusone,plotdata_meso$shannon, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(0.5,3), pch=16, frame=F ) mtext(bquote("'taxon' Shannon"), side=2, line=2.5, cex=1) points(plotdata_meso$log_worm_bmass_plusone[plotdata_meso$worms==0],plotdata_meso$shannon[plotdata_meso$worms==0], col="#9CEE00", pch=21) points(plotdata_meso$log_worm_bmass_plusone[plotdata_meso$worms==1],plotdata_meso$shannon[plotdata_meso$worms==1], col="#7A9E35", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod5.5_me_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod5.5_me_w)$tTable[2,5]<=0.05){lines(x,y, col="#9CEE00", lwd=2.5)} if(summary(mod5.5_me_w)$tTable[2,5]>=0.05){lines(x,y, col="#9CEE00", lwd=2.5, lty=2)} ### Nema ### plotdata_nema$log_shannon <- log10(plotdata_nema$shannon) mod5.5_mi_w <- lme(shannon~log_worm_bmass_plusone, random=~1|forest, data=plotdata_nema, na.action=na.omit) summary(mod5.5_mi_w) plot(plotdata_nema$log_worm_bmass_plusone,plotdata_nema$shannon, xlab="", ylab="", type="n", xlim=c(0,2.5), ylim=c(1.5,3.5), pch=16, frame=F ) mtext(bquote("genus Shannon"), side=2, line=2.5, cex=1) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==0],plotdata_nema$shannon[plotdata_nema$worms==0], col="#9CEE00", pch=21) points(plotdata_nema$log_worm_bmass_plusone[plotdata_nema$worms==1],plotdata_nema$shannon[plotdata_nema$worms==1], col="#7A9E35", pch=16) x <- seq(min(log_worm_bmass_plusone),max(log_worm_bmass_plusone),length.out=length(log_worm_bmass_plusone)) y <- predict(mod5.5_mi_w, newdata=list(log_worm_bmass_plusone=x), level=0) if(summary(mod5.5_mi_w)$tTable[2,5]<=0.05){lines(x,y, col="#9CEE00", lwd=2.5)} if(summary(mod5.5_mi_w)$tTable[2,5]>=0.05){lines(x,y, col="#9CEE00", lwd=2.5, lty=2)} #dev.off() ################################################################################################################################## # Extracting model parameters for WORM BIOMASS models (n, slope estimate, confint, p-value, R2) model_names <- c("macro_abundance", "macro_biomass", "macro_richness", "macro_shannon","meso_abundance", "meso_biomass", "meso_richness", "meso_shannon","micro_abundance", "micro_biomass", "micro_richness", "micro_shannon") #### #sample size model_ns <- c(summary(mod1.5_ma_w)$dims$N, summary(mod2.5_ma_w)$dims$N, summary(mod3.5_ma_w)$dims$N, summary(mod5.5_ma_w)$dims$N, summary(mod1.5_me_w)$dims$N, summary(mod2.5_me_w)$dims$N, summary(mod3.5_me_w)$dims$N, summary(mod5.5_me_w)$dims$N, summary(mod1.5_mi_w)$dims$N, summary(mod2.5_mi_w)$dims$N, summary(mod3.5_mi_w)$dims$N, summary(mod5.5_mi_w)$dims$N) #### #slopes model_slopes <- c(coef(mod1.5_ma_w)[1,2], coef(mod2.5_ma_w)[1,2], coef(mod3.5_ma_w)[1,2], coef(mod5.5_ma_w)[1,2], coef(mod1.5_me_w)[1,2], coef(mod2.5_me_w)[1,2], coef(mod3.5_me_w)[1,2], coef(mod5.5_me_w)[1,2], coef(mod1.5_mi_w)[1,2], coef(mod2.5_mi_w)[1,2], coef(mod3.5_mi_w)[1,2], coef(mod5.5_mi_w)[1,2]) #### #confidence intervals (using intervals() approximation as it is lme's) model_intlow <- c(intervals(mod1.5_ma_w)$fixed[2,1], intervals(mod2.5_ma_w)$fixed[2,1], intervals(mod3.5_ma_w)$fixed[2,1], intervals(mod5.5_ma_w)$fixed[2,1], intervals(mod1.5_me_w)$fixed[2,1], intervals(mod2.5_me_w,which=c("fixed"))$fixed[2,1], intervals(mod3.5_me_w,which=c("fixed"))$fixed[2,1], intervals(mod5.5_me_w)$fixed[2,1], # some weird issue otherwise... intervals(mod1.5_mi_w)$fixed[2,1], intervals(mod2.5_mi_w)$fixed[2,1], intervals(mod3.5_mi_w)$fixed[2,1], intervals(mod5.5_mi_w,which=c("fixed"))$fixed[2,1]) model_intupp <- c(intervals(mod1.5_ma_w)$fixed[2,3], intervals(mod2.5_ma_w)$fixed[2,3], intervals(mod3.5_ma_w)$fixed[2,3], intervals(mod5.5_ma_w)$fixed[2,3], intervals(mod1.5_me_w)$fixed[2,3], intervals(mod2.5_me_w,which=c("fixed"))$fixed[2,3], intervals(mod3.5_me_w,which=c("fixed"))$fixed[2,3], intervals(mod5.5_me_w)$fixed[2,3], # some weird issue otherwise... intervals(mod1.5_mi_w)$fixed[2,3], intervals(mod2.5_mi_w)$fixed[2,3], intervals(mod3.5_mi_w)$fixed[2,3], intervals(mod5.5_mi_w,which=c("fixed"))$fixed[2,3]) #### #F-values model_Fvals <- c(anova(mod1.5_ma_w)[2,3] , anova(mod2.5_ma_w)[2,3] , anova(mod3.5_ma_w)[2,3] , anova(mod5.5_ma_w)[2,3] , anova(mod1.5_me_w)[2,3] , anova(mod2.5_me_w)[2,3] , anova(mod3.5_me_w)[2,3] , anova(mod5.5_me_w)[2,3] , anova(mod1.5_mi_w)[2,3] , anova(mod2.5_mi_w)[2,3] , anova(mod3.5_mi_w)[2,3] , anova(mod5.5_mi_w)[2,3] ) #### #p-values model_pvals <- c(summary(mod1.5_ma_w)$tTable[2,5], summary(mod2.5_ma_w)$tTable[2,5], summary(mod3.5_ma_w)$tTable[2,5], summary(mod5.5_ma_w)$tTable[2,5], summary(mod1.5_me_w)$tTable[2,5], summary(mod2.5_me_w)$tTable[2,5], summary(mod3.5_me_w)$tTable[2,5], summary(mod5.5_me_w)$tTable[2,5], summary(mod1.5_mi_w)$tTable[2,5], summary(mod2.5_mi_w)$tTable[2,5], summary(mod3.5_mi_w)$tTable[2,5], summary(mod5.5_mi_w)$tTable[2,5]) #### library(piecewiseSEM) model_mar_r2 <- c(rsquared(mod1.5_ma_w)[1,5], rsquared(mod2.5_ma_w)[1,5], rsquared(mod3.5_ma_w)[1,5], rsquared(mod5.5_ma_w)[1,5], rsquared(mod1.5_me_w)[1,5], rsquared(mod2.5_me_w)[1,5], rsquared(mod3.5_me_w)[1,5], rsquared(mod5.5_me_w)[1,5], rsquared(mod1.5_mi_w)[1,5], rsquared(mod2.5_mi_w)[1,5], rsquared(mod3.5_mi_w)[1,5], rsquared(mod5.5_mi_w)[1,5]) model_con_r2 <- c(rsquared(mod1.5_ma_w)[1,6], rsquared(mod2.5_ma_w)[1,6], rsquared(mod3.5_ma_w)[1,6], rsquared(mod5.5_ma_w)[1,6], rsquared(mod1.5_me_w)[1,6], rsquared(mod2.5_me_w)[1,6], rsquared(mod3.5_me_w)[1,6], rsquared(mod5.5_me_w)[1,6], rsquared(mod1.5_mi_w)[1,6], rsquared(mod2.5_mi_w)[1,6], rsquared(mod3.5_mi_w)[1,6], rsquared(mod5.5_mi_w)[1,6]) model_transf <- c("log10", "log10", "log10", "none", "log10", "log10", "none", "none", "log10", "log10", "none", "none") #### put all together model_params_wormbmass <- data.frame(model=model_names, n=model_ns, resp_transf=model_transf, slope=round(model_slopes,2), low=round(model_intlow,2), upp=round(model_intupp,2), Fval=round(model_Fvals,2), pval=round(model_pvals,2), pval_holm=round(p.adjust(model_pvals, method="holm"),2), # Holm's adjustment to p-values when running multiple tests pseudoR2mar=round(model_mar_r2,2), pseudoR2con=round(model_con_r2,2)) ############# Table 2 model_params_wormbmass ####################################################################################################################################### ############### V. Run a model-selection procedure to test if biodiversity facet (community property) identity and size group affect the response of soil fauna to ewarthworm invasion ################### ####################################################################################################################################### # Please note that what is called "community property" or "commprop" in this section is referred to as "facet" or "biodiversity facet" in the paper. # Both terms refer to: richness, abundance, biomass, shannon index. rstandardise <- function(x){(x-mean(x))/(sd(x))} # here, they are standardized to zero mean and divided by their standard deviation # Range-standardize Macrofauna data plotdata_macro_ranged <- data.frame( sizegroup = rep("macro",80), plot = plotdata_macro$plot, forest = plotdata_macro$forest, worms = plotdata_macro$worms, std_abundance = rstandardise(plotdata_macro$abundance.m2), # log10??? std_biomass = rstandardise(plotdata_macro$biomass.g.m2), std_richness = rstandardise(plotdata_macro$richness), std_shannon = rstandardise(plotdata_macro$shannon) ) # Range-standardize Mesofauna data plotdata_meso_ranged <- data.frame( sizegroup = rep("meso",80), plot = plotdata_meso$plot, forest = plotdata_meso$forest, worms = plotdata_meso$worms, std_abundance = rstandardise(plotdata_meso$abundance.m2), # log10??? std_biomass = rstandardise(plotdata_meso$biomass.m2), std_richness = rstandardise(plotdata_meso$richness), std_shannon = rstandardise(plotdata_meso$shannon) ) plotdata_nema2 <- na.omit(plotdata_nema) # Range-standardize nemafauna data plotdata_nema_ranged <- data.frame( sizegroup = rep("nema",length(plotdata_nema2$plot)), # 75 because of NA-removal... plot = plotdata_nema2$plot, forest = plotdata_nema2$forest, worms = plotdata_nema2$worms, std_abundance = rstandardise(plotdata_nema2$abundance), # log10??? std_biomass = rstandardise(plotdata_nema2$biomass.g.per.m2), std_richness = rstandardise(plotdata_nema2$richness), std_shannon = rstandardise(plotdata_nema2$shannon) ) ############################################ ### combine three frames with four different response columns for test if sizegroup is important! modeled per response variable plotdata_ranged <- rbind(plotdata_meso_ranged,plotdata_macro_ranged,plotdata_nema_ranged) #################### ### AIC comparison based on dredge AICc library(MuMIn) ############################################ ### combine three frames with four different response columns for test if community property is important! done for each size group ################# macro plotdata_ranged_macro_commprop <- data.frame(plot=rep(plotdata_macro_ranged$plot,4), forest=rep(plotdata_macro_ranged$forest,4), worms=as.factor(rep(plotdata_macro_ranged$worms,4)), commprop=as.factor(c(rep("abundance",80),rep("biomass",80),rep("richness",80),rep("shannon",80))), std_response=c(plotdata_macro_ranged$std_abundance, plotdata_macro_ranged$std_biomass, plotdata_macro_ranged$std_richness, plotdata_macro_ranged$std_shannon)) ################# meso plotdata_ranged_meso_commprop <- data.frame(plot=rep(plotdata_meso_ranged$plot,4), forest=rep(plotdata_meso_ranged$forest,4), worms=as.factor(rep(plotdata_meso_ranged$worms,4)), commprop=as.factor(c(rep("abundance",80),rep("biomass",80),rep("richness",80),rep("shannon",80))), std_response=c(plotdata_meso_ranged$std_abundance, plotdata_meso_ranged$std_biomass, plotdata_meso_ranged$std_richness, plotdata_meso_ranged$std_shannon)) ################# nema plotdata_ranged_nema_commprop <- data.frame(plot=rep(plotdata_nema_ranged$plot,4), forest=rep(plotdata_nema_ranged$forest,4), worms=as.factor(rep(plotdata_nema_ranged$worms,4)), commprop=as.factor(c(rep("abundance",78),rep("biomass",78),rep("richness",78),rep("shannon",78))), std_response=c(plotdata_nema_ranged$std_abundance, plotdata_nema_ranged$std_biomass, plotdata_nema_ranged$std_richness, plotdata_nema_ranged$std_shannon)) ############################################################################################################################ ## just for fun: get all data in one df (all size groups, all commprops), run dredge on model response~worms*sizegroup*commprop # get sizegroup column included in singe size group data frames: plotdata_ranged_macro_commprop$sizegroup <- as.factor(rep("macro",80)) plotdata_ranged_meso_commprop$sizegroup <- as.factor(rep("meso",80)) plotdata_ranged_nema_commprop$sizegroup <- as.factor(rep("nema",78)) plotdata_ranged_all <- rbind(plotdata_ranged_macro_commprop,plotdata_ranged_meso_commprop,plotdata_ranged_nema_commprop) mod_r_all <- lme(std_response~worms*sizegroup*commprop, random=~1|forest, method="ML", data=plotdata_ranged_all) # ML for AIC comparison summary(mod_r_all) modsset <- dredge(mod_r_all) ######################### Supporting Information Table S6 modsset