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, 3Fixed 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: [EMAIL PROTECTED]
tel.: +261 32 40 165 53 (cell)
+261 20 22 665 36 ext. 225 (work)
+261 20 22 494 37 (home)______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
