grpVars <- tapply(inverts\$PLD, inverts\$tft, var) summary(grpVars) grpMeans <- tapply(inverts\$PLD,inverts\$tft, mean) summary(grpMeans) # Quasi-Poisson lm1 <- lm(grpVars~grpMeans-1) phi.fit <- coef(lm1) # The -1 specifies a model with the intercept set to zero # Negative binomial lm2 <- lm(grpVars ~ I(grpMeans^2) + offset(grpMeans)-1) k.fit <- 1/coef(lm2) # Offset() used to specify that we want group means added as a term with its coefficient fixed to 1 # Non-parametric loess fit Lfit <- loess(grpVars~grpMeans) plot(grpVars ~ grpMeans, xlab="group means", ylab="group variances" ) abline(a=0,b=1, lty=2) text(60,200,"Poisson") curve(phi.fit*x, col=2,add=TRUE) # bquote() is used to substitute numeric values in equations with symbols text(60,800, bquote(paste("QP: ",sigma^2==.(round(phi.fit,1))*mu)),col=2) curve(x*(1+x/k.fit),col=4,add=TRUE) text(60,1600,paste("NB: k=",round(k.fit,1),sep=""),col=4) mvec <- 0:120 lines(mvec,predict(Lfit,mvec),col=5) text(50,1300,"loess",col=5) # Poisson GLMM mp1 <- glmer(PLD ~ temp*feeding.type + (1|taxon), data=inverts, family="poisson") overdisp_fun(mp1) # ratio is significantly >1 # NB GLMM mnb1 <- glmer.nb(PLD ~ temp*feeding.type + (1|taxon), data=inverts) overdisp_fun(mnb1) # Looks good!