[R] How to use quasibinomial?

Henric Nilsson hennil at statisticon.se
Thu Jul 3 11:33:38 CEST 2003


Dear all,

I've got some questions, probably due to misunderstandings on my behalf, related
to fitting overdispersed binomial data using glm().

1. I can't seem to get the correct p-values from anova.glm() for the F-tests when
supplying the dispersion argument and having fitted the model using
family=quasibinomial. Actually the p-values for the F-tests seems identical to the
p-values for the Chi-squared tests. When not supplying the dispersion argument,
i.e. when anova.glm() uses the default scaled Pearson statistic from
family=quasibinomial, both tests returns the p-values I'd expect. What am I doing
wrong here and how can I make it work?

> fit.1<-glm(y/n~host*variety,family=quasibinomial,weights=n)
> dscale<-sum(residuals(fit.1,type="deviance")^2/fit.1$df.residual)
> dscale
[1] 1.957517

> anova(fit.1,test="F",dispersion=dscale)
Analysis of Deviance Table
Model: quasibinomial, link: logit
Response: y/n
Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev       F    Pr(>F)
NULL                            20     98.719
host          1   55.969        19     42.751 28.5916 8.937e-08
variety       1    3.065        18     39.686  1.5657    0.2108
host:variety  1    6.408        17     33.278  3.2736    0.0704

I expected:
> 1-pf(3.2736,1,17)
[1] 0.08812074

> anova(fit.1,test="Chisq",dispersion=dscale)
Analysis of Deviance Table
Model: quasibinomial, link: logit
Response: y/n
Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL                            20     98.719
host          1   55.969        19     42.751 8.937e-08
variety       1    3.065        18     39.686     0.211
host:variety  1    6.408        17     33.278     0.070

As expected:
> 1-pchisq(6.408/dscale,1)
[1] 0.07040576

2. When using summary.glm() on a glm object fitted using family=quasibinomial the
reported tests are t-tests. Why?

Thanks,
Henric




More information about the R-help mailing list