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

[Ferenci Tamas]

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!

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