Extracting the coefficients from corStruct objects

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.

It would be interesting to know what the transformation is or how to find out. Is it tanh(x/2)?

Best regards,
Fredrik

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>

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