[R] accuracy of GLM dispersion parameters
Ben Bolker
bbolker at gmail.com
Tue Nov 30 16:30:09 CET 2010
<Timothy_Handley <at> nps.gov> writes:
> I'm confused as to the trustworthiness of the dispersion parameters
> reported by glm. Any help or advice would be greatly appreciated.
>
[snip]
> Specifics: The summary function says that my fitted GLM has a dispersion
> parameter=15.8. On the other hand, the gamma.dispersion function (MASS)
> says that the GLM uses a dispersion parameter of 1.86. I could understand
> some modest difference, as the help for gamma.shape() says that the MASS
> functions return a more accurate dispersion value than summary(). However,
> these two numbers differ by a factor of 8, which is quite a lot. Is this
> normal? Would you folks expect such a large difference? Which value should
> I trust?
>
> R terminal excerpt:
>
>
> > summary(tempglm_g2)
>
> Call:
> glm(formula = precip_sbi ~ precip_oxx + precip_oxx_sq, family = Gamma(link
> = identity),
> data = w.combo, start = c(0.1, 0.4, 0.02))
>
> Deviance Residuals:
> Min 1Q Median 3Q Max
> -2.99999 -1.63183 -1.00720 0.04878 8.93461
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.09236 0.04834 1.911 0.0583 .
> precip_oxx 0.26848 0.35891 0.748 0.4558
> precip_oxx_sq 0.05138 0.13418 0.383 0.7024
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> (Dispersion parameter for Gamma family taken to be 15.78978)
>
> Null deviance: 528.73 on 130 degrees of freedom
> Residual deviance: 305.81 on 128 degrees of freedom
> AIC: -100.33
>
> Number of Fisher Scoring iterations: 5
>
> > library(MASS)
> > gamma.shape(tempglm_g2)
>
> Alpha: 0.53807358
> SE: 0.05526108
> > gamma.dispersion(tempglm_g2)
> [1] 1.858482
>
I would trust the gamma.dispersion() result more, although I
agree that the difference is worrisome. The way to look at this
further would be to profile the dispersion parameter. As I recall
there isn't such a built in option in MASS (profile.glm only
profiles the coefficients), but you may be able to do it
*approximately* like this:
library(bbmle)
m1 <- mle2(precip_sbi ~ dgamma(shape=a,scale=mu/a),
parameters=list(mu~precip_oxx+precip_oxx_sq),
data=w.combo, start=list(mu=0.1,a=2))
p1 <- profile(m1)
plot(p1)
More information about the R-help
mailing list