bug or numerical problem in nlme's corCAR1

The alarming part is the lack of any warning.

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