[R] Specification of factorial random-effects model

[Renaud Lancelot]

Berton Gunter a écrit :
> 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.
> 
> -- Bert Gunter
> Genentech Non-Clinical Statistics
> South San Francisco, CA
>  
> "The business of the statistician is to catalyze the scientific learning
> process."  - George E. P. Box
>  
>  
[snip]

You can use the package lme4 to fit models with crossed random effects. 
However, it looks like you have to explicitly create the interaction term:

 > library(lme4)
Loading required package: Matrix
Loading required package: latticeExtra
 >
 > 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)
 > Data <- data.frame(a=a, b=b, ab=paste(a, b, sep = ""), c=c, y=y)
 > rm(a, b, c, y)
 > fm1 <- lme(y ~ c, random = ~ 1 | a * b, data = Data)
 > fm1
Linear mixed-effects model
Fixed: y ~ c
  Data: Data
  log-restricted-likelihood:  -183.6605
Random effects:
  Groups   Name        Variance Std.Dev.
  b        (Intercept) 286.5391 16.9275
  a        (Intercept)  86.3823  9.2942
  Residual               4.0039  2.0010
# of obs: 81, groups: b, 3; a, 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
 >
 > fm2 <- lme(y ~ c, random = ~ 1 | a + b + ab, data = Data)
 > fm2
Linear mixed-effects model
Fixed: y ~ c
  Data: Data
  log-restricted-likelihood:  -181.0074
Random effects:
  Groups   Name        Variance Std.Dev.
  ab       (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: ab, 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

#######

Beware: if you loaded nlme before, you have to start a new session to 
use lme4 which conflicts with nlme.

Best,

Renaud


-- 
Dr Renaud Lancelot, vétérinaire
C/0 Ambassade de France - SCAC
BP 834 Antananarivo 101 - Madagascar

e-mail: renaud.lancelot at cirad.fr
tel.:   +261 32 40 165 53 (cell)
         +261 20 22 665 36 ext. 225 (work)
         +261 20 22 494 37 (home)