[R] AR1 covariance structure for lme object and R/SAS differences in model output

[Ben Bolker]

anord <andreas.nord <at> biol.lu.se> writes:

> 


  [snip snip]

> We are working on a data set in which we have measured repeatedly a
> physiological response variable (y) 
> every 20 min for 12 h (time variable; 'x') in subjects ('id') beloning to
> one of five groups ('group'; 'A' to 'E'). Data are located at:
> https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0
> 
> We are interested to model if the response in y differences 
> with time (i.e.
> 'x') for the two groups. Thus:
> require(nlme)
> m1<-lme(y~group*x+group*I(x^2),random=~x|id,
>   data=data.df,na.action=na.omit)
> 
> But because data are collected repeatedly over 
> short time intervals for each
> subject, it seemed prudent to consider an autoregressive covariance
> structure. Thus:
> m2<-update(m1,~.,corr=corCAR1(form=~x|id))
> 
> AIC values indicate the latter (i.e. m2) as more appropriate:
>   anova(m1,m2)
> #   Model df      AIC      BIC       logLik        Test  L.Ratio     
> p-value
> #m1     1 19 2155.996 2260.767 -1058.9981                        
> #m2     2 20 2021.944 2132.229  -990.9718 1 vs 2 136.0525  <.0001
> 
> Fixed effects and test statistics differ between models.
>  A look at marginal
> ANOVA tables suggest inference might differ somewhat between models:
> 
> anova.lme(m1,type="m")
> #              numDF denDF  F-value p-value
> #(Intercept)      1  1789 63384.80  <.0001
> #group             4    45      1.29  0.2893
> #x                   1  1789     0.05  0.8226
> #I(x^2)            1  1789     4.02  0.0451
> #group:x          4  1789     2.61  0.0341
> #group:I(x^2)   4  1789     4.37  0.0016
> 
> anova.lme(m2,type="m")
> #             numDF denDF  F-value p-value
> #(Intercept)      1  1789 59395.79  <.0001
> #group             4    45      1.33  0.2725
> #x                   1  1789     0.04  0.8379
> #I(x^2)            1  1789     2.28  0.1312
> #group:x          4  1789     2.09  0.0802
> #group:I(x^2)   4  1789     2.81  0.0244
> 
> Now, this is all well. But: my colleagues have been running 
> the same data
> set using PROC MIXED in SAS and come up with 
> substantially different results
> when comparing SAS default covariance structure (variance components) and
> AR1. Specifically, there is virtually no change 
> in either test statistics or
> fitted values when using AR1 instead of Variance Components in SAS, which
> fits the observation that AIC values (in SAS) indicate both covariance
> structures fit data equally well. 
> 
> This is not very satisfactory to me, and I would be 
> interesting to know what
> is happening here. Realizing
> this might not be the correct forum for this question, I would like to ask
> you all if anyone would have any
> input as to what is going on here, e.g. am I setting up my model
> erroneously, etc.? 
> 
> N.b. I have no desire to replicate SAS results, but I would most certainly
> be interested to know what could possibly explain  
> such a large discrepancy
> between the two platforms. Any suggestions greatly welcomed.
> 
> (Data are located at:
> https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0)

  Could you repost this on r-sig-mixed-models at r-project.org ?

It would be useful to see some of the comparisons (fixed effects,
RE variance-covariance, correlation parameter) between SAS and
R; are they actually fitting the same model?  (e.g., does SAS
allow for covariance between the slope and intercept random effects?)
Maybe they're getting the same estimates but computing df/p-values
in different ways?

  I thought I would try this with orthogonal polynomials in case
some of the fits were unstable ...

data.df <- read.csv2("ar1_data.csv")
library("nlme")
m1 <- lme(y~group*x+group*I(x^2),random=~x|id,
          data=data.df,na.action=na.omit,method="ML")
## use method="ML" so we can compare orthog. and non-orthog. polynomials
m1B <- update(m1,.~group*poly(x,2))
m2 <- update(m1,corr=corCAR1(form=~x|id))
m2B <- update(m1B,corr=corCAR1(form=~x|id))
AIC(m1,m1B,m2,m2B)