Getting mixed answers on when to include random slopes into crossed effects linear mixed model

Question

I've tried reading through the other threads here on this subject but I still don't really understand the specifics to the answer I'm looking for. I have a crossed effects model that looks something like this, where items are just the words they have to read:

lmer(Word_Reading ~ XS_1 * XS_2 * XI + (1|Subject) + (1|Item), data=data)

The XS predictors here are those that vary for each subject. Like IQ, which is an internal characteristic of the subject, if their XS_1 or XS_2 is higher, their chances of reading accurately should also be higher. Similarly, XI is a predictor that is specific to the items. Lets say for example XI is vowel count. Vowel count should certainly vary for items. For inclusion of slopes, my understanding was that it should go like this, following the logic above.

lmer(Word_Reading ~ XS_1 * XS_2 * XI + (1+XS_1+XS_2|Subject) + (1+XI|Item), data=data)

My brain doesn't handle formulas well, but I believe the structure is similar to this:

$$ \begin{aligned} \text{Word Reading}_{i} &\sim N \left(\mu, \sigma^2 \right) \\ \mu &=\alpha_{ij,ik} + \beta_{1ik}(\text{XS_1}) + \beta_{2}(\text{XS_2}) + \beta_{3ij}(\text{XI}) + \beta_{4}(\text{XS_1} \times \text{XS_2}) + \beta_{5}(\text{XI} \times \text{XS_1}) + \beta_{6}(\text{XI} \times \text{XS_2}) + \beta_{7}(\text{XI} \times \text{XS_1} \times \text{XS_2}) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{3j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\mu_{\alpha_{j}} \\ &\mu_{\beta_{3j}} \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{j}} & \rho_{\alpha_{j}\beta_{3j}} \\ \rho_{\beta_{3j}\alpha_{j}} & \sigma^2_{\beta_{3j}} \end{array} \right) \right) \text{, for Item j = 1,} \dots \text{,J} \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{k} \\ &\beta_{1k} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\mu_{\alpha_{k}} \\ &\mu_{\beta_{1k}} \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{k}} & \rho_{\alpha_{k}\beta_{1k}} \\ \rho_{\beta_{1k}\alpha_{k}} & \sigma^2_{\beta_{1k}} \end{array} \right) \right) \text{, for Subject k = 1,} \dots \text{,K} \end{aligned} $$

The logic is that XS_1 and XS_2 vary within subjects, so you are accounting for this variance among the clusters made by subjects. Likewise, XI varies by item, so the slope included in this model should specify that among items it influences the random effect. For example, if we reduced to just XS_1 and Subject models, I see that looking like this below from Harrison et al. 2018:

However, I was referred to an article that specifies that it should instead be the opposite (Brauer & Curtin, 2017), which gives this example:

Consider a study in which subjects evaluate 20 well-known politicians (10 Democrats and 10 Republicans, “party affiliation,” a dichotomous predictor). The researchers also measure subjects’ level of openness to experience (“openness,” a continuous predictor). The fixed effects structure is relatively easy to determine: it contains the overall intercept, the effect for party affiliation, the effect for openness, and the interaction between the two predictors (four effects in total). To find the maximal random-effects structure, one can use the rules in Table 11. There is nonindependence due to subjects and items, and we thus want to include a by-subject random intercept and a by-item random intercept (first rule). Party affiliation varies within-subjects but between-items. We thus need to specify a by-subject random slope for party affiliation (second rule). Openness to experience varies between-subjects but within items. The design thus calls for a by-item random slope for openness (second rule). It is not the case that both predictors vary within-subjects, and it is also not the case that they both vary within-items. We thus do not have to include any random effects for the interaction term (third rule). The maximal random effects structure contains four random effects: a by-subject random intercept, a by-subject random slope for party affiliation, a by-item random intercept, and a by-item random slope for openness to experience.

It is this part specifically that doesn't make sense to me. How does openness to experience vary between subjects but also within items for this example? The logic for their model is shown directly with this syntax:

model_11 <- lmer(like ~ 1 + affilC * openC + (1 + affilC|subject.ID) + (1 + openC|item.ID), data = d)

Are my model random slopes correct, or am I completely missing what is going on here?

Answer 1

You've answered your own question. I'd like to add a few comments.

As you point out, in this case there are only group-level predictors: Each subject is a group characterized by XS_1 and XS_2; each item is a group characterized by XI. There are no observation-level predictors. If you take two observations for the same subject, they share the XS_i values. If you take two observations for the same item, they share the XI value.

