# R code of statistical analysis of the paper entitled 'Climate-induced phenological shift of apple trees has diverse effects on pollinators, herbivores and natural enemies' # Authors: Kőrösi Á, Markó V, Kovács-Hostyánszki A, Somay L, Varga Á, Elek Z, Boreux V, Klein A-M, Földesi R, Báldi A # Date: February 2018 # Data can be found in apple_phenology_mismatch.xls; each sheet corresponds to a csv file library(nlme) library(lme4) library(MASS) library(MuMIn) library(mgcv) library(glmmADMB) library(vegan) #setwd("/home/korosi/Dokumentumok/research_projects/alma/2013/data/repository") # Levels of the variable "treatment2": gh1=Advanced1, gh2=Advanced2, control=Control, ch1=Delayed1, ch2=Delayed2 ####################################### # NUMBER OF FLOWERS ####################################### # datasheet for pollinators and flower numbers pollan<-read.csv("poll_anal.csv",sep=",",dec=".") # Apimel: Apis mellifera (honey bee) abundance; Bombus_t: abundance of bumblebees; Wildbee: abundance of wild bees; # Wildbee_t: abundance of wild bees including bumblebees; block2: a combination of site and block # mbloss: the mean no. flowers of the two sampling events pollan$block2=as.factor(pollan$block2) # a factor with eight levels, used as a random factor later bllme<-lme(sqrt(mbloss)~relevel(treatment2,ref="control"),random=~1|block2,data=pollan,na.action=na.pass) summary(bllme) # treatment sigf VarCorr(bllme) # plot(bllme) EP<-resid(bllme,type="pearson") boxplot(EP~pollan[["block2"]],notch=T) boxplot(EP~pollan[["treatment2"]],notch=T) qqnorm(EP) ####################################### # LEAF SIZE ####################################### leaf<-read.csv2("leaf_size.csv",sep=",",dec=".") leaf$newblock<-as.factor(paste(leaf[["site"]],leaf[["block"]])) leaf_lme<-lme(Area~relevel(treatment,ref="k"),random=~1|newblock/tree,data=leaf) summary(leaf_lme) plot(leaf_lme) qqnorm(resid(leaf_lme)) VarCorr(leaf_lme) ####################################### # SHOOT GROWTH ####################################### sg<-read.csv2("sg.csv",sep=",",dec=".") sglme<-lme(sg~relevel(treat,ref="k"),random=~1|newblock/tid,data=sg) summary(sglme) plot(sglme) qqnorm(resid(sglme)) VarCorr(sglme) # one outlier hist(sg$sg) sg2=sg[c(1:368,370:546),] # omit the outlier sg=9 cm sglme<-lme(sg~relevel(treat,ref="k"),random=~1|newblock/tid,data=sg2) summary(sglme) plot(sglme) qqnorm(resid(sglme)) VarCorr(sglme) # by omitting one outlier the estimations hardly change, but the diagnostic plots look better # model estimates including all data are reported in the paper ####################################### # FRUIT SET ####################################### succ<-read.csv("poll_succ.csv",sep=";",dec=",") boxplot(succ$ps~succ$site) # here 'mbloss' is the max number of flowers found on each tree in the two sampling events # ripe apples were found in only one orchard ("Szabó") mydata<-subset(succ,site=="Szabó") mydata[["block1"]]<-as.factor(mydata[["block1"]]) # mixed effects model with normal error distribution pslme=lme(ps~relevel(treatment2,ref="control"),random=~1|block1,data=mydata) summary(pslme) # VarCorr(pslme) plot(pslme) EP<-resid(pslme,type="pearson") boxplot(EP~mydata[["block1"]],notch=T) boxplot(EP~mydata[["treatment2"]],notch=T) qqnorm(EP) # residuals are not normally distributed # mixed effects model with a so-called 'quasipoisson' error structure suc2<-glmmPQL(ps~relevel(treatment2,ref="control"),random=~1|block1,family=quasipoisson,data=mydata) summary(suc2) VarCorr(suc2) plot(suc2) EP1<-resid(suc2,type="pearson") EP2<-resid(suc2,type="response") boxplot(EP1~mydata$block1,notch=T) boxplot(EP1~mydata$treatment2,notch=T) boxplot(EP2~mydata$block1,notch=T) boxplot(EP2~mydata$treatment2,notch=T) ####################################### ############ POLLINATORS ############## ####################################### # datasheet for pollinators and flower numbers pollan<-read.csv("poll_anal.csv",sep=",",dec=".") # Apimel: Apis mellifera (honey bee) abundance; Bombus_t: abundance of bumblebees; Wildbee: abundance of wild bees; # Wildbee_t: abundance of wild bees including bumblebees; block2: a combination of site and block # mbloss: the mean no. flowers of the two sampling events pollan$block2=as.factor(pollan$block2) # factor with eight levels pollan$cmbloss<-sqrt(pollan[["mbloss"]])-mean(sqrt(pollan[["mbloss"]])) # square root of flower number is centered # it is used as a covariate # model selection procedure is not shown, only the best model ####################################### # ABUNDANCE OF HONEY BEES ####################################### apis.nbinom3<-glmmadmb(Apimel~relevel(treatment2,ref="control")+cmbloss+(1|block2),zeroInflation=FALSE,family="nbinom1",data=pollan) summary(apis.nbinom3)# AIC: 1018.1 var(residuals(apis.nbinom3)) # 0.64 ####################################### # ABUNDANCE OF HOVERFLIES ####################################### p3<-read.csv("poll_anal3.csv",sep=",",dec=".") p3[["block2"]]<-as.factor(p3[["block2"]]) p3$cmbloss<-sqrt(p3[["mbloss"]])-mean(sqrt(p3[["mbloss"]])) # square root of flower number is centered hnb12<-glmmadmb(hoverfly~relevel(treatment,ref="k")+cmbloss+(1|block2),zeroInflation=FALSE,family="nbinom1",data=p3) summary(hnb12) # AIC: 518.7 ch sigf nagyobb, cmbloss sigf pozitiv var(residuals(hnb12)) # 0.76 ####################################### # ABUNDANCE OF WILD BEES (including bumblebees) ####################################### nbinom2c<-glmmadmb(Wildbee_t~relevel(treatment2,ref="control")+cmbloss+(1|block2),zeroInflation=FALSE,family="nbinom1",data=pollan) summary(nbinom2c) # AIC: 782.3 var(residuals(nbinom2c))# 1.1 ####################################### ############ HERBIVORES ############## ####################################### # model selection procedure is not shown, only the best model ######################################## # OCCURRENCE OF A. pomorum ######################################## anthonomus=read.csv(file="anthonomus_pomorum.csv",sep=",",dec=".") # this table contains data of only one orchard ('Donka'), because A. pomorum occurred there only # one level of treatment ('ch1': delayed1) was omitted, because no A. pomorum was observed on that trees # treatment codes: gh1=advanced1, gh2=advanced2, k=control, ch2=delayed2 anthonomus$block=as.factor(anthonomus$block) anpob2<-glmer(anpobin~relevel(treatment,ref="k")+(1|block),family="binomial",data=anthonomus) summary(anpob2) # ch2 sigf kisebb, gh marg non-sigf kisebb # the random factor had only two levels, because there were only two blocks in this orchard ######################################## # NUMBER OF CAPPED BUDS ######################################## bud<-read.csv2("capped_buds.csv",header=TRUE,sep=",",dec=".") # data of capped buds counted on 34 trees in one orchard ('Donka') # since there were only two blocks in this orchard, we did not use mixed-effect models # capped buds were found only on advanced and control trees (three levels of treatment) # Poisson model on the number of capped buds showed overdispersion rgp=glm(rg~relevel(treatment,ref="k"),family=poisson,data=bud) summary(rgp) bud$lrg<-log(bud[["rg"]]+1) # log-transformation # Linear model with normal error distribution on the log-transformed number of capped buds rgl=lm(lrg~relevel(treatment,ref="k"),data=bud) summary(rgl) plot(rgl) # non-normal error distribution # GLM with quasipoisson error structure lrgq<-glm(lrg~relevel(treatment,ref="k"),family=quasipoisson,data=bud) summary(lrgq) # no overdispersion ############################################ # APHID ABUNDANCE ############################################ aphids=read.csv(file="aphid_data.csv") # response variable m3atotal: the mean of estimated aphid abundance in three samplings for each tree # treatment codes: gh1=advanced1, gh2=advanced2, k=control, ch1=delayed1, ch2=delayed2 # cm3gp: centred mean proportion of growing shoots of three samplings for each tree ma_lme<-lme(m3atotal~relevel(treatment,ref="k")*cm3gp,random=~1|newblock,data=aphids,na.action=na.pass) summary(ma_lme) # AIC: 158.44 dredge(ma_lme) best_ma<-lme(m3atotal~relevel(treatment,ref="k")*cm3gp,random=~1|newblock,data=aphids) summary(best_ma) plot(best_ma) qqnorm(resid(best_ma)) plot(resid(best_ma)~a614[["treatment"]],notch=TRUE) VarCorr(best_ma) ########################################### # ABUNDANCE OF Stephanitis pyri ########################################### nzspyri=read.csv(file="stephanitis_pyri.csv") # response variable ln: log-transfomed non-zero abundance # treatment codes: gh1=advanced1, gh2=advanced2, k=control, ch1=delayed1, ch2=delayed2 # newblock: factor with 8 levels qnzsp3<-glmmPQL(ln~relevel(treatment,ref="k"),random=~1|newblock,family="quasipoisson",data=nzspyri) summary(qnzsp3) plot(qnzsp3) qnzEP<-resid(qnzsp3,type="pearson") plot(qnzEP~nzsp3[["treatment"]],notch=TRUE) plot(qnzEP~nzsp3[["newblock"]],notch=TRUE) VarCorr(qnzsp3) ################################################# # ABUNDANCE OF PHYTOPHAGOUS BUGS (excl. S. pyri) ################################################# fb_sm<-read.csv("fitobug_fajhely.csv",sep=",",header=TRUE) # species matrix fb_sm=fb_sm[,-c(10,28)] # omit columns Coreidae and Pentatomidae fb_sm$abund=apply(fb_sm[,5:33],1,sum) # calculating total abundance fb_sm$newblock=as.factor(paste(fb_sm[["site"]],fb_sm[["block"]])) # factor for blocks # response variable abund: abundance of phytophagous bugs # treatment codes: gh1=advanced1, gh2=advanced2, k=control, ch1=delayed1, ch2=delayed2 nb1=glmer.nb(abund~relevel(treatment,ref="k")+(1|newblock),data=fb_sm) summary(nb1) # ############################################### # NATURAL ENEMIES ############################################### ##################################### # ABUNDANCE OF APHIDOPHAGOUS BEETLES ##################################### beetle_aphid=read.csv(file="ladybug_aphid.csv") # m3atotal is the log-transformed mean abundance of aphids (see above) # only the best model is shown here adp3best<-glmer(n~relevel(treatment.x,ref="k")+m3atotal+(1|newblock.x),family="poisson",data=beetle_aphid,na.action=na.pass) summary(adp3best) vc=VarCorr(adp3best) print(vc,comp=c("Variance","Std.Dev."),digits=3) var(resid(adp3best)) ################################### # ABUNDANCE OF ZOOPHAGOUS BUGS ################################### zb_aphid=read.csv(file="zoobug_aphid.csv") zbnb3<-glmmadmb(n~relevel(treatment.x,ref="k")+(1|newblock.x),family="nbinom",zeroInflation=FALSE,data=zb_aphid) summary(zbnb3) # var(resid(zbnb3)) ################################### # ABUNDANCE OF SPIDERS ################################### spiders=read.csv(file="spider_spyri.csv") pokcspnb<-glmmadmb(n.x~relevel(treatment.x,ref="k")*cspln+(1|newblock.x),family="nbinom",zeroInflation=FALSE,data=spiders) summary(pokcspnb) # var(resid(pokcspnb)) ############################################# # ANALYSIS OF SPECIES COMPOSITION ############################################# ############################## # WILDBEES ############################## wb_sm<-read.csv("wildbee_fajhely.csv",sep=",",header=TRUE) # Redundancy analysis std_wb<-decostand(wb_sm[,5:43],"hellinger") rda1<-rda(std_wb~wb_sm[["treatment"]]) rda1 permutest(rda1,permu=10000) # sigf, p<0.0001, PseudoF(4,177)=6.187; 12.3% rda2<-rda(std_wb~wb_sm[["site"]]) rda2 permutest(rda2,permu=10000) # sigf, p<0.0001, PseudoF(2,179)=3.17; 3.4% # Renyi diversity profile kezeles<-aggregate(wb_sm[,5:43],by=list(treat=wb_sm[["treatment"]]),sum) renyis<-renyi(kezeles[,2:40]) renyis<-as.matrix(renyis) min(renyis) # 0.52 max(renyis) # 2.89 xl<-c(1:11) plot(renyis[1,]~xl,ylim=c(0.5,3),las=1,bty="l",type="l",xlab="Scale parameter",ylab="Diversity",col="blue",lwd=2,axes=FALSE,cex.lab=1.4,lty=3,main="") # ch1 lines(renyis[2,]~xl,col="blue",lty=2,lwd=2) # ch2 lines(renyis[3,]~xl,col="red",lty=2,lwd=2) # gh1 lines(renyis[4,]~xl,col="red",lty=3,lwd=2) # gh2 lines(renyis[5,]~xl,col="black",lty=1,lwd=2) # k legend(8,3,c("Advanced1","Advanced2","Control","Delayed1","Delayed2"),lty=c(2,3,1,3,2),col=c("red","red","black","blue","blue"),bty="n",lwd=2,cex=1.2) axis(1,at=xl,labels=c("0","0.25","0.5","1","2","4","8","16","32","64","Inf"),cex.axis=1.4) axis(2,at=c(0.5,1,1.5,2,2.5,3),las=1,cex.axis=1.4) mtext(c("A","Wild bees"),3,1,at=c(1,6),font=c(2,2),cex=1.4) ############################## # Phytophagous bugs ############################## fb_sm<-read.csv("fitobug_fajhely.csv",sep=",",header=TRUE) # species matrix fb_sm=fb_sm[,-c(10,28)] # omit Coreidae and Pentatomidae std_fb<-decostand(fb_sm[,5:33],"hellinger") # Redundancy analysis rda1<-rda(std_fb~fb_sm[["treatment"]]) rda1 permutest(rda1,permu=10000) # marg non-sigf, p=0.067; 3.3% Pseudo-F(4,177): 1.49 rda2<-rda(std_fb~fb_sm[["site"]]) rda2 permutest(rda2,permu=10000) # sigf, p=0.00001 11.8% Pseudo-F(2,179): 12.02 # Renyi diversity profile kezeles<-aggregate(fb_sm[,5:33],by=list(treat=fb_sm[["treatment"]]),sum) renyis<-renyi(kezeles[,2:29]) renyis<-as.matrix(renyis) min(renyis) # 0.525 max(renyis) # 2.77 # plot xl<-c(1:11) plot(renyis[1,]~xl,ylim=c(0.5,3),las=1,bty="l",type="l",xlab="Scale parameter",ylab="Diversity",col="blue",lwd=2,axes=FALSE,cex.lab=1.4,lty=3,main="") # ch1 lines(renyis[2,]~xl,col="blue",lty=2,lwd=2) # ch2 lines(renyis[3,]~xl,col="red",lty=2,lwd=2) # gh1 lines(renyis[4,]~xl,col="red",lty=3,lwd=2) # gh2 lines(renyis[5,]~xl,col="black",lty=1,lwd=2) # k legend(8,3,c("Advanced1","Advanced2","Control","Delayed1","Delayed2"),lty=c(2,3,1,3,2),col=c("red","red","black","blue","blue"),bty="n",lwd=2,cex=1.2) axis(1,at=xl,labels=c("0","0.25","0.5","1","2","4","8","16","32","64","Inf"),cex.axis=1.4) axis(2,at=c(0.5,1,1.5,2,2.5,3),las=1,cex.axis=1.4) mtext(c("B","Phytophagous bugs"),3,1,at=c(1,6),font=c(2,2),cex=1.4) #################################### # APHIDOPHAGOUS BEETLES #################################### ab_sm<-read.csv("afidofag_beetle_fajhely.csv",sep=",",header=TRUE) # Redundancy analysis std_ab<-decostand(ab_sm[,5:16],"hellinger") rda1<-rda(std_ab~ab_sm[["treatment"]]) rda1 permutest(rda1,permu=10000) # sigf p=0.025 Pseudo-F(4,177): 1.91, 18.3% rda2<-rda(std_ab~ab_sm[["site"]]) rda2 permutest(rda2,permu=10000) # non-sigf p=0.126 Pseudo-F(2,179): 1.6, 0.8% # Renyi diversity profile kezeles<-aggregate(ab_sm[,5:16],by=list(treat=ab_sm[["treatment"]]),sum) renyis<-renyi(kezeles[,2:13]) renyis<-as.matrix(renyis) min(renyis) # 0.33 max(renyis) # 2.2 xl<-c(1:11) plot(renyis[1,]~xl,ylim=c(0,2.5),las=1,bty="l",type="l",xlab="Scale parameter",ylab="Diversity",col="blue",lwd=2,axes=FALSE,cex.lab=1.4,lty=3,main="") # ch1 lines(renyis[2,]~xl,col="blue",lty=2,lwd=2) # ch2 lines(renyis[3,]~xl,col="red",lty=2,lwd=2) # gh1 lines(renyis[4,]~xl,col="red",lty=3,lwd=2) # gh2 lines(renyis[5,]~xl,col="black",lty=1,lwd=2) # k legend(8,2.5,c("Advanced1","Advanced2","Control","Delayed1","Delayed2"),lty=c(2,3,1,3,2),col=c("red","red","black","blue","blue"),bty="n",lwd=2,cex=1.2) axis(1,at=xl,labels=c("0","0.25","0.5","1","2","4","8","16","32","64","Inf"),cex.axis=1.4) axis(2,at=c(0,0.5,1,1.5,2,2.5),las=1,cex.axis=1.4) mtext(c("C","Aphidophagous beetles"),3,1,at=c(1,6),font=c(2,2),cex=1.4) #################################### # SPIDERS #################################### pok_sm<-read.csv("spider_fajhely.csv",sep=",",header=TRUE) std_pok<-decostand(pok_sm[,5:44],"hellinger") # Redundancy analysis rda1<-rda(std_pok~pok_sm[["treatment"]]) rda1 permutest(rda1,permu=10000) # non-sigf, p=0.059; 2.95% Pseudo-F(4,177): 1.34 rda2<-rda(std_pok~pok_sm[["site"]]) rda2 permutest(rda2,permu=10000) # sigf, p<0.001, 3.52% Pseudo-F(2,179): 3.27 # Renyi diversity profile kezeles<-aggregate(pok_sm[,5:44],by=list(treat=pok_sm[["treatment"]]),sum) renyis<-renyi(kezeles[,2:41]) renyis<-as.matrix(renyis) min(renyis) # 1.196 max(renyis) # 3.14 xl<-c(1:11) plot(renyis[1,]~xl,ylim=c(1,3.5),las=1,bty="l",type="l",xlab="Scale parameter",ylab="Diversity",col="blue",lwd=2,axes=FALSE,cex.lab=1.4,lty=3,main="") # ch1 lines(renyis[2,]~xl,col="blue",lty=2,lwd=2) # ch2 lines(renyis[3,]~xl,col="red",lty=2,lwd=2) # gh1 lines(renyis[4,]~xl,col="red",lty=3,lwd=2) # gh2 lines(renyis[5,]~xl,col="black",lty=1,lwd=2) # k legend(8,3.5,c("Advanced1","Advanced2","Control","Delayed1","Delayed2"),lty=c(2,3,1,3,2),col=c("red","red","black","blue","blue"),bty="n",lwd=2,cex=1.2) axis(1,at=xl,labels=c("0","0.25","0.5","1","2","4","8","16","32","64","Inf"),cex.axis=1.4) axis(2,at=c(1,1.5,2,2.5,3,3.5),las=1,cex.axis=1.4) mtext(c("D","Spiders"),3,1,at=c(1,6),font=c(2,2),cex=1.4)