# R-Code accompanying Radny & Meyer: # The role of biotic factors during plant establishment in novel # communities assessed with an agent-based simulation model # Aim: Analyze and plot local-scale simulation model output on local establishment # in terms of survival of non-native plant species with respect # to competition and herbivory as well as competitive and defensive traits # No warranty for anything # Last update: June 28, 2018 # Read simulation model output of last time step (end of first year) radny <- read.table("radny.txt", header=TRUE) # Turn Variables comp and def into factor radny$comp <- as.factor(radny$comp) radny$def <- as.factor(radny$def) # Calculate and add survival rate radny$survivalrate <- radny$finalEx / radny$initialEx # Variables str(radny) # $ run : Simulation run (replicate) # $ invasion.level : Initial number of non-native species # $ herbivore.density: Number of herbivores # $ a.symmetry : Competition asymmetry parameter theta # $ density : Density of native species # $ mixture : Mixture of native species # $ initialEx : Initial number of non-natives # $ finalEx : Final number of non-natives at the end of the first year # $ comp : Competitive strength trait # $ def : Defensive strength trait # $ survivalrate : Proportion of initial non-native plants # that survived until the end of the simulated year # --------- # Analysis # --------- # Constructing the y for a binomial model for non-native species number.of.survivors <- radny$finalEx number.of.deaths <- radny$initialEx - radny$finalEx survival.binomial <- cbind(number.of.survivors, number.of.deaths) # Analyze survival of non-native plants with glm # with main effects and all two-way interactions model.survival <- glm(survival.binomial ~ comp * def + initialEx + herbivore.density + a.symmetry + density + mixture + comp:initialEx + comp:herbivore.density + comp:a.symmetry + comp:density + comp:mixture + def:initialEx + def:herbivore.density + def:a.symmetry + def:density + def:mixture + initialEx:herbivore.density + initialEx:a.symmetry + initialEx:density + initialEx:mixture + herbivore.density:a.symmetry + herbivore.density:density + herbivore.density:mixture + a.symmetry:density + a.symmetry:mixture + density:mixture, data=radny, family=quasibinomial) summary(model.survival) # Model simplification by hand: # Always removing the term with the highest p-value # starting with the interactions # until only significant terms are retained model.survival.1 <- update(model.survival, . ~ . -density:mixture) summary(model.survival.1) model.survival.2 <- update(model.survival.1, . ~ . -def:mixture) summary(model.survival.2) model.survival.3 <- update(model.survival.2, . ~ . -comp:def) summary(model.survival.3) # (Note that non-significant main effects are retained # when they are part of a significant interaction) # Check model fit pchisq(model.survival.3$null.deviance - model.survival.3$deviance, model.survival.3$df.null - model.survival.3$df.residual, lower=F) # -> p = 0, so model is significantly better than null model. # Check explicitly for over/underdispersion deviance(model.survival.3) / df.residual(model.survival.3) # -> 2.05 -> slight overdispersion -> family=quasibinomial accounts for that # -------------------------------- # Plotting figures for publication # -------------------------------- # Load library library(ggplot2) # Recode individual plant profiles attach(radny) radny$profile2[comp == 1 & def == 1] <- "HiAll" radny$profile2[comp == 1 & def == 0] <- "HiComp" radny$profile2[comp == 0 & def == 1] <- "HiDef" radny$profile2[comp == 0 & def == 0] <- "LowAll" detach(radny) # Define mytheme for publishable graphics mytheme<-theme(panel.grid.minor = element_line(colour = "white"), title = element_text(size = rel(1.5), hjust = 0.5, vjust = rel(0.5)), axis.text.x = element_text(size = rel(1.5), colour = "black"), axis.text.y = element_text(size = rel(1.5), colour = "black"), axis.line = element_line(colour = "black"), panel.background = element_rect(fill = F), strip.text.x = element_text(size = rel(1.5), colour = "black", hjust = 0.0), strip.text.y = element_text(size = rel(1.5), colour = "black"), strip.background = element_blank(), legend.text = element_text(size = rel(1)), legend.title = element_text(size = rel (1.1)) ) # Figure 4: Plot survival of non-natives with respect to trait profiles traits = ggplot(radny, aes(x=profile2, y=survivalrate, fill=profile2)) + geom_boxplot() + scale_y_continuous(name="Survival rate of non-native plants") + scale_x_discrete(name="Profile") + scale_fill_discrete(name="Profile") + mytheme + guides(fill=FALSE) traits # Figure 5: Plot survival of non-natives with respect to initial number of non-natives, # community density and trait profiles: # First relabel density high -> A, low -> B radny$density <- factor(radny$density, labels=c("A", "B")) # then plot traits_invasionlevel = ggplot(radny, aes(x=as.factor(initialEx), y=survivalrate, fill=profile2)) + geom_boxplot() + scale_y_continuous(name="Survival rate of non-native plants") + scale_x_discrete(name="Initial number of non-native plants") + scale_fill_discrete(name="Profile") + facet_grid(.~density) + mytheme + guides(fill=FALSE) traits_invasionlevel # Label back to original radny$density <- factor(radny$density, labels=c("high", "low")) # Figure 6: Plot survival of non-natives with respect to intensity of competition (theta), # community density, community mixture, and trait profiles # 'even mixture': # First relabel density high -> A, low -> B radny$density <- factor(radny$density, labels=c("A", "B")) # then plot for 'even' mixture traits_asymmetry = ggplot(radny[radny$mixture == "even", ], aes(x=as.factor(a.symmetry), y=survivalrate, fill=profile2)) + geom_boxplot() + scale_y_continuous(name="Survival rate of non-native plants") + scale_x_discrete(name="Intensity of competition (theta)") + scale_fill_discrete(name="Profile") + mytheme + facet_grid(.~density) + guides(fill=FALSE) traits_asymmetry # 'seedbank mixture': # First relabel density high -> C, low -> D radny$density <- factor(radny$density, labels=c("C", "D")) # then plot for 'seedbank' mixture traits_asymmetry = ggplot(radny[radny$mixture == "seedbank", ], aes(x=as.factor(a.symmetry), y=survivalrate, fill=profile2)) + geom_boxplot() + scale_y_continuous(name="Survival rate of non-native plants") + scale_x_discrete(name="Intensity of competition (theta)") + scale_fill_discrete(name="Profile") + mytheme + facet_grid(.~density) + guides(fill=FALSE) traits_asymmetry # Label back to original radny$density <- factor(radny$density, labels=c("high", "low")) # Figure 7: Plot survival of non-natives with respect to trait profile # and herbivore density traits_herbivores = ggplot(radny, aes(x=as.factor(profile2), y=survivalrate, fill=as.factor(herbivore.density))) + geom_boxplot() + scale_y_continuous(name="Survival rate of non-native plants") + scale_x_discrete(name="Profile") + scale_fill_manual(values=c("white", "lightgrey", "darkgrey")) + mytheme + guides(fill=FALSE) traits_herbivores # ------------------------------------------------------------------- # Sensitivity analysis # ------------------------------------------------------------------- # Note that simulation model output was too large (800MB) to attach # Therefore, this script is for illustrational purpose only # Read simulation model output sensi1 <- read.table("sensi1.csv", header = TRUE, sep=";", dec = ".") sensi2 <- read.table("sensi2.csv", header = TRUE, sep=";", dec = ".") sensi3 <- read.table("sensi3.csv", header = TRUE, sep=";", dec = ".") sensi4 <- read.table("sensi4.csv", header = TRUE, sep=";", dec = ".") sensi5 <- read.table("sensi5.csv", header = TRUE, sep=";", dec = ".") sensi6 <- read.table("sensi6.csv", header = TRUE, sep=";", dec = ".") sensi7 <- read.table("sensi7.csv", header = TRUE, sep=";", dec = ".") sensi8 <- read.table("sensi8.csv", header = TRUE, sep=";", dec = ".") sensi9 <- read.table("sensi9.csv", header = TRUE, sep=";", dec = ".") sensi10 <- read.table("sensi10.csv", header = TRUE, sep=";", dec = ".") sensi11 <- read.table("sensi11.csv", header = TRUE, sep=";", dec = ".") sensi12 <- read.table("sensi12.csv", header = TRUE, sep=";", dec = ".") sensi13 <- read.table("sensi13.csv", header = TRUE, sep=";", dec = ".") sensi14 <- read.table("sensi14.csv", header = TRUE, sep=";", dec = ".") sensi15 <- read.table("sensi15.csv", header = TRUE, sep=";", dec = ".") # Bind all results into columns of a data.frame sensitivity <- rbind(sensi1, sensi2, sensi3, sensi4, sensi5, sensi6, sensi7, sensi8, sensi9, sensi10, sensi11, sensi12, sensi13, sensi14, sensi15) # Calculate survival rate of non-natives and attach it to the data.frame sensitivity$survivalrate.exo <- sensitivity$finalEx / sensitivity$initialEx # Only consider data of the maximum time step (This is usually 75, # but might be less when non-natives go extinct earlier) survival <- setNames(aggregate(sensitivity$step, by = list(sensitivity$run, sensitivity$invasion.level, sensitivity$herbivore_density, sensitivity$density, sensitivity$a.symmetry, sensitivity$seedweight.nat, sensitivity$biomass_max.nat, sensitivity$repulse.nat, sensitivity$germprob.nat, sensitivity$seedweight.exo, sensitivity$biomass_max.exo, sensitivity$repulse.exo, sensitivity$germprob.exo, sensitivity$survivalrate, sensitivity$survivalrate.exo, sensitivity$initialEx, sensitivity$finalEx), which.max), c("run", "invasion.level", "herbivore_density", "density", "a.symmetry", "seedweight.nat", "biomass_max.nat", "repulse.nat", "germprob.nat", "seedweight.exo", "biomass_max.exo", "repulse.exo", "germprob.exo", "survivalrate.nat", "survivalrate.exo", "initial.exo", "final.exo", "timestep")) # -------- # Analysis # -------- # Analyze survival of non-native plants with glm with main effects and all two-way interactions # Constructing the y for a binomial model for non-native species number.of.survivors <- survival$final.exo number.of.deaths <- survival$initial.exo - survival$final.exo survival.binomial <- cbind(number.of.survivors, number.of.deaths) # Binomial model model.full.twoway <- glm(survival.binomial ~ (invasion.level + herbivore_density + density + a.symmetry + seedweight.nat + biomass_max.nat + repulse.nat + germprob.nat + seedweight.exo + biomass_max.exo + repulse.exo + germprob.exo) ^ 2, data=survival, family=quasibinomial) # quasibinomial due to overdispersion summary(model.full.twoway) # Model simplification by hand: # always removing the least significant term # starting with the interactions and keeping main terms # that are part of a significant interaction model.full.twoway1 <- update(model.full.twoway, . ~ . -repulse.exo:germprob.exo) summary(model.full.twoway1) model.full.twoway2 <- update(model.full.twoway1, . ~ . -biomass_max.exo:germprob.exo) summary(model.full.twoway2) model.full.twoway3 <- update(model.full.twoway2, . ~ . -biomass_max.exo:repulse.exo) summary(model.full.twoway3) model.full.twoway4 <- update(model.full.twoway3, . ~ . -seedweight.exo:germprob.exo ) summary(model.full.twoway4) model.full.twoway5 <- update(model.full.twoway4, . ~ . -seedweight.exo:repulse.exo ) summary(model.full.twoway5) model.full.twoway6 <- update(model.full.twoway5, . ~ . -seedweight.exo:biomass_max.exo ) summary(model.full.twoway6) model.full.twoway7 <- update(model.full.twoway6, . ~ . -germprob.nat:germprob.exo) summary(model.full.twoway7) model.full.twoway8 <- update(model.full.twoway7, . ~ . -germprob.nat:repulse.exo) summary(model.full.twoway8) model.full.twoway9 <- update(model.full.twoway8, . ~ . -germprob.nat:biomass_max.exo) summary(model.full.twoway9) model.full.twoway10 <- update(model.full.twoway9, . ~ . -germprob.nat:seedweight.exo ) summary(model.full.twoway10) model.full.twoway11 <- update(model.full.twoway10, . ~ . -repulse.nat:germprob.exo ) summary(model.full.twoway11) model.full.twoway12 <- update(model.full.twoway11, . ~ . -repulse.nat:repulse.exo) summary(model.full.twoway12) model.full.twoway13 <- update(model.full.twoway12, . ~ . -repulse.nat:biomass_max.exo) summary(model.full.twoway13) model.full.twoway14 <- update(model.full.twoway13, . ~ . -repulse.nat:seedweight.exo) summary(model.full.twoway14) model.full.twoway15 <- update(model.full.twoway14, . ~ . -repulse.nat:germprob.nat) summary(model.full.twoway15) model.full.twoway16 <- update(model.full.twoway15, . ~ . -biomass_max.nat:germprob.exo) summary(model.full.twoway16) model.full.twoway17 <- update(model.full.twoway16, . ~ . -biomass_max.nat:repulse.exo) summary(model.full.twoway17) model.full.twoway18 <- update(model.full.twoway17, . ~ . -biomass_max.nat:biomass_max.exo) summary(model.full.twoway18) model.full.twoway19 <- update(model.full.twoway18, . ~ . -biomass_max.nat:seedweight.exo) summary(model.full.twoway19) model.full.twoway20 <- update(model.full.twoway19, . ~ . -biomass_max.nat:germprob.nat) summary(model.full.twoway20) model.full.twoway21 <- update(model.full.twoway20, . ~ . -biomass_max.nat:repulse.nat) summary(model.full.twoway21) model.full.twoway22 <- update(model.full.twoway21, . ~ . -seedweight.nat:germprob.exo) summary(model.full.twoway22) model.full.twoway23 <- update(model.full.twoway22, . ~ . -density:repulse.exo) # from here on removing interactions with highest p-values summary(model.full.twoway23) model.full.twoway24 <- update(model.full.twoway23, . ~ . -herbivore_density:a.symmetry) summary(model.full.twoway24) model.full.twoway25 <- update(model.full.twoway24, . ~ . -herbivore_density:seedweight.exo) summary(model.full.twoway25) model.full.twoway26 <- update(model.full.twoway25, . ~ . -herbivore_density:biomass_max.exo) summary(model.full.twoway26) model.full.twoway27 <- update(model.full.twoway26, . ~ . -a.symmetry:germprob.nat) summary(model.full.twoway27) model.full.twoway28 <- update(model.full.twoway27, . ~ . -invasion.level:seedweight.nat) summary(model.full.twoway28) model.full.twoway29 <- update(model.full.twoway28, . ~ . -herbivore_density:biomass_max.nat) summary(model.full.twoway29) model.full.twoway30 <- update(model.full.twoway29, . ~ . -herbivore_density:repulse.nat) summary(model.full.twoway30) model.full.twoway31 <- update(model.full.twoway30, . ~ . -invasion.level:herbivore_density) summary(model.full.twoway31) model.full.twoway32 <- update(model.full.twoway31, . ~ . -invasion.level:biomass_max.nat ) summary(model.full.twoway32) model.full.twoway33 <- update(model.full.twoway32, . ~ . -density:germprob.exo) summary(model.full.twoway33) model.full.twoway34 <- update(model.full.twoway33, . ~ . -herbivore_density:repulse.exo) summary(model.full.twoway34) model.full.twoway35 <- update(model.full.twoway34, . ~ . -seedweight.nat:biomass_max.nat) summary(model.full.twoway35) model.full.twoway36 <- update(model.full.twoway35, . ~ . -herbivore_density:seedweight.nat) summary(model.full.twoway36) model.full.twoway37 <- update(model.full.twoway36, . ~ . -invasion.level:density) summary(model.full.twoway37) model.full.twoway38 <- update(model.full.twoway37, . ~ . -invasion.level:biomass_max.exo) summary(model.full.twoway38) model.full.twoway39 <- update(model.full.twoway38, . ~ . -density:repulse.nat) summary(model.full.twoway39) model.full.twoway40 <- update(model.full.twoway39, . ~ . -density:biomass_max.nat) summary(model.full.twoway40) model.full.twoway41 <- update(model.full.twoway40, . ~ . -invasion.level:repulse.nat) summary(model.full.twoway41) model.full.twoway42 <- update(model.full.twoway41, . ~ . -invasion.level:repulse.exo) summary(model.full.twoway42) model.full.twoway43 <- update(model.full.twoway42, . ~ . -invasion.level:a.symmetry) summary(model.full.twoway43) model.full.twoway44 <- update(model.full.twoway43, . ~ . -a.symmetry:biomass_max.exo) summary(model.full.twoway44) model.full.twoway45 <- update(model.full.twoway44, . ~ . -density:a.symmetry) summary(model.full.twoway45) model.full.twoway46 <- update(model.full.twoway45, . ~ . -herbivore_density:germprob.nat) summary(model.full.twoway46) # Call: # glm(formula = survival.binomial ~ invasion.level + herbivore_density + # density + a.symmetry + seedweight.nat + biomass_max.nat + # repulse.nat + germprob.nat + seedweight.exo + biomass_max.exo + # repulse.exo + germprob.exo + invasion.level:germprob.nat + # invasion.level:seedweight.exo + invasion.level:germprob.exo + # herbivore_density:density + herbivore_density:germprob.exo + # density:seedweight.nat + density:germprob.nat + density:seedweight.exo + # density:biomass_max.exo + a.symmetry:seedweight.nat + a.symmetry:biomass_max.nat + # a.symmetry:repulse.nat + a.symmetry:seedweight.exo + a.symmetry:repulse.exo + # a.symmetry:germprob.exo + seedweight.nat:repulse.nat + seedweight.nat:germprob.nat + # seedweight.nat:seedweight.exo + seedweight.nat:biomass_max.exo + # seedweight.nat:repulse.exo, family = quasibinomial, data = survival) # # Deviance Residuals: # Min 1Q Median 3Q Max # -21.485 -5.728 -1.193 4.873 27.290 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 3.569e+00 1.323e-01 26.967 < 2e-16 *** # invasion.level -1.025e-01 8.335e-03 -12.295 < 2e-16 *** # herbivore_density 1.617e-02 1.022e-02 1.582 0.113682 # density"low" 3.166e-01 2.565e-02 12.342 < 2e-16 *** # a.symmetry -1.507e+00 7.017e-02 -21.480 < 2e-16 *** # seedweight.nat -1.006e+00 4.406e-02 -22.830 < 2e-16 *** # biomass_max.nat -1.049e-04 2.630e-06 -39.898 < 2e-16 *** # repulse.nat -1.083e+00 8.729e-02 -12.408 < 2e-16 *** # germprob.nat -4.624e+00 1.672e-01 -27.659 < 2e-16 *** # seedweight.exo -9.003e-01 3.446e-02 -26.130 < 2e-16 *** # biomass_max.exo 3.375e-04 8.928e-06 37.796 < 2e-16 *** # repulse.exo -2.556e+00 1.254e-01 -20.375 < 2e-16 *** # germprob.exo -1.951e-01 7.987e-02 -2.442 0.014597 * # invasion.level:germprob.nat 2.145e-01 1.729e-02 12.404 < 2e-16 *** # invasion.level:seedweight.exo -1.678e-02 2.235e-03 -7.506 6.11e-14 *** # invasion.level:germprob.exo -5.748e-02 1.695e-02 -3.391 0.000698 *** # herbivore_density:density"low" -7.922e-02 1.107e-02 -7.157 8.29e-13 *** # herbivore_density:germprob.exo -1.442e-01 2.168e-02 -6.652 2.90e-11 *** # density"low":seedweight.nat 1.040e-01 7.583e-03 13.716 < 2e-16 *** # density"low":germprob.nat 1.460e-01 3.688e-02 3.960 7.51e-05 *** # density"low":seedweight.exo -7.399e-02 5.685e-03 -13.014 < 2e-16 *** # density"low":biomass_max.exo -1.350e-05 1.572e-06 -8.586 < 2e-16 *** # a.symmetry:seedweight.nat -4.240e-01 1.076e-02 -39.398 < 2e-16 *** # a.symmetry:biomass_max.nat -3.230e-05 2.318e-06 -13.935 < 2e-16 *** # a.symmetry:repulse.nat 5.775e-01 5.074e-02 11.383 < 2e-16 *** # a.symmetry:seedweight.exo 3.672e-01 8.734e-03 42.037 < 2e-16 *** # a.symmetry:repulse.exo 1.599e-01 5.409e-02 2.957 0.003112 ** # a.symmetry:germprob.exo 4.843e-01 4.936e-02 9.812 < 2e-16 *** # seedweight.nat:repulse.nat 1.232e+00 5.003e-02 24.622 < 2e-16 *** # seedweight.nat:germprob.nat 2.487e+00 7.484e-02 33.238 < 2e-16 *** # seedweight.nat:seedweight.exo 4.159e-01 1.176e-02 35.367 < 2e-16 *** # seedweight.nat:biomass_max.exo -1.097e-04 3.306e-06 -33.169 < 2e-16 *** # seedweight.nat:repulse.exo -4.263e-01 7.745e-02 -5.504 3.72e-08 *** # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # # (Dispersion parameter for quasibinomial family taken to be 52.59424) # # Null deviance: 12811887 on 161289 degrees of freedom # Residual deviance: 10904051 on 161257 degrees of freedom # AIC: NA # # Number of Fisher Scoring iterations: 6 # Rename resulting model model.twoway.final <- model.full.twoway46 # Check model fit pchisq(model.twoway.final$null.deviance - model.twoway.final$deviance, model.twoway.final$df.null - model.twoway.final$df.residual, lower=F) # -> p = 0, so model fit is ok. # Check explicitly for over/underdispersion deviance(model.twoway.final) / df.residual(model.twoway.final) # -> 67.62 -> Overdispersion -> but family=quasibinomial already accounted for that # Standardized regression coefficients: # Estimate / Standard Error (excluding the Intercept) coef.standardized <- summary(model.twoway.final)$coefficients[-1, 1] / summary(model.twoway.final)$coefficients[-1, 2] # Normalized absolute standardized regression coefficients: sensitivities <- abs(coef.standardized) / sum(abs(coef.standardized)) # Sort and turn into data.frame for plotting ii <- order(sensitivities) names(sensitivities)[ii] <- c("DensHerb", "tGermProb2", "Theta:tRepuls2", "InitNN:tGermProb2", "DensComm:tGermProb1", "tSeedWeight1:tRepuls2", "DensHerb:tGermProb2", "DensHerb:DensComm", "InitNN:tSeedWeight2", "DensComm:tMaxBiomass2", "Theta:tGermProb2", "Theta:tRepuls1", "InitNN", "DensComm", "InitNN:tGermProb1", "tRepuls1", "DensComm:tSeedWeight2", "DensComm:tSeedWeight1", "Theta:tMaxBiomass1", "tRepuls2", "Theta", "tSeedweight1", "tSeedWeight1:tRepuls1", "tSeedweight2", "tGermProb1", "tSeedWeight1:tMaxBiomass2", "tSeedWeight1:tGermProb1", "tSeedWeight1:tSeedWeight2", "tMaxBiomass2", "Theta:tSeedWeight1", "tMaxBiomass1", "Theta:tSeedWeight2") sensi.data.frame <- data.frame(Parameters=names(sensitivities[ii]), Sensitivities=sensitivities[ii], row.names=NULL) # Plotting sensitivities library(ggplot2) # Define mytheme for publishable graphics mytheme<-theme(panel.grid.minor = element_line(colour = "white"), title = element_text(size = rel(1.5), hjust = 0.5, vjust = rel(0.5)), axis.text.x = element_text(size = rel(1.5), colour = "black"), axis.text.y = element_text(size = rel(0.8), colour = "black"), axis.line = element_line(colour = "black"), panel.background = element_rect(fill = F), strip.text.x = element_text(size = rel(1.5), colour = "black", hjust = 0.0), strip.text.y = element_text(size = rel(1.5), colour = "black"), strip.background = element_blank(), legend.text = element_text(size = rel(1)), legend.title = element_text(size = rel (1.1)) ) # Figure 8: Plot sorted sensitivities sensi <- ggplot(sensi.data.frame, aes(x=Parameters, y=Sensitivities)) + geom_col(color="black", fill="black", width=0.75) + scale_y_continuous(name="Standardized sensitivities") + scale_x_discrete(name="Model parameters", limits=sensi.data.frame$Parameters) + coord_flip() + mytheme sensi