**4. Multilevel model methodology**

To propose a multilevel model, it is necessary to determine which variables will be in the fixed part, which in the random part, and the cross-level interactions. This task can be complex, so we need a strategy to build the model.

### **4.1 Multilevel model building strategy**

In this section, we show an adaptation of the bottom-up strategy to specify a threelevel multilevel model. The bottom-up methodology is used in the frequentist approach [3]. Our adaptation proposal tries to use the Bayesian LOO-CV from Step 2 to Step 7 for model comparison.

*Step 1.* Fitting the intercept-only model:

$$\mathcal{Y}\_{ijk} = \mathfrak{xi}\_{000} + \left(\nu\_{00k} + \mu\_{0jk} + \mathfrak{e}\_{ijk}\right) \tag{37}$$

This model gives a basal line to compare with the next models. *Step 2.* Add all the level-one independent variables fixed:

$$\mathcal{Y}\_{ijk} = \xi\_{000} + \sum\_{p=1}^{P} a\_p \mathbf{x}\_{p\bar{i}k} + \sum\_{q=P+1}^{P+Q} \xi\_{q00} \mathbf{x}\_{q\bar{i}k} + \left(\nu\_{00k} + \mu\_{0\bar{j}k} + \mathfrak{e}\_{\bar{j}k}\right) \tag{38}$$

It must be determined which level-one variable has a significant effect on *y*. We will assume that all the *P* level-one variables are statistically significant. Models 38 and 37 must be compared.

*Step 3.* Add the level-two independent variables fixed:

$$y\_{qjk} = \xi\_{000} + \sum\_{p=1}^{P} a\_p \mathbf{x}\_{pjk} + \sum\_{q=P+1}^{P+Q} \xi\_{q00} \mathbf{x}\_{qjk} + \sum\_{m=1}^{M} \xi\_{0m0} w\_{mjk} + \left(\nu\_{00k} + u\_{0jk} + e\_{jk}\right) \tag{39}$$

It must be determined which level-two variable has a significant effect on *y*. If the variables *wm* explain the variability of *y*, Model 39 should be superior to Model 38. Again, we assume that all the *M* level-two independent variables are statistically significant.

*Step 4.* Add the level-three independent variables fixed:

$$\begin{split} y\_{ijk} &= \xi\_{000} + \sum\_{p=1}^{P} a\_p \mathbf{x}\_{p\bar{i}k} + \sum\_{q=P+1}^{P+Q} \xi\_{q00} \mathbf{x}\_{q\bar{i}k} + \sum\_{m=1}^{M} \xi\_{0m0} w\_{m\bar{j}k} + \sum\_{l=1}^{L} \xi\_{00l} \mathbf{z}\_{l\bar{k}} \\ &+ \left( \nu\_{00k} + \mu\_{0\bar{j}k} + \mathfrak{e}\_{\bar{j}\bar{k}} \right) \end{split} \tag{40}$$

It must be determined which level-three variable has a significant effect on *y*. If the variables *zl* explain the variability of *y*, Model 40 should be superior to Model 39.

Steps 1–3 consider the specification of the fixed part of the three-level multilevel model. Now we will specify the random part of the model.

*Step 5.* Assessing whether any of the slopes of the independent variables at level one has a significative variance component between groups at level two or level three.

$$\begin{split} y\_{ijk} &= \xi\_{000} + \sum\_{p=1}^{P} a\_p \mathbf{x}\_{pjk} + \sum\_{q=P+1}^{P+Q} \xi\_{q00} \mathbf{x}\_{qjk} + \sum\_{m=1}^{M} \xi\_{0m0} w\_{mjk} + \sum\_{l=1}^{L} \xi\_{00l} \mathbf{z}\_{lk} \\ &+ \left( \nu\_{00k} + u\_{0jk} + \sum\_{q=P+1}^{P+Q} \nu\_{q0k} \mathbf{x}\_{qjk} + \sum\_{q=P+1}^{P+Q} u\_{qjk} \mathbf{x}\_{qjik} + \mathbf{e}\_{ijk} \right) \end{split} \tag{41}$$

where *uqjk* are the level-two residuals of the slopes of the level-one independent variable *xq*, and *ν*'s are the level-three residuals of the slopes of the level-two independent variable *wm*.

