# This script prepares the raw data for analysis # Journal: Journal of Ecology # Title: Negative density dependence in the mortality and growth of tropical tree seedlings is strong, and primarily caused by fungal pathogens # Authors: Kirstie Hazelwood, Harald Beck, C. E. Timothy Paine* # * Corresponding Author: C. E. T. Paine, Environmental and Rural Science, University of New England, Armidale, New South Wales 2350 Australia. Email: timothy.paine@une.edu.au #### Load some required libraries #### library(lme4) library(ggplot2) library(pbkrtest) library(effects) library(blmeco) library(ggpubr) library(DHARMa) library(doSNOW) library(MuMIn) library(RColorBrewer) source("code/Helper functions.R") #Start a cluster ncores <- parallel:::detectCores(all.tests = T)-1 clust <- makeCluster(ncores, type = 'SOCK') registerDoSNOW(clust) clusterEvalQ(clust, library(lme4, pbkrtest)) #stopCluster(clust) clust #### import data #### dat <- read.csv('data/Processed data.csv') divsdl <- read.csv('data/Diversity data.csv') all.equal(levels(divsdl$Treatment), levels(dat$Treatment)) # check to make sure that teh two datasets are set up similarly: they are. #### Adjust Dataset #### ## Remove ferns, herbs and Wendlandiellas, unidentified and those identified only to family, and compute summary statistics ## dat <- dat[dat$ht<101,] # remove saplings (they remain in the analysis until they exceed 1m and are not counted as dead) dat <- dat[dat$genus!="Wendlandiella",] dat$ht[dat$family =='Arecaceae'] <- NA # Palms are super hard to measure the height of. So don't. ONly keep them for survival length(unique(dat$tag)) # 10557 length(unique(dat$Name)) # 638 spp or morphospecies dat <- dat[dat$ID.to!="Unknown",] dat <- dat[dat$ID.to!="Family",] length(unique(dat$tag)) # 9240 (1317 individuals were unidentified, or ided only to family) ## Remove very rare species, those with < 10 individuals ## abund <- tapply(dat$tag, dat$Name, function(x){length(unique(x))}) # hist(sort(abund), breaks = 30) dat <- dat[dat$Name %in% names(abund[abund >= 10]),] # Keep only species with >= 10 individuals length(unique(dat$tag)) # 8018 length(unique(dat$Name)) # 149 spp or morphospecies round(100 * table(dat$ID.to)/nrow(dat), 1)# Family 2.3%, Genus 20.8%, Species 76.8% # Add some columns dat$consp.crowd.l <- log(1+dat$consp.crowd) dat$hetsp.crowd.l <- log(1+dat$hetsp.crowd) dat$nei.crowd.l <- log(1+dat$nei.crowd) dat$ht.l <- log(dat$ht) dat$offset_time <- log(dat$difftime/365.25) dat$plot <- factor(dat$plot) dat$trans <- factor(dat$trans) # make colors for use in plotting COL <- c(brewer.pal(6, "Spectral")) COL <- c(COL[-4], "#846386") names(COL) <- levels(factor(dat$Treatment))[c(1, 2, 5, 6, 7, 3)] COL <- c(COL, `Large Mammals` = "#CD3333", `Small Mammals` = "#EE3A8C", Control = '#BEBEBE') COL_band <- paste(COL, 30, sep = "") names(COL_band) <- names(COL) pie(rep(1, 9), col = COL, labels = names(COL)) #### Figure 1: NDD effects on mortality and RGR in all control & no-treatment plots #### # Filter dataset # table(dat$Treatment) dat1 <- dat[dat$Treatment %in% c("Control", "No Treatment"),] # non-treatment plots before the start of the experimental study table(dat1$Treatment) # some indivioduals had extremen growth or shrinkage, usyally because they fell over, or sprang back up. dat1[dat1$rgr > 2 & !is.na(dat1$rgr),] dat1[dat1$rgr < -2 & !is.na(dat1$rgr),] dim(dat1) temp1 <- na.omit(dat1[dat1$rgr > -2 & dat1$rgr < 2 & !is.na(dat1$rgr), c("rgr", "consp.crowd.l", "hetsp.crowd.l", "trans", "plot", "Census", "Name", 'ht.l'),]) # cut down dat to the key variables, and exclude observations where size changes unreasonably fast. Most of these are seedlings that either fell down, or had fallen down and now stood back up. dim(temp1) # Mortality # Compare random effect structures m1.mort.0 <- glmer(event ~ consp.crowd.l + hetsp.crowd.l + ht.l + (1|trans/plot) + (1|Census), offset=offset_time, family=binomial(link=cloglog), data=dat1, verbose = 2) m1.mort.1 <- update(m1.mort.0, . ~ . + (1|Name)) m1.mort.2a <- update(m1.mort.0, . ~ . + (consp.crowd.l|Name)) m1.mort.2b <- update(m1.mort.0, . ~ . + (hetsp.crowd.l|Name)) m1.mort.2c <- update(m1.mort.0, . ~ . + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m1.mort.0, m1.mort.1, m1.mort.2a, m1.mort.2b, m1.mort.2c) # strong evidence for different slopes among spp for both types of crowding m1.mort.3a <- update(m1.mort.2c, . ~ . + consp.crowd.l*ht.l) m1.mort.3b <- update(m1.mort.2c, . ~ . + hetsp.crowd.l*ht.l) m1.mort.3c <- update(m1.mort.2c, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m1.mort.2c, m1.mort.3a, m1.mort.3b, m1.mort.3c) # NO sign that height affects responses to either type of crowding summary(m1.mort.2c) r.squaredGLMM(m1.mort.2c) # R2m R2c, theoretical 0.2383904 0.4725918, delta 0.1698800 0.3367749 pval_calc(m1.mort.2c, nsim = 1000) # X.Intercept. consp.crowd.l hetsp.crowd.l ht.l # 0.000 0.000 0.102 0.000 inverse_cloglog(predict(m1.mort.2c, newdata = data.frame(ht.l = median(dat1$ht.l, na.rm = T), consp.crowd.l = 0, hetsp.crowd.l = 0), re.form = NA)) # 25.3% mortality per year for median-sized seedlings without neighbors inverse_cloglog(predict(m1.mort.2c, newdata = data.frame(ht.l = median(dat1$ht.l, na.rm = T), consp.crowd.l = c(log(1+1), log(10+1)), hetsp.crowd.l = 0), re.form = NA)) # Mortality increased from 27.1% to 32.1% with an order of magnitude increase in conspecific crowding # Examine the residuals for the model par(mfrow = c(1,1)) simulationOutput <- simulateResiduals(fittedModel = m1.mort.2c, plot = T) # looks good par(mfrow = c(2, 2)) plot(residuals(simulationOutput), dat1$consp.crowd.l) plot(residuals(simulationOutput), dat1$hetsp.crowd.l) plot(residuals(simulationOutput), dat1$ht.l) dispersion_glmer(m1.mort.2c) #0.9962439, no prob # All residuals look OK, also when plotted against predictor variables # Growth m1.grow.0 <- lmer(rgr ~ consp.crowd.l + hetsp.crowd.l + (1|trans/plot) + (1|Census), temp1, REML=FALSE, verbose = 1) m1.grow.1 <- update(m1.grow.0, . ~ . + (1|Name)) m1.grow.2a <- update(m1.grow.0, . ~ . + (consp.crowd.l|Name)) m1.grow.2b <- update(m1.grow.0, . ~ . + (hetsp.crowd.l|Name)) m1.grow.2c <- update(m1.grow.0, . ~ . + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m1.grow.0, m1.grow.1, m1.grow.2a, m1.grow.2b, m1.grow.2c) # strong evidence for diff slopes among spp for both types of crowding on rgr summary(grow.2c) m1.grow.3a <- update(m1.grow.2c, . ~ . + consp.crowd.l*ht.l) m1.grow.3b <- update(m1.grow.2c, . ~ . + hetsp.crowd.l*ht.l) m1.grow.3c <- update(m1.grow.2c, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m1.grow.2c, m1.grow.3a, m1.grow.3b, m1.grow.3c) # Height strongly affects responses to both types of crowding m1.grow.3c.ss <- getME(m1.grow.3c, c("theta","fixef")) m1.grow.3c <- update(m1.grow.3c, start=grow.3c.ss, control=lmerControl(optCtrl=list(maxfun=2e4))) summary(m1.grow.3c) predict(m1.grow.3c, newdata = data.frame(ht.l = median(dat1$ht.l, na.rm = T), consp.crowd.l = 0, hetsp.crowd.l = 0), re.form = NA) # 0.061 cm/cm/year for median-sized seedlings without neighbors r.squaredGLMM(m1.grow.3c) # R2m R2c: 0.06571776 0.140407 pval_calc(m1.grow.3c, nsim = 10000) # X.Intercept. consp.crowd.l hetsp.crowd.l ht.l consp.crowd.l.ht.l hetsp.crowd.l.ht.l # 0.0000 0.0000 0.0000 0.0016 0.0000 0.0000 # examine model residuals par(mfrow = c(1, 1)); simulationOutput <- simulateResiduals(fittedModel = grow.2c, plot = T) # bad QQ plot testZeroInflation(simulationOutput) # predict results for plotting # range(dat1$consp.crowd.l, na.rm = T) consp.nd <- expand.grid( consp.crowd.l = seq(0,4.3,length = 101), hetsp.crowd.l = median(dat1$hetsp.crowd.l, na.rm = T), ht.l = log(c(10, 20, 50, 100)), type = 'consp') # range(dat1$hetsp.crowd.l, na.rm = T) hetsp.nd <- expand.grid( consp.crowd.l = median(temp1$consp.crowd.l, na.rm = T), hetsp.crowd.l = seq(0, 5.8,length = 101), ht.l = log(c(10, 20, 50, 100)), type = 'hetsp') nd <- rbind(consp.nd, hetsp.nd) nd$pred_mort <- inverse_cloglog(predict(m1.mort.2c, newdata = nd, re.form = NA)) nd$pred_grow <- predict(m1.grow.3c, newdata = nd, re.form = NA) boot_mort <- data.frame(bootMer(m1.mort.2c, FUN = function(x, newdata){inverse_cloglog(predict(x, newdata = nd, re.form = NA))}, nsim = 1000, parallel = 'multicore', ncpus = ncores, verbose = T)) boot_grow <- data.frame(bootMer(m1.grow.3c, FUN = function(x, newdata){predict(x, newdata = nd, re.form = NA)}, nsim = 1000, parallel = 'multicore', ncpus = ncores, verbose = T)) nd$CI_mort <- t(apply(boot_mort, 2, FUN = quantile, c(0.025, 0.975))) nd$CI_grow <- t(apply(boot_grow, 2, FUN = quantile, c(0.025, 0.975))) nd.consp <- nd[nd$type == 'consp',] nd.hetsp <- nd[nd$type == 'hetsp',] pdf("figs/Figure 1 Overview of NDD.pdf", paper = 'a4', useDingbats = F) BREAKS_con <- seq(0, max(dat1$consp.crowd.l, na.rm = T), length = 20) BREAKS_het <- seq(0, max(dat1$hetsp.crowd.l, na.rm = T), length = 20) Ymax <- 0.65 Hts <- unique(nd$ht.l)[c(1,2,3)] par(mfrow = c(2, 2), las = 1, bty = "L", mar = c(2, 2, 1.5, 0.5), oma = c(3, 3, 1, 3)) # mortality - conspecific NDD plot(event ~ consp.crowd.l, data = dat1, xlab = "", ylab = "", type = 'n', ylim = c(0, Ymax), xaxt = "n") axis(1, labels = NA) title(ylab = expression(Probability~~of~~mortality~~(y^-1)), xpd = NA) top_hist <- hist(dat1$consp.crowd.l[dat1$event == 1], breaks = BREAKS_con, plot = F) low_hist <- hist(dat1$consp.crowd.l[dat1$event == 0], breaks = BREAKS_con, plot = F) rect(BREAKS_con[-20]+0.02, Ymax, BREAKS_con[-1]-0.02, Ymax - 0.3*top_hist$counts/6152, col= 'gray', border = NA) rect(BREAKS_con[-20]+0.02, 0, BREAKS_con[-1]-0.02, 0.3*low_hist$counts/6152, col= 'gray', border = NA) axis(4, at = Ymax - 0.3*(0:2*2000)/6152, labels = NA, xpd = NA, las = 1, line = 0, cex.axis = 0.6, xpd = NA) axis(4, at = 0.3*(0:2*2000)/6152, labels = NA, xpd = NA, las = 1, line = 0, cex.axis = 0.6, xpd = NA) for(i in 1:3){ nd.i <- nd.consp[nd.consp$ht.l == Hts[i],] lines(nd.i$consp.crowd.l, nd.i$pred_mort, lty = i, lwd = 2) } mtext("A)", adj = 0.05, line = 0) # mortality - heterospecific NDD plot(event ~ hetsp.crowd.l, data = dat1, xlab = "", ylab = "", type = 'n', ylim = c(0, Ymax), xaxt = "n", yaxt = 'n') axis(1, labels = NA) axis(2, labels = NA) top_hist <- hist(dat1$hetsp.crowd.l[dat1$event == 1], breaks = BREAKS_het, plot = F) low_hist <- hist(dat1$hetsp.crowd.l[dat1$event == 0], breaks = BREAKS_het, plot = F) rect(BREAKS_het[-20]+0.02, Ymax, BREAKS_het[-1]-0.02, Ymax - 0.3*top_hist$counts/6152, col= 'gray', border = NA) rect(BREAKS_het[-20]+0.02, 0, BREAKS_het[-1]-0.02, 0.3 * low_hist$counts/6152, col= 'gray', border = NA) axis(4, at = Ymax - 0.3*(0:2*2000)/6152, labels = 0:2*2000, xpd = NA, las = 1, line = 0, cex.axis = 0.6, xpd = NA) axis(4, at = 0.3*(0:2*2000)/6152, labels = 0:2*2000, xpd = NA, las = 1, line = 0, cex.axis = 0.6, xpd = NA) mtext(side = 4, at = Ymax/2, 'Number of seedlings dying or surviving\neach census interval', line = 2, cex = 0.8, srt = 90, xpd = NA) for(i in 1:3){ nd.i <- nd.hetsp[nd.hetsp$ht.l == Hts[i],] lines(nd.i$hetsp.crowd.l, nd.i$pred_mort, lty = i, lwd = 2) } mtext("B)", adj = 0.05, line = 0) # growth - conspecific NDD plot(rgr ~ consp.crowd.l, data = temp1, xlab = "", ylab = "", pch = 16, col = "#00000020", cex = 0.5, ylim = c(-0.3, 0.2)) title(xlab = "Conspecific neighbourhood crowding index", xpd = NA) title(ylab = expression(Relative~~growth~~rate~~(cm %.% cm %.% y^-1)), xpd = NA) for(i in 1:3){ nd.i <- nd.consp[nd.consp$ht.l == Hts[i],] lines(nd.i$consp.crowd.l, nd.i$pred_grow, lty = i, lwd = 2) } mtext("C)", adj = 0.05, line = 0) legend("top", lty = 1:3, col = "black", lwd = 2, title = "Seedling height (cm)", legend = c(10, 20, 50), ncol = 3, xpd = NA, inset = -0.16, bg = 'white') # growth - heterospecific NDD plot(rgr ~ hetsp.crowd.l, data = temp1, xlab = "", ylab = "", pch = 16, col = "#00000020", cex = 0.5, ylim = c(-0.3, 0.2), yaxt = 'n') axis(2, labels = NA) title(xlab = "Heterospecific neighbourhood crowding index", xpd = NA) for(i in 1:3){ nd.i <- nd.hetsp[nd.hetsp$ht.l == Hts[i],] lines(nd.i$hetsp.crowd.l, nd.i$pred_grow, lty = i, lwd = 2) } mtext("D)", adj = 0.05, line = 0) mtext("Figure 1", family = 'Times', outer = T, adj = 0.05) dev.off() #### Figure 2: Effects of taxon exclusion on MORTALITY #### dat2 <- dat[dat$Treatment %in% c("All Mammals", "Control", "Fungicide", "Insecticide", "Large Mammals"),] dat2$Treatment <- relevel(factor(dat2$Treatment), "Control") table(dat2$Treatment) contrasts(dat2$Treatment) contrasts(dat2$Treatment) <- (contrast.matrix1) contrasts(dat2$Treatment) dat2$consp.crowd.class <- round(dat2$consp.crowd.l) dat2$consp.crowd.class[dat2$consp.crowd.class == 4] <- 3 dat2$hetsp.crowd.class <- round(dat2$hetsp.crowd.l) dat2$hetsp.crowd.class[dat2$hetsp.crowd.class > 4] <- 4 dat3 <- dat[dat$Treatment %in% c("Control", "Fungicide", "Insecticide", "Fungicide & Insecticide"),] dat3$Insecticide_treatment <- ifelse(dat3$Treatment %in% c("Insecticide", "Fungicide & Insecticide"), "Insecticide", "Control") dat3$Fungicide_treatment <- ifelse(dat3$Treatment %in% c("Fungicide", "Fungicide & Insecticide"), "Fungicide", "Control") table(dat3$Insecticide_treatment, dat3$Fungicide_treatment) # nicely balanced dat3$consp.crowd.class <- round(dat3$consp.crowd.l) dat3$consp.crowd.class[dat3$consp.crowd.class == 4] <- 3 dat3$hetsp.crowd.class <- round(dat3$hetsp.crowd.l) dat3$hetsp.crowd.class[dat3$hetsp.crowd.class > 4] <- 4 dat4 <- dat[dat$Treatment %in% c("Control", "Fungicide", "All Mammals", "All Mammals & Fungicide"),] dat4$Mammal_treatment <- ifelse(dat4$Treatment %in% c("All Mammals", "All Mammals & Fungicide"), "All Mammals", "Control") dat4$Fungicide_treatment <- ifelse(dat4$Treatment %in% c("Fungicide", "All Mammals & Fungicide"), "Fungicide", "Control") table(dat4$Mammal_treatment, dat4$Fungicide_treatment) # nicely balanced dat4$consp.crowd.class <- round(dat4$consp.crowd.l) dat4$consp.crowd.class[dat4$consp.crowd.class == 4] <- 3 dat4$hetsp.crowd.class <- round(dat4$hetsp.crowd.l) dat4$hetsp.crowd.class[dat4$hetsp.crowd.class > 4] <- 4 dat5 <- dat[dat$Treatment %in% c("Control", "All Mammals", "Insecticide", "All Mammals & Insecticide"),] dat5$Mammal_treatment <- ifelse(dat5$Treatment %in% c("All Mammals", "All Mammals & Insecticide"), "All Mammals", "Control") dat5$Insecticide_treatment <- ifelse(dat5$Treatment %in% c("Insecticide", "All Mammals & Insecticide"), "Insecticide", "Control") table(dat5$Mammal_treatment, dat5$Insecticide_treatment) # nicely balanced dat5$consp.crowd.class <- round(dat5$consp.crowd.l) dat5$consp.crowd.class[dat5$consp.crowd.class == 4] <- 3 dat5$hetsp.crowd.class <- round(dat5$hetsp.crowd.l) dat5$hetsp.crowd.class[dat5$hetsp.crowd.class > 4] <- 4 # Observe the real patterns in the data COL2 <- COL[c(9, 1, 3, 5, 7)] COL3 <- COL[c(9, 3, 5, 4)] COL4 <- COL[c(1, 2, 9, 3)] COL5 <- COL[c(1, 6, 9, 5)] PCH <- c(1, 3, 4, 8) pdf("figs/NOT IN MANUSCRIPT - Real main effects on mortality.pdf", paper = 'a4', width = 8, height = 11, useDingbats = F) par(mfrow = c(4, 2), mar = c(3, 3 ,0, 0), oma = c(1, 1, 1, 1), xpd = NA, bty = "L", las = 1) means <- data.frame(tapply(dat2$event, list(dat2$consp.crowd.class, dat2$Treatment), mean)) matplot(means, type = "l", lty = 1, lwd = 2, ylab = "Change in mortality rate", col = COL2) legend('topleft', col = COL2, pch = PCH, legend = names(COL2), ncol = 2, lwd = 2) means <- data.frame(tapply(dat2$event, list(dat2$hetsp.crowd.class, dat2$Treatment), mean)) matplot(means, type = "l", lty = 1, lwd = 2, ylab = "", col = COL2) means <- aggregate(dat3$event, list(crowding = dat3$consp.crowd.class, dat3$Insecticide_treatment, dat3$Fungicide_treatment), mean, drop = F) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL3[trt], pch = PCH[trt], ylab = "Mortality rate", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL3[i], lwd = 2)} legend('topleft', col = COL3, pch = PCH, legend = names(COL3), ncol = 2, lwd = 2) means <- aggregate(dat3$event, list(crowding = dat3$hetsp.crowd.class, dat3$Insecticide_treatment, dat3$Fungicide_treatment), mean, drop = F) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL3[trt], pch = PCH[trt], ylab = "", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL3[i], lwd = 2)} means <- aggregate(dat4$event, list(crowding = dat4$consp.crowd.class, dat4$Mammal_treatment, dat4$Fungicide_treatment), mean, drop = F) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL4[trt], pch = PCH[trt], ylab = "Mortality rate", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL4[i], lwd = 2)} legend('topleft', col = COL4, pch = PCH, legend = names(COL4), ncol = 2, lwd = 2) means <- aggregate(dat4$event, list(crowding = dat4$hetsp.crowd.class, dat4$Mammal_treatment, dat4$Fungicide_treatment), mean, drop = F) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL4[trt], pch = PCH[trt], ylab = "", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL4[i], lwd = 2)} means <- aggregate(dat5$event, list(crowding = dat5$consp.crowd.class, dat5$Mammal_treatment, dat5$Insecticide_treatment), mean, drop = F) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL5[trt], pch = PCH[trt], ylab = "Mortality rate", xlab = "Conspecific neighborhood crowding index") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL5[i], lwd = 2)} legend('topleft', col = COL5, pch = PCH, legend = names(COL5), ncol = 2, lwd = 2) means <- aggregate(dat5$event, list(crowding = dat5$hetsp.crowd.class, dat5$Mammal_treatment, dat5$Insecticide_treatment), mean, drop = F) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL5[trt], pch = PCH[trt], ylab = "", xlab = "Heterospecific neighborhood crowding index") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL5[i], lwd = 2)} dev.off() # main effects m2.mort.0 <- glmer(event ~ (consp.crowd.l + hetsp.crowd.l)*Treatment + ht.l + (1|trans/plot) + (1|Census), offset= offset_time, family=binomial(link=cloglog), data=dat2, verbose = 1) m2.mort.1 <-update(m2.mort.0, . ~ . + (1|Name)) anova(m2.mort.0, m2.mort.1) # with the RE for species is MUCH better summary(m2.mort.1) anova(m2.mort.1) m2.mort.2a <- update(m2.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name)) m2.mort.2b <- update(m2.mort.1, . ~ . - (1|Name) + (hetsp.crowd.l|Name)) m2.mort.2c <- update(m2.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m2.mort.1, m2.mort.2a, m2.mort.2b, m2.mort.2c) # No evidence of different effects of crowding on different species in this dataset (less power to detect than in Figure 1). So use the simpler model. m2.mort.3a <- update(m2.mort.1, . ~ . + consp.crowd.l*ht.l) m2.mort.3b <- update(m2.mort.1, . ~ . + hetsp.crowd.l*ht.l) m2.mort.3c <- update(m2.mort.1, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m2.mort.1, m2.mort.3a, m2.mort.3b, m2.mort.3c) # Also, no sign that height affects responses to either type of crowding. Again, stick with the simpler model # significance tests of each predictor summary(m2.mort.1) r.squaredGLMM(m2.mort.1) # R2m R2c: theoretical 0.2718192 0.5779930, delta 0.1297976 0.2760001 pval_calc(m2.mort.1, nsim = 10000) # X.Intercept. consp.crowd.l hetsp.crowd.l # 0.0030 0.1347 0.2513 # TreatmentFungi TreatmentInsects TreatmentSmall.Mammals TreatmentLarge.Mammals # 0.0005 0.3760 0.1252 0.0429 # ht.l # 0.0000 # consp.crowd.l.TreatmentFungi consp.crowd.l.TrtInsects consp.crowd.l.TrtSmall.Mamm consp.crowd.l.TrtLarge.Mam # 0.0476 0.3085 0.1663 0.1289 # hetsp.crowd.l.TreatmentFungi hetsp.crowd.l.TrtInsects hetsp.crowd.l.TrtSmall.Mam hetsp.crowd.l.TrtLarge.Mammals # 0.0001 0.4715 0.1342 0.0461 nd_fungal_effect <- data.frame(ht.l = median(dat1$ht.l, na.rm = T), consp.crowd.l = log(10+1), hetsp.crowd.l = 0, Treatment = factor(levels(dat2$Treatment))) nd_fungal_effect$pred <- inverse_cloglog(predict(m2.mort.1, newdata = nd_fungal_effect, re.form = NA)) # Mortality increased from 27.1% to 32.1% with an order of magnitude increase in conspecific crowding # Insecticide * Fungicide interaction # m3.mort.0 <- glmer(event ~ (consp.crowd.l + hetsp.crowd.l)*Insecticide_treatment*Fungicide_treatment + ht.l + (1|trans/plot) + (1|Census), offset= offset_time, family=binomial(link=cloglog), data=dat3, verbose = 1) m3.mort.1 <-update(m3.mort.0, . ~ . + (1|Name)) anova(m3.mort.0, m3.mort.1) # with the RE for species is MUCH better summary(m3.mort.1) anova(m3.mort.1) m3.mort.2a <- update(m3.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name)) m3.mort.2b <- update(m3.mort.1, . ~ . - (1|Name) + (hetsp.crowd.l|Name)) m3.mort.2c <- update(m3.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m3.mort.1, m3.mort.2a, m3.mort.2b, m3.mort.2c) # No evidence of different effects of crowding on different species in this dataset (less power to detect than in Figure 1). So use the simpler model. m3.mort.3 <- update(m3.mort.1, . ~ . + consp.crowd.l*ht.l) m3.mort.3b <- update(m3.mort.1, . ~ . + hetsp.crowd.l*ht.l) m3.mort.3c <- update(m3.mort.1, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m3.mort.1, m3.mort.3, m3.mort.3b, m3.mort.3c) # Also, no sign that height affects responses to either type of crowding. Again, stick with the simpler model # significance tests of each predictor summary(m3.mort.1) r.squaredGLMM(m3.mort.1) # R2m R2c: theoretical 0.2284872 0.5170969, delta 0.0956022 0.2163605 pval_calc(m3.mort.1, nsim = 10000) # X.Intercept. # 0.1735 # consp.crowd.l hetsp.crowd.l # 0.1002 0.0016 # Insecticide_treatmentInsecticide Fungicide_treatmentFungicide # 0.1187 0.0003 # ht.l # 0.0000 # consp.crowd.l.Insecticide_treatmentInsecticide hetsp.crowd.l.Insecticide_treatmentInsecticide # 0.1187 0.0514 # consp.crowd.l.Fungicide_treatmentFungicide hetsp.crowd.l.Fungicide_treatmentFungicide # 0.0321 0.0000 # Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide # 0.0217 # consp.crowd.l.Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide # 0.3074 # hetsp.crowd.l.Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide # 0.0212 # Fungicide * Mammal interaction # m4.mort.0 <- glmer(event ~ (consp.crowd.l + hetsp.crowd.l)*Mammal_treatment*Fungicide_treatment + ht.l + (1|trans/plot) + (1|Census), offset= offset_time, family=binomial(link=cloglog), data=dat4, verbose = 1) m4.mort.1 <-update(m4.mort.0, . ~ . + (1|Name)) anova(m4.mort.0, m4.mort.1) # with the RE for species is MUCH better summary(m4.mort.1) anova(m4.mort.1) m4.mort.2a <- update(m4.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name)) m4.mort.2b <- update(m4.mort.1, . ~ . - (1|Name) + (hetsp.crowd.l|Name)) m4.mort.2c <- update(m4.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m4.mort.1, m4.mort.2a, m4.mort.2b, m4.mort.2c) # No evidence of different effects of crowding on different species in this dataset (less power to detect than in Figure 1). So use the simpler model. m4.mort.3a <- update(m4.mort.1, . ~ . + consp.crowd.l*ht.l) m4.mort.3b <- update(m4.mort.1, . ~ . + hetsp.crowd.l*ht.l) m4.mort.3c <- update(m4.mort.1, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m4.mort.1, m4.mort.3a, m4.mort.3b, m4.mort.3c) # Also, no sign that height affects responses to either type of crowding. Again, stick with the simpler model # significance tests of each predictor summary(m4.mort.1) r.squaredGLMM(m4.mort.1) # R2m R2c: theoretical 0.21716120 0.5220259, delta 0.09377857 0.2254309 pval_calc(m4.mort.1, nsim = 10000) # NEEDS MORE ITER! # X.Intercept. # 0.0936 # consp.crowd.l hetsp.crowd.l # 0.1320 0.0732 # Mammal_treatmentControl Fungicide_treatmentFungicide # 0.1611 0.1299 # ht.l # 0.0000 # consp.crowd.l.Mammal_treatmentControl hetsp.crowd.l.Mammal_treatmentControl # 0.4314 0.1063 # consp.crowd.l.Fungicide_treatmentFungicide hetsp.crowd.l.Fungicide_treatmentFungicide # 0.1828 0.2003 # Mammal_treatmentControl.Fungicide_treatmentFungicide # 0.0002 # consp.crowd.l.Mammal_treatmentControl.Fungicide_treatmentFungicide # 0.3248 # hetsp.crowd.l.Mammal_treatmentControl.Fungicide_treatmentFungicide # 0.0001 # Mammal * Insecticide interaction # m5.mort.0 <- glmer(event ~ (consp.crowd.l + hetsp.crowd.l)*Mammal_treatment*Insecticide_treatment + ht.l + (1|trans/plot) + (1|Census), offset= offset_time, family=binomial(link=cloglog), data=dat5, verbose = 1) m5.mort.1 <-update(m5.mort.0, . ~ . + (1|Name)) anova(m5.mort.0, m5.mort.1) # with the RE for species is MUCH better summary(m5.mort.1) anova(m5.mort.1) m5.mort.2a <- update(m5.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name)) m5.mort.2b <- update(m5.mort.1, . ~ . - (1|Name) + (hetsp.crowd.l|Name)) m5.mort.2c <- update(m5.mort.1, . ~ . - (1|Name) + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m5.mort.1, m5.mort.2a, m5.mort.2b, m5.mort.2c) # No evidence of different effects of crowding on different species in this dataset (less power to detect than in Figure 1). So use the simpler model. m5.mort.3a <- update(m5.mort.1, . ~ . + consp.crowd.l*ht.l) m5.mort.3b <- update(m5.mort.1, . ~ . + hetsp.crowd.l*ht.l) m5.mort.3c <- update(m5.mort.1, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m5.mort.1, m5.mort.3a, m5.mort.3b, m5.mort.3c) # Also, no sign that height affects responses to either type of crowding. Again, stick with the simpler model # significance tests of each predictor summary(m5.mort.1) r.squaredGLMM(m5.mort.1) # R2m R2c: theoretical 0.2461192 0.5198019, delta 0.1023426 0.2161468 pval_calc(m5.mort.1, nsim = 10000) # X.Intercept. # 0.1599 # consp.crowd.l hetsp.crowd.l # 0.2620 0.0933 # Mammal_treatmentControl Insecticide_treatmentInsecticide # 0.1627 0.3172 # ht.l # 0.0088 # consp.crowd.l.Mammal_treatmentControl hetsp.crowd.l.Mammal_treatmentControl # 0.3442 0.1207 # consp.crowd.l.Insecticide_treatmentInsecticide hetsp.crowd.l.Insecticide_treatmentInsecticide # 0.1139 0.2145 # Mammal_treatmentControl.Insecticide_treatmentInsecticide # 0.3012 # consp.crowd.l.Mammal_treatmentControl.Insecticide_treatmentInsecticide # 0.4080 # hetsp.crowd.l.Mammal_treatmentControl.Insecticide_treatmentInsecticide # 0.2510 # Prepare for graphing: Mortality and Growth affected by crowding and main effects # eff_mort2 <- mort_calculator(m2.mort.1, nsim = 1000) saveRDS(eff_mort2, file = "results/eff_mort2.rds") eff_mort2 <- readRDS("results/eff_mort2.rds") eff_mort3 <- interaction_mort_calculator(m3.mort.1, "Insecticide_treatment", "Fungicide_treatment", dat3, npred = 101, nsim = 1000) saveRDS(eff_mort3, file = "results/eff_mort3.rds") eff_mort3 <- readRDS("results/eff_mort3.rds") eff_mort4 <- interaction_mort_calculator(m4.mort.1, "Mammal_treatment", "Fungicide_treatment", dat4, npred = 101, nsim = 1000) saveRDS(eff_mort4, file = "results/eff_mort4.rds") eff_mort4 <- readRDS("results/eff_mort4.rds") eff_mort5 <- interaction_mort_calculator(m5.mort.1, "Mammal_treatment", "Insecticide_treatment", dat5, npred = 101, nsim = 1000) saveRDS(eff_mort5, file = "results/eff_mort5.rds") eff_mort5 <- readRDS("results/eff_mort5.rds") COL2 <- COL[c(3, 5, 7, 8)] COL3 <- COL[3:5] COL4 <- COL[3:1] COL5 <- COL[c(5, 6, 1)] COL_band2 <- COL_band[c(3, 5, 7, 8)] COL_band3 <- COL_band[3:5] COL_band4 <- COL_band[3:1] COL_band5 <- COL_band[c(5, 6, 1)] pdf("figs/Figure 2 - Effects on mortality.pdf", paper = 'a4', useDingbats = F, width = 8, height = 12) par(mfrow = c(4, 2), las = 1, bty = "L", mar = c(2, 2, 0.5, 0.5), oma = c(3, 3, 1, 1)) # mortality - conspecific NDD matplot(x = eff_mort2$mean$consp.crowd.l, y = eff_mort2$mean[,1:4], type = "n", ylab = "", xlab = "", ylim = c(-0.5, 0.5), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:4){ polygon(c(eff_mort2$mean$consp.crowd.l, rev(eff_mort2$mean$consp.crowd.l)), c(eff_mort2$CI_low[,i], rev(eff_mort2$CI_high[,i])), col = COL_band2[i], border = NA) } matlines(eff_mort2$mean$consp.crowd.l, eff_mort2$mean[,1:4], col = COL2, lwd = 2, lty = 1) mtext("A)", adj = 0.05, line = -1) legend("bottomleft", legend = c("Fungicide: 0.0476", "Insecticide: ns", "No large mammals: ns", "No small mammals: ns"), lwd = 2, col = COL2, cex = 0.8, ncol = 2, xpd = NA) # mortality - heterospecific NDD matplot(x = eff_mort2$mean$hetsp.crowd.l, y = eff_mort2$mean[,7:10], type = "n", ylab = '', xlab = "", ylim = c(-0.5, 0.5), axes = F) axis(1, labels = NA) axis(2, labels = NA) box() abline(h = 0, col = 'gray', lty = 2) for(i in 1:4){ polygon(c(eff_mort2$mean$hetsp.crowd.l, rev(eff_mort2$mean$hetsp.crowd.l)), c(eff_mort2$CI_low[,i+6], rev(eff_mort2$CI_high[,i+6])), col = COL_band2[i], border = NA) } matlines(eff_mort2$mean$hetsp.crowd.l, eff_mort2$mean[,7:10], col = COL2, lwd = 2, lty = 1) legend("bottom", legend = c("Fungicide: 0.0001", "Insecticide: ns", "No large mammals: ns", "No small mammals: 0.0461"), lwd = 2, col = COL2, cex = 0.8, ncol = 2, xpd = NA) mtext("B)", adj = 0.05, line = -1) # mortality - conspecific NDD matplot(x = eff_mort3$mean$consp.crowd.l, y = eff_mort3$mean[,1:3], type = "n", ylab = "", xlab = "", ylim = c(-0.5, 0.5), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_mort3$mean$consp.crowd.l, rev(eff_mort3$mean$consp.crowd.l)), c(eff_mort3$CI_low[,i], rev(eff_mort3$CI_high[,i])), col = COL_band3[i], border = NA) } matlines(eff_mort3$mean$consp.crowd.l, eff_mort3$mean[,1:3], col = COL3, lwd = 2, lty = 1) mtext("C)", adj = 0.05, line = -1) legend("bottomleft", legend = c("Fungicide: 0.0321", "F x I: ns", "Insecticide: ns"), lwd = 2, col = COL3, cex = 0.8, ncol = 1, xpd = NA) # mortality - heterospecific NDD matplot(x = eff_mort3$mean$hetsp.crowd.l, y = eff_mort3$mean[,6:8], type = "n", ylab = '', xlab = "", ylim = c(-0.5, 0.5), axes = F) axis(1, labels = NA) axis(2, labels = NA) box() abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_mort3$mean$hetsp.crowd.l, rev(eff_mort3$mean$hetsp.crowd.l)), c(eff_mort3$CI_low[,i+5], rev(eff_mort3$CI_high[,i+5])), col = COL_band3[i], border = NA) } matlines(eff_mort3$mean$hetsp.crowd.l, eff_mort3$mean[,6:8], col = COL3, lwd = 2, lty = 1) legend("bottom", legend = c("Fungicide: <0.0001", "F x I: 0.0212", "Insecticide: ns"), lwd = 2, col = COL3, cex = 0.8, ncol = 1, xpd = NA) mtext("D)", adj = 0.05, line = -1) matplot(x = eff_mort4$mean$consp.crowd.l, y = eff_mort4$mean[,1:3], type = "n", ylab = "", xlab = "", ylim = c(-0.5, 0.5), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_mort4$mean$consp.crowd.l, rev(eff_mort4$mean$consp.crowd.l)), c(eff_mort4$CI_low[,i], rev(eff_mort4$CI_high[,i])), col = COL_band4[i], border = NA) } matlines(eff_mort4$mean$consp.crowd.l, eff_mort4$mean[,1:3], col = COL4, lwd = 2, lty = 1) mtext("E)", adj = 0.05, line = -1) legend("bottomleft", legend = c("Fungicide: ns", "F x M: ns", "No mammals: ns"), lwd = 2, col = COL4, cex = 0.8, ncol = 1, xpd = NA) # mortality - heterospecific NDD matplot(x = eff_mort4$mean$hetsp.crowd.l, y = eff_mort4$mean[,6:8], type = "n", ylab = '', xlab = "", ylim = c(-0.5, 0.5), axes = F) axis(1, labels = NA) axis(2, labels = NA) box() abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_mort4$mean$hetsp.crowd.l, rev(eff_mort4$mean$hetsp.crowd.l)), c(eff_mort4$CI_low[,i+5], rev(eff_mort4$CI_high[,i+5])), col = COL_band4[i], border = NA) } matlines(eff_mort4$mean$hetsp.crowd.l, eff_mort4$mean[,6:8], col = COL4, lwd = 2, lty = 1) legend("bottom", legend = c("Fungicide: 0.1828", "F x M: 0.0001", "No mammals: ns"), lwd = 2, col = COL4, cex = 0.8, ncol = 1, xpd = NA) mtext("F)", adj = 0.05, line = -1) # mortality - conspecific NDD matplot(x = eff_mort5$mean$consp.crowd.l, y = eff_mort5$mean[,1:3], type = "n", ylab = "", xlab = "", ylim = c(-0.5, 0.5)) abline(h = 0, col = 'gray', lty = 2) title(xlab = "Conspecific neighborhood crowding index", xpd = NA) for(i in 1:3){ polygon(c(eff_mort5$mean$consp.crowd.l, rev(eff_mort5$mean$consp.crowd.l)), c(eff_mort5$CI_low[,i], rev(eff_mort5$CI_high[,i])), col = COL_band5[i], border = NA) } matlines(eff_mort5$mean$consp.crowd.l, eff_mort5$mean[,1:3], col = COL5, lwd = 2, lty = 1) legend("bottomleft", legend = c("Insecticide: ns", "I x M: ns", "No mammals: ns"), lwd = 2, col = COL5, cex = 0.8, ncol = 1, xpd = NA) mtext("G)", adj = 0.05, line = -1) matplot(x = eff_mort5$mean$hetsp.crowd.l, y = eff_mort5$mean[,6:8], type = "n", ylab = '', xlab = "", ylim = c(-0.5, 0.5), yaxt = 'n') axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) title(xlab = "Heterospecific neighborhood crowding index", xpd = NA) for(i in 1:3){ polygon(c(eff_mort5$mean$hetsp.crowd.l, rev(eff_mort5$mean$hetsp.crowd.l)), c(eff_mort5$CI_low[,i+5], rev(eff_mort5$CI_high[,i+5])), col = COL_band5[i], border = NA) } matlines(eff_mort5$mean$hetsp.crowd.l, eff_mort5$mean[,6:8], col = COL5, lwd = 2, lty = 1) legend("bottom", legend = c("Insecticide: ns", "I x M: ns", "No mammals: ns"), lwd = 2, col = COL5, cex = 0.8, ncol = 1, xpd = NA) mtext("H)", adj = 0.05, line = -1) title(ylab = expression(Effect~~of~~each~~taxon~~on~~probability~~of~~mortality~~(y^-1)), outer = T, line = -0.0) mtext("Figure 2", family = 'Times', outer = T, adj = 0.05) dev.off() #### Figure 3: Effects of taxon exclusion on GROWTH #### dat2 <- dat[dat$Treatment %in% c("All Mammals", "Control", "Fungicide", "Insecticide", "Large Mammals"),] dat2$Treatment <- relevel(factor(dat2$Treatment), "Control") table(dat2$Treatment) contrasts(dat2$Treatment) contrasts(dat2$Treatment) <- (contrast.matrix1) contrasts(dat2$Treatment) dat2$consp.crowd.class <- round(dat2$consp.crowd.l) dat2$consp.crowd.class[dat2$consp.crowd.class == 4] <- 3 dat2$hetsp.crowd.class <- round(dat2$hetsp.crowd.l) dat2$hetsp.crowd.class[dat2$hetsp.crowd.class > 4] <- 4 temp2 <- na.omit(dat2[dat2$rgr > -2 & dat2$rgr < 2 & !is.na(dat2$rgr), c("rgr", "consp.crowd.l", "hetsp.crowd.l", "trans", "plot", "Census", "Name", 'ht.l', 'Treatment', 'consp.crowd.class', 'hetsp.crowd.class'),]) # cut down dat to the key variables, and exclude observations where size changes unreasonably fast. Most of these are seedlings that either fell down, or had fallen down and now stood back up. dim(temp2) dat3 <- dat[dat$Treatment %in% c("Control", "Fungicide", "Insecticide", "Fungicide & Insecticide"),] dat3$Insecticide_treatment <- ifelse(dat3$Treatment %in% c("Insecticide", "Fungicide & Insecticide"), "Insecticide", "Control") dat3$Fungicide_treatment <- ifelse(dat3$Treatment %in% c("Fungicide", "Fungicide & Insecticide"), "Fungicide", "Control") table(dat3$Insecticide_treatment, dat3$Fungicide_treatment) # nicely balanced dat3$consp.crowd.class <- round(dat3$consp.crowd.l) dat3$consp.crowd.class[dat3$consp.crowd.class == 4] <- 3 dat3$hetsp.crowd.class <- round(dat3$hetsp.crowd.l) dat3$hetsp.crowd.class[dat3$hetsp.crowd.class > 4] <- 4 temp3 <- na.omit(dat3[dat3$rgr > -2 & dat3$rgr < 2 & !is.na(dat3$rgr), c("rgr", "consp.crowd.l", "hetsp.crowd.l", "trans", "plot", "Census", "Name", 'ht.l', 'Insecticide_treatment', 'Fungicide_treatment', 'consp.crowd.class', 'hetsp.crowd.class'),]) # cut down dat to the key variables, and exclude observations where size changes unreasonably fast. Most of these are seedlings that either fell down, or had fallen down and now stood back up. dim(temp3) dat4 <- dat[dat$Treatment %in% c("Control", "Fungicide", "All Mammals", "All Mammals & Fungicide"),] dat4$Mammal_treatment <- ifelse(dat4$Treatment %in% c("All Mammals", "All Mammals & Fungicide"), "All Mammals", "Control") dat4$Fungicide_treatment <- ifelse(dat4$Treatment %in% c("Fungicide", "All Mammals & Fungicide"), "Fungicide", "Control") table(dat4$Mammal_treatment, dat4$Fungicide_treatment) # nicely balanced dat4$consp.crowd.class <- round(dat4$consp.crowd.l) dat4$consp.crowd.class[dat4$consp.crowd.class == 4] <- 3 dat4$hetsp.crowd.class <- round(dat4$hetsp.crowd.l) dat4$hetsp.crowd.class[dat4$hetsp.crowd.class > 4] <- 4 temp4 <- na.omit(dat4[dat4$rgr > -2 & dat4$rgr < 2 & !is.na(dat4$rgr), c("rgr", "consp.crowd.l", "hetsp.crowd.l", "trans", "plot", "Census", "Name", 'ht.l', 'Mammal_treatment', 'Fungicide_treatment', 'consp.crowd.class', 'hetsp.crowd.class'),]) # cut down dat to the key variables, and exclude observations where size changes unreasonably fast. Most of these are seedlings that either fell down, or had fallen down and now stood back up. dim(temp4) dat5 <- dat[dat$Treatment %in% c("Control", "All Mammals", "Insecticide", "All Mammals & Insecticide"),] dat5$Mammal_treatment <- ifelse(dat5$Treatment %in% c("All Mammals", "All Mammals & Insecticide"), "All Mammals", "Control") dat5$Insecticide_treatment <- ifelse(dat5$Treatment %in% c("Insecticide", "All Mammals & Insecticide"), "Insecticide", "Control") table(dat5$Mammal_treatment, dat5$Insecticide_treatment) # nicely balanced dat5$consp.crowd.class <- round(dat5$consp.crowd.l) dat5$consp.crowd.class[dat5$consp.crowd.class == 4] <- 3 dat5$hetsp.crowd.class <- round(dat5$hetsp.crowd.l) dat5$hetsp.crowd.class[dat5$hetsp.crowd.class > 4] <- 4 temp5 <- na.omit(dat5[dat5$rgr > -2 & dat5$rgr < 2 & !is.na(dat5$rgr), c("rgr", "consp.crowd.l", "hetsp.crowd.l", "trans", "plot", "Census", "Name", 'ht.l', 'Mammal_treatment', 'Insecticide_treatment', 'consp.crowd.class', 'hetsp.crowd.class'),]) # cut down dat to the key variables, and exclude observations where size changes unreasonably fast. Most of these are seedlings that either fell down, or had fallen down and now stood back up. dim(temp5) # Observe the real patterns in the data pdf("figs/NOT IN PAPER - Real effects on growth.pdf", paper = 'a4', height = 12, width = 8, useDingbats = F) par(mfrow = c(4, 2), mar = c(3, 3 ,0, 0), oma = c(1, 1, 1, 1), xpd = NA, bty = "L", las = 1) COL2 <- COL[c(9, 1, 3, 5, 7)] COL3 <- COL[c(9, 3, 5, 4)] COL4 <- COL[c(1, 2, 9, 3)] COL5 <- COL[c(1, 6, 9, 5)] PCH <- c(3, 8, 1, 4) means <- data.frame(tapply(temp2$rgr, list(temp2$consp.crowd.class, temp2$Treatment), mean, na.rm = T)) matplot(means, type = "l", lty = 1, lwd = 2, ylab = "", xlab = "", col = COL2) legend("topleft", legend = names(COL2), lwd = 1, col = COL2, ncol = 2, cex = 0.8) means <- data.frame(tapply(temp2$rgr, list(temp2$hetsp.crowd.class, temp2$Treatment), mean)) matplot(means, type = "l", lty = 1, lwd = 2, ylab = "", xlab = "", col = COL2) means <- aggregate(temp3$rgr, list(crowding = temp3$consp.crowd.class, temp3$Insecticide_treatment, temp3$Fungicide_treatment), mean, drop = F, na.rm = T) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL3[trt], pch = PCH[trt], ylab = "", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL3[i], lwd = 2)} legend("topleft", legend = names(COL3), lwd = 1, col = COL3, ncol = 2, cex = 0.8) means <- aggregate(temp3$rgr, list(crowding = temp3$hetsp.crowd.class, temp3$Insecticide_treatment, temp3$Fungicide_treatment), mean, drop = F, na.rm = T) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL3[trt], pch = PCH[trt], ylab = "", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL3[i], lwd = 2)} # It does not look like either treatment affected the impacts of heterospecific crowding. Insecticide seems to have INCREASED mortality rates at high levels of conspeficic crowding (which is not predicted), and fungicide INCREASED rgr at high conspecific crowding (which is expected) means <- aggregate(temp4$rgr, list(crowding = temp4$consp.crowd.class, temp4$Mammal_treatment, temp4$Fungicide_treatment), mean, drop = F, na.rm = T) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL4[trt], pch = PCH[trt], ylab = "", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL4[i], lwd = 2)} legend("topleft", legend = names(COL4), lwd = 1, col = COL4, ncol = 2, cex = 0.8) means <- aggregate(temp4$rgr, list(crowding = temp4$hetsp.crowd.class, temp4$Mammal_treatment, temp4$Fungicide_treatment), mean, drop = F, na.rm = T) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL4[trt], pch = PCH[trt], ylab = "", xlab = "") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL4[i], lwd = 2)} # It does not look like either treatment affected the impacts of heterospecific crowding. Insecticide seems to have INCREASED mortality rates at high levels of conspeficic crowding (which is not predicted), and fungicide INCREASED rgr at high conspecific crowding (which is expected) means <- aggregate(temp5$rgr, list(crowding = temp5$consp.crowd.class, temp5$Mammal_treatment, temp5$Insecticide_treatment), mean, drop = F, na.rm = T) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL5[trt], pch = PCH[trt], ylab = "", xlab = "Conspecific crowding") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL5[i], lwd = 2)} legend("topleft", legend = names(COL5), lwd = 1, col = COL5, ncol = 2, cex = 0.8) means <- aggregate(temp5$rgr, list(crowding = temp5$hetsp.crowd.class, temp5$Mammal_treatment, temp5$Insecticide_treatment), mean, drop = F, na.rm = T) means$trt <- factor(paste(means$Group.2, means$Group.3)) plot(x ~ crowding, data = means, type = 'p', col = COL5[trt], pch = PCH[trt], ylab = "", xlab = "Heterospecific crowding") for(i in 1:4){lines(x ~ crowding, data = means[as.numeric(means$trt) == i,], col = COL5[i], lwd = 2)} # It does not look like either treatment affected the impacts of heterospecific crowding. Insecticide seems to have INCREASED mortality rates at high levels of conspeficic crowding (which is not predicted), and fungicide INCREASED rgr at high conspecific crowding (which is expected) title(ylab = "Relative growth rate (cm/cm/yr)", outer = T, xpd = NA, line = -0.5) dev.off() # Main effects # m2.grow.0 <- lmer(rgr ~ (consp.crowd.l + hetsp.crowd.l)*Treatment + (1|trans/plot) + (1|Census), temp2, REML=FALSE, verbose = 1) summary(grow.0) m2.grow.1 <- update(m2.grow.0, . ~ . + (1|Name)) m2.grow.2a <- update(m2.grow.0, . ~ . + (consp.crowd.l|Name)) m2.grow.2b <- update(m2.grow.0, . ~ . + (hetsp.crowd.l|Name)) m2.grow.2c <- update(m2.grow.0, . ~ . + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m2.grow.0, m2.grow.1, m2.grow.2a, m2.grow.2b, m2.grow.2c) # NO evidence for diff slopes among spp for either type of crowding on rgr summary(m2.grow.1) m2.grow.3a <- update(m2.grow.1, . ~ . + consp.crowd.l*ht.l) m2.grow.3b <- update(m2.grow.1, . ~ . + hetsp.crowd.l*ht.l) m2.grow.3c <- update(m2.grow.1, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m2.grow.1, m2.grow.3a, m2.grow.3b, m2.grow.3c) # Height strongly affects responses to crowding from conspecifics m2.grow.3a.ss <- getME(m2.grow.3a, c("theta","fixef")) m2.grow.3a <- update(m2.grow.3a, start=grow.3a.ss, control=lmerControl(optCtrl=list(maxfun=2e4))) summary(m2.grow.3a) r.squaredGLMM(m2.grow.3a) # R2m R2c: 0.04461328 0.1018713 pval_calc(m2.grow.3a, nsim = 10000) # X.Intercept. consp.crowd.l # 0.0017 0.0537 # hetsp.crowd.l TreatmentFungi # 0.1602 0.0890 # TreatmentInsects TreatmentSmall.Mammals # 0.2908 0.4683 # TreatmentLarge.Mammals ht.l # 0.4344 0.0037 # consp.crowd.l.TreatmentFungi consp.crowd.l.TreatmentInsects # 0.0172 0.0000 # consp.crowd.l.TreatmentSmall.Mammals consp.crowd.l.TreatmentLarge.Mammals # 0.1953 0.0615 # hetsp.crowd.l.TreatmentFungi hetsp.crowd.l.TreatmentInsects # 0.1283 0.1940 # hetsp.crowd.l.TreatmentSmall.Mammals hetsp.crowd.l.TreatmentLarge.Mammals # 0.2542 0.2642 # consp.crowd.l.ht.l # 0.0514 # Fungicide * Insecticide interactions # m3.grow.0 <- lmer(rgr ~ (consp.crowd.l + hetsp.crowd.l)*Insecticide_treatment*Fungicide_treatment + (1|trans/plot) + (1|Census), temp3, REML=FALSE, verbose = 1) summary(m3.grow.0) m3.grow.1 <- update(m3.grow.0, . ~ . + (1|Name)) m3.grow.2a <- update(m3.grow.0, . ~ . + (consp.crowd.l|Name)) m3.grow.2b <- update(m3.grow.0, . ~ . + (hetsp.crowd.l|Name)) m3.grow.2c <- update(m3.grow.0, . ~ . + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m3.grow.0, m3.grow.1, m3.grow.2a, m3.grow.2b, m3.grow.2c) # NO evidence for diff intercepts or slopes among spp for either type of crowding on rgr summary(m3.grow.0) m3.grow.3 <- update(m3.grow.0, . ~ . + consp.crowd.l*ht.l) m3.grow.3b <- update(m3.grow.0, . ~ . + hetsp.crowd.l*ht.l) m3.grow.3c <- update(m3.grow.0, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m3.grow.0, m3.grow.3, m3.grow.3b, m3.grow.3c) # No sign that height affects responses to crowding from conspecifics m3.grow.0.ss <- getME(m3.grow.0, c("theta","fixef")) m3.grow.0 <- update(m3.grow.0, start=m3.grow.0.ss, control=lmerControl(optCtrl=list(maxfun=2e4))) summary(m3.grow.0) r.squaredGLMM(m3.grow.0) # R2m R2c: 0.03257394 0.05741377 pval_calc(m3.grow.0, nsim = 10000) # X.Intercept. # 0.0603 # consp.crowd.l hetsp.crowd.l # 0.3526 0.3136 # Insecticide_treatmentInsecticide Fungicide_treatmentFungicide # 0.3775 0.0474 # consp.crowd.l.Insecticide_treatmentInsecticide hetsp.crowd.l.Insecticide_treatmentInsecticide # 0.0118 0.4765 # consp.crowd.l.Fungicide_treatmentFungicide hetsp.crowd.l.Fungicide_treatmentFungicide # 0.1113 0.0694 # Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide # 0.1827 # consp.crowd.l.Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide # 0.0091 # hetsp.crowd.l.Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide # 0.4534 # Fungicide * Mammal interaction # m4.grow.0 <- lmer(rgr ~ (consp.crowd.l + hetsp.crowd.l)*Mammal_treatment*Fungicide_treatment + (1|trans/plot) + (1|Census), temp4, REML=FALSE, verbose = 1) summary(m4.grow.0) m4.grow.1 <- update(m4.grow.0, . ~ . + (1|Name)) m4.grow.2a <- update(m4.grow.0, . ~ . + (consp.crowd.l|Name)) m4.grow.2b <- update(m4.grow.0, . ~ . + (hetsp.crowd.l|Name)) m4.grow.2c <- update(m4.grow.0, . ~ . + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m4.grow.0, m4.grow.1, m4.grow.2a, m4.grow.2b, m4.grow.2c) # NO evidence for diff intercepts or slopes among spp for either type of crowding on rgr summary(m4.grow.0) m4.grow.3a <- update(m4.grow.0, . ~ . + consp.crowd.l*ht.l) m4.grow.3b <- update(m4.grow.0, . ~ . + hetsp.crowd.l*ht.l) m4.grow.3c <- update(m4.grow.0, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m4.grow.0, m4.grow.3a, m4.grow.3b, m4.grow.3c) # No height affects responses to crowding from conspecifics ONLY m4.grow.3a.ss <- getME(m4.grow.3a, c("theta","fixef")) m4.grow.3a <- update(m4.grow.3a, start=m4.grow.3a.ss, control=lmerControl(optCtrl=list(maxfun=2e4))) summary(m4.grow.3a) r.squaredGLMM(m4.grow.3a) # R2m R2c: 0.03377294 0.04704284 pval_calc(m4.grow.3a, nsim = 10000) # X.Intercept. # 0.0078 # consp.crowd.l hetsp.crowd.l # 0.3696 0.4854 # Mammal_treatmentControl Fungicide_treatmentFungicide # 0.2162 0.0311 # ht.l # 0.0058 # consp.crowd.l.Mammal_treatmentControl hetsp.crowd.l.Mammal_treatmentControl # 0.1337 0.1225 # consp.crowd.l.Fungicide_treatmentFungicide hetsp.crowd.l.Fungicide_treatmentFungicide # 0.0838 0.0800 # Mammal_treatmentControl.Fungicide_treatmentFungicide # 0.3584 # consp.crowd.l.ht.l # 0.1797 # consp.crowd.l.Mammal_treatmentControl.Fungicide_treatmentFungicide # 0.0250 # hetsp.crowd.l.Mammal_treatmentControl.Fungicide_treatmentFungicide # 0.4268 # Mammal * Insecticide interaction # m5.grow.0 <- lmer(rgr ~ (consp.crowd.l + hetsp.crowd.l)*Mammal_treatment*Insecticide_treatment + (1|trans/plot) + (1|Census), temp5, REML=FALSE, verbose = 1) summary(m5.grow.0) m5.grow.1 <- update(m5.grow.0, . ~ . + (1|Name)) m5.grow.2a <- update(m5.grow.0, . ~ . + (consp.crowd.l|Name)) m5.grow.2b <- update(m5.grow.0, . ~ . + (hetsp.crowd.l|Name)) m5.grow.2c <- update(m5.grow.0, . ~ . + (consp.crowd.l|Name) + (hetsp.crowd.l|Name)) anova(m5.grow.0, m5.grow.1, m5.grow.2a, m5.grow.2b, m5.grow.2c) # NO evidence for diff intercepts or slopes among spp for either type of crowding on rgr summary(m5.grow.0) m5.grow.3a <- update(m5.grow.0, . ~ . + consp.crowd.l*ht.l) m5.grow.3b <- update(m5.grow.0, . ~ . + hetsp.crowd.l*ht.l) m5.grow.3c <- update(m5.grow.0, . ~ . + consp.crowd.l*ht.l + hetsp.crowd.l*ht.l) anova(m5.grow.0, m5.grow.3a, m5.grow.3b, m5.grow.3c) # No height affects responses to crowding from conspecifics ONLY m5.grow.3a.ss <- getME(m5.grow.3a, c("theta","fixef")) m5.grow.3a <- update(m5.grow.3a, start=m5.grow.3a.ss, control=lmerControl(optCtrl=list(maxfun=2e4))) summary(m5.grow.3a) r.squaredGLMM(m5.grow.3a) # R2m R2c: 0.0289217 0.06202766 pval_calc(m5.grow.3a, nsim = 10000) # X.Intercept. # 0.0092 # consp.crowd.l hetsp.crowd.l # 0.4277 0.4052 # Mammal_treatmentControl Insecticide_treatmentInsecticide # 0.2505 0.4142 # ht.l # 0.0073 # consp.crowd.l.Mammal_treatmentControl hetsp.crowd.l.Mammal_treatmentControl # 0.1360 0.1625 # consp.crowd.l.Insecticide_treatmentInsecticide hetsp.crowd.l.Insecticide_treatmentInsecticide # 0.0731 0.2670 # Mammal_treatmentControl.Insecticide_treatmentInsecticide # 0.3371 # consp.crowd.l.ht.l # 0.3767 # consp.crowd.l.Mammal_treatmentControl.Insecticide_treatmentInsecticide # 0.4680 # hetsp.crowd.l.Mammal_treatmentControl.Insecticide_treatmentInsecticide # 0.2714 # Prepare for graphing: Growth affected by crowding and main effects # eff_grow2 <- grow_calculator(m2.grow.3a, nsim = 10000) saveRDS(eff_grow2, "results/eff_grow2.rds") eff_grow2 <- readRDS("results/eff_grow2.rds") eff_grow3 <- interaction_grow_calculator(m3.grow.0, "Insecticide_treatment", "Fungicide_treatment", temp3, npred = 101, nsim = 10000) saveRDS(eff_grow3, "results/eff_grow3.rds") eff_grow3 <- readRDS("results/eff_grow3.rds") eff_grow4 <- interaction_grow_calculator(m4.grow.3a, "Mammal_treatment", "Fungicide_treatment", temp4, npred = 101, nsim = 10000) saveRDS(eff_grow4, "results/eff_grow4.rds") eff_grow4 <- readRDS("results/eff_grow4.rds") eff_grow5 <- interaction_grow_calculator(m5.grow.3a, "Mammal_treatment", "Insecticide_treatment", temp5, npred = 101, nsim = 10000) saveRDS(eff_grow5, "results/eff_grow5.rds") eff_grow5 <- readRDS("results/eff_grow5.rds") COL2 <- COL[c(3, 5, 7, 8)] COL3 <- COL[3:5] COL4 <- COL[3:1] COL5 <- COL[c(5, 6, 1)] COL_band2 <- COL_band[c(3, 5, 7, 8)] COL_band3 <- COL_band[3:5] COL_band4 <- COL_band[3:1] COL_band5 <- COL_band[c(5, 6, 1)] pdf("figs/Figure 3 - Effects on growth.pdf", paper = 'a4', height = 11, width = 8, useDingbats = F) par(mfrow = c(4, 2), las = 1, bty = "L", mar = c(2, 2, 0.5, 0.5), oma = c(3, 3, 1, 1)) # Growth - conspecific NDD matplot(x = eff_grow2$mean$consp.crowd.l, y = eff_grow2$mean[,1:4], type = "n", ylab = "", xlab = "", ylim = c(-0.5, 0.5), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:4){ polygon(c(eff_grow2$mean$consp.crowd.l, rev(eff_grow2$mean$consp.crowd.l)), c(eff_grow2$CI_low[,i], rev(eff_grow2$CI_high[,i])), col = COL_band2[i], border = NA) } matlines(eff_grow2$mean$consp.crowd.l, eff_grow2$mean[,1:4], col = COL2, lwd = 2, lty = 1) legend("top", legend = c("Fungicide: 0.0172", "Insecticide: <0.0001", "Large mammals: ns", "Small mammals: ns"), lwd = 2, col = COL2, ncol = 2, cex = 0.8) mtext("A)", adj = 0.05, line = -1) # Growth - heterospecific NDD matplot(x = eff_grow2$mean$hetsp.crowd.l, y = eff_grow2$mean[,7:10], type = "n", ylab = '', xlab = "", ylim = c(-0.5, 0.5), axes = F) axis(1, labels = NA) axis(2, labels = NA) box() abline(h = 0, col = 'gray', lty = 2) for(i in 1:4){ polygon(c(eff_grow2$mean$hetsp.crowd.l, rev(eff_grow2$mean$hetsp.crowd.l)), c(eff_grow2$CI_low[,i+6], rev(eff_grow2$CI_high[,i+6])), col = COL_band2[i], border = NA) } matlines(eff_grow2$mean$hetsp.crowd.l, eff_grow2$mean[,7:10], col = COL2, lwd = 2, lty = 1) legend("top", legend = c("Fungicide: ns", "Insecticide: ns", "Large mammals: ns", "Small mammals: ns"), lwd = 2, col = COL2, ncol = 2, cex = 0.8) mtext("B)", adj = 0.05, line = -1) # Growth - conspecific NDD matplot(x = eff_grow3$mean$consp.crowd.l, y = eff_grow3$mean[,1:3], type = "n", ylab = "", xlab = "", ylim = c(-0.5, 0.5), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_grow3$mean$consp.crowd.l, rev(eff_grow3$mean$consp.crowd.l)), c(eff_grow3$CI_low[,i], rev(eff_grow3$CI_high[,i])), col = COL_band3[i], border = NA) } matlines(eff_grow3$mean$consp.crowd.l, eff_grow3$mean[,1:3], col = COL3, lwd = 2, lty = 1) legend("top", legend = c("Fungicide: ns", "F x I: 0.0091", "Insecticide: ns"), lwd = 2, col = COL3, ncol = 1, cex = 0.8) mtext("C)", adj = 0.05, line = -1) # Growth - heterospecific NDD matplot(x = eff_grow3$mean$hetsp.crowd.l, y = eff_grow3$mean[,6:8], type = "n", ylab = '', xlab = "", ylim =c(-0.5, 0.5), axes = F) axis(1, labels = NA) axis(2, labels = NA) box() abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_grow3$mean$hetsp.crowd.l, rev(eff_grow3$mean$hetsp.crowd.l)), c(eff_grow3$CI_low[,i+5], rev(eff_grow3$CI_high[,i+5])), col = COL_band3[i], border = NA) } matlines(eff_grow3$mean$hetsp.crowd.l, eff_grow3$mean[,6:8], col = COL3, lwd = 2, lty = 1) legend("top", legend = c("Fungicide: ns", "F x I: ns", "Insecticide: ns"), lwd = 2, col = COL3, ncol = 1, cex = 0.8) mtext("D)", adj = 0.05, line = -1) # Growth - conspecific NDD matplot(x = eff_grow4$mean$consp.crowd.l, y = eff_grow4$mean[,1:3], type = "n",ylab = "", xlab = "", ylim = c(-0.5, 0.5), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_grow4$mean$consp.crowd.l, rev(eff_grow4$mean$consp.crowd.l)), c(eff_grow4$CI_low[,i], rev(eff_grow4$CI_high[,i])), col = COL_band4[i], border = NA) } matlines(eff_grow4$mean$consp.crowd.l, eff_grow4$mean[,1:3], col = COL4, lwd = 2, lty = 1) legend("top", legend = c("Fungicide: ns", "F x M: 0.0250", "All mammals: ns"), lwd = 2, col = COL4, ncol = 1, cex = 0.8) mtext("E)", adj = 0.05, line = -1) # Growth - heterospecific NDD matplot(x = eff_grow4$mean$hetsp.crowd.l, y = eff_grow4$mean[,6:8], type = "n", ylab = '', xlab = "", ylim = c(-0.5, 0.5), axes = F) axis(1, labels = NA) axis(2, labels = NA) box() abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_grow4$mean$hetsp.crowd.l, rev(eff_grow4$mean$hetsp.crowd.l)), c(eff_grow4$CI_low[,i+5], rev(eff_grow4$CI_high[,i+5])), col = COL_band4[i], border = NA) } matlines(eff_grow4$mean$hetsp.crowd.l, eff_grow4$mean[,6:8], col = COL4, lwd = 2, lty = 1) legend("top", legend = c("Fungicide: ns", "F x M: ns", "All mammals: ns"), lwd = 2, col = COL4, ncol = 1, cex = 0.8) mtext("F)", adj = 0.05, line = -1) matplot(x = eff_grow5$mean$consp.crowd.l, y = eff_grow5$mean[,1:3], type = "l=n", ylab = "", xlab = "", ylim = c(-0.5, 0.5)) title(xlab = "Conspecific neighborhood crowding index", xpd = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_grow5$mean$consp.crowd.l, rev(eff_grow5$mean$consp.crowd.l)), c(eff_grow5$CI_low[,i], rev(eff_grow5$CI_high[,i])), col = COL_band5[i], border = NA) } matlines(eff_grow5$mean$consp.crowd.l, eff_grow5$mean[,1:3], col = COL5, lwd = 2, lty = 1) legend("top", legend = c("Insecticide: ns", "M x I: ns", "All mammals: ns"), lwd = 2, col = COL5, ncol = 1, cex = 0.8) mtext("G)", adj = 0.05, line = -1) # Growth - heterospecific NDD matplot(x = eff_grow5$mean$hetsp.crowd.l, y = eff_grow5$mean[,6:8], type = "n", ylab = '', xlab = "", ylim = c(-0.5, 0.5), yaxt = 'n') title(xlab = "Heterospecific neighborhood crowding index", xpd = NA) axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ polygon(c(eff_grow5$mean$hetsp.crowd.l, rev(eff_grow5$mean$hetsp.crowd.l)), c(eff_grow5$CI_low[,i+5], rev(eff_grow5$CI_high[,i+5])), col = COL_band5[i], border = NA) } matlines(eff_grow5$mean$hetsp.crowd.l, eff_grow5$mean[,6:8], col = COL5, lwd = 2, lty = 1) legend("top", legend = c("Insecticide: ns", "M x I: ns", "All mammals: ns"), lwd = 2, col = COL5, ncol = 1, cex = 0.8) mtext("H)", adj = 0.05, line = -1) title(ylab = expression(Effect~~of~~each~~taxon~~on~~relative~~growth~~rate~~(cm%.%cm^-1%.%y^-1)), outer = T, line = 0.5) mtext("Figure 3", family = 'Times', outer = T, adj = 0.05) dev.off() #### Figure 4 - Changes in Shannon diversity #### table(divsdl$Treatment) levels(divsdl$Treatment) <- levels(divsdl$Treatment)[c(1:8, 4)] table(divsdl$Treatment) temp <- tapply(as.Date(dat$obs.date), dat$plot.census, mean) divsdl$obs.date <- as.Date("1970-01-01") + temp[match(divsdl$plot.census, names(temp))] temp <- tapply(divsdl$obs.date, divsdl$Census, mean, na.rm = T) divsdl$obs.date.simple <- as.Date("1970-01-01") + temp[match(divsdl$Census, names(temp))] hist(divsdl$obs.date.simple, breaks = 1000) divsdl <- data.frame(divsdl, matrix(as.numeric(unlist(strsplit(as.character(divsdl$location), "_"))), ncol = 2, byrow = T)) names(divsdl)[24:25] <- c("trans", 'plot') divsdl$trans <- factor(divsdl$trans) divsdl$plot <- factor(divsdl$plot) head(divsdl) means <- tapply(divsdl$diversity.shannon.exp, list(divsdl$obs.date.simple, divsdl$Treatment), mean, na.rm = T) sds <- tapply(divsdl$diversity.shannon.exp, list(divsdl$obs.date.simple, divsdl$Treatment), sd, na.rm = T) ns <- tapply(divsdl$diversity.shannon.exp, list(divsdl$obs.date.simple, divsdl$Treatment), function(x){sum(!is.na(x))}) ses <- sds/sqrt(ns) xpos <- outer(as.Date(row.names(means)), 10*1:ncol(means), "+") matplot(x = xpos, y = means, lty = 1, type = 'l') for(i in 1:ncol(means)){segments(xpos[,i], means[,i] - ses[,i], xpos[,i], means[,i] + ses[,i], col = i)} # rescale the diversity for each plot, so that at Census 9, every plot has diversity 100. plots <- sort(unique(divsdl$location[divsdl$Census %in% 8:9])) for(i in 1:length(plots)){ divsdl.i <- divsdl[divsdl$location == plots[i],] if(9 %in% divsdl.i$Census){ divsdl$diversity.shannon.exp.scaled[divsdl$location == plots[i]] <- 100 * divsdl$diversity.shannon.exp[divsdl$location == plots[i]]/divsdl$diversity.shannon.exp[divsdl$location == plots[i] & divsdl$Census == 9] } else { divsdl$diversity.shannon.exp.scaled[divsdl$location == plots[i]] <- 100 * divsdl$diversity.shannon.exp[divsdl$location == plots[i]]/divsdl$diversity.shannon.exp[divsdl$location == plots[i] & divsdl$Census == 8] } } summary( divsdl$diversity.shannon.exp.scaled) # THe NAs are from plost that were not censused at Census 9 or 8, mostly because they were crushed under a treefall or something like that. means <- data.frame(tapply(divsdl$diversity.shannon.exp.scaled, list(divsdl$Census, divsdl$Treatment), mean, na.rm = T)) sds <- data.frame(tapply(divsdl$diversity.shannon.exp.scaled, list(divsdl$Census, divsdl$Treatment), sd, na.rm = T)) ns <- data.frame(tapply(divsdl$diversity.shannon.exp.scaled, list(divsdl$Census, divsdl$Treatment), function(x){sum(!is.na(x))})) ses <- sds/sqrt(ns) xpos <- outer(as.numeric(row.names(means)), -3.5:3.5/40, "+") names(means) <- levels(divsdl$Treatment) matplot(x = xpos, y = means, lty = 1, type = 'l') for(i in 1:ncol(means)){segments(xpos[,i], means[,i] - ses[,i], xpos[,i], means[,i] + ses[,i], col = i)} mains <- c('All Mammals', "Fungicide", "Insecticide", 'Large Mammals', 'Control') dim(divsdl) divsdl1 <- divsdl[divsdl$Census > 8 & divsdl$Treatment %in% c('All Mammals', "Fungicide", "Insecticide", 'Large Mammals', 'Control'),] divsdl1$Treatment <- factor(divsdl1$Treatment) contrasts(divsdl1$Treatment) contrasts(divsdl1$Treatment) <- contrast.matrix1 contrasts(divsdl1$Treatment) divsdl1$Census <- factor(divsdl1$Census) dim(divsdl1) divsdl2 <- divsdl[divsdl$Census > 8 & divsdl$Treatment %in% c("Control", "Fungicide", "Insecticide", "Fungicide & Insecticide"),] divsdl2$Census <- factor(divsdl2$Census) divsdl2$Insecticide_treatment <- factor(ifelse(divsdl2$Treatment %in% c("Insecticide", "Fungicide & Insecticide"), "Insecticide", "Control")) divsdl2$Fungicide_treatment <- factor(ifelse(divsdl2$Treatment %in% c("Fungicide", "Fungicide & Insecticide"), "Fungicide", "Control")) table(divsdl2$Insecticide_treatment, divsdl2$Fungicide_treatment) divsdl3 <- divsdl[divsdl$Census > 8 & divsdl$Treatment %in% c("Control", "Fungicide", "All Mammals", "All Mammals & Fungicide"),] divsdl3$Census <- factor(divsdl3$Census) divsdl3$Mammal_treatment <- factor(ifelse(divsdl3$Treatment %in% c("All Mammals", "All Mammals & Fungicide"), "All Mammals", "Control")) divsdl3$Fungicide_treatment <- factor(ifelse(divsdl3$Treatment %in% c("Fungicide", "All Mammals & Fungicide"), "Fungicide", "Control")) table(divsdl3$Mammal_treatment, divsdl3$Fungicide_treatment) divsdl4 <- divsdl[divsdl$Census > 8 & divsdl$Treatment %in% c("Control", "Insecticide", "All Mammals", "All Mammals & Insecticide"),] divsdl4$Census <- factor(divsdl4$Census) divsdl4$Mammal_treatment <- factor(ifelse(divsdl4$Treatment %in% c("All Mammals", "All Mammals & Insecticide"), "All Mammals", "Control")) divsdl4$Insecticide_treatment <- factor(ifelse(divsdl4$Treatment %in% c("Insecticide", "All Mammals & Insecticide"), "Insecticide", "Control")) table(divsdl4$Mammal_treatment, divsdl4$Insecticide_treatment) COL_main <- COL[names(COL) %in% mains] COL_main <- COL_main[c(1, 5, 2:4)] means_main <- means[,names(means) %in% mains] ses_main <- ses[,names(means) %in% mains] xpos_main <-outer(1:nrow(means_main), -2:2/40, "+") par(mfrow = c(1, 1), las = 1, bty = "L") matplot(x = xpos_main, y = means_main, lty = 1, type = 'l', col= COL_main, lwd = 2, ylab = "Shannon diversity (scaled so that census 8 or 9 = 100)") legend('top', ncol = 2, legend = names(COL_main), col = COL_main, lwd = 2) m1.div.0 <- lmer(diversity.shannon.exp ~ Treatment*Census + (1|location), divsdl1, REML=FALSE, verbose = 1) summary(m1.div.0) anova(m1.div.0) r.squaredGLMM(m1.div.0) # R2m R2c: 0.04890621 0.8612194 pval1 <- pval_calc(m1.div.0, nsim = 10000) # X.Intercept. TreatmentFungi TreatmentInsects # 0.0000 0.0379 0.4876 # TreatmentSmall.Mammals TreatmentLarge.Mammals Census10 # 0.2183 0.3574 0.0001 # Census11 TreatmentFungi.Census10 TreatmentInsects.Census10 # 0.0000 0.2702 0.4183 # TreatmentSmall.Mammals.Census10 TreatmentLarge.Mammals.Census10 TreatmentFungi.Census11 # 0.2753 0.4431 0.0027 # TreatmentInsects.Census11 TreatmentSmall.Mammals.Census11 TreatmentLarge.Mammals.Census11 # 0.4768 0.0039 0.4653 m2.div.0 <- lmer(diversity.shannon.exp ~ Insecticide_treatment*Fungicide_treatment*Census + (1|location), divsdl2, REML=FALSE, verbose = 1) summary(m2.div.0) r.squaredGLMM(m2.div.0) # R2m R2c: 0.01665264 0.602411 pval2 <- pval_calc(m2.div.0, nsim = 10000) # X.Intercept. # 0.0000 # Insecticide_treatmentInsecticide # 0.3166 # Fungicide_treatmentFungicide # 0.0168 # Census10 # 0.0105 # Census11 # 0.2793 # Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide # 0.1573 # Insecticide_treatmentInsecticide.Census10 # 0.3393 # Insecticide_treatmentInsecticide.Census11 # 0.0358 # Fungicide_treatmentFungicide.Census10 # 0.4612 # Fungicide_treatmentFungicide.Census11 # 0.4900 # Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide.Census10 # 0.3725 # Insecticide_treatmentInsecticide.Fungicide_treatmentFungicide.Census11 # 0.3820 m3.div.0 <- lmer(diversity.shannon.exp ~ Mammal_treatment*Fungicide_treatment*Census + (1|location), divsdl3, REML=FALSE, verbose = 1) summary(m3.div.0) r.squaredGLMM(m3.div.0) # R2m R2c: 0.02556795 0.6063025 pval3 <- pval_calc(m3.div.0, nsim = 10000) # X.Intercept. # 0.0000 # Mammal_treatmentControl # 0.2564 # Fungicide_treatmentFungicide # 0.1137 # Census10 # 0.0130 # Census11 # 0.0000 # Mammal_treatmentControl.Fungicide_treatmentFungicide # 0.3789 # Mammal_treatmentControl.Census10 # 0.1586 # Mammal_treatmentControl.Census11 # 0.0007 # Fungicide_treatmentFungicide.Census10 # 0.4558 # Fungicide_treatmentFungicide.Census11 # 0.2371 # Mammal_treatmentControl.Fungicide_treatmentFungicide.Census10 # 0.4376 # Mammal_treatmentControl.Fungicide_treatmentFungicide.Census11 # 0.3031 m4.div.0 <- lmer(diversity.shannon.exp ~ Insecticide_treatment*Mammal_treatment*Census + (1|location), divsdl4, REML=FALSE, verbose = 1) summary(m4.div.0) r.squaredGLMM(m4.div.0) # R2m R2c: 0.02083459 0.6112499 pval4 <- pval_calc(m4.div.0, nsim = 10000) # X.Intercept. # 0.0000 # Insecticide_treatmentInsecticide # 0.1428 # Mammal_treatmentControl # 0.2618 # Census10 # 0.0064 # Census11 # 0.0000 # Insecticide_treatmentInsecticide.Mammal_treatmentControl # 0.2762 # Insecticide_treatmentInsecticide.Census10 # 0.3448 # Insecticide_treatmentInsecticide.Census11 # 0.2450 # Mammal_treatmentControl.Census10 # 0.1358 # Mammal_treatmentControl.Census11 # 0.0001 # Insecticide_treatmentInsecticide.Mammal_treatmentControl.Census10 # 0.4864 # Insecticide_treatmentInsecticide.Mammal_treatmentControl.Census11 # 0.2013 eff_div1 <- div_calculator(m1.div.0, nsim = 10000) COL_div1 <- COL[c(3, 5, 7, 8)] eff_div1$Pval # Applying fungicide significantly REDUCES diversity. Insecticide has no effect. Large mammals and Small mamals do, but only at Census 11, and they both INCREASE diversity, whcih is hared to explain. eff_div1$obs.date <- divsdl1$obs.date.simple[match(eff_div1$Census, divsdl1$Census)] eff_div1$xpos <- eff_div1$obs.date + rep((-1.5:1.5)*10, each = 3) saveRDS(eff_div1, "results/eff_div1.rds") eff_div1 <- readRDS("results/eff_div1.rds") eff_div2 <- div_interaction_calculator(m2.div.0, trt1 = "Insecticide_treatment", trt2 = "Fungicide_treatment", dat_pred = divsdl2, nsim = 10000) COL_div2 <- COL[5:3] eff_div2$obs.date <- divsdl2$obs.date.simple[match(eff_div2$Census, divsdl2$Census)] eff_div2$xpos <- eff_div2$obs.date+rep((-1:1)*10, each = 3) eff_div2$Pval # BOth trts significantly REDUCE diversity, but fungicide is stronger. the interaction is simply between them (no sig interaction) saveRDS(eff_div2, "results/eff_div2.rds") eff_div2 <- readRDS("results/eff_div2.rds") eff_div3 <- div_interaction_calculator(m3.div.0, trt1 = "Mammal_treatment", trt2 = "Fungicide_treatment", dat_pred = divsdl3, nsim = 10000) COL_div3 <- COL[1:3] eff_div3$obs.date <- divsdl3$obs.date.simple[match(eff_div3$Census, divsdl3$Census)] eff_div3$xpos <- eff_div3$obs.date+rep((-1:1)*10, each = 3) eff_div3$Pval # Applying fungicide significantly REDUCES diversity. Insecticide has no effectNor does large mammals, and Small mammals saveRDS(eff_div3, "results/eff_div3.rds") eff_div3 <- readRDS("results/eff_div3.rds") eff_div4 <- div_interaction_calculator(m4.div.0, trt1 = "Insecticide_treatment", trt2 = "Mammal_treatment", dat_pred = divsdl4, nsim = 10000) COL_div4 <- COL[c(5,6,1)] eff_div4$obs.date <- divsdl4$obs.date.simple[match(eff_div4$Census, divsdl4$Census)] eff_div4$xpos <- eff_div4$obs.date+rep((-1:1)*10, each = 3) eff_div4$Pval # Applying fungicide significantly REDUCES diversity. Insecticide has no effectNor does large mammals, and Small mamals saveRDS(eff_div4, "results/eff_div4.rds") eff_div4 <- readRDS("results/eff_div4.rds") # Now plot up the effects on diversity pdf("figs/Figure 4 = Effects on diversity.pdf", paper = 'a4', useDingbats = F) YLIM <- c(-6, 3) par(mfrow = c(2, 2), las = 1, bty = "L", mar = c(2, 2, 0.5, 0.5), oma = c(3, 3, 1, 1)) Trts <- unique(eff_div1$Treatment) plot(mean ~ xpos, data = eff_div1, col = rep(COL_div1, each = 3), xlab = "", ylab = "", ylim = YLIM, lwd = 1, pch = ifelse(Pval<0.05, 16, 1), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:4){ eff_div1.i <- eff_div1[eff_div1$Treatment == Trts[i],] lines(eff_div1.i$xpos, eff_div1.i$mean, col = COL_div1[i]) segments(eff_div1.i$xpos, eff_div1.i$CI_low, eff_div1.i$xpos, eff_div1.i$CI_high, col = COL_div1[i]) } legend('bottomleft', ncol = 2, cex = 0.8, lwd = 2, legend = names(COL_div1), col = COL_div1, inset = 0.0) mtext("A)", adj = 0.05, line = -1) Trts <- unique(eff_div2$Treatment) plot(mean ~ xpos, data = eff_div2, col = rep(COL_div2, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1, xaxt = 'n', yaxt = 'n') axis(1, labels = NA) axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_div2.i <- eff_div2[eff_div2$Treatment == Trts[i],] lines(eff_div2.i$xpos, eff_div2.i$mean, col = COL_div2[i]) segments(eff_div2.i$xpos, eff_div2.i$CI_low, eff_div2.i$xpos, eff_div2.i$CI_high, col = COL_div2[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div2), col = COL_div2, inset = 0.00) mtext("B)", adj = 0.05, line = -1) Trts <- unique(eff_div3$Treatment) plot(mean ~ xpos, data = eff_div3, col = rep(COL_div3, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_div3.i <- eff_div3[eff_div3$Treatment == Trts[i],] lines(eff_div3.i$xpos, eff_div3.i$mean, col = COL_div3[i]) segments(eff_div3.i$xpos, eff_div3.i$CI_low, eff_div3.i$xpos, eff_div3.i$CI_high, col = COL_div3[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div3), col = COL_div3, inset = 0.00) mtext("C)", adj = 0.05, line = -1) Trts <- unique(eff_div4$Treatment) plot(mean ~ xpos, data = eff_div4, col = rep(COL_div4, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1, yaxt = 'n') axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_div4.i <- eff_div4[eff_div4$Treatment == Trts[i],] lines(eff_div4.i$xpos, eff_div4.i$mean, col = COL_div4[i]) segments(eff_div4.i$xpos, eff_div4.i$CI_low, eff_div4.i$xpos, eff_div4.i$CI_high, col = COL_div4[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div4), col = COL_div4, inset = 0.00) mtext("D)", adj = 0.05, line = -1) mtext("Figure 4", family = 'Times', outer = T, adj = 0.05) title (ylab = expression(Effect~~of~~each~~taxon~~on~~Shannon~~diversity~~(e^("H'"))), xlab = "Observation date", outer = T, xpd = NA, line = 0) dev.off() #### Supplemental Figure 1 - Effects on richness & evenness #### m1.even.0 <- lmer(evenness ~ Treatment*Census + (1|location), divsdl1, REML=FALSE, verbose = 1) summary(m1.even.0) anova(m1.even.0) m2.even.0 <- lmer(evenness ~ Insecticide_treatment*Fungicide_treatment*Census + (1|location), divsdl2, REML=FALSE, verbose = 1) summary(m2.even.0) m3.even.0 <- lmer(evenness ~ Mammal_treatment*Fungicide_treatment*Census + (1|location), divsdl3, REML=FALSE, verbose = 1) summary(m3.even.0) m4.even.0 <- lmer(evenness ~ Insecticide_treatment*Mammal_treatment*Census + (1|location), divsdl4, REML=FALSE, verbose = 1) summary(m4.even.0) eff_even1 <- div_calculator(m1.even.0, nsim = 10000) COL_even1 <- COL[c(3, 5, 7, 8)] eff_even1$obs.date <- divsdl1$obs.date.simple[match(eff_even1$Census, divsdl1$Census)] eff_even1$xpos <- eff_even1$obs.date + rep((-1.5:1.5)*10, each = 3) saveRDS(eff_even1, "results/eff_even1.rds") eff_even1 <- readRDS("results/eff_even1.rds") eff_even2 <- div_interaction_calculator(m2.even.0, trt1 = "Insecticide_treatment", trt2 = "Fungicide_treatment", dat_pred = divsdl2, nsim = 10000) COL_div2 <- COL[5:3] eff_even2$obs.date <- divsdl2$obs.date.simple[match(eff_even2$Census, divsdl2$Census)] eff_even2$xpos <- eff_even2$obs.date+rep((-1:1)*10, each = 3) saveRDS(eff_even2, "results/eff_even2.rds") eff_even2 <- readRDS("results/eff_even2.rds") eff_even3 <- div_interaction_calculator(m3.even.0, trt1 = "Mammal_treatment", trt2 = "Fungicide_treatment", dat_pred = divsdl3, nsim = 10000) COL_div3 <- COL[1:3] eff_even3$obs.date <- divsdl3$obs.date.simple[match(eff_even3$Census, divsdl3$Census)] eff_even3$xpos <- eff_even3$obs.date+rep((-1:1)*10, each = 3) saveRDS(eff_even3, "results/eff_even3.rds") eff_even3 <- readRDS("results/eff_even3.rds") eff_even4 <- div_interaction_calculator(m4.even.0, trt1 = "Insecticide_treatment", trt2 = "Mammal_treatment", dat_pred = divsdl4, nsim = 10000) COL_div4 <- COL[c(5,6,1)] eff_even4$obs.date <- divsdl4$obs.date.simple[match(eff_even4$Census, divsdl4$Census)] eff_even4$xpos <- eff_even4$obs.date+rep((-1:1)*10, each = 3) saveRDS(eff_even4, "results/eff_even4.rds") eff_even4 <- readRDS("results/eff_even4.rds") ## Effects on richness m1.rich.0 <- lmer(richness ~ Treatment*Census + (1|location), divsdl1, REML=FALSE, verbose = 1) summary(m1.rich.0) anova(m1.rich.0) m2.rich.0 <- lmer(richness ~ Insecticide_treatment*Fungicide_treatment*Census + (1|location), divsdl2, REML=FALSE, verbose = 1) summary(m2.rich.0) m3.rich.0 <- lmer(richness ~ Mammal_treatment*Fungicide_treatment*Census + (1|location), divsdl3, REML=FALSE, verbose = 1) summary(m3.rich.0) m4.rich.0 <- lmer(richness ~ Insecticide_treatment*Mammal_treatment*Census + (1|location), divsdl4, REML=FALSE, verbose = 1) eff_rich1 <- div_calculator(m1.rich.0, nsim = 10000) COL_rich1 <- COL[c(3, 5, 7, 8)] eff_rich1$obs.date <- divsdl1$obs.date.simple[match(eff_rich1$Census, divsdl1$Census)] eff_rich1$xpos <- eff_rich1$obs.date + rep((-1.5:1.5)*10, each = 3) saveRDS(eff_rich1, "results/eff_rich1.rds") eff_rich1 <- readRDS("results/eff_rich1.rds") eff_rich2 <- div_interaction_calculator(m2.rich.0, trt1 = "Insecticide_treatment", trt2 = "Fungicide_treatment", dat_pred = divsdl2, nsim = 10000) COL_div2 <- COL[5:3] eff_rich2$obs.date <- divsdl2$obs.date.simple[match(eff_rich2$Census, divsdl2$Census)] eff_rich2$xpos <- eff_rich2$obs.date+rep((-1:1)*10, each = 3) saveRDS(eff_rich2, "results/eff_rich2.rds") eff_rich2 <- readRDS("results/eff_rich2.rds") eff_rich3 <- div_interaction_calculator(m3.rich.0, trt1 = "Mammal_treatment", trt2 = "Fungicide_treatment", dat_pred = divsdl3, nsim = 10000) COL_div3 <- COL[1:3] eff_rich3$obs.date <- divsdl3$obs.date.simple[match(eff_rich3$Census, divsdl3$Census)] eff_rich3$xpos <- eff_rich3$obs.date+rep((-1:1)*10, each = 3) saveRDS(eff_rich3, "results/eff_rich3.rds") eff_rich3 <- readRDS("results/eff_rich3.rds") eff_rich4 <- div_interaction_calculator(m4.rich.0, trt1 = "Insecticide_treatment", trt2 = "Mammal_treatment", dat_pred = divsdl4, nsim = 10000) COL_div4 <- COL[c(5,6,1)] eff_rich4$obs.date <- divsdl4$obs.date.simple[match(eff_rich4$Census, divsdl4$Census)] eff_rich4$xpos <- eff_rich4$obs.date+rep((-1:1)*10, each = 3) saveRDS(eff_rich4, "results/eff_rich4.rds") eff_rich4 <- readRDS("results/eff_rich4.rds") pdf("figs/Supplemental Figure 1 - Effects on richness & evenness.pdf", paper = 'a4', useDingbats = F, width = 8, height = 10) YLIM <- c(-8, 5) par(mfrow = c(4, 2), las = 1, bty = "L", mar = c(2, 2, 0.5, 0.5), oma = c(3, 3, 1, 1)) Trts <- unique(eff_rich1$Treatment) plot(mean ~ xpos, data = eff_rich1, col = rep(COL_div1, each = 3), xlab = "", ylab = "", ylim = YLIM, lwd = 1, pch = ifelse(Pval<0.05, 16, 1), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:4){ eff_rich1.i <- eff_rich1[eff_rich1$Treatment == Trts[i],] lines(eff_rich1.i$xpos, eff_rich1.i$mean, col = COL_rich1[i]) segments(eff_rich1.i$xpos, eff_rich1.i$CI_low, eff_rich1.i$xpos, eff_rich1.i$CI_high, col = COL_rich1[i]) } legend('bottomleft', ncol = 2, cex = 0.8, lwd = 2, legend = names(COL_div1), col = COL_div1, inset = 0.0) mtext("A)", adj = 0.05, line = -1) Trts <- unique(eff_rich2$Treatment) plot(mean ~ xpos, data = eff_rich2, col = rep(COL_div2, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1, xaxt = 'n', yaxt = 'n') axis(1, labels = NA) axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_rich2.i <- eff_rich2[eff_rich2$Treatment == Trts[i],] lines(eff_rich2.i$xpos, eff_rich2.i$mean, col = COL_div2[i]) segments(eff_rich2.i$xpos, eff_rich2.i$CI_low, eff_rich2.i$xpos, eff_rich2.i$CI_high, col = COL_div2[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div2), col = COL_div2, inset = 0.00) mtext("B)", adj = 0.05, line = -1) Trts <- unique(eff_rich3$Treatment) plot(mean ~ xpos, data = eff_rich3, col = rep(COL_div3, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_rich3.i <- eff_rich3[eff_rich3$Treatment == Trts[i],] lines(eff_rich3.i$xpos, eff_rich3.i$mean, col = COL_div3[i]) segments(eff_rich3.i$xpos, eff_rich3.i$CI_low, eff_rich3.i$xpos, eff_rich3.i$CI_high, col = COL_div3[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div3), col = COL_div3, inset = 0.00) mtext("C)", adj = 0.05, line = -1) Trts <- unique(eff_rich4$Treatment) plot(mean ~ xpos, data = eff_rich4, col = rep(COL_div4, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1, yaxt = 'n') axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_rich4.i <- eff_rich4[eff_rich4$Treatment == Trts[i],] lines(eff_rich4.i$xpos, eff_rich4.i$mean, col = COL_div4[i]) segments(eff_rich4.i$xpos, eff_rich4.i$CI_low, eff_rich4.i$xpos, eff_rich4.i$CI_high, col = COL_div4[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div4), col = COL_div4, inset = 0.00) mtext("D)", adj = 0.05, line = -1) mtext(side = 2, text = "Effect of each taxon on richness", outer = T, xpd = NA, line = 0, adj = 0.9, las = 3) YLIM <- c(-0.1, 0.1) Trts <- unique(eff_even1$Treatment) plot(mean ~ xpos, data = eff_even1, col = rep(COL_div1, each = 3), xlab = "", ylab = "", ylim = YLIM, lwd = 1, pch = ifelse(Pval<0.05, 16, 1), xaxt = 'n') axis(1, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:4){ eff_even1.i <- eff_even1[eff_even1$Treatment == Trts[i],] lines(eff_even1.i$xpos, eff_even1.i$mean, col = COL_even1[i]) segments(eff_even1.i$xpos, eff_even1.i$CI_low, eff_even1.i$xpos, eff_even1.i$CI_high, col = COL_even1[i]) } legend('bottomleft', ncol = 2, cex = 0.8, lwd = 2, legend = names(COL_div1), col = COL_div1, inset = 0.0) mtext("E)", adj = 0.05, line = -1) Trts <- unique(eff_even2$Treatment) plot(mean ~ xpos, data = eff_even2, col = rep(COL_div2, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1, xaxt = 'n', yaxt = 'n') axis(1, labels = NA) axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_even2.i <- eff_even2[eff_even2$Treatment == Trts[i],] lines(eff_even2.i$xpos, eff_even2.i$mean, col = COL_div2[i]) segments(eff_even2.i$xpos, eff_even2.i$CI_low, eff_even2.i$xpos, eff_even2.i$CI_high, col = COL_div2[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div2), col = COL_div2, inset = 0.00) mtext("F)", adj = 0.05, line = -1) Trts <- unique(eff_even3$Treatment) plot(mean ~ xpos, data = eff_even3, col = rep(COL_div3, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_even3.i <- eff_even3[eff_even3$Treatment == Trts[i],] lines(eff_even3.i$xpos, eff_even3.i$mean, col = COL_div3[i]) segments(eff_even3.i$xpos, eff_even3.i$CI_low, eff_even3.i$xpos, eff_even3.i$CI_high, col = COL_div3[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div3), col = COL_div3, inset = 0.00) mtext("G)", adj = 0.05, line = -1) Trts <- unique(eff_even4$Treatment) plot(mean ~ xpos, data = eff_even4, col = rep(COL_div4, each = 3), pch = ifelse(Pval<0.05, 16, 1), xlab = "", ylab = "", ylim = YLIM, lwd = 1, yaxt = 'n') axis(2, labels = NA) abline(h = 0, col = 'gray', lty = 2) for(i in 1:3){ eff_even4.i <- eff_even4[eff_even4$Treatment == Trts[i],] lines(eff_even4.i$xpos, eff_even4.i$mean, col = COL_div4[i]) segments(eff_even4.i$xpos, eff_even4.i$CI_low, eff_even4.i$xpos, eff_even4.i$CI_high, col = COL_div4[i]) } legend('bottomleft', ncol = 1, cex = 0.8, lwd = 2, legend = names(COL_div4), col = COL_div4, inset = 0.00) mtext("H)", adj = 0.05, line = -1) mtext("Supplemental Figure 1", family = 'Times', outer = T, adj = 0.05) mtext(side = 2, text = "Effect of each taxon on evenness", outer = T, xpd = NA, line = 0, adj = 0.1, las = 3) mtext(side = 1, text = "Observation date", outer = T, xpd = NA, line = 0) dev.off()