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

[Douglas Bates]

On 5/10/07, Nilsson Fredrik X <Fredrik.X.Nilsson at skane.se> wrote:
> Dear Markus,

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

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.

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