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

[Douglas Bates]

On 5/11/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.

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

Just look at the code for corAR1. The unconstrained parameter is
log((1+r)/(1-r)).

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