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

[Ferenci Tamas]

Dear Ben,

In addition to that, if you run that script with age, it will be
correct... *regardless* of the starting value!

> The alarming part is the lack of any warning.

Yes, what I found suspicious was that it was _exactly_ 0.2... and it
was only accidental that I remembered that it is the default value in
corCAR1, so I realized what is going on. But I can imagine that it
might go unnoticed...

Tamas


2018. január 25., 18:57:25, írtad:

> The alarming part is the lack of any warning.  I don't know whether
> there's any easy way to detect this case, though ...  FWIW the
> starting value has to be above about 0.62 before it works

> ff <- function(rhostart) {
>   m <- lme(distance ~ agemos + factor(Sex),random = ~ 1 | Subject,
>     cor=corCAR1(rhostart,form=~agemos|Subject),
>     data = Orthodont)
>   return(coef(m$modelStruct$corStruct,unconstrained=FALSE))
> }

> startvec <- seq(0.02,0.98,by=0.02)
> estvec <- sapply(startvec,ff)

> plot(startvec,estvec)

> On Thu, Jan 25, 2018 at 12:40 PM, Fox, John <jfox at mcmaster.ca> wrote:
>> 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
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models