#### Description: This script runs the forward model selection approach advocated by Blanchet et al. (2008, Ecology) to identify which willow traits best predict herbivore community responses. After identifying these traits, I then determined whether willow genotype still had a significant effect. These results are reported in Table 3 of the manuscript entitled, "Multiple plant traits shape the genetic basis of herbivore community assembly." # Reference: Blanchet, F. G., Legendre, P. & Borcard, D. (2008) Forward selection of explanatory variables. Ecology 89, 2623–2632. ## load required scripts and libraries for analyses source('~/Documents/Genotype_Willow_Community/datasets_&_Rscripts/arthropods_&_traits/herb_trait_manage.r') library(rockchalk) # calculate delta R2 for dropping variables from the model. library(vegan) library(MASS) library(car) ## Function for generating multiple regression data tables (mr_table) mr_table <- function(name, model){ require(rockchalk) data.frame(Response = c(name, rep("", length(get(paste("model"))$coefficients)-2)), Variable = names(get(paste("model"))$coefficients)[-1], coef.mean = round(summary(get(paste("model")))$coefficients[-1,-3][,1],2), coef.se = round(summary(get(paste("model")))$coefficients[-1,-3][,2],3), P = round(summary(get(paste("model")))$coefficients[-1,-3][,3],3), partial.r2 = round(getDeltaRsquare(get(paste("model"))),3)) } ## Note that for all univariate responses (e.g., herbivore abundance, species and guild responses), I manually implemented Blanchet et al.'s (2008, Ecology) approach, whereas for the community composition analyses, I used the built in function 'ordiR2step' in the package 'vegan'. I also coded out all residual plots to make the script easier to run. ## Herbivore abundance model herb.null <- lm(herb.com.trait$log_herbabund ~ 1, herb.traits) herb.full <- lm(herb.com.trait$log_herbabund ~ ., herb.traits) summary(herb.full) # adj R2 = 0.3071 stepAIC(herb.null, scope=formula(herb.full), direction = "forward") herb.update <- update(herb.null, .~. + log_size + flavonOLES.PC2 + log_trich) # adj R2 = 0.1992, 0.2415 summary(herb.update) #plot(herb.update) herb.mr.tab <- mr_table(name = "log(Herbivore abundance)", model = herb.update) herb.update.Geno <- update(herb.update, .~. + herb.com.trait$Genotype) #to examine whether Genotype still has an effect anova(herb.update.Geno) #Genotype has only a marginally significant effect. ## determining whether my approach to reduce multicollinearity was effect. vif(herb.full) # variance inflation factors for each variable in the full trait model were all less than 1.8. Note that although I only ran this for 'herbivore abundance', VIF only depends on the predictor variables, so this holds for all analyses. ## Herbivore community composition model. Currently, I'm thinking its better to include all of the traits and interpret the marginal contributions of each variable. rda0 <- rda(decostand(herb.com.trait[ ,herbivore.names], method="norm") ~ 1, herb.traits) rda.trait <- rda(decostand(herb.com.trait[ ,herbivore.names], method = "norm") ~ ., herb.traits) RsquareAdj(rda.trait) anova(rda.trait, by="margin") rda.R2step <- ordiR2step(rda0, scope=formula(rda.trait)) RsquareAdj(rda.R2step) plot(rda.R2step, display=c("sp","bp")) rda.R2step.anova <- anova(rda.R2step, by="margin", step = 1000) anova(rda.R2step, step = 1000) rda.R2step.tab <- data.frame(Response = c("Community composition",rep("",4)), Variable = c("cinn.PC1","cinn.PC2","log_size","flavanonOLES.PC1","flavonOLES.PC2"), coef.mean = rep("",5), coef.se = rep("",5), P = rda.R2step.anova$"Pr(>F)"[-6], partial.r2 = round((rda.R2step.anova$Var/sum(rda.R2step.anova$Var))[-6],3)) rda.R2step.geno <- update(rda.R2step, .~. + herb.com.trait$Genotype) anova(rda.R2step.geno, by="term", step = 1000) # Genotype still has a significant effect. ## Leaf damage model. ld.null <- lm(log(herb.com.trait$avg.LeafDam.score) ~ 1, herb.traits) ld.full <- lm(log(herb.com.trait$avg.LeafDam.score) ~ ., herb.traits) summary(ld.full) # full model not significant #plot(ld.full) ## Gall abundance model. gall.null <- lm(herb.com.trait$log_gall.abund ~ 1, herb.traits) gall.full <- lm(herb.com.trait$log_gall.abund ~ ., herb.traits) plot(gall.full) summary(gall.full) # Global model not significant. ## Leaf mine abundance model leafmine.null <- lm(herb.com.trait$log_leafminer.abund ~ 1, herb.traits) leafmine.full <- lm(herb.com.trait$log_leafminer.abund ~ ., herb.traits) summary(leafmine.full) # adj R2 = 0.3201 stepAIC(leafmine.null, scope = formula(leafmine.full), direction="forward") leafmine.update <- update(leafmine.null, .~. + log_size + density_resid + cinn.PC2 + cinn.PC1 + height_resid) summary(leafmine.update) # adj R2 = 0.1133, 0.1921, 0.2362, 0.2768, 0.3028 #plot(leafmine.update) leafmine.mr.tab <- mr_table(name = "log(leaf miner abundance+1)", model = leafmine.update) # is it really log(leaf miner + 1)? leafmine.update.Geno <- update(leafmine.update, .~. + herb.com.trait$Genotype) anova(leafmine.update.Geno) # marginally sig. effect #plot(leafmine.update.Geno) ## Leaf chewer abundance model. leafchew.null <- lm(herb.com.trait$log_leafchew.abund ~ 1, herb.traits) leafchew.full <- lm(herb.com.trait$log_leafchew.abund ~ ., herb.traits) summary(leafchew.full) # adj R2 = 0.1968 stepAIC(leafchew.null, scope = formula(leafchew.full), direction="forward") leafchew.update <- update(leafchew.null, .~. + log_size) summary(leafchew.update) # sal_tannin.PC1 not included in the model. #leafchew.mr.tab <- mr_table(name = "log(leaf chewer abundance = 1)", model = leafchew.update) # not working because there is only one significant trait effect. leafchew.update.Geno <- update(leafchew.update, .~. + herb.com.trait$Genotype) anova(leafchew.update.Geno) ## Xylem sucker abundance model. These results are not reported in the manuscript, because xylem sucker abundance only varied marginally among willow genotypes. xylem.null <- lm(herb.com.trait$log_xylem.abund ~ 1, herb.traits) xylem.full <- lm(herb.com.trait$log_xylem.abund ~ ., herb.traits) summary(xylem.full) # adj R2 = 0.1299 stepAIC(xylem.null, scope = formula(xylem.full), direction="forward") xylem.update <- update(xylem.null, .~. + log_size + flavonOLES.PC2) summary(xylem.update) # adj R2 = 0.0666, 0.1275 #plot(xylem.update) xylem.mr.tab <- mr_table(name = "log(Xylem sucker abundance+1)", model = xylem.update) xylem.update.Geno <- update(xylem.update, .~. + herb.com.trait$Genotype) anova(xylem.update.Geno) #plot(xylem.update.Geno) # Phloem sucker abundance model. phloem.null = lm(herb.com.trait$log_phloem.abund ~ 1, herb.traits) phloem.full = lm(herb.com.trait$log_phloem.abund ~ ., herb.traits) summary(phloem.full) # adj R2 = 0.1869 stepAIC(phloem.null, scope = formula(phloem.full), direction = "forward") phloem.update = update(phloem.null, .~. + log_size + sla_resid + cinn.PC1 + height_resid) summary(phloem.update) #plot(phloem.update) # only a little wonky phloem.mr.tab <- mr_table(name = "log(Phloem sucker abundance+1)", model = phloem.update) phloem.update.Geno <- update(phloem.update, .~. + herb.com.trait$Genotype) anova(phloem.update.Geno) ## generate table herb.trait_reg_table <- rbind.data.frame(herb.mr.tab, rda.R2step.tab, xylem.mr.tab, phloem.mr.tab, leafmine.mr.tab) # Need to add leafchewer results manually. gall abundance global model was not significant rownames(herb.trait_reg_table) <- NULL ## Tachyerges salicicis (Curc1) abundance curc1.null <- lm(herb.com.trait$log_curc1 ~ 1, herb.traits) curc1.full <- lm(herb.com.trait$log_curc1 ~ ., herb.traits) summary(curc1.full) # adj. R2 = 0.2022 curc1.step <- stepAIC(curc1.null, scope = formula(curc1.full), direction = "forward") curc1.update <- update(curc1.null, .~. + cinn.PC1 + density_resid + log_size + water_content) # flavonOLES.PC2 not included (P > 0.05) summary(curc1.update) #plot(curc1.update) curc1.mr.tab <- mr_table(name = "log(Tachyerges salicis + 1)", model = curc1.update) curc1.update.Geno <- update(curc1.update, .~. + herb.com.trait$Genotype) #to examine whether Genotype still has an effect anova(curc1.update.Geno) # P = 0.04960 ## Caloptilia (LTF) abundance model LTF.null <- lm(herb.com.trait$sqrt_LTF ~ 1, herb.traits) LTF.full <- lm(herb.com.trait$sqrt_LTF ~ ., herb.traits) summary(LTF.full) # adj. R2 = 0.2319 LTF.step <- stepAIC(LTF.null, scope = formula(LTF.full), direction = "forward") LTF.update <- update(LTF.null, .~. + log_size + flavonOLES.PC1) # border line, but excluding flavanonOLES.PC1 because including brings adj R2 = 0.2361 summary(LTF.update) #plot(LTF.update) LTF.mr.tab <- mr_table(name = "sqrt(LTF)", model = LTF.update) LTF.update.Geno <- update(LTF.update, .~. + herb.com.trait$Genotype) #to examine whether Genotype still has an effect anova(LTF.update.Geno) # P = 0.0368 # Iteomyia salicisverruca (vLG) abundance model vLG.null <- lm(herb.com.trait$log_vLG ~ 1, herb.traits) vLG.full <- lm(herb.com.trait$log_vLG ~ ., herb.traits) summary(vLG.full) # adj. R2 = 0.1187 vLG.step <- stepAIC(vLG.null, scope = formula(vLG.full), direction = "forward") vLG.update <- update(vLG.null, .~. + log_size + C_N_imputed) # water content not significant at alpha < 0.05 summary(vLG.update) #plot(vLG.update) vLG.mr.tab <- mr_table(name = "log(Iteomyia salicisverruca + 1)", model = vLG.update) vLG.update.Geno <- update(vLG.update, .~. + herb.com.trait$Genotype) #to examine whether Genotype still has an effect anova(vLG.update.Geno) # still a strong effect of genotype that is not being explained. # Aculus tetanothrix (mLG) abundance model mLG.null <- lm(herb.com.trait$log_mLG ~ 1, herb.traits) mLG.full <- lm(herb.com.trait$log_mLG ~ ., herb.traits) summary(mLG.full) # adj. R2 = 0.1523 mLG.step <- stepAIC(mLG.null, scope = formula(mLG.full), direction = "forward") mLG.update <- update(mLG.null, .~. + cinn.PC2) # although trichomes are significant, it brings it above the adj. R2 cut off. summary(mLG.update) #plot(mLG.update) # residuals don't look good; however, model is still significant after running negative binomial glm so I'm going to stick with it. mLG.glm <- glm.nb(herb.com.trait$mLG ~ cinn.PC2, herb.traits) summary(mLG.glm) #mLG.mr.tab <- mr_table(name = "log(Aculus tetanothrix + 1)", model = mLG.update) # not working because there is only one significant trait effect. mLG.update.Geno <- update(mLG.update, .~. + herb.com.trait$Genotype) # to examine whether Genotype still has an effect anova(mLG.update.Geno) # Genotype no longer has an effect. ## Cicadellidae nymph sp. 1 abundance model. Not reported in manuscript because it only exhibited a marginally significant response to willow genotype. cicnymph.null <- lm(herb.com.trait$log_cicnymph ~ 1, herb.traits) cicnymph.full <- lm(herb.com.trait$log_cicnymph ~ ., herb.traits) summary(cicnymph.full) # not significant, therefore I can't follow up with model selection #plot(cicnymph.full) ## Generate species table. Added mLG results manually because there was only one significant trait effect, so it was difficult to use the 'mr_table' function. species.trait_reg_table <- rbind.data.frame(curc1.mr.tab, LTF.mr.tab, vLG.mr.tab) # Need to add mLG.mr.tab manually. rownames(species.trait_reg_table) <- NULL