[R-sig-ME] [R] lme nesting/interaction advice

[Kingsford Jones]

Hi Rolf,

On Tue, May 13, 2008 at 1:59 PM, Rolf Turner <r.turner at auckland.ac.nz> wrote:

< in response to Doug Bates' useful tutorial...>

>  Thanks very much for your long, detailed, patient, and lucid
>  response to my cri de coeur.  That helps a *great* deal.
>

Hear Hear!

<snip>
>  One point that I'd like to spell out very explicitly, to make sure
>  that I'm starting from the right square:
>
>  The model that you start with in the Machines/Workers example is
>
>
>
> >
> > > fm1 <- lmer(score ~ Machine + (1|Worker) + (1|Machine:Worker), Machines)
> > >
> >
>
>
>  My understanding is that this is the ``only'' model that could be fitted
>  by the Package-Which-Must-Not-Be-Named.  I.e. this package *could not* fit
>  the second, more complex model:
>
>
>
> >
> > > fm2 <- lmer(score ~ Machine + (Machine|Worker), Machines)
> > >
> >
>
>  (At least not directly.)  Can you (or someone) confirm that I've got that
>  right?

Compare:

## R
> m4
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.64098 4.07934
          MachineB 74.39558 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

## "The-Package"

proc mixed data = machines;
class worker machine;
model score = machine  / solution;
random machine / subject = worker type = un;
run;

 Covariance Parameter Estimates

Cov Parm     Subject    Estimate

UN(1,1)      Worker      16.6405
UN(2,1)      Worker      28.2447
UN(2,2)      Worker      74.3956
UN(3,1)      Worker      11.1465
UN(3,2)      Worker      29.1841
UN(3,3)      Worker      19.2675
Residual                  0.9246


The two outputs report essentially the same thing.
Note that e.g, UN(2,1) = 28.2477 approx= .803*4.07934*8.62529
(And, as usual, the fixed effects estimates match too once the
contrasts and 'types' of SS for an ANOVA table are set up)

UN is short for 'unstructured' - a term Doug has pointed out is not
particularly fitting because the covariance matrix is symmetric
positive definite.

>
>  It seems to me to be the key to why I've had such trouble in the past
>  in grappling with mixed models in R.  I.e. I've been thinking like
>  the Package-Which-Must-Not-Be-Named --- that the simple,
> everything-independent
>  model was the only possible model.
>

Although this may well not apply to you, another area of confusion
arises not so much from differences between stats packages but by
differences between methods. I'm not an expert in the estimation
methods but, as I understand it, classic texts describe fitting mixed
models in terms of ANOVA in the OLS framework, calculating method of
moments estimators for the variances of the random effects by equating
observed and expected mean squares (I believe using aov and lm with an
'Error' term would fall into this category, and proc anova and proc
glm would also).  Starting in the 90's these methods started falling
out of fashion in favor of ML/REML/GLS methods (likelihood based),
which offer more flexibility in structuring both the error and random
effects covariance matrices, will not produce negative variance
estimates, and have other nice properties that someone more 'mathy'
than me could explain.  Tools like lme, lmer, proc mixed and proc
glimmix fall into this category.

hoping this helps,

Kingsford Jones







>  Thanks again.
>
>         cheers,
>
>                 Rolf
>
>
>
>  ######################################################################
>  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.
>