Level-one independent variables that do not have a significant slope may have a significant random slope. This step and the next should be carefully performed, because the model can easily become overparameterized and/or have problems such as non-convergence or extremely slow calculations. It is advisable to assess significance of the slopes variable by variable. Next, the model is formulated with all the variables with significant random slopes. If Model 41 is not better than Model 40, the procedure for specifying a three-level multilevel model stops.

*Step 6.* Assessing whether any of the slopes of the level-two independent variable has a significant variance component among level-three groups.

$$\begin{split} y\_{\bar{q}\bar{k}} &= \xi\_{000} + \sum\_{p=1}^{P} a\_{p} \mathbf{x}\_{p\bar{j}k} + \sum\_{q=P+1}^{P+Q} \xi\_{q00} \mathbf{x}\_{q\bar{j}k} + \sum\_{m=1}^{M} \xi\_{0m0} w\_{m\bar{j}k} + \sum\_{l=1}^{L} \xi\_{0l\bar{l}} \mathbf{z}\_{l\bar{k}} \\ &+ \left( \nu\_{00k} + u\_{0\bar{j}k} + \sum\_{q=P+1}^{P+Q} \nu\_{q0k} \mathbf{x}\_{q\bar{j}k} + \sum\_{q=P+1}^{P+Q} u\_{q\bar{j}k} \mathbf{x}\_{q\bar{j}k} + \sum\_{m=1}^{M} \nu\_{0mk} w\_{m\bar{j}k} + e\_{\bar{q}k} \right) \end{split} \tag{42}$$

where the *ν*'s are the level-three residuals of the slopes of the level-two independent variable *wm*.

The assessment of random slopes of the level-two variables should be performed variable by variable, and then, all these variables should be included into a model to assess the improvement of the model with respect to Model 41.

*Step 7.* Adding interactions between level-three independent variables and the level-one and level-two independent variables that have a significant slope variance in Steps 5 and 6. This produces the full model:

*yijk* <sup>¼</sup> *<sup>ξ</sup>*<sup>000</sup> <sup>þ</sup><sup>X</sup> *P p*¼1 *αpxpijk* þ *P* X þ*Q q*¼*P*þ1 *<sup>ξ</sup>q*00*xqijk* <sup>þ</sup><sup>X</sup> *M m*¼1 *<sup>ξ</sup>*0*m*0*wmjk* <sup>þ</sup><sup>X</sup> *L l*¼1 *ξ*00*lzlk* þ X *M m*¼1 X *L l*¼1 *ξ*0*mlzlkwmjk* þ *P* X þ*Q q*¼*P*þ1 X *L l*¼1 *ξ<sup>q</sup>*0*lzlkxqijk* þ *P* X þ*Q q*¼*P*þ1 X *M m*¼1 *ξqm*0*wmjkxqijk* þ *P* X þ*Q q*¼*P*þ1 X *M m*¼1 X *L l*¼1 *ξqmlzlkwmjkxqijk* þ *ν*00*<sup>k</sup>* þ *u*0*jk* þ *P* X þ*Q q*¼*P*þ1 *<sup>ν</sup>q*0*kxqijk* þ *P* X þ*Q q*¼*P*þ1 *uqjkxqijk* <sup>þ</sup><sup>X</sup> *M m*¼1 *ν*0*mkwmjk* þ *P* X þ*Q q*¼*P*þ1 X *M m*¼1 *<sup>ν</sup>qmkwmjkxqijk* <sup>þ</sup> *eijk*! *:* (43)

When explaining the variances of the random slopes in terms of contextual variables, the model automatically includes interaction terms between levels that compose the fixed part of the model. It is recommended to add variables that explain the variance of the random slope coefficients one by one and not as shown in this step (this was done here to avoid specifying more equations).

When it comes to an MGLM, the methodology changes slightly; that is, instead of defining models in terms of *y*, models are defined in terms of *<sup>g</sup> <sup>μ</sup>ijk* � � <sup>¼</sup> *gEYijk*j*ν*, *<sup>u</sup>* � � � � , and the residual errors at the individual level are no longer specified. An example of this methodology for an MGLM is illustrated below.
