[R] anova of lme objects (model1, model2) gives different results depending on order of models

[Albyn Jones]

No, both yield the same result: reject the null hypothesis,
which always corresponds to the restricted (smaller) model.

albyn

On Thu, May 31, 2012 at 12:47:30PM +0100, Chris Beeley wrote:
> Hello-
> 
> I understand that it's convention, when comparing two models using
> the anova function anova(model1, model2), to put the more
> "complicated" (for want of a better word) model as the second model.
> However, I'm using lme in the nlme package and I've found that the
> order of the models actually gives opposite results. I'm not sure if
> this is supposed to be the case or if I have missed something
> important, and I can't find anything in the Pinheiro and Bates book
> or in ?anova, or in Google for that matter which unfortunately only
> returns results about ANOVA which isn't much help. I'm using the
> latest version of R and nlme, just checked both.
> 
> Here is the code and output:
> 
> > PHQmodel1=lme(PHQ~Age+Gender+Date*Treatment, data=compfinal,
> random=~1|Case, na.action=na.omit)
> >
> > PHQmodel2=lme(PHQ~Age+Gender+Date*Treatment, data=compfinal,
> random=~1|Case, na.action=na.omit,
> +              correlation=corAR1(form=~Date|Case))
> 
> > anova(PHQmodel1, PHQmodel2) # accept model 2
>                     Model df      AIC      BIC    logLik   Test
> L.Ratio p-value
> PHQmodel1     1  8 48784.57 48840.43 -24384.28
> PHQmodel2     2  9 48284.68 48347.51 -24133.34 1 vs 2 501.8926 <.0001
> 
> > PHQmodel1=lme(PHQ~Age+Gender+Date*Treatment, data=compfinal,
> random=~1|Case, na.action=na.omit,
> +              correlation=corAR1(form=~Date|Case))
> >
> > PHQmodel2=lme(PHQ~Age+Gender+Date*Treatment, data=compfinal,
> random=~1|Case, na.action=na.omit)
> 
> > anova(PHQmodel1, PHQmodel2) # accept model 2
>                      Model df      AIC      BIC    logLik   Test
> L.Ratio p-value
> PHQmodel1     1  9 48284.68 48347.51 -24133.34
> PHQmodel2     2  8 48784.57 48840.43 -24384.28 1 vs 2 501.8926 <.0001
> 
> In both cases I am led to accept model 2 even though they are
> opposite models. Is it really just that you have to put them in the
> right order? It just seems like if there were say four models you
> wouldn't necessarily be able to determine the correct order.
> 
> Many thanks,
> Chris Beeley, Institute of Mental Health, UK
> 
> ...session info follows
> 
> > sessionInfo()
> R version 2.15.0 (2012-03-30)
> Platform: i386-pc-mingw32/i386 (32-bit)
> 
> locale:
> [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
> Kingdom.1252
> [3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
> [5] LC_TIME=English_United Kingdom.1252
> 
> attached base packages:
> [1] grid      stats     graphics  grDevices utils     datasets
> methods   base
> 
> other attached packages:
>  [1] gridExtra_0.9      RColorBrewer_1.0-5 car_2.0-12
> nnet_7.3-1         MASS_7.3-17
>  [6] xtable_1.7-0       psych_1.2.4        languageR_1.4
> nlme_3.1-104       ggplot2_0.9.1
> 
> loaded via a namespace (and not attached):
>  [1] colorspace_1.1-1 dichromat_1.2-4  digest_0.5.2     labeling_0.1
> lattice_0.20-6   memoise_0.1
>  [7] munsell_0.3      plyr_1.7.1       proto_0.3-9.2
> reshape2_1.2.1   scales_0.2.1     stringr_0.6
> [13] tools_2.15.0
> 
> > packageDescription("nlme")
> Package: nlme
> Version: 3.1-104
> Date: 2012-05-21
> Priority: recommended
> Title: Linear and Nonlinear Mixed Effects Models
> Authors at R: c(person("Jose", "Pinheiro", comment = "S version"),
> person("Douglas", "Bates", comment =
>            "up to 2007"), person("Saikat", "DebRoy", comment = "up
> to 2002"), person("Deepayan",
>            "Sarkar", comment = "up to 2005"), person("R-core", email
> = "R-core at R-project.org", role =
>            c("aut", "cre")))
> Author: Jose Pinheiro (S version), Douglas Bates (up to 2007),
> Saikat DebRoy (up to 2002), Deepayan
>            Sarkar (up to 2005), the R Core team.
> Maintainer: R-core <R-core at R-project.org>
> Description: Fit and compare Gaussian linear and nonlinear
> mixed-effects models.
> Depends: graphics, stats, R (>= 2.13)
> Imports: lattice
> Suggests: Hmisc, MASS
> LazyLoad: yes
> LazyData: yes
> License: GPL (>= 2)
> BugReports: http://bugs.r-project.org
> Packaged: 2012-05-23 07:28:59 UTC; ripley
> Repository: CRAN
> Date/Publication: 2012-05-23 07:37:45
> Built: R 2.15.0; x86_64-pc-mingw32; 2012-05-29 12:36:01 UTC; windows
> 
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
> 

-- 
Albyn Jones
Reed College
jones at reed.edu