[R-sig-ME] bug or numerical problem in nlme's corCAR1

[Fox, John]

Dear Tamas Farenci,

You're mistaken: For the continuous first-order AR process, the unit of time matters. Measuring time in months implies a much larger autocorrelation at lag 1 than measuring time in years. As it turns out, 0.2 is a poor start values for the former:

> library(nlme)
> lme(distance ~ age + factor(Sex),random = ~ 1 | Subject,
+     cor=corCAR1(form=~age|Subject),data = Orthodont)
Linear mixed-effects model fit by REML
  Data: Orthodont 
  Log-restricted-likelihood: -218.6984
  Fixed: distance ~ age + factor(Sex) 
      (Intercept)               age factor(Sex)Female 
       17.7214161         0.6594049        -2.3274848 

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    1.788899 1.454494

Correlation Structure: Continuous AR(1)
 Formula: ~age | Subject 
 Parameter estimate(s):
      Phi 
0.2418536 
Number of Observations: 108
Number of Groups: 27 

> Orthodont$agemos <- Orthodont$age*12
> lme(distance ~ agemos + factor(Sex),random = ~ 1 | Subject,
+     cor=corCAR1(0.8, form=~agemos|Subject),data = Orthodont)
Linear mixed-effects model fit by REML
  Data: Orthodont 
  Log-restricted-likelihood: -221.1833
  Fixed: distance ~ agemos + factor(Sex) 
      (Intercept)            agemos factor(Sex)Female 
      17.72141619        0.05495041       -2.32748480 

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    1.788899 1.454494

Correlation Structure: Continuous AR(1)
 Formula: ~agemos | Subject 
 Parameter estimate(s):
      Phi 
0.8884427 
Number of Observations: 108
Number of Groups: 27 

> .8884427^12
[1] 0.2418539

Thus, the two solutions agree. I agree, however, that it's slightly disconcerting that lme() can't overcome the poor start value.

I hope this helps,
 John

-----------------------------
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
Web: socialsciences.mcmaster.ca/jfox/




> -----Original Message-----
> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-
> project.org] On Behalf Of Ferenci Tamas
> Sent: Thursday, January 25, 2018 12:05 PM
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] bug or numerical problem in nlme's corCAR1
> 
> Dear list members,
> 
> Consider the following example (yes, I haven't centered age, corCAR1 is not
> necessarily needed here etc., so it is perhaps not the most meaningful model,
> I use it just to show the problem):
> 
> lme(distance ~ age + factor(Sex),random = ~ 1 | Subject,
>    cor=corCAR1(form=~age|Subject),data = Orthodont)
> 
> Everything works perfectly.
> 
> Let's now multiply age, say, we measure it in months:
> 
> Orthodont$agemos <- Orthodont$age*12
> lme(distance ~ agemos + factor(Sex),random = ~ 1 | Subject,
>     cor=corCAR1(form=~agemos|Subject),data = Orthodont)
> 
> It shouldn't make any difference, but check the autocorrelation coefficient!
> It changed from 0.2418536 to 0.2... i.e. it stuck at its default value, as if it
> was not optimized at all! (One can verify that it is indeed the case, and 0.2
> was not a coincidence by calling lme(distance ~ agemos +
> factor(Sex),random = ~ 1 | Subject,
> cor=corCAR1(0.1234,form=~agemos|Subject),data = Orthodont) which will
> return a coefficient of 0.1234 and so on.)
> 
> What's going on...?
> 
> Thank you in advance,
> Tamas Ferenci
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models