Skip to content

Extracting the coefficients from corStruct objects

6 messages · Markus Jäntti, Nilsson Fredrik X, Douglas Bates

Original
Markus Jäntti
# 
Dear All:

I have been trying to get the coefficients from corAR1 objects in lme objects in 
then lme package (3.1-80), but coef(obj$modelStruct$corStruct) returns something 
rather different that summary(object). I have tried to follow the various 
summary. methods in the source code, but have kind of gotten lost. Would be 
greatful if someone could point to a simple solution.

It would be helpful to manually extract the AR(1) coefficients, because I am 
data sets consisting of very large cohorts of individuals. I extract the 
information I need from each estimated object and store them in small matrices. 
I have not been able to figure out a way to extract the AR(1) parameters, however.

I include example code below, which produces the following results:
 > summary(fm2)
<snip>
Correlation Structure: AR(1)
  Formula: ~year | ind
  Parameter estimate(s):
   Phi
0.424
<snip>

but

 > coef(fm2$modelStruct$corStruct)
[1] 0.906

Version info for nlme:

 > packageDescription("nlme")
Package: nlme
Version: 3.1-80
Date: 2007-03-30
Priority: recommended
Title: Linear and Nonlinear Mixed Effects Models
Author: Jose Pinheiro <Jose.Pinheiro at pharma.novartis.com>, Douglas
         Bates <bates at stat.wisc.edu>, Saikat DebRoy
         <saikat at stat.wisc.edu>, and Deepayan Sarkar
         <deepayan at stat.wisc.edu>
Maintainer: R-core <R-core at R-project.org>
Description: Fit and compare Gaussian linear and nonlinear
         mixed-effects models.
Depends: graphics, stats, R (>= 2.3.0)
Imports: lattice
LazyLoad: yes
LazyData: yes
License: GPL version 2 or later
Packaged: Fri Mar 30 07:37:02 2007; ripley
Built: R 2.5.0; i486-pc-linux-gnu; 2007-04-25 03:13:00; unix

Example code:
<code>
library(nlme)
## number of time periods
years <- 5
## individuals
n <- 100

## the random effects
a <- as.vector(sapply(rnorm(n), rep, years))

## indicators for individuals
ind <- as.vector(sapply(1:n, rep, years))
year <- rep(1:years, n)

## the white noise
v <- rnorm(n*years)
## the AR(1) parameter
rho <- .5
## generate the AR1 erroer
## randomly draw the starting value of the AR(1) error for each unit
u1 <- rnorm(n)
## store values here:
u <- numeric(n*years)
## a counter:
k <- 1
## there has to be a much more elegant way to do this!
for(i in 1:n)
   {
     u[k] <- u1[i]
     for(j in 2:years)
       {
         k <- k+1
         u[k] <- rho*u[k-1] + v[k]
       }
     k <- k+1
   }

## the covariate
x <- runif(n*years)
## the intercept and slope
b0 <- 10
b1 <- 4
## and the dependent variable
y <- b0 + b1*x + a + u

test.d <- data.frame(y=y, x=x, ind=ind, year=year)
rm("x", "y", "ind", "years")

fm1 <- lme(y ~ x, random = ~ 1 | ind,
            data=test.d)
summary(fm1)
fm2 <- lme(y ~ x, random = ~ 1 | ind, correlation=corAR1(form=~year),
            data=test.d)
summary(fm2)
coef(fm2$modelStruct$corStruct)
</code>
Nilsson Fredrik X
# 
Dear Markus,

I don't realize why one gets these strange results and I would like to understand this too!

But here's a solution while waiting from the enlightenment:

Look at:
intervals(ModelwAR1.lme)$corStruct

And then 

intervals(ModelwAR1.lme)$corStruct[2]

is what you look for.

Cheers,

Fredrik

-----Ursprungligt meddelande-----
Fr?n: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] F?r Markus J?ntti
Skickat: den 9 maj 2007 23:49
Till: R-SIG-Mixed-Model
?mne: [R-sig-ME] Extracting the coefficients from corStruct objects

Dear All:

