[R-sig-ME] Extracting the coefficients from corStruct objects
Markus Jäntti
markus.jantti at iki.fi
Wed May 9 23:49:21 CEST 2007
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
More information about the R-sig-mixed-models
mailing list