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

Chris Beeley chris.beeley at gmail.com
Thu May 31 13:47:30 CEST 2012


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



More information about the R-help mailing list