A possible solution would be to see whether the estimated
autoregressive parameter is stuck at the start value and then try
something like a bisection search. That's pretty crude and I bet there's a smarter way to do it.
-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com]
Sent: Thursday, January 25, 2018 12:57 PM
To: Fox, John <jfox at mcmaster.ca>
Cc: Ferenci Tamas <tamas.ferenci at medstat.hu>; r-sig-mixed-models at r-
project.org
Subject: Re: [R-sig-ME] 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
[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
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
return a coefficient of 0.1234 and so on.)
What's going on...?
Thank you in advance,
Tamas Ferenci