[R-sig-ME] bug or numerical problem in nlme's corCAR1
Ben Bolker
bbolker at gmail.com
Thu Jan 25 18:57:25 CET 2018
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
More information about the R-sig-mixed-models
mailing list