lme nesting/interaction advice
I've been looking back over this discussion.
Another model that one can fit using lme is:
> lme(score~Machine, random=list(Worker=pdIdent(~0+Machine)),
+ weights=varIdent(form=~1|Machine), data=Machines)
Linear mixed-effects model fit by REML
Data: Machines
Log-restricted-likelihood: -108.9547
Fixed: score ~ Machine
(Intercept) MachineB MachineC
52.355556 7.966667 13.916667
Random effects:
Formula: ~0 + Machine | Worker
Structure: Multiple of an Identity
MachineA MachineB MachineC Residual
StdDev: 6.06209 6.06209 6.06209 1.148581
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | Machine
Parameter estimates:
A B C
1.0000000 0.8713263 0.5859709
Number of Observations: 54
Number of Groups: 6
This insists (I think) that conditional on the random effect for the
worker,
the worker/machine random effects be independent,
but allows them to have different variances. I am wondering whether
it is possible to fit such a model using lmer().
[In this example the large estimated correlations suggest that it is not
a sensible model.]
John Maindonald email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473 fax : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.
On 14 May 2008, at 2:49 AM, Douglas Bates wrote:
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
######################################################################
<Machines.pdf>_______________________________________________ R-sig-mixed-models at r-project.org mailing list