I found an answer to my question on this thread: Repeated measures ANOVA with lme in R for two within-subject factors (somehow this thread was already one of my favorites, I must have forgotten about it). The specification is a little unhandy.
m6 <- lme(mean ~ condition*group*problem*topic,
random = list(code=pdBlocked(list(~1, pdIdent(~problem-1), pdIdent(~topic-1)))), data = d)
anova(m6)
However, the denominator dfs are still wrong, as noted in the thread and apparent in comparisons between the ANOVA and lme dfs.
data.frame(effect = rownames(anova(m6)), denDf= anova(m6)$denDF)
m4$ANOVA[,c("Effect", "DFd")]
As long as there are no other ideas, I think I will need to do the analysis in lme4, for which I wil need to post another question.
I believe that lmer won't be able to duplicate what comes out of aov because it does not have the capability of restricting the variance-covariance matrix of the random components to compound symmetry as done in aov. However, you can still try something like
require(lme4)
# assuming a simple symmetric positive-definite structure of variance-covariance matrix
anova(m2 <- lmer(mean ~ condition*group*problem*topic + (0+problem | code) + (0+topic | code), data = d))
or a simple model
anova(m3 <- lmer(mean ~ condition*group*problem*topic + (1|code), data = d))
Then you can compare the two models:
anova(m2, m3)
Models:
m3: mean ~ condition * group * problem * topic + (1 | code)
m2: mean ~ condition * group * problem * topic + (0 + problem | code) +
m2: (0 + topic | code)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
m3 98 24985 25572 -12395
m2 117 24899 25599 -12332 124.4 19 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The above result indicates the complicated model (the 1st one) is much better. In terms of model complexity, m1 with aov in your OP fits between m2 and m3.
To obtain p-values and confidence intervals for specific effects, do this
require(languageR)
set.seed(101)
# It will take a long time to run the MCMC simulations due to the huge number of effects in the model
mcmc2 <- pvals.fnc(m2, nsim=10000, withMCMC=TRUE)
mcmc2$fixed
The last line will show you the confidence levels and two p-values (one from MCMC simulations, and one from fitted t-statistic) for all the possible fixed effects specified in the model. I am not copying the result here because it's a long table.
If you want to know if the those groups have different variance-covariance structure, you may try
anova(m20 <- lmer(mean ~ condition*group*problem*topic + (0+problem|code/group) + (0+condition|code/group), data = d))
or,
anova(m21 <- lmer(mean ~ condition*group*problem*topic + (0+problem|code/condition) + (0+topic|code/condition), data = d))
And then compare the above two models with the simple one:
anova(m2, m20, m21)
Best Answer
When you are specifying random effects in an
lme4::lmermodel, the random factors go on the left of the pipe and the non-independence grouping variables go on the right, so the fully specified model in your question would very likely be:I took some time to explore the difference between a random effect on the left of the pipe to a random effect on the right side of the pipe but it made a better post on it's own than as an answer to your particular question.
RPubs doc
Gist code
The most noticeable difference between the following two models...
...is that the latter estimates a random correlation parameter between random intercepts and random slopes. If you remove that random correlation parameter...
...the two models return the exact same fixed effects (parameter estimates and associated errors), though I would guess that similarity depends on the particular design of the study.