honeydew<-read.table("Vosteen et al. 2016_Fig. 4b honeydew quantity.txt",header=T) attach(honeydew) str(honeydew) nr<-as.factor(nr) is.factor(nr) honeydew_dropletts<-as.factor(honeydew_dropletts) is.factor(honeydew_dropletts) boxplot(eggs~honeydew_dropletts,ylab="aphids") #### generalized linear mixed effect model with poisson distribution #### library(lme4) model0<-glmer(eggs~1+(1|nr),family=poisson(link=log)) model1<-glmer(eggs~honeydew_dropletts+(1|nr),family=poisson(link=log)) anova(model1,model0) summary(model1) #### test for overdispersion #### ## number of fixed effects = 3, number of random effects = 1 ## ## number of observations = 30 ## ## residual degrees of freedom: 30-4 ## E<-resid(model1,type="pearson") (sum(E)^2))/(26) shapiro.test(unlist(ranef(model1)$nr)) levels(honeydew_dropletts) honeydew_dropletts2<-(honeydew_dropletts) levels(honeydew_dropletts2)[c(1,2)]<-"6 + 12" levels(honeydew_dropletts2) model2<-glmer(eggs~honeydew_dropletts2+(1|nr),family=poisson(link=log)) anova(model1,model2) detach(honeydew)