I have been trying to get the coefficients from corAR1 objects in lme objects in 
then lme package (3.1-80), but coef(obj$modelStruct$corStruct) returns something 
rather different that summary(object). I have tried to follow the various 
summary. methods in the source code, but have kind of gotten lost. Would be 
greatful if someone could point to a simple solution.

It would be helpful to manually extract the AR(1) coefficients, because I am 
data sets consisting of very large cohorts of individuals. I extract the 
information I need from each estimated object and store them in small matrices. 
I have not been able to figure out a way to extract the AR(1) parameters, however.

I include example code below, which produces the following results:
 > summary(fm2)
<snip>
Correlation Structure: AR(1)
  Formula: ~year | ind
  Parameter estimate(s):
   Phi
0.424
<snip>

but

 > coef(fm2$modelStruct$corStruct)
[1] 0.906

Version info for nlme:

 > packageDescription("nlme")
Package: nlme
Version: 3.1-80
Date: 2007-03-30
Priority: recommended
Title: Linear and Nonlinear Mixed Effects Models
Author: Jose Pinheiro <Jose.Pinheiro at pharma.novartis.com>, Douglas
         Bates <bates at stat.wisc.edu>, Saikat DebRoy
         <saikat at stat.wisc.edu>, and Deepayan Sarkar
         <deepayan at stat.wisc.edu>
Maintainer: R-core <R-core at R-project.org>
Description: Fit and compare Gaussian linear and nonlinear
         mixed-effects models.
Depends: graphics, stats, R (>= 2.3.0)
Imports: lattice
LazyLoad: yes
LazyData: yes
License: GPL version 2 or later
Packaged: Fri Mar 30 07:37:02 2007; ripley
Built: R 2.5.0; i486-pc-linux-gnu; 2007-04-25 03:13:00; unix

Example code:
<code>
library(nlme)
## number of time periods
years <- 5
## individuals
n <- 100

## the random effects
a <- as.vector(sapply(rnorm(n), rep, years))

## indicators for individuals
ind <- as.vector(sapply(1:n, rep, years))
year <- rep(1:years, n)

## the white noise
v <- rnorm(n*years)
## the AR(1) parameter
rho <- .5
## generate the AR1 erroer
## randomly draw the starting value of the AR(1) error for each unit
u1 <- rnorm(n)
## store values here:
u <- numeric(n*years)
## a counter:
k <- 1
## there has to be a much more elegant way to do this!
for(i in 1:n)
   {
     u[k] <- u1[i]
     for(j in 2:years)
       {
         k <- k+1
         u[k] <- rho*u[k-1] + v[k]
       }
     k <- k+1
   }

## the covariate
x <- runif(n*years)
## the intercept and slope
b0 <- 10
b1 <- 4
## and the dependent variable
y <- b0 + b1*x + a + u

test.d <- data.frame(y=y, x=x, ind=ind, year=year)
rm("x", "y", "ind", "years")

fm1 <- lme(y ~ x, random = ~ 1 | ind,
            data=test.d)
summary(fm1)
fm2 <- lme(y ~ x, random = ~ 1 | ind, correlation=corAR1(form=~year),
            data=test.d)
summary(fm2)
coef(fm2$modelStruct$corStruct)
</code>
1 day later
Douglas Bates
#
On 5/10/07, Nilsson Fredrik X <Fredrik.X.Nilsson at skane.se> wrote:
The reason for the strange parameterization is that the nlme package
always uses unconstrained parameters for the optimization.  When the
parameters have natural constraints it is necessary to transform them
to produce unconstrained versions.
Nilsson Fredrik X
# 
It would be interesting to know what the transformation is or how to find out. Is it tanh(x/2)?

Best regards,
Fredrik
Douglas Bates
#
On 5/11/07, Nilsson Fredrik X <Fredrik.X.Nilsson at skane.se> wrote:
Just look at the code for corAR1. The unconstrained parameter is
log((1+r)/(1-r)).
Markus Jäntti
#
Douglas Bates wrote:
So for the AR1, the transform is  "Fishers z-tranform"?

I was trying to figure this out by looking at the summary methods. Embarasingly, 
it never occurred to me to look in the corAR1 code itself.

Thanks for the help to both Fredrik and Doug.

Markus