[R-sig-ME] Extracting the coefficients from corStruct objects

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

-- 
Markus Jantti
Abo Akademi University
markus.jantti at iki.fi
http://www.iki.fi/~mjantti