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

[Markus Jäntti]

Douglas Bates wrote:
> 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)).

So for the AR1, the transform is  "Fishers z-tranform"?

I was trying to figure this out by looking at the summary methods. Embarasingly, 
it never occurred to me to look in the corAR1 code itself.

Thanks for the help to both Fredrik and Doug.

Markus

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


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