[R] PROC MIXED user trying to use (n)lme...

[Douglas Bates]

Søren Højsgaard <sorenh at mail.dk> writes:

> Coming from a proc mixed (SAS) background I am trying to get into
> the use of (n)lme.

You may want to check the package SASmixed which provides lme analyses
for several of the examples in the book "SAS System for Mixed Models"
by Littel, Milliken, Stroup and Wolfinger.  Studying those examples
done in both PROC MIXED and lme can ease the transition.

> In this connection, I have some (presumably stupid) questions
> which I am sure someone out there can answer:
> 
> 1) With proc mixed it is easy to get a hold on the estimated
> variance parameters as they can be put out into a SAS data set.
> How do I do the same with lme-objects? For example, I can see the
> estimated variance matrix when typing the name of the lme-object,
> but that does not give me a "handle" on the matrix.

The VarCorr function returns a table of variance estimates, the
corresponding standard deviations, and the estimated correlations.
The value of VarCorr is a character matrix, not a numeric matrix.  If
you want the numeric matrix you need to do a bit of digging.  The
fitted model contains a "modelStruct" component which, in turn,
contains an "reStruct" component.  Each element of the reStruct can be
converted to a scaled variance-covariance matrix by as.matrix.  To
remove the scaling you need to multiply by sigma^2.  This is shown below.

 > library(nlme)
 Loading required package: nls 
 > data(Oxboys)
 > fm1 <- lme(height ~ age, data = Oxboys)
 > VarCorr(fm1)
 Subject = pdSymm(age) 
             Variance   StdDev   Corr  
 (Intercept) 65.3038110 8.081077 (Intr)
 age          2.8248080 1.680717 0.641 
 Residual     0.4354534 0.659889       
 > str(VarCorr(fm1))
  chr [1:3, 1:3] "65.3038110" " 2.8248080" " 0.4354534" "8.081077" ...
  - attr(*, "dimnames")=List of 2
   ..$ : chr [1:3] "(Intercept)" "age" "Residual"
   ..$ : chr [1:3] "Variance" "StdDev" "Corr"
  - attr(*, "title")= chr "Subject = pdSymm(age)"
  - attr(*, "class")= chr "VarCorr.lme"
 > lapply(lapply(fm1$modelStruct$reStruct, as.matrix), function(x) x * fm1$sigma^2)
 $Subject
             (Intercept)      age
 (Intercept)   65.303811 8.709806
 age            8.709806 2.824808

> 2) I am trying to fit a mixed model, say
> 
>     y_{ij} = X_{ij}\beta + U_i + e_{ij}
> 
> where Y_{ij}, U_i and e_{ij} all are D-dimensional, and U_i ~
> N_D(0,G), and e_{ij} ~ N(0,R). This is easy to do in proc mixed by
> including a factor keeping track of how the elements of y_{ij} are
> "stacked" on top of each other. This factor is then used in the
> model, the random and the repeated statements. In the book by
> Pinheiro and Bates (p.28 ff) there is a similar example. Yet, I
> don't quite get the "generality" of the example. Can anyone help
> with that?

I don't happen to have our book handy (I'm writing this from Belgium
where I am attending a conference) so I can't reply regarding the
example on p. 28.  In general the random statement in PROC MIXED is
easy to translate into nlme.  The repeated statement isn't.

The basic approach in lme is to define one or more factors that
represent the grouping of the observations.  In the example above the
data contain the information that the grouping is by the "Subject"
factor

 > formula(Oxboys)
 height ~ age | Subject

A full specification of the model in fm1 is
 > fm2 <- lme(fixed = height ~ age, data = Oxboys, random = ~ age | Subject)
 > VarCorr(fm2)
 Subject = pdLogChol(age) 
             Variance   StdDev   Corr  
 (Intercept) 65.3038159 8.081078 (Intr)
 age          2.8248081 1.680717 0.641 
 Residual     0.4354534 0.659889       

so the index i (or perhaps ij?) in your notation corresponds to the
levels of Subject.  We tend to express the model in the Laird-Ware
notation of a vector of responses y_i for the i'th subject along with
model matrices X_i and Z_i and a random within-subject noise vector
eps_i.  This gives y_i = X_i \beta + Z_i b_i + eps_i I'm not sure how
that correponds to your U_i in that I don't see a model matrix
associated with U_i.

It happens that I am giving a presentation on "Mixed-effects Models in
Practice" here in Belgium.  I put the slides in PDF format on
http://nlme.stat.wisc.edu/BSS2001.pdf (A4 sized) and as "4-up" copies
(USletter sized) on http://nlme.stat.wisc.edu/BSS2001_4up.pdf.  Those
slides contain examples of the model matrices associated with this
example.

> 3) With proc mixed it is easy to specify that some covariance
> parameters are known (fixed) in advance. How can I do the same
> with lme? (To be specific, I would like to specify that the
> elements of R in the example above are known).

You need to check for the available correlation structures in the
corStruct family and what their arguments are.  If your correlation
structure is not among those listed you would need to code it.

I hope this helps.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._