# How to analyze this dataset with 2 between-, two within- factors using lme4?

This post is a follow up on my previous post (which was interested in lme) and uses the same dataset. Now I would like to know how to analyze it using lme4.

## The data

The data is from a behavioral experiment in which participants in 6 groups (based on two crossed factors) worked on 16 trials (two crossed 4-level factors). That is, we have a dataset d with two between-subject factors, group and condition, and two within-subject factors (i.e., repeated-measures factors), topic and problem (I uploaded the data to pastebin, so everybody should be able to obtain it), the participant id is code, the dv is mean:

> d <- read.table(url("http://pastebin.com/raw.php?i=4hRFyaRj"), colClasses = c(rep("factor", 6), "numeric"))
> str(d)
'data.frame':   2928 obs. of  6 variables:
$code : Factor w/ 183 levels "A03U","A08C",..: 1 1 1 1 1 1 1 1 1 1 ...$ group    : Factor w/ 2 levels "control","experimental": 2 2 2 2 2 2 2 2 2 2 ...
$condition: Factor w/ 3 levels "alternatives",..: 3 3 3 3 3 3 3 3 3 3 ...$ topic    : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 2 2 2 2 3 3 ...
$problem : Factor w/ 4 levels "AC","DA","MP",..: 3 4 1 2 3 4 1 2 3 4 ...$ mean     : num  94.5 94.5 86.5 84.5 80 46.5 73.5 43.5 51 39 ...


The usual way to analyze this data would to fit an ANOVA on this data (note how the error term is constructed for the within-subject factors):

m1 <- aov(mean ~ (condition*group*problem*topic) + Error(code/(problem*topic)), d)


## The Question

My main interest in the data is the following:

• Is there an effect of the group factor on any level (i.e., main effect or interaction)? I hope there is not.
• Is there an interaction of condition with problem ? Or even an interaction of condition with problem and topic?

I have two questions regarding the analysis in lme4:

1. How can I specify these questions using lme4?
2. As lme4 does not provide p-values, how do I determine whether a variable (e.g., group) has any effect (I imagine using some kind of likelihood ratio test) and what is the critical value above which I need to accept effect to be 'significant'?

As is probably obvious from the above description I am no expert in lme4 neither a statistician, so both Venables & Ripley and the lme4 Book by Bates gave me a hard time. Leaving me kind of clueless as before.

-
Have you read Baayen, Davidson, Bates (2008), "Mixed-effects modeling with crossed random effects for subjects and items"? It can be found at ualberta.ca/~baayen/publications/baayenDavidsonBates.pdf –  Aaron Feb 29 '12 at 17:05
To be honest, I am reading it right now. But so far I have problems understanding what type of slopes I need or not. I hope with the help of the text and here, I can sort this out. –  Henrik Feb 29 '12 at 17:10
You might try to get help in-person too; it's very difficult to fully understand the details of the study over email. Perhaps portal.uni-freiburg.de/imbi/StatConsult/details might be a resource? –  Aaron Feb 29 '12 at 17:22

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)

-
Thanks for your very interesting answer. And I know, that I cannot duplciate the results from aov (there is a reason I want to go away from it). What I want is a lme4 model that mimics the general structure (i.e., the nesting of problem and topic within code). I just do not understand how to probably convert the structure with the two within-subject factors to lme4 syntax. But from your answer I get that it is not that easy to decide on the "correct" random terms. –  Henrik Feb 29 '12 at 22:08
Can you explain what is the difference between models m2, m20, and m21 regarding the within-subject factors? And is there a way to specify something like (problem*topic)/code (cf. the the aov call) that would make sense? –  Henrik Feb 29 '12 at 22:51
@Henrik: What do you mean by converting "the structure with the two within-subject factors to lme4 syntax"? The m3 model in my response exactly handles the possible correlation structure among those levels for each of those two within-subject factors. The only difference here is that, your aov command assumes spherecity which may or may not be reasonable while the m3 model does not put such a strong constraint on the correlation structure. You may argue that m3 might be over-fitting, but, with enough amount of data, this should not be a big issue. –  bluepole Feb 29 '12 at 22:53
Thanks for the clear clarification. It improves my understanding drastically. m3 actually is what I wanted. So what is the difference between m3 and m2? Different random factors or different slopes (i.e., a different slope for topic within code and problem within code with code remaining the only random effect)? Btw, I would upvote more if I could. –  Henrik Feb 29 '12 at 23:09
@Henrik: Sorry I meant m2 instead of m3 in my previous comment. m3 has a simple random intercept while m2 models the pair-wise correlation among those levels of each within-subject factor. With a continuous variable, that would be a model with a random slope. If you feel that my response answers your OP, you can accept it by clicking the check mark. –  bluepole Feb 29 '12 at 23:48