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

[Ben Bolker]

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