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

[Fox, John]

Hi Ben,

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.

Best,
 John

> -----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
> >
> >> .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