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, 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-16For 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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Nicholas Galwey
Sent: Wednesday, January 26, 2005 1:45 PM
To: [email protected]
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)
______________________________________________ [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
