[R] Mixed model question.

From: Rolf Turner <r.turner_at_auckland.ac.nz>
Date: Mon, 28 Jul 2008 14:06:13 +1200

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. How can I (is there any way that I can) tell lmer() to fit the most general possible covariance structure?

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.

Thanks.

                cheers,

                        Rolf Turner

######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

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. Received on Mon 28 Jul 2008 - 02:13:42 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 28 Jul 2008 - 23:32:49 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive