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)
>
> 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)