[R] [R-sig-ME] lme nesting/interaction advice
John Maindonald
john.maindonald at anu.edu.au
Fri May 16 06:28:59 CEST 2008
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
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-help
mailing list