[R] extracting estimated covariance parameters from lme fit

[Kingsford Jones]

The VarCorr function will extract the components of the random effects
covariance matrix, but note the quirk that it returns values as
characters:

library(nlme)
f1 <- lme(distance ~ age, data = Orthodont, random = ~1 + age|Subject)
(vc <- VarCorr(f1))
# Subject = pdLogChol(1 + age)
#             Variance   StdDev    Corr
# (Intercept) 5.41508724 2.3270340 (Intr)
# age         0.05126955 0.2264278 -0.609
# Residual    1.71620401 1.3100397

str(vc)
#  'VarCorr.lme' chr [1:3, 1:3] "5.41508724" "0.05126955" "1.71620401" ...
# - attr(*, "dimnames")=List of 2
#   ..$ : chr [1:3] "(Intercept)" "age" "Residual"
#   ..$ : chr [1:3] "Variance" "StdDev" "Corr"
#  - attr(*, "title")= chr "Subject = pdLogChol(1 + age)"

(sigma2.age <- as.numeric(vc[2, 1]))
# [1] 0.05126955


hth,

Kingsford Jones

On Mon, Nov 16, 2009 at 9:25 AM, Green, Gerwyn (greeng6)
<g.green8 at lancaster.ac.uk> wrote:
> Dear all
>
> Apologies in advance as this seems like a trivial question. Nonetheless,
> a question I haven't been able to resolve myself !. Within a single
> repetition of a simulation (to be repeated many times) I am fitting the
> following linear mixed model using lme...
>
> Y_{gtr} = \mu + U_{g} +  W_{gt} + Z_{gtr}
>
> U_{g} ~ N(0,\gamma^{2}), W_{gt} ~ N(0,\kappa^{2}), Z_{gtr} ~
> N(0,\tau^{2})
>
> g = 1,...,G
> t = 1,...,T
> r= 1,...,R
>
>
> ...by doing
>
>> Model.fit <- lme(Y ~ 1, data=data, random= ~1|gene/treatment)
>
> I would like to be able to extract the estimated covariance parameters
> contained within the lme object. I know if I type...
>
>> Model.fit$sigma
>
>
> ...then I get the estimated residual variance, i.e. within the context
> of the above model, the estimate for \tau. But I would also like to
> extract the estimates for \gamma and \kappa by doing
> Model.fit$"something". I am aware that I can view the output using the
> extractor function "summary", but within a single repetition of my
> simulation routine I want to be able to code something like
>
> gamma <- Model.fit$.....
> kappa <- Model.fit$.....
>
>
> and then plug `gamma' and `kappa' into some formulae. This process of
> fitting and extracting will be repeated many times, which is why I wish
> to automate everything.
>
> Again, any help would be greatly appreciated
>
> Best
>
> Gerwyn Green
> School of Health and Medicine
> Lancaster University
>
>
>
>
>
>
> Any help would be greatly appreciated
>
> Best
>
>
> Gerwyn Green
> School of Health and Medicine
> Lancaster Uinversity
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>