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

[Chris Beeley]

Well that's that cleared up then. Thanks to all.

Chris B.

On 31/05/2012 17:51, Albyn Jones wrote:
> 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.
>>