[R] Difference between drop1() vs. anova() for Gaussian glm models
peter dalgaard
pdalgd at gmail.com
Mon Jul 20 19:03:20 CEST 2015
I am somewhat surprised that _anything_ sensible comes out of anova.glm(l, test="Chisq"). I think it is mostly expected that you use F tests for that case.
What does seem to come out is the same as for drop1(l, test="Rao"), which gives the scaled score test, which would seem to be equivalent to scaled deviance in this case.
drop1.glm(l, test="Chisq") appears to be calculating the "real" likelihood ratio test, evaluated in its asymptotic chi-square distribution:
> 2*(logLik(l) - logLik(update(l,.~1)))
'log Lik.' 6.944717 (df=3)
(Apologies for the daft output there... Why does "-" not either subtract the df or unclass the whole thing?)
Notice that the scaled tests basically assume that the scale is known, even if it is estimated, so in that sense, the real LRT should be superior. However, in that case it is well known that the asymptotic approximation can be improved by transforming the LRT to the F statistic, whose exact distribution is known.
The remaining part of the riddle is why anova.glm doesn't do likelihood differences in the same fashion as drop1.glm. My best guess is that it tries to be consistent with anova.lm and anova.lm tries not to have to refit the sequence of submodels.
> On 20 Jul 2015, at 14:59 , Karl Ove Hufthammer <karl at huftis.org> wrote:
>
> Dear list members,
>
> I’m having some problems understanding why drop1() and anova() gives
> different results for *Gaussian* glm models. Here’s a simple example:
>
> d = data.frame(x=1:6,
> group=factor(c(rep("A",2), rep("B", 4))))
> l = glm(x~group, data=d)
>
> Running the following code gives *three* different p-values. (I would expect
> it to give two different p-values.)
>
> anova(l, test="F") # p = 0.04179
> anova(l, test="Chisq") # p = 0.00313
> drop1(l, test="Chisq") # p = 0.00841
>
> I’m used to anova() and drop1() giving identical results for the same ‘test’
> argument. However, it looks like the first two tests above use the F-
> statistic as a test statistic, while the last one uses a ‘scaled deviance’
> statistic:
>
> 1-pf(8.7273, 1, 4) # F-statistic
> 1-pchisq(8.7273, 1) # F-statistic
> 1-pchisq(6.9447, 1) # Scaled deviance
>
> I couldn’t find any documentation on this difference. The help page for
> drop1() does say:
>
> The F tests for the "glm" methods are based on analysis of
> deviance tests, so if the dispersion is estimated it is based
> on the residual deviance, unlike the F tests of anova.glm.
>
> But here it’s talking about *F* tests. And drop1() with test="F" actually
> gives the *same* p-value as anova() with test="F":
>
> drop1(l, test="F") # p = 0.04179
>
> Any ideas why anova() and drop(1) uses different test statistics for the
> same ‘test’ arguments? And why the help page implies (?) that the results
> should differ for F-tests (while not mentioning chi-squared test), but here
> they do not (and the chi-squared tests do)?
>
> $ sessionInfo()
> R version 3.2.1 (2015-06-18)
> Platform: x86_64-suse-linux-gnu (64-bit)
> Running under: openSUSE 20150714 (Tumbleweed) (x86_64)
>
> locale:
> [1] LC_CTYPE=nn_NO.UTF-8 LC_NUMERIC=C
> [3] LC_TIME=nn_NO.UTF-8 LC_COLLATE=nn_NO.UTF-8
> [5] LC_MONETARY=nn_NO.UTF-8 LC_MESSAGES=nn_NO.UTF-8
> [7] LC_PAPER=nn_NO.UTF-8 LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=nn_NO.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets
> [6] methods base
>
> loaded via a namespace (and not attached):
> [1] tools_3.2.1
>
> --
> Karl Ove Hufthammer
> E-mail: karl at huftis.org
> Jabber: huftis at jabber.no
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
More information about the R-help
mailing list