--- title: "Tielensetal2020_EcoEnto_code" author: "Tielens" date: "6/24/2020" output: html_document --- ```{r setup} rm(list=ls()) library(ggplot2) library(scales) library(lattice) library(vegan) library(devtools) library(MuMIn) library(iNEXT) library(plyr) library(mvabund) library(reshape2) library(lme4) # Utilizes phytophagous insect data sampled in 2014 and 2015, both in long and in wide format. # full dataset saved as: data_Tielensetal2020_EcoEnto ``` ```{r Sampling effort } # sampling summary: sum(rowSums(hem)) mean(ChaoRichness(t(hem[env$sampleyear=="2014",]))[,2]) mean(ChaoRichness(t(hem[env$sampleyear=="2015",]))[,2]) # species richness and shannnon diversity across samples: rich<-ChaoRichness(as.data.frame(t(hem.2014) ))[,2] shan<-exp(ChaoShannon(as.data.frame(t(hem.2014)) )[,2]) ``` ```{r (i) How do phytophagous insect abundance and diversity vary across host plant phenotype? } # Generating figure 2: tiff(filename = "fig_richandabun_morph.tiff", width =8, height = 6, units = "in", compression = c("lzw"), res = 300) par(mfrow=c(1,2)) par(oma=c(1,0.5,1,0.5)) plot(rowSums(hem)~factor(env$morphotype),names = c("Glabrous" , "Intermediate","Pubescent"), xlab="", cex.lab=0.8, cex.axis=0.7, ylab="") mtext(text = "Abundance", side = 2, line = 2) mtext(paste0("(a)"), side = 3, adj = 0.05, line = -1.3) plot(ChaoRichness(as.data.frame(t(hem)), datatype="abundance")[,2]~factor(env$morphotype), ylab="", xlab="",names = c("Glabrous" , "Intermediate","Pubescent"), cex.axis=0.7, cex.lab=0.8 ) mtext(text = "Effective Species Richness", side = 2, line = 2) mtext(paste0("(b)"), side = 3, adj = 0.05, line = -1.3) mtext(paste0("a"), side = 3, adj = 0.19, line = -5.5) mtext(paste0("b"), side = 3, adj = 0.5, line = -5.5) mtext(paste0("a"), side = 3, adj = 0.81, line = -5.5) dev.off() # Center and scale leaf traits: hem.2014.long$percN.standard <-scale(hem.2014.long$percN, center=T, scale=T) hem.2014.long$percP.standard <-scale(hem.2014.long$percP, center=T, scale=T) hem.2014.long$meanSLA.standard <-scale(hem.2014.long$meanSLA, center=T, scale=T) hem.2014.long$meanWaterstandard <-scale(hem.2014.long$meanWatercontent, center=T, scale=T) # Models: # Does species richness differ by morphotype? mm<-lmer(ChaoRichness(t(hem), datatype="abundance")[,2]~env$morphotype+(1|env$site)) summary(mm) mnull<-lmer(ChaoRichness(t(hem), datatype="abundance")[,2]~(1|env$site)) anova(mnull,mm) # Does abundance differ by morphotype? mabun<-lmer(rowSums(hem)~env$morphotype+(1|env$site)) summary(mabun) mabunnull<-lmer(rowSums(hem)~(1|env$site)) anova(mabunnull,mabun) # Does shannon diversity vary with plant traits? mod1<-lmer(shan~env.2014$morphotype+leaftraits$percN+leaftraits$percP+leaftraits$meanSLA+leaftraits$meanWatercontent+(1|env.2014$site)+offset(log(env.2014$DW))) shan.null<-lmer(shan~(1|env.2014$site)+offset(log(env.2014$DW))) anova(mod1,shan.null) # null model has lowest AIC qqnorm(resid(shan.null)) plot(shan.null) mod2<-lmer(shan~env.2014$morphotype+(1|env.2014$site)+offset(log(env.2014$DW))) mod3<-lmer(shan~leaftraits$percN+(1|env.2014$site)+offset(log(env.2014$DW))) mod4<-lmer(shan~leaftraits$percP+(1|env.2014$site)+offset(log(env.2014$DW))) mod5<-lmer(shan~leaftraits$meanSLA+(1|env.2014$site)+offset(log(env.2014$DW))) mod6<-lmer(shan~leaftraits$meanWatercontent+(1|env.2014$site)+offset(log(env.2014$DW))) mod7<-lmer(shan~leaftraits$percN+leaftraits$morphotype+(1|env.2014$site)+offset(log(env.2014$DW))) anova(shan.null,mod7) anova(shan.null,mod1) anova(shan.null,mod2) anova(shan.null,mod6) summary(mod3) # Conclusion: null model for shannon diversity # Does species richness differ with leaf traits? modnull<-lmer(rich~(1|env.2014$site)+offset(log(env.2014$DW))) mod2<-lmer(rich~env.2014$morphotype+(1|env.2014$site)+offset(log(env.2014$DW))) mod3<-lmer(rich~leaftraits$percN+(1|env.2014$site)+offset(log(env.2014$DW))) mod4<-lmer(rich~leaftraits$percP+(1|env.2014$site)+offset(log(env.2014$DW))) mod5<-lmer(rich~leaftraits$meanSLA+(1|env.2014$site)+offset(log(env.2014$DW))) mod6<-lmer(rich~leaftraits$meanWatercontent+(1|env.2014$site)+offset(log(env.2014$DW))) mod7<-lmer(rich~leaftraits$percN+leaftraits$morphotype+(1|env.2014$site)+offset(log(env.2014$DW))) anova(mod5,modnull) # mod 4 summary(mod4) # Conclusion: model including foliar phosphorus for species richness ``` ```{r (ii) Which host plant traits drive insect community composition? } #################### ##### Model selection #################### # First pass: what type of model mod.nb<-glmer.nb(data=hem.long, value~(1|variable)+(1|site)+offset(log(DW)), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 10000000), tol=0.0002)) mod.poisson =glmer(value ~ (1|variable)+(1|site)+offset(log(DW)), data =hem.long, family =poisson) anova(mod.poisson, mod.nb) # check for overdispersion: pcsq<-sum(resid(mod.nb, type="pearson")^2) rdf<-df.residual(mod.nb) pchisq(pcsq,rdf,lower.tail=T) pcsq/rdf # Diagnostics using Dharma: simulationOutputM <- simulateResiduals(fittedModel = mod.nb , n = 1000) hist(simulationOutputM$scaledResiduals) plot(simulationOutputM) testDispersion(simulationOutput = simulationOutputM) # dispersion looks good testZeroInflation(simulationOutput = simulationOutputM) # does not indicate being zero inflated testTemporalAutocorrelation(simulationOutput = simulationOutputM) testSpatialAutocorrelation(simulationOutput = simulationOutputM) # some spatial autocorrrelation, NS # Negative binomial model. Test predictors using dredge: mod3<-glmer.nb(data=hem.long, value~percN.standard +morphotype+meanWaterstandard +meanSLA.standard+ percP.standard+(1|variable) +(1|site)+offset(log(DW)), na.action="na.fail", glmerControl(optimizer= "bobyqa", optCtrl = list(maxfun = 10000000), tol=0.0002)) null<-glmer.nb(data=hem.long, value~(1|variable)+(1|site)+offset(log(DW)), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 10000000), tol=0.0002)) anova(mod3, null) dredge(mod3, beta = c("none"), evaluate = TRUE, rank = "AICc") #################### ##### NMDS visualization of differences between morphotypes within sites #################### # Subset data for easy plotting: hem.ero<-hem[env$site=="ERO",] hem.ali<-hem[env$site=="ALI",] hem.kaiy<-hem[env$site=="KAIY",] hem.olaa<-hem[env$site=="OLAA",] # Plot: momentKAIY<-metaMDS(hem.kaiy, k=3) mds.fig <- ordiplot(momentKAIY, type = "none", main="Kaiholena") points(mds.fig, "sites", pch = 17, cex=1.5,col = "blue", select = env.kaiy$morphotype == "G" & env.kaiy$sampleyear=="2014") points(mds.fig, "sites", pch = 2, cex=1.5,col = "blue", select = env.kaiy$morphotype == "G" & env.kaiy$sampleyear=="2015") points(mds.fig, "sites", pch = 19, cex=1.5,col = "red", select = env.kaiy$morphotype == "P" & env.kaiy$sampleyear=="2014") points(mds.fig, "sites", pch = 1, cex=1.5,col = "red", select = env.kaiy$morphotype == "P" & env.kaiy$sampleyear=="2015") #text(mds.fig, "sites",adj=c(1,-1.2),labels=env.kaiy$sampleyear) Lab<-c( "Glabrous - 2014","Glabrous - 2015", "Pubescent - 2014", "Pubescent - 2015") legend( "topright", legend=Lab,pch=c(17,2,19,1),col=c("blue","blue","red", "red"), pt.cex=0.9, cex=0.7, xpd=TRUE) mtext(paste0("(a)"), side = 3, adj = 0.05, line = -1.3) momentOLAA<-metaMDS(hem.olaa, k=3) mds.fig <- ordiplot(momentOLAA, type = "none", main="HAVO Olaa") points(mds.fig, "sites", pch = 17,cex=1.5, col = "blue", select = env[env$site=="OLAA",]$morphotype == "G" & env[env$site=="OLAA",]$sampleyear == "2014") points(mds.fig, "sites", pch = 2,cex=1.5, col = "blue", select = env[env$site=="OLAA",]$morphotype == "G" & env[env$site=="OLAA",]$sampleyear == "2015") points(mds.fig, "sites", pch = 15,cex=1.5, col = "purple", select = env[env$site=="OLAA",]$morphotype == "H" & env[env$site=="OLAA",]$sampleyear == "2014") points(mds.fig, "sites", pch = 0,cex=1.5, col = "purple", select = env[env$site=="OLAA",]$morphotype == "H" & env[env$site=="OLAA",]$sampleyear == "2015") #text(mds.fig, "sites",adj=c(1,-1.2),labels=env[env$site=="OLAA",]$sampleyear) Lab<-c("Glabrous - 2014","Glabrous - 2015","Intermediate - 2014", "Intermediate - 2015") color<-c("blue", "purple") legend( "topright", legend=Lab,pch=c(17,2,15, 0),col=c("blue", "blue", "purple", "purple"),pt.cex=0.9,cex=0.7,xpd=TRUE) mtext(paste0("(b)"), side = 3, adj = 0.05, line = -1.3) # Test: adonis(hem.ali~env[env$site=="ALI",]$morphotype, by="margin") adonis(hem.kaiy~env[env$site=="KAIY",]$morphotype, by="margin") adonis(hem.ero~env[env$site=="ERO",]$morphotype, by="margin") adonis(hem.olaa~env[env$site=="OLAA",]$morphotype, by="margin") ``` ```{r (iii) Does the interaction between insect traits and host plant traits further explain insect community composition? } # Function to plot figure nicely envtraitplot<-function(modelname, xname, yname){ a = max( abs(modelname$fourth.corner),na.rm=T ) colort = colorRampPalette(c("blue","white","yellow")) plot.4th = levelplot(t(as.matrix(modelname$fourth.corner)), xlab=xname, ylab=yname, col.regions=colort(100), at=seq(-a, a, length=100), scales = list( x= list(rot = 25), fontface="bold")) print(plot.4th) } # Confirm that datasets are ordered correctly leaftraits<-leaftraits[order(rownames(leaftraits)),] hem.2014<-hem.2014[order(rownames(hem.2014)),] # Trait * Env Interaction Models: # resampling based testing of the interaction term ftSmall=traitglm(hem.2014,leaftraits,traits, method="manyglm" ,family="negative.binomial") anova.traitglm(ftSmall, nBoot=999, show.time=T) envtraitplot(ftSmall) # model including LASSO penalty model4th<-traitglm(hem.sub, leaftraits, traits, method="glm1path", family="negative.binomial") model4th$fourth envtraitplot(model4th) ```