HD<-read.table("Vosteen et al. 2016_Fig. 3 honeydew attractivity.txt",header=T) attach(HD) names(HD) View(HD) str(HD) ## treat: influence of previous infestation (induced) and honeydew nr<-as.factor(nr) is.factor(nr) library(lme4) # generalized linar mixed effect model with poisson distribution #### # full model mg1<-glmer(hoverfly_eggs~treat+plant+aphid_pres+ treat:plant+treat:aphid_pres+plant:aphid_pres+ treat:plant:aphid_pres+ (1|nr),family=poisson(link=log)) # error, model too complex, leave out 3-way interaction # new model mg2<-glmer(hoverfly_eggs~treat+plant+aphid_pres+ treat:plant+treat:aphid_pres+plant:aphid_pres+ (1|nr),family=poisson(link=log)) # plant:aphid_pres mg3<-glmer(hoverfly_eggs~treat+plant+aphid_pres+ treat:plant+treat:aphid_pres+ (1|nr),family=poisson(link=log)) anova(mg2,mg3) #Chi=1.3405 df=1 p=0.247 # treat:plant mg4<-glmer(hoverfly_eggs~treat+plant+aphid_pres+ treat:aphid_pres+ (1|nr),family=poisson(link=log)) anova(mg3,mg4) #Chi=9.2615 df=2 p=0.009748 ** # treat:aphid_pres mg5<-glmer(hoverfly_eggs~treat+plant+aphid_pres+ (1|nr),family=poisson(link=log)) anova(mg4,mg5) #Chi=21.272 df=2 p=2.403e-05 *** #aphid_pres mg6<-glmer(hoverfly_eggs~treat+plant+ (1|nr),family=poisson(link=log)) anova(mg5,mg6) #Chi=42.489 df=1 p=7.108e-11 *** # plant mg7<-glmer(hoverfly_eggs~treat+ (1|nr),family=poisson(link=log)) anova(mg6,mg7) #Chi=18.886 df=1 p=1.388e-05 *** # treat mg8<-glmer(hoverfly_eggs~1+ (1|nr),family=poisson(link=log)) anova(mg7,mg8) #Chi=67.052 df=2 p=2.754e-15 *** # thus the minimal model contains treat+plant+aphid_pres+treat:aphid_pres+treat:plant = mg3 summary(mg3) #### test overdispersion mg3 #### ## number of fixed effects = 9, number of random effects = 1 ## ## number of observations = 119 ## ## residual degrees of freedom: 119-10=109 ## (sum(residuals(mg3)^2))/(109) # overdispersed #### negative binomial #### # full model nbg1<-glmer.nb(hoverfly_eggs~treat+plant+aphid_pres+ treat:plant+treat:aphid_pres+plant:aphid_pres+ treat:plant:aphid_pres+(1|nr)) # error, model too complex, leave out 3-way interaction # new model nbg2<-glmer.nb(hoverfly_eggs~treat+plant+aphid_pres+ treat:plant+treat:aphid_pres+plant:aphid_pres+(1|nr)) #plant:aphid_pres nbg3<-glmer.nb(hoverfly_eggs~treat+plant+aphid_pres+ treat:plant+treat:aphid_pres+(1|nr)) anova(nbg2,nbg3) # Chi=0.5523 df=1 p=0.4574 #treat:plant nbg4<-glmer.nb(hoverfly_eggs~treat+plant+aphid_pres+ treat:aphid_pres+(1|nr)) anova(nbg3,nbg4) # Chi=5.0977 df=2 p=0.07817 . # treat:aphid_pres nbg5<-glmer.nb(hoverfly_eggs~treat+plant+aphid_pres+(1|nr)) anova(nbg4,nbg5) # Chi=12.787 df=2 p=0.001673 ** # plant nbg6<-glmer.nb(hoverfly_eggs~treat+aphid_pres+(1|nr)) anova(nbg5,nbg6) # Chi=8.0123 df=1 p=0.004646 ** # aphid_pres nbg7<-glmer.nb(hoverfly_eggs~treat+(1|nr)) anova(nbg6,nbg7) # Chi=16.379 df=1 p=5.186e-05 *** # treat nbg8<-glmer.nb(hoverfly_eggs~1+(1|nr)) anova(nbg7,nbg8) # Chi=18.683 df=2 p=8.769e-05 *** # minimal model contains treat+plant+aphid_pres+treat:aphid_pres = nbg4 summary(nbg4) ##### test overdispersion nbg4 #### ## number of fixed effects = 7, number of random effects = 1 ## ## number of observation observations = 119 ## ## residual degrees of freedom: 119-8=111 ## (sum(residuals(nbg4)^2))/(111) # no overdispersion #### post-hoc test #### require(multcomp) # for treat eggs.multcomp1 <- glht(nbg5, linfct=mcp(treat="Tukey")) summary(eggs.multcomp1) # for treat:aphid_pres interaction eggs.multcomp2 <- glht(nbg9, linfct=mcp(treatment="Tukey")) summary(eggs.multcomp2) detach(HD)