(Brauer & Curtin, 2017) says this in a more sophisticated (but harder to understand?) way:

Party affiliation varies within-subjects but between-items.

This means that two evaluations (observations) by the same study participant of two different politicians have the same "openness" but (possibly) different party affiliation.

Openness to experience varies between-subjects but within items.

This means that two evaluations (observations) of the same politician by different participants have the same party affiliation but different openness.

I also suggest a different way to write the model. It's less mathematically rigorous but easier to read.

I introduce notation first. Let $j[i]$ and $k[i]$ denote the subject and the item of observation $Y_i$, respectively. Since there are no observation-level predictors, the model for $Y_i$ is determined by $j[i]$ and $k[i]$.

For simplicity, I assume there is one predictor for subjects, $x_1$, and one predictor for items, $x_2$, with an interaction. It's straightforward to add more group predictors.

$$ \begin{aligned} Y_i &= \left[\beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2\right] \\ &+ \left[\color{blue}{u_{0,j[i]}} + \color{blue}{u_{1,j[i]}}x_1\right] \\ &+ \left[\color{blue}{v_{0,k[i]}} + \color{blue}{v_{1,k[i]}}x_2\right] \\ &+ \epsilon \end{aligned} $$ where (in blue) the $u_0$s and $v_0$s are random intercepts for subjects and items, respectively; and the $u_1$s and $v_1$s are random slopes. The notation is a bit unwieldy but indexing is important in a multilevel model. To recap, $u_{0,k[i]}$ is the random intercept of the $k$th subject who produced the $i$th observation.

This formula can be re-written to add the fixed main effects to the random components. The result is an equivalent re-parameterized specification of the same model.

$$ \begin{aligned} Y_i &= \left[\beta_0 + \beta_3x_1x_2\right] \\ &+ \left[\color{blue}{u_{0,j[i]}} + \color{blue}{\left(\beta_1 + u_{1,j[i]}\right)}x_1\right] \\ &+ \left[\color{blue}{v_{0,k[i]}} + \color{blue}{\left(\beta_2 + v_{1,k[i]}\right)}x_2\right] \\ &+ \epsilon \\ &= \left[\beta_0 + \beta_3x_1x_2\right] \\ &+ \left[\color{blue}{u_{0,j[i]}} + \color{blue}{u^*_{1,j[i]}}x_1\right] \\ &+ \left[\color{blue}{v_{0,k[i]}} + \color{blue}{v^*_{1,k[i]}}x_2\right] \\ &+ \epsilon \end{aligned} $$

The second formulation highlights that the model allows the main effects of $x_1$ and $x_2$ to vary by subject and by item but the subject-item interaction is fixed (it has no random component). This is perhaps not a very natural assumption, though a random interaction effect can only be estimated if there are repeated observations for each subject-item pair; in that case, add (1 | Subject:Item).

You may want to consider Bayesian hierarchical modeling as well. The following introduction describes how to use the brms package to analyze phonetic data.

L. Nalborczyk, C. Batailler, H. Lœvenbruck, A. Vilain, P.-C. Bürkner. An introduction to Bayesian multilevel models using brms: A case study of gender effects on vowel variability in Standard Indonesian. Journal of Speech, Language, and Hearing Research, 62(5):1225–1242, 2019.

Answer 2

I think I understand the distinction now. The issue is how many data points you would have for a random effects model with crossed effects. For the subject level data point, there would only be one per person, effectively meaning one X score but multiple Y scores. For item it would be the same. I tried simulating this idea with a very small dataset to illustrate for myself what this meant, using just the subject level factor of IQ and the outcome of reading the words.

set.seed(123)
subject <- 1:5
iq <- round(rnorm(n=5,
            mean=120,
            sd=15))
word.read <- rbinom(n=15,
                    size=1,
                    prob = .9)
df <- data.frame(subject,
                 iq)
df2 <- expand_grid(df,
                   word.read)
df2$item <- rep(1:25,3)

Grouping the points by subject yields no logistic regression lines, as each subject has only one measure of IQ.

df2 %>% 
  ggplot(aes(x=iq,
             y=word.read,
             color=as.factor(subject)))+
  geom_point()+
  stat_smooth(method="glm",
              se=F,
              method.args = list(family=binomial))

Grouping them by item, however, gives per item logistic regression fits, as there are multiple points for each X and Y per item.

df2 %>% 
  ggplot(aes(x=iq,
             y=word.read,
             color=as.factor(item)))+
  geom_point()+
  stat_smooth(method="glm",
              se=F,
              method.args = list(family=binomial))