data(wolbachia)
# Binomial fit
model1 <- glm(cbind(y, m-y) ~ treat, family=binomial, data=wolbachia)
anova(model1, test="Chisq")
# Quasi-binomial fit
model2 <- glm(cbind(y, m-y) ~ treat, family=quasibinomial, data=wolbachia)
summary(model2)
anova(model2,test="F")
## half normal plots
par(mfrow=c(1,2),cex=1.2, cex.main=1.1)
hnp(model1, print=TRUE, pch=4, main="(a) Binomial",
xlab="Half-normal scores", ylab="Deviance residuals")
hnp(model2, print=TRUE, pch=4, main="(b) Quasi-binomial",
xlab="Half-normal scores", ylab='')
if (FALSE) {
# Beta-binomial fit
### using package aods3
require(aods3)
model3 <- aodml(cbind(y, m-y) ~ treat, family='bb', data=wolbachia)
hnp(model3, print=TRUE, pch=4,
xlab="Half-normal scores", ylab='Deviance residuals')
### using package VGAM
require(VGAM)
model3a <- vglm(cbind(y, m-y) ~ treat, family=betabinomial,
data=wolbachia)
summary(model3a)
1/(1+exp(-coef(model3a)[2])) # phi estimate
hnp(model3a, data=wolbachia)
}
## for discussion on the analysis of this data set,
## see Demetrio et al. (2014)
Run the code above in your browser using DataLab