I have a crossed-classified (Hox, 2010) mixed effects model—also known as crossed random factors (West, Welch, & Galecki, 2015), but I am struggling with how to write the formula for an interaction.

I have observations ($i$) nested within items ($j$) and respondents ($k$). items and respondents are not nested; they are both the units of analysis at "level two." I measured a characteristic of the item, called $Z_j$, and a characteristic of the person, called $X_k$. Thus, both of these are "level two" predictors, but are crossed—not measured at the same unit of analysis.

The primary hypothesis is: Does the relationship between $X_k$ and the outcome, $y_{i(jk)}$, depend on $Z_j$?

The R code seems harmless:

y ~ x * z + (1 + x | j) + (1 + z | k)

But I hit a snag when writing down the formula.

The "level one" equation is simple enough, where we let the intercept vary by items $j$ and respondents $k$, which are crossed (hence in parentheses together, showing they are both at the same "level"):

$y_{i(jk)} = \beta_{0(jk)} + \epsilon_{i(jk)}$

And now the main effects come in when they predict the intercept:

$\beta_{0(jk)} = \gamma_{00} + \gamma_{01j}X_k + \gamma_{02k}Z_j + u_{0j} + v_{0k}$

This shows that the relationship between $X$ and $y$ varies by $j$; and it shows that the relationship between $Z$ and $y$ varies by $k$.

Now, my hypothesis is specified as the $y \sim X$ relationship depending on $Z$. I could write this as:

$\gamma_{01j} = \gamma_{010} + \gamma_{011}Z_j + u_{01j}$

When we substitute everything in, we get two main effects, an interaction, and random effects around intercepts at both $j$ and $k$ as well as random slopes for the main effects.

However, if I add this equation:

$\gamma_{02k} = \gamma_{020} + \gamma_{021}X_k + v_{01k}$

Then we get two identical interaction effects, both $\gamma_{011}$ and $\gamma_{021}$. But since interactions are multiplicative terms, then the way that I formulate my hypothesis shouldn't matter. This is where my confusion arises. Can I include both? It appears not, since they model the same thing. If not, are the two equivalent? Thus, my questions are:

  • Am I specifying this model correctly? (That is, before adding $\gamma_{02k}$)

  • If so, then does that mean: $\gamma_{01j} = \gamma_{010} + \gamma_{011}Z_j + u_{01j}$ and $\gamma_{02k} = \gamma_{020} + \gamma_{021}X_k + v_{01k}$ are equivalent?


Hox, J. (2010). Multilevel Analysis, 2nd Edition.

West, B., Welch, K., & Galecki, A. (2015). Linear Mixed Models, 2nd Edition.

  • 2
    $\begingroup$ Your claim is: $$ (Y\mid \mathcal B =\mathbf b, \mathcal G = \mathbf g) \sim \mathcal{N}(\mathbf{C\beta} + \mathbf{Vb} + \mathbf{Tg}, \sigma) \\ \text{where } \mathbf C = [1, X, Z, X\times Z],\ \mathbf V = [1, X], \\ \quad\mathbf T = [1, Z],\ \mathcal B \sim \mathcal N_2(\mathbf{0, \Sigma_0}),\ \mathcal G \sim \mathcal N_2(\mathbf{0, \Sigma_1}) $$ $\endgroup$ Jan 29, 2019 at 15:19
  • 1
    $\begingroup$ @HeteroskedasticJim could you relate this to the vector notation I use above? $\endgroup$
    – Mark White
    Jan 29, 2019 at 16:32
  • $\begingroup$ I don't think it'd be any different from the answer you already have below. That's why I left it as a comment. $\endgroup$ Jan 30, 2019 at 0:15

1 Answer 1


You could write the model as follows:

$$\left \{ \begin{array}{l} y_{ijk} = \beta_0 + \beta_1 X_{ijk} + \beta_2 Z_{ijk} + \beta_3 \{X_{ijk} \times Z_{ijk} \}+ \\ \quad\quad\quad\quad b_{0k} + b_{1k} X_{ijk} + u_{0j} + u_{1j} Z_{ijk} + \varepsilon_{ijk},\\\\ b_k = (b_{0k}, b_{1k}) \sim \mathcal N(0, D), \\ u_j = (u_{0j}, u_{1j}) \sim \mathcal N(0, V), \\ \varepsilon_{ijk} \sim \mathcal N(0, \sigma^2), \end{array} \right.$$

where $D$ and $V$ denote the variance-covariance matrices for the two sets of random effects, $b_k$ and $u_j$, respectively.

  • $\begingroup$ What about the interaction? $\endgroup$
    – Mark White
    Jan 29, 2019 at 13:50
  • $\begingroup$ Yes, I'll update it. $\endgroup$ Jan 29, 2019 at 14:05
  • 1
    $\begingroup$ Shouldn't the lowercase $b_1$ and $u_1$ terms swap subscripts? How can the relationship between $X_k$ and $y$ vary by $k$? it seems like the relationship would vary by $j$. $\endgroup$
    – Mark White
    Jan 29, 2019 at 14:41
  • 1
    $\begingroup$ @MarkWhite actually, you can permit the sex effect to be different across persons. See page 16 of this JStatSoft paper using lme4 for a variety of IRT jstatsoft.org/v39/i12 $\endgroup$ Jan 30, 2019 at 0:41
  • 1
    $\begingroup$ @MarkWhite I assumed you were familiar with IRT since you were modeling item-response data using similar techniques to the one in the paper. I just think of it simply as allowing the effect of sex to be different by person. If you use indicator coding with male (1) female (0), not all females are zero and not all males are whatever the male coefficient is. $\endgroup$ Jan 30, 2019 at 1:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.