How to choose random- and fixed-effects structure in linear mixed models?

Question

Consider the following data from a two-way within subjects design:

df <- "http://personality-project.org/r/datasets/R.appendix4.data"
df <- read.table(df,header=T)
head(df)

Observation Subject Task Valence Recall
1           1     Jim Free     Neg      8
2           2     Jim Free     Neu      9
3           3     Jim Free     Pos      5
4           4     Jim Cued     Neg      7
5           5     Jim Cued     Neu      9
6           6     Jim Cued     Pos     10

I would like to analyze this using mixed-linear models. Considering all possible fixed- and random-effects there are multiple possible models:

# different fixed effects with random-intercept
a0 <- lmer(Recall~1 + (1|Subject), REML=F,df)
a1 <- lmer(Recall~Task + (1|Subject), REML=F,df)
a2 <- lmer(Recall~Valence + (1|Subject), REML=F,df)
a3 <- lmer(Recall~Task+Valence + (1|Subject), REML=F,df)
a4 <- lmer(Recall~Task*Valence + (1|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope
b0 <- lmer(Recall~1 + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b1 <- lmer(Recall~Task + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b2 <- lmer(Recall~Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b3 <- lmer(Recall~Task+Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b4 <- lmer(Recall~Task*Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope including variance-covariance matrix
c0 <- lmer(Recall~1 + (1 + Valence + Task|Subject), REML=F,df)
c1 <- lmer(Recall~Task + (1 + Valence + Task|Subject), REML=F,df)
c2 <- lmer(Recall~Valence + (1 + Valence + Task|Subject), REML=F,df)
c3 <- lmer(Recall~Task+Valence + (1 + Valence + Task|Subject), REML=F,df)
c4 <- lmer(Recall~Task*Valence + (1 + Valence + Task|Subject), REML=F,df)

1. What is the recommended way to select the best fitting model in this context? When using log-likelihood ratio tests what is the recommended procedure? Generating models upwards (from null model to most complex model) or downwards (from most complex model to null model)? Stepwise inclusion or exclusion? Or is it recommended to put all models in one log-likelihood ratio test and select the model with the lowest p-value? How to compare models that are not nested?

2. Is it recommended to first find the appropriate fixed-effects structure and then the appropriate random-effects structure or the other way round (I have found references for both options...)?

3. What is the recommended way of reporting results? Reporting the p-value from the log-likelihood ratio test comparing the full mixed-model (with the effect in question) to reduced model (without the effect in question). Or is it better to use log-likelihood ratio test to find the best fitting model and then use lmerTest to report p-values from the effects in the best fitting model?

Answer 1

I'm not sure there's really a canonical answer to this, but I'll give it a shot.

What is the recommended way to select the best fitting model in this context? When using log-likelihood ratio tests what is the recommended procedure? Generating models upwards (from null model to most complex model) or downwards (from most complex model to null model)? Stepwise inclusion or exclusion? Or is it recommended to put all models in one log-likelihood ratio test and select the model with the lowest p-value? How to compare models that are not nested?

It depends what your goals are.

• In general you should be very, very careful about model selection (see e.g. this answer, or this post, or just Google "Harrell stepwise" ...).
• If you are interested in having your p-values be meaningful (i.e. you are doing confirmatory hypothesis testing), you should not do model selection. However: it's not so clear to me whether model selection procedures are quite as bad if you are doing model selection on non-focal parts of the model, e.g. doing model selection on the random effects if your primary interest is inference on the fixed effects.
• There's no such thing as "putting all the models in one likelihood ratio test" -- likelihood ratio testing is a pairwise procedure. If you wanted to do model selection (e.g.) on the random effects, I would probably recommend an "all at once" approach using information criteria as in this example -- that at least avoids some of the problems of stepwise approaches (but not of model selection more generally).
• Barr et al. 2013 "Keep it maximal" Journal of Memory and Language (doi:10.1016/j.jml.2012.11.001) would recommend using the maximal model (only).
• Shravan Vasishth disagrees, arguing that such models are going to be underpowered and hence problematic unless the data set is very large (and the signal-to-noise ratio is high)
• Another reasonably defensible approach is to fit a large but reasonable model and then, if the fit is singular, remove terms until it isn't any more
• With some caveats (enumerated in the GLMM FAQ), you can use information criteria to compare non-nested models with differing random effects (although Brian Ripley disagrees: see bottom of p. 6 here)

Is it recommended to first find the appropriate fixed-effects structure and then the appropriate random-effects structure or the other way round (I have found references for both options...)?

I don't think anyone knows. See previous answer about model selection more generally. If you could define your goals sufficiently clearly (which few people do), the question might be answerable. If you have references for both options, it would be useful to edit your question to include them ... (For what it's worth, this example (already quoted above) uses information criteria to select the random effects part, then eschews selection on the fixed-effect part of the model.

What is the recommended way of reporting results? Reporting the p-value from the log-likelihood ratio test comparing the full mixed-model (with the effect in question) to reduced model (without the effect in question). Or is it better to use log-likelihood ratio test to find the best fitting model and then use lmerTest to report p-values from the effects in the best fitting model?

This is (alas) another difficult question. If you report the marginal effects as reported by lmerTest, you have to worry about marginality (e.g., whether the estimates of the main effects of A and B are meaningful when there is an A-by-B interaction in the model); this is a huge can of worms, but is somewhat mitigated if you use contrasts="sum" as recommend by afex::mixed(). Balanced designs help a little bit too. If you really want to paper over all these cracks, I think I would recommend afex::mixed, which gives you output similar to lmerTest, but tries to deal with these issues.

Answer 2

Update May 2017: As it turns out, a lof of what I have written here is kind of wrongish. Some updates are made throughout the post.

---

I agree a lot with what has been said by Ben Bolker already (thanks for the shout-out to afex::mixed()) but let me add a few more general and specific thoughts on this issue.

Focus on fixed versus random effects and how to report results

For the type of experimental research that is represented in the example data set from Jonathan Baron you use the important question is usually whether or not a manipulated factor has an overall effect. For example, do we find an overall main effect or interaction of Task? An important point is that in those data sets usually all factors are under complete experimental control and randomly assigned. Consequently, the focus of interest is usually on the fixed effects.
In contrast, the random effects components can be seen as "nuisance" parameters that capture systematic variance (i.e., inter-individual differences in the size of the effect) that are not necessarily important for the main question. From this point of view the suggestion of using the maximal random effects structure as advocated by Barr et al. follows somewhat naturally. It is easy to imagine that an experimental manipulation does not affect all individuals in the exact same way and we want to control for this. On the other hand, the number of factors or levels is usually not too large so that the danger of overfitting seems comparatively small.

Consequently, I would follow the suggestion of Barr et al. and specify a maximal random effects structure and report tests of the fixed effects as my main results. To test the fixed effects I would also suggest to use afex::mixed() as it reports tests of effects or factors (instead of test of parameters) and calculates those tests in a somewhat sensible way (e.g., uses the same random effects structure for all models in which a single effect is removed, uses sum-to-zero-contrasts, offers different methods to calculate p-values, ...).

What about the example data

The problem with the example data you gave is that for this dataset the maximal random effects structure leads to an oversaturated model as there is only one data point per cell of the design:

> with(df, table(Valence, Subject, Task))
, , Task = Cued

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

, , Task = Free

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

Consequently, lmer chokes on the maximal random effects structure:

> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable

Unfortunately, there is to my knowledge no agreed upon way to deal with this problem. But let me sketch and discuss some:

1. A first solution could be to remove the highest random slope and test the effects for this model:

   require(afex)
   mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 6.56   1 4.00      1.00     .06
   2      Valence 0.80   2 3.00      0.75     .53
   3 Task:Valence 0.42   2 8.00      1.00     .67
   

   However, this solution is a little ad-hoc and not overly motivated.

   Update May 2017: This is the approach I am currently endorsing. See this blog post and the draft of the chapter I am co-authoring, section "Random Effects Structures for Traditional ANOVA Designs".

2. An alternative solution (and one that could be seen as advocated by Barr et al.'s discussion) could be to always remove the random slopes for the smallest effect. This has two problems though: (1) Which random effects structure do we use to find out what the smallest effect is and (2) R is reluctant to remove a lower-order effect such as a main effect if higher order effects such as an interaction of this effect is present (see here). As a consequence one would need to set up this random effects structure by hand and pass the so constructed model matrix to the lmer call.

3. A third solution could be to use an alternative parametrization of the random effects part, namely one that corresponds to the RM-ANOVA model for this data. Unfortunately (?), lmer doesn't allow "negative variances" so this parameterization doesn't exactly correspond to the RM-ANOVA for all data sets, see discussion here and elsewhere (e.g. here and here). The "lmer-ANOVA" for these data would be:

   > mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 7.35   1 4.00      1.00     .05
   2      Valence 1.46   2 8.00      1.00     .29
   3 Task:Valence 0.29   2 8.00      1.00     .76
   

Given all this problems I simply wouldn't use lmer for fitting data sets for which there is only one data point per cell of the design unless a more agreed upon solution for the problem of the maximal random effects structure is available.

4. Instead, I would One also could still use the classical ANOVA. Using one of the wrappers to car::Anova() in afex the results would be:

   > aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
           Effect         df  MSE      F  ges   p
   1      Valence 1.44, 5.75 4.67   1.46  .02 .29
   2         Task       1, 4 4.08 7.35 +  .07 .05
   3 Valence:Task 1.63, 6.52 2.96   0.29 .003 .71
   

   Note that afex now also allows to return the model fitted with aov which can be passed to lsmeans for post-hoc tests (but for test of effects the ones reported by car::Anova are still more reasonable):

   > require(lsmeans)
   > m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
   > lsmeans(m, ~Task+Valence)
    Task Valence lsmean       SE   df lower.CL upper.CL
    Cued Neg       11.8 1.852026 5.52  7.17157 16.42843
    Free Neg       10.2 1.852026 5.52  5.57157 14.82843
    Cued Neu       13.0 1.852026 5.52  8.37157 17.62843
    Free Neu       11.2 1.852026 5.52  6.57157 15.82843
    Cued Pos       13.6 1.852026 5.52  8.97157 18.22843
    Free Pos       11.0 1.852026 5.52  6.37157 15.62843
   
   Confidence level used: 0.95