[6565] in s-news-athena
Re: Deviance tests for glm objects
daemon@ATHENA.MIT.EDU (Trevor Hastie)
Sat Mar 4 14:33:38 1995
From: trevor@mallet.Stanford.EDU (Trevor Hastie)
To: S-news@utstat.toronto.edu
Date: Sat, 4 Mar 1995 11:16:51 -0800 (PST)
In-Reply-To: <Pine.HPP.3.90.950304101511.17657A@baloun.entu.cas.cz> from "Petr Smilauer" at Mar 4, 95 10:26:19 am
Petr Smilauer writes:
>
> Dear colleagues,
> I was wondering whether someone could give a hint with what seems to
> be a problem to me.
> Using S-Plus version 3.2, MS Windows brand, say, on sample data
> corn.yield and corn.rain:
> >glm.0<-glm(corn.yield~+1,poisson)
> >glm.1<-glm(corn.yield~corn.rain,poisson)
> Now, if I do:
> >anova(glm.0,glm.1,test="F")
> I get:
> ...
> Terms Resid.Df Resid.Dev Test Df Deviance F Value Pr(F)
> 1 37 23.22
> corn.rain 36 19.68 1 3.544 6.483 0.0153
>
> Note, that I've rounded the figures.
> Now, if I do:
> >anova(glm.1,test="F")
> [where I suppose the object's method anova.glm() is called]
> I get something like:
> ....
> Df Deviance Resid.Df Resid.Dev F value Pr(F)
> NULL 37 23.22
> corn.rain 1 3.544 36 19.68 6.682 0.0139
>
> NOTE, THAT the resid.Dev, deviance difference between the models, res.d.f.
> etc are the same but not F value and, hence, Pr(F). Tracing the source
> for anova.glm() down, I find it calls (at the end) stat.anova() method,
> with third parameter ('scale') being set with expression
> deviance.lm(object)/object$df.resid . It seems to me that 'deviance.lm'
> (which calculates only sum of squares of residuals) seems to be the
> problem.
> With the presumption of innocence, I guess I should be wrong, using
> an inappropriate function. Sight, using F test here might not be best
> pick (well, is NOT) here, but still, the program shouldn't behave so..
> Any help is welcome.
> Thanks!
> Petr Smilauer
> Univers.of South Bohemia
> Czech Republic
> <petrsm@entu.cas.cz>
>
This does seem to be a bug in deviance.glmlist. When an "F" test is
requested, the underlying model is either an over or under dispersed
binomial or else a quasi likelihood model all of which require a
dispersion parameter. The McCullagh and Nelder recommended estimate is
the scaled residual Chi-square, as delivered by deviance.lm(). So its
seems anova.glm() is correct and anova.glmlist() is in error. The fix
is to change the line (near the end):
stat.anova(aod, test, dev[o[1]]/dfres[o[1]], dfres[o[1]], n)
to
stat.anova(aod, test,deviance.lm(object[[o[1]]])/dfres[o[1]], dfres[o[1]], n)
--
Trevor Hastie trevor@playfair.stanford.edu
Phone: 415-725-2231 Fax: 415-725-8977
Statistics Department, Sequoia Hall, Stanford University, Ca94305
