[R] Mixed model question.

[Douglas Bates]

On Sun, Jul 27, 2008 at 9:06 PM, Rolf Turner <r.turner at auckland.ac.nz> wrote:

> I continue to struggle with mixed models.  The square zero version
> of the problem that I am trying to deal with is as follows:
>
> A number (240) of students are measured (tested; for reading comprehension)
> on 6 separate occasions.  Initially (square zero) I want to treat the
> test time as a factor (with 6 levels).  The students are of course
> ``random effects''.  Later I want to look at reduced models, with the
> means at the 6 times following patterns:
>
>        (mu1, mu2, mu2, mu3, mu3, mu4)
>
> or
>
>        (mu1, mu1+theta, mu1+theta, mu1+2*theta, mu1+2*theta, mu1+3*theta)
>
> the idea begin that the tests are given ``in pairs'' at the beginning and
> end of the school year.  Progress is expected over the school year, no
> progress over the two intervening summers.  The summers fall between
> times 2 and 3 and between times 4 and 5.  The second of the two models
> assumes a constant amount of progress (``theta'') in each of three
> school years in which the students are observed.
>
> But the square zero model is just a trivial repeated measures model.
>
> The estimate of the means at the 6 observation times is just
> mu-hat = xbar <- apply(X,2,mean) where X is the ``data matrix'',
> and the estimate of the covariance structure is just given by
> Sigma-hat = S <- var(X).
>
> No problem so far.  I can also fit the two reduced models (I'm pretty
> sure) via maximum likelihood, using optim(), for example.  Assuming
> normal data.  Rash assumption, but that's not the issue here.
>
> But the thing is, I also want (later!) to include such things as
> a school effect (there are 6 different schools), a sex effect, and
> an ethnicity effect.
>
> Things start to get complicated --- sounds like a job for lmer().
>
> So I'd just like to get a toe in the door by fitting the trivial
> (square 0) model with lmer() --- and then if I can get my head
> around that, move on to the two reduced models.  I.e. I'd like to
> reproduce my simple-minded computations using lmer() --- which would
> give me a little bit of confidence that I'm driving lmer() correctly.
>
> *Can* the trivial model be fitted in lmer()?  I tried using
>
>        fit <- lmer(y ~ tstnum + (1|stdnt), data=schooldat)
>
> and got estimates of the coefficients for tstnum as follows:
>
> Fixed effects:
>            Estimate Std. Error t value
> (Intercept)  3.22917    0.09743   33.14
> tstnum2      0.46667    0.08461    5.52
> tstnum3      0.50000    0.08461    5.91
> tstnum4      0.83750    0.08461    9.90
> tstnum5      0.47083    0.08461    5.56
> tstnum6      0.97500    0.08461   11.52
>
> The mean of (the columns of) the data matrix is
>
> 3.229167 3.695833 3.729167 4.066667 3.700000 4.204167
>
> which is in exact agreement with the lmer() results when converted to
> the same parameterization (mu_i = mu + alpha_i, with alpha_1 = 0).
>
> (Notice the surprizing, depressing, and so far unexplained *drop*
> in the response over the second summer.)
>
> What I *don't* understand is the correlation structure of the estimates
> produced by lmer(), which is:
>
> Correlation of Fixed Effects:
>        (Intr) tstnm2 tstnm3 tstnm4 tstnm5
> tstnum2 -0.434
> tstnum3 -0.434  0.500
> tstnum4 -0.434  0.500  0.500
> tstnum5 -0.434  0.500  0.500  0.500
> tstnum6 -0.434  0.500  0.500  0.500  0.500
>
> So apparently the way I called lmer() places substantial constraints
> on the covariance structure.

That's the correlation matrix of the fixed-effects parameters.  You
should have separately gotten estimates of the variance-covariance of
the random effects, which you coyly did not show us :-).  Because you
are allowing only a simple, scalar random effect per student there
will be an estimate of the variance of this random effect and an
estimate of the residual variance.

> How can I (is there any way that I can)
> tell lmer() to fit the most general possible covariance structure?

It sounds like you want a model formula of

lmer(y ~ tstnum + (0 + tstnum|stdnt), data=schooldat)

but that model will have 21 variance-covariance terms to estimate (22
if you count the residual variance but that one gets profiled out of
the optimization).  I would not be surprised if the estimated
variance-covariance matrix for the random effects turns out to be
singular.

> As usual, advice, insight, tutelage humbly appreciated.

> If anyone wishes to experiment with the real data set (it's a bit
> too big to post here) I can make it available to them via email.

Generally I would jump at the chance but not with my "To Do" list in
its current, sadly over-committed, state.