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

Fox, John jfox at mcmaster.ca
Thu Jan 25 23:43:46 CET 2018 

Dear Tamas,

> -----Original Message-----
> From: Ferenci Tamas [mailto:tamas.ferenci at medstat.hu]
> Sent: Thursday, January 25, 2018 4:18 PM
> To: Fox, John <jfox at mcmaster.ca>; Ben Bolker <bbolker at gmail.com>
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] bug or numerical problem in nlme's corCAR1
> 
> Dear John,
> 
> That's a possible solution IF you realize something is wrong, but what I find
> problematic is that one might not realize at all that such measures are needed!

I was suggesting to Ben that lme() might do it.

John

> 
> Tamas
> 
> 
> 2018. január 25., 19:18:40, írtad:
> 
> > 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

More information about the R-sig-mixed-models mailing list