[R] [R-sig-ME] lme nesting/interaction advice
Douglas Bates
bates at stat.wisc.edu
Tue May 13 18:49:55 CEST 2008
On Mon, May 12, 2008 at 5:39 PM, Rolf Turner <r.turner at auckland.ac.nz> wrote:
>
> On 13/05/2008, at 4:09 AM, Douglas Bates wrote:
>
>> I'm entering this discussion late so I may be discussing issues that
>> have already been addressed.
>>
>> As I understand it, Federico, you began by describing a model for data
>> in which two factors have a fixed set of levels and one factor has an
>> extensible, or "random", set of levels and you wanted to fit a model
>> that you described as
>>
>> y ~ effect1 * effect2 * effect3
>>
>> The problem is that this specification is not complete.
>
> At *last* (as Owl said to Rabbit) we're getting somewhere!!!
>
> I always knew that there was some basic fundamental point
> about this business that I (and I believe many others) were
> simply missing. But I could not for the life of me get anyone
> to explain to me what that point was. Or to put it another
> way, I was never able to frame a question that would illuminate
> just what it was that I wasn't getting.
>
> I now may be at a stage where I can start asking the right
> questions.
>
>> An interaction of factors with fixed levels and a factor with random
>> levels can mean, in the lmer specification,
>>
>> lmer(y ~ effect1 * effect2 + (1| effect3) + (1|effect1:effect2:effect3),
>> ...)
>>
>> or
>>
>> lmer(y ~ effect1 * effect2 + (effect1*effect2 | effect3), ...)
>>
>> or other variations. When you specify a random effect or an random
>> interaction term you must, either explicitly or implicitly, specify
>> the form of the variance-covariance matrix associated with those
>> random effects.
>>
>> The "advantage" that other software may provide for you is that it
>> chooses the model for you but that, of course, means that you only
>> have the one choice.
>
> Now may I start asking what I hope are questions that will lift
> the fog a bit?
>
> Let us for specificity consider a three-way model with two
> fixed effects and one random effect from the good old Rothamstead
> style
> agricultural experiment context: Suppose we have a number of
> species/breeds of wheat (say) and a number of fertilizers.
> These are fixed effects. And we have a number of fields (blocks)
> --- a random effect. Each breed-fertilizer combination is
> applied a number of times in each field. We ***assume*** that
> that the field or block effect is homogeneous throughout. This
> may or may not be a ``good'' assumption, but it's not completely
> ridiculous and would often be made in practice. And probably
> *was* made at Rothamstead. The response would be something like
> yield in bushels per acre.
>
> The way that I would write the ``full'' model for this setting,
> in mathematical style is:
>
> Y_ijkl = mu + alpha_i + beta_j + (alpha.beta)_ij + C_k + (alpha.C)_ik
> + (beta.C)_jk + (alpha.beta.C)_ijk + E_ijkl
>
> The alpha_i and beta_j are parameters corresponding to breed and
> fertilizer
> respectively; the C_k are random effects corresponding to fields or
> blocks.
> Any effect ``involving'' C is also random.
>
> The assumptions made by the Package-Which-Must-Not-Be-Named are (I
> think)
> that
>
> C_k ~ N(0,sigma_C^2)
> (alpha.C)_ik ~ N(0,sigma_aC^2)
> (beta.C)jk ~ N(0,sigma_bC^2)
> (alpha.beta.C)_ijk ~ N(0,sigma_abC^2)
> E_ijkl ~ N(0,sigma^2)
>
> and these random variables are *all independent*.
>
> Ahhhhhhhh ... perhaps I'm on the way to answering my own question.
> Is
> it this assumption of ``all independent'' which is questionable? It
> seemed innocent enough when I first learned about this stuff, lo
> these
> many years ago. But .... mmmmmaybe not!
>
> To start with: What would be the lmer syntax to fit the foregoing
> (possibly naive) model? I am sorry, but I really cannot get my head
> around the syntax of lmer model specification, and I've tried. I
> really have. Hard. I know I must be starting from the wrong place,
> but I haven't a clue as to what the *right* place to start from is.
> And if I'm in that boat, I will wager Euros to pretzels that there
> are others in it. I know that I'm not the brightest bulb in the
> chandelier, but I'm not the dullest either.
Thanks for the questions, Rolf. I completely agree that mixed model
specification can be an extremely confusing area.
Let's consider a set of models for the Machines data reproduced (from
Milliken and Johnson) in Pinheiro and Bates and available in the MEMSS
package.
> library(lme4)
Loading required package: Matrix
Loading required package: lattice
Attaching package: 'Matrix'
The following object(s) are masked from package:stats :
xtabs
> data("Machines", package = "MEMSS")
> str(Machines)
'data.frame': 54 obs. of 3 variables:
$ Worker : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
$ Machine: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
$ score : num 52 52.8 53.1 51.8 52.8 53.1 60 60.2 58.4 51.1 ...
We consider the Machine factor to have a fixed set of levels in that
we only consider these three machines. The levels of the Worker
factor represent a sample from the set of potential operators. As you
might imagine from this description, I now think of the distinction
between "fixed" and "random" as being associated with the factor, not
necessarily the "effects".
If you plot these data
> dotplot(reorder(Worker, score) ~ score, Machines, groups = Machine, type = c("g", "p", "a"), xlab = "Efficiency score", ylab = "Worker", auto.key = list(columns = 3, lines = TRUE))
(resulting PDF enclosed) you will see evidence of an interaction.
That is, some workers have noticeably different score patterns on the
three machines than do others. Worker 6 on machine B is the most
striking example.
One way to model this interaction is to say that there is a random
effect for each worker and a separate random effect for each
worker/machine combination. If the random effects for the
worker/machine combinations are assumed to be independent with
constant variance then one expresses the model as
> print(fm1 <- lmer(score ~ Machine + (1|Worker) + (1|Machine:Worker), Machines), corr = FALSE)
Linear mixed model fit by REML
Formula: score ~ Machine + (1 | Worker) + (1 | Machine:Worker)
Data: Machines
AIC BIC logLik deviance REMLdev
227.7 239.6 -107.8 225.5 215.7
Random effects:
Groups Name Variance Std.Dev.
Machine:Worker (Intercept) 13.90963 3.72956
Worker (Intercept) 22.85528 4.78072
Residual 0.92464 0.96158
Number of obs: 54, groups: Machine:Worker, 18; Worker, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 52.356 2.486 21.062
MachineB 7.967 2.177 3.659
MachineC 13.917 2.177 6.393
An equivalent formulation is
> print(fm1 <- lmer(score ~ Machine + (1|Worker/Machine), Machines), corr = FALSE)
Linear mixed model fit by REML
Formula: score ~ Machine + (1 | Worker/Machine)
Data: Machines
AIC BIC logLik deviance REMLdev
227.7 239.6 -107.8 225.5 215.7
Random effects:
Groups Name Variance Std.Dev.
Machine:Worker (Intercept) 13.90963 3.72956
Worker (Intercept) 22.85528 4.78072
Residual 0.92464 0.96158
Number of obs: 54, groups: Machine:Worker, 18; Worker, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 52.356 2.486 21.062
MachineB 7.967 2.177 3.659
MachineC 13.917 2.177 6.393
The expression (1|Worker/Machine) is just "syntactic sugar". It is
expanded to (1|Worker) + (1|Machine:Worker) before the model matrices
are created.
If you want to start with the formula and see what that means for the
model then use these rules:
- a term including the '|' operator is a random effects term
- if the left-hand side of the '|' operator is 1 then the random
effects are scalar random effects, one for each level of the factor on
the right of the '|'
- random effects associated with different terms are independent
- random effects associated with different levels of the factor within
a term are independent
- the variance of the random effects within the same term is constant
However, there is another mixed-effects model that could make sense
for these data. Suppose I consider the variations associated with
each worker as a vector of length 3 (Machines A, B and C) with a
symmetric, positive semidefinite 3 by 3 variance-covariance matrix. I
fit that model as
> print(fm2 <- lmer(score ~ Machine + (Machine|Worker), Machines), corr = FALSE)
Linear mixed model fit by REML
Formula: score ~ Machine + (Machine | Worker)
Data: Machines
AIC BIC logLik deviance REMLdev
228.3 248.2 -104.2 216.6 208.3
Random effects:
Groups Name Variance Std.Dev. Corr
Worker (Intercept) 16.64051 4.07928
MachineB 34.54670 5.87764 0.484
MachineC 13.61398 3.68971 -0.365 0.297
Residual 0.92463 0.96158
Number of obs: 54, groups: Worker, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 52.356 1.681 31.151
MachineB 7.967 2.421 3.291
MachineC 13.917 1.540 9.037
It may be more meaningful to write it as
> print(fm3 <- lmer(score ~ Machine + (0+Machine|Worker), Machines), corr = FALSE)
Linear mixed model fit by REML
Formula: score ~ Machine + (0 + Machine | Worker)
Data: Machines
AIC BIC logLik deviance REMLdev
228.3 248.2 -104.2 216.6 208.3
Random effects:
Groups Name Variance Std.Dev. Corr
Worker MachineA 16.64097 4.07933
MachineB 74.39557 8.62529 0.803
MachineC 19.26646 4.38936 0.623 0.771
Residual 0.92463 0.96158
Number of obs: 54, groups: Worker, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 52.356 1.681 31.150
MachineB 7.967 2.421 3.291
MachineC 13.917 1.540 9.037
Now we are fitting 3 variances and 3 covariances for the random
effects instead of the previous models which had two variances. The
difference in the models is exactly what made you pause - the simpler
model assumes that, conditional on the random effect for the worker,
the worker/machine random effects are independent and have constant
variance. In the more general models the worker/machine interactions
are allowed to be correlated within worker.
It is more common to allow this kind of correlation within subject in
models for longitudinal data (the Laird-Ware formulation) where each
subject has a random effect for the intercept and a random effect for
the slope with respect to time and these can be correlated. However,
this type of representation can make sense with a factor on the left
hand side of the '|' operator, like the Machine factor here. If that
factor has a large number of levels then the model quickly becomes
unwieldy because the number of variance-covariance parameters to
estimate is quadratic in the number of levels of the factor on the
lhs.
I hope this helps.
> Having got there: Presuming that I'm more-or-less on the right track
> in my foregoing conjecture that it's the over-simple dependence
> structure
> that is the problem with what's delivered by the
> Package-Which-Must-Not-Be-Named,
> how might one go about being less simple-minded? I.e. what might be
> some
> more realistic dependence structures, and how would one specify these
> in lmer?
> And how would one assess whether the assumed dependence structure
> gives
> a reasonable fit to the data?
>
>> If you can describe how many variance components you think should be
>> estimated in your model and what they would represent then I think it
>> will be easier to describe how to fit the model.
>
> How does this fit in with my conjecture (above) about what I've been
> missing all these years? Does it fit? How many variance components
> are there in the ``naive'' model? It looks like 5 to me ... but
> maybe
> I'm totally out to lunch in what I think I'm understanding at this
> stage.
> (And besides --- there are three sorts of statistician; those who can
> count,
> and those who can't.)
>
> Thank you for your indulgence.
>
> cheers,
>
> Rolf Turner
>
> ######################################################################
> Attention:This e-mail message is privileged and confidential. If you are not
> theintended recipient please delete the message and notify the sender.Any
> views or opinions presented are solely those of the author.
>
> This e-mail has been scanned and cleared by
> MailMarshalwww.marshalsoftware.com
> ######################################################################
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Machines.pdf
Type: application/pdf
Size: 38675 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20080513/ae7a3c02/attachment.pdf>
More information about the R-help
mailing list