[R-sig-ME] bug or numerical problem in nlme's corCAR1
Ben Bolker
bbolker at gmail.com
Fri Jan 26 01:19:30 CET 2018
Yes (although I'd like to remind everyone that I am *not* the
maintainer of nlme; R-core is).
cheers
Ben
On 18-01-25 05:43 PM, Fox, John wrote:
> 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