bug or numerical problem in nlme's corCAR1

7 messages · Ferenci Tamas, John Fox, Ben Bolker

Original

1

7

John Fox

Thu, Jan 25, 2018 9:40 AM #

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:

+     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

+     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

[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/

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Ben Bolker

Thu, Jan 25, 2018 9:57 AM #

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

John Fox

Thu, Jan 25, 2018 10:18 AM #

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

Thu, Jan 25, 2018 1:16 PM #

Dear Ben,

In addition to that, if you run that script with age, it will be
correct... *regardless* of the starting value!

Yes, what I found suspicious was that it was _exactly_ 0.2... and it
was only accidental that I remembered that it is the default value in
corCAR1, so I realized what is going on. But I can imagine that it
might go unnoticed...

Tamas


2018. janu?r 25., 18:57:25, ?rtad:

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

Thu, Jan 25, 2018 1:17 PM #

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

John Fox

Thu, Jan 25, 2018 2:43 PM #

Dear Tamas,

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

Ben Bolker

Thu, Jan 25, 2018 4:19 PM #

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