[R] Multilevel model ("lme") question

[Douglas Bates]

On 10/21/06, Lukas Rode <lukas.rode at googlemail.com> wrote:
> Dear list,
>
> I'm trying to fit a multilevel (mixed-effects) model using the lme function
> (package nlme) in R 2.4.0. As a mixed-effects newbie I'm neither sure about
> the modeling nor the correct R syntax.

You may also want to consider using the lmer function from package
lme4.  It's a later version of lme for R and is generally faster than
lme.  I'll show the syntax for both.

> My data is structured as follows: For each subject, a quantity Y is measured
> at a number (>= 2) of time points. Moreover, at time point 0 ("baseline"), a
> quantity X is measured for each subject (I am interested to see whether X,
> or log(X), is a predictor of Y). I saw in some examples that time-invariant
> predictors should be repeated for all rows of a subject, which is why I
> copied the baseline value of X also to time points > 0.

That's correct.  For both lme and lmer the data need to be in the
'long' form, also called the 'person-period' form.  You can reshape
the data by hand, as you have done, or you can use the 'reshape'
function or the more general capabilities in Hadley Wickham's reshape
package.

> The resulting data frame looks like this:

> Grouped Data: Y ~ TIME | Subject
> Y         TIME        Subject         X.Baseline
> 9           0.0             1                1084
> 7           3.7             1                1084
> 11         0.0             2                7150
> 8           9.2             2                7150

> Intra-subject trajectories of Y very close to linear. I'd like to check
> whether slope (and maybe also offset) of this line are (in part) predicted
> by X.baseline.

I would say that this question would be addressed as part of the
fixed-effects specification.  You would be comparing a model with main
effects for TIME and X.Baseline and their interaction to various
submodels.  In the S language formula notation these fixed effects
specifications are

 Y ~ TIME * X.Baseline  # main effects for TIME and X.Baseline and
their interaction
 Y ~ TIME + X.Baseline + TIME:X.Baseline # equivalent to first model
 Y ~ TIME + X.Baseline # main effects but no interaction.  That is
X.Baseline does not affect slope
 Y ~ TIME + TIME:X.Baseline # X.Baseline affects slope (wrt TIME) but not offset
 Y ~ TIME # X.Baseline has no effect on Y

> Thus, I think the multilevel model specification should be as follows (i =
> subject, j=measurement):
> y_ij = \beta_i  + b_i * TIME_ij + \epsilon_ij,
> with
> b_i = \zeta_i0 + \zeta_i1 * X.Baseline
> Is this correct?

Well that model doesn't have a main effect for X.Baseline so you can't
test for an effect for X.Baseline on the slope.

> Now, I am completely unsure how to "translate" this into the syntax needed
> by lme.
> Is there any standard procedure on how to get from e.g. the Laird&Ware'82
> matrix model notation to the lme input?

I think that is more of a question of how to use the model formula
language than a question about lme.  In lme the right hand side of the
formula argument generates Laird&Ware's matrix X and the random
specification generates the matrix Z.

> And, in my case, is the correct model as follows, or am I wrong?
> my.model <- lme( Y ~ TIME, my.data.frame, random=~X.Baseline | Subject )
> [in case this is correct, it is by accident, since I made it up from some
> examples I found without really understanding the translation from the model
> formulation]

Actually this model specification should throw an error but it doesn't.

Estimates for the random effects in this model are not well defined
because X.Baseline is constant within Subject.  Thus within Subject
the Intercept term and the X.Baseline term are confounded.  I would
begin with the model

lme(Y ~ TIME * X.Baseline, my.data.frame, random = ~ 1|Subject)

which is equivalent to

lmer(Y ~ TIME * X.Baseline + (1|Subject), my.data.frame)

Then examine the estimates for the fixed-effects coefficients and
their standard errors to see what reductions of this model may be
appropriate.

> I'd greatly appreciate some advice or help.

I hope this helps.