bug or numerical problem in nlme's corCAR1
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