Mixed Model – How to Analyze a Dataset with 2 Between- and 2 Within-Factors Using lme4 in Repeated Measures

lme4-nlmemixed modelrepeated measures

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.

Best Answer

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)