[R] Specification of factorial random-effects model

[Douglas Bates]

Berton Gunter wrote:
> If you read the Help file for lme (!), you'll see that ~1|a*b is certainly
> incorrect.
> 
> Briefly, the issue has been discussed before on this list: the current
> version of lme() follows the original Laird/Ware formulation for **nested**
> random effects. Specifying **crossed** random effects is possible but
> difficult, and the fitting algorithm is not optimized for this. See p. 163
> in Bates and Pinheiro for an example.

The development package lme4 has a version of a linear mixed model 
function that does handle crossed random effects.  In lme4_0.8-1 and 
later the new version of lme, called lmer (which could mean either "lme 
revised" or "lme for R"), has a different syntax for specifying mixed 
models.  A random effects specification is indicated by a `|' character 
which  separates a linear model expression on the left side from the 
grouping factor on the right side.  Because the | operator has very low 
precedence, such terms usually must be enclosed in parentheses.

The same type of specification is used for nested or crossed or 
partially crossed grouping factors.  The only restriction is that the 
grouping factor must have a unique level for each group, which is to say 
that you must explicitly create nested factors - you cannot specify them 
implicitly.

This example could be fit as

 > (fm1 <- lmer(y ~ c + (1|a) +(1|b) + (1|a:b)))
Linear mixed-effects model fit by REML
Formula: y ~ c + (1 | a) + (1 | b) + (1 | a:b)
       AIC      BIC    logLik MLdeviance REMLdeviance
  376.0148 392.7759 -181.0074   369.2869     362.0148
Random effects:
  Groups   Name        Variance Std.Dev.
  a:b      (Intercept)   1.1118  1.0544
  b        (Intercept) 286.8433 16.9364
  a        (Intercept)  86.2138  9.2851
  Residual               3.4626  1.8608
# of obs: 81, groups: a:b, 9; b, 3; a, 3

Fixed effects:
              Estimate Std. Error DF  t value  Pr(>|t|)
(Intercept)  65.91259   11.16262 78   5.9048 8.707e-08
c2           -9.47000    0.50645 78 -18.6989 < 2.2e-16
c3          -10.88259    0.50645 78 -21.4881 < 2.2e-16

For the random effects the Variance column is the estimate of the 
variance component.  The Std.Dev. column is simply the square root of 
the estimated variance.  I find it easier to think in terms of standard 
deviations rather than variances because I can compare the standard 
deviations to the scale of the data.  Note that this column is *not* a 
standard error of the estimated variance component (and purposely so 
because I feel that such quantities are often nonsensical).

A test of, say, whether the variance component for the interaction could 
be zero is performed by fitting the reduced model and using the anova 
function to compare the fitted models.  The p-value quoted for this test 
is conservative because the null hypothesis is on the boundary of the 
parameter space.


 > (fm2 <- lmer(y ~ c + (1|a) +(1|b)))
Linear mixed-effects model fit by REML
Formula: y ~ c + (1 | a) + (1 | b)
       AIC      BIC    logLik MLdeviance REMLdeviance
  379.3209 393.6876 -183.6605   374.8822     367.3209
Random effects:
  Groups   Name        Variance Std.Dev.
  a        (Intercept)  86.3823  9.2942
  b        (Intercept) 286.5391 16.9275
  Residual               4.0039  2.0010
# of obs: 81, groups: a, 3; b, 3

Fixed effects:
             Estimate Std. Error DF  t value  Pr(>|t|)
(Intercept)  65.9126    11.1560 78   5.9083  8.58e-08
c2           -9.4700     0.5446 78 -17.3890 < 2.2e-16
c3          -10.8826     0.5446 78 -19.9829 < 2.2e-16
Warning message:
optim returned message ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
  in: "LMEoptimize<-"(`*tmp*`, value = list(maxIter = 50, msMaxIter = 50,
 > anova(fm1, fm2)
Data:

Models:
fm2: y ~ c + (1 | a) + (1 | b)
fm1: y ~ c + (1 | a) + (1 | b) + (1 | a:b)

     Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
fm2  6  386.88  401.25 -187.44
fm1  7  383.29  400.05 -184.64 5.5953      1    0.01801

>>-----Original Message-----
>>From: r-help-bounces at stat.math.ethz.ch 
>>[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Nicholas Galwey
>>Sent: Wednesday, January 26, 2005 1:45 PM
>>To: r-help at stat.math.ethz.ch
>>Subject: [R] Specification of factorial random-effects model
>>
>>I want to specify two factors and their interaction as random 
>>effects using
>>the function lme().   This works okay when I specify these 
>>terms using the
>>function Error() within the function aov(), but I can't get 
>>the same model
>>fitted using lme().   The code below illustrates the problem.
>>
>> 
>>
>>a <- factor(rep(c(1:3), each = 27))
>>
>>b <- factor(rep(rep(c(1:3), each = 9), times = 3))
>>
>>c <- factor(rep(rep(c(1:3), each = 3), times = 9))
>>
>>y <- c(74.59,75.63,76.7,63.48,63.17,65.99,64,66.35,64.5,
>>
>>   46.57,44.16,47.96,35.09,36.14,35.16,36.4,34.72,34.58,
>>
>>   41.82,47.35,45.74,33.33,36.8,33.38,34.13,34.39,34.48,
>>
>>   89.73,85.24,90.86,82.5,79.44,81.65,77.74,77.02,81.62,
>>
>>   59.32,62.29,60.7,55.42,55.5,51.17,50.54,53.54,51.85,
>>
>>   64.5,63.6,65.19,55.07,50.26,53.73,54.57,47.8,48.8,91.56,
>>
>>   94.49,92.17,82.14,83.16,81.31,83.58,78.63,77.08,60.53,
>>
>>   60.79,58.57,51.28,52.9,51.54,49.15,48.97,51.61,59.44,
>>
>>   60.07,60.07,51.94,52.2,50.2,49.45,50.75,49.56)
>>
>>anovamodel <- aov(y ~ c + Error(a*b))
>>
>>summary(anovamodel)
>>
>>lmemodel <- lme(y ~ c, random = ~ 1|a*b)
>>
>>anova(lmemodel)