Edit: Just adding a relevant blog post that discusses checking if a random effect should be included or not, but my question is more specifically based on deciding if the intercepts and slopes of random effects should be included or not (particularly in the case of a crossed design with 2+ random effects).
Note: I'm using the same sample data as my previous (different) question here. I am aware of previous related questions such as this and this but I don't think the question was asked/answered in a general sense.
Consider the following data (R code):
sample <- read.table(text = " Session ID Response Predictor
1: 1 1 60.47 0.012
2: 2 1 42.53 -0.50
3: 3 1 64.45 0.01
4: 1 2 67.64 0.01
5: 2 2 51.77 0.04
6: 3 2 68.84 0.09
7: 1 3 79.80 -0.05
8: 2 3 46.95 0.43
9: 3 3 83.3 -0.05 ", h = T)
sample$Session <- factor(sample$Session)
sample$ID <- factor(sample$ID)
Where Session and ID are crossed random effects and we are examining the relationship between Predictor and Response. Now considering the possible combinations of random slopes and intercepts I believe (based on this thread) that there are a number of ways to model this situation.
model.1 <- lmer(Response ~ Predictor + (1 + Predictor| ID) +
(1 + Predictor| Session), data=sample,
REML = FALSE)
model.2 <- lmer(Response ~ Predictor + (1 | ID) +
(1 + Predictor| Session), data=sample,
REML = FALSE)
model.3 <- lmer(Response ~ Predictor + (1 + Predictor| ID) +
(1 | Session), data=sample, REML = FALSE)
model.4 <- lmer(Response ~ Predictor + (0 + Predictor| ID) +
(0 + Predictor| Session), data=sample,
REML = FALSE)
model.5 <- lmer(Response ~ Predictor + (0 + Predictor | ID) +
(1 + Predictor| Session), data=sample,
REML = FALSE)
model.6 <- lmer(Response ~ Predictor + (1 + Predictor| ID) +
(0 + Predictor | Session), data=sample,
REML = FALSE)
This doesn't even include all of the possibilities where we find one of the random effects to not be significant e.g.
model.7 <- lmer(Response ~ Predictor + (1 + Predictor| ID) ,
data=sample, REML = FALSE)
...
...
etc
So the question is: What is a systematic way of deciding which of these models should be used and thus should different random intercepts and slopes be included?
Some possibilities I've considered are:
- Using the
anova
function and choosing the lowest AIC - Comparing to the null model e.g.
model.null <- lmer(Response ~ 1 + (1 | ID) + (1 | Session), data=sample, REML = FALSE)
- Some version of backwards/forwards stepwise elimination based on likelihood ratio tests
- Rejecting models that have a standard deviation of 0 for any random effect
(1 + Predictor | ID:Session)
. If all these specifications make sense and you have sufficient information in the data to estimate them, you can attempt the more flexible specifications. In your example above, you can't do this as the interaction of ID and Session is a single data point. $\endgroup$