Recovering Crossed Random Effects in Mixed-Effects Models Using Model Averaging

Model Averaging With Akaike Weights

Model averaging with Akaike weights uses the AIC index to compute an average estimate of the target parameters (e.g., Akaike, 1979; Burnham & Anderson, 2002; Claeskens & Hjort, 2008; Kishino et al., 1991; Steele et al., 2014). This procedure involves fitting the competing models, M₁, M₂, …, M_r, and computing the so-called Akaike weights using the AIC index (e.g., Akaike, 1974) for each competing model as relative evidence in favor of each one among all the competing models (Burnham & Anderson, 2002). As explained, the AIC indices of all R models are ranked as follows:

where ∆_r represents the difference between the AIC of each r competing model and that of the best fitting model (AIC_min). Then, Akaike weights are then calculated as:

where ω_r is the resulting Akaike weight for each competing model based on ∆_r of all R competing models (Burnham & Anderson, 2002). Next, a weighted estimation of the target parameters is obtained by weighting the estimations of each model r using its ω_r. These averaged estimations have been found to generate estimations of SEs of fixed effects of MEMs-CR with small bias, being a less risky approach than model selection strategies (Martínez-Huertas et al., 2022). In the present study, model averaging was applied on the estimation of the random effects of MEMs-CR.

Bayesian Model Averaging With BIC Posterior Probabilities

BMA shares the same rationale than model averaging with Akaike weights about the computation of averaged estimates for target parameters, but using a Bayesian framework (e.g., Chatfield, 1995; Draper, 1995; Fragoso et al., 2018; Kaplan & Lee, 2018; Hinne et al., 2020; Steel, 2020). It is worth mentioning that this approach works with posterior probabilities and that Bayes factor is one of the most accurate approximations to them, but in this study we are going to use BIC as an approximation to the posterior probability of the competing models due to its computational simplicity and effective performance (see Schwarz, 1978, and Neath & Cavanaugh, 2012 for a demonstration of the derivation of prior distributions based on BIC). Once different competing models M₁, M₂, …, M_r, have been fitted to the data set D, we can extract their BIC indices (BIC₁, BIC₂, …, BIC_r). Again, it is necessary to rank all the R competing models using:

where ∆_r represents the difference between the BIC of each r competing model and that of the best fitting model (BIC_min). Then, it is possible to approximate the posterior probability of each r model given data D as:

where $Pr\left(M_j\right|D)$ is an approximation to the posterior probability of each competing model, and shares an important equivalence with Akaike weights as they are based on ∆_r of all R competing models. Similarly, a weighted estimation of the target parameters is obtained weighting the estimations of each model r using its posterior probability.

The Present Study

Random effects should be understood as informative parameters of the target processes under study, making them a relevant substantive part of the results report (Barr, 2013 presented a similar rationale on the use of random effects from a confirmatory perspective). Large intercept variances for subjects or items would indicate that there are important differences in the mean (e.g., some participants could be slower than others, or some items could be more difficult than others). Similarly, large slope variances for subjects or items would indicate that there are important differences in how the fixed effect and the dependent variable are related (e.g., some participants or items could be more sensitive to the experimental conditions than others). Thus, random effects should not be understood as of secondary interest, but as relevant and substantive information about the modeled processes. That is one of the main reasons to study the performance of model averaging, which aims to deal with model uncertainty, and to compare it with model selection all-or-nothing decision making about the inclusion of such random effects. In this line, it is known that the estimation of fixed effects does not present bias regardless of whether the random effects of the fitted model are correct or not, and that they are very similar across different random structures (Hoffman, 2015; Hox et al., 2018; Meyers & Beretvas, 2006). Previous research found that the estimation of the standard errors of fixed effects was more efficient for model averaging than for model selection (Martínez-Huertas et al., 2022). We think that such efficiency is related with the consideration of random effects.

In this line, there is a handicap when researchers want to recover the random structure of MEMs-CR using model averaging. This procedure uses all the available information from the competing models but, unfortunately, some models do not include all the parameters of interest. This is the case of random effects. Researchers could try to average random intercepts or random slopes when some of the competing models do not include such parameters by weighting null information (a zero) in model averaging. Thus, the usefulness of model averaging could be compromised recovering different random structures. The present study aims to evaluate the bias of averaged crossed random effects in simulation scenarios where the random structure of various competing models is incomplete. For this purpose, we are going to compare the average estimates of the two model averaging strategies introduced previously (Akaike weights and BMA with BIC posterior probabilities), with model selection of AIC and BIC indices as benchmarks. We think that there will be relevant differences between both perspectives because model averaging is supposed to deal with model uncertainty while model selection always does an all-or-nothing decision making. Given that all the MEMs-CR of the study are nested, the maximal MEM-CR (Barr et al., 2013; which is a MEM-CR with random intercepts and random slopes for both subjects and items in this study) was used as an additional benchmark to compare with the averaged estimations. All these conditions were studied in two different simulation scenarios under the presence of symmetric and asymmetric variances of random intercepts and slopes of subjects and items, that is, when both random effects present the same or different variability. Although the simulation scenario of symmetric variances of random effects is more artificial, it is useful to analyze if the quality of the estimations of variances of random intercepts and random slopes tends to be similar. The simulation scenario of asymmetric variances of random effects is more ecologic because of usually the random effects present different variabilities in real data sets.

Simulation Study

An experimental design in psycholinguistics was simulated using Equation 1 as the generating model. Specifically, response time (ms) of multiple items by different participants were analyzed in two different conditions (control vs. experimental). Thus, a within-subject effect was simulated here. Two different sources of variability were simulated: participants and items. These random effects were fully crossed due to all the participants answered all the items. In our simulation study, we used similar fixed and random effects population parameters as those in previous related work (Baayen et al., 2008; Barr et al., 2013; Matuschek et al., 2017; Martínez-Huertas et al., 2022). Table 1 presents the simulation parameters. The within-subject effect was set to 0 in all simulation conditions (fixed effects do not present bias as has been evidenced previously: Hoffman, 2015; Hox et al., 2018; Meyers & Beretvas, 2006). The mean of the control condition (intercept) was set to 2,000 ms. The residual level-1 variance ${\mathrm{\sigma }}_{\mathrm{e}}^2$ was 90,000, so a 300 standard deviation should be understood as the residual level-1 error. Given that the focus of this study was the study of random effects, almost all the random structure of this model was manipulated in the simulation. The cluster-specific random effects for subjects and items intercepts are represented by U_0s0 and U_00i, respectively, and their variances ${\mathrm{\sigma }}_{\mathrm{0s0}}^2$ and ${\mathrm{\sigma }}_{\mathrm{00i}}^2$, respectively) were set to 0, 10,000 and 40,000 ms (standard deviations equal to 0, 100 and 200 ms). The subject and item cluster-specific random slopes are represented by U_1s0 and U_10i. Two different variances ${\mathrm{\sigma }}_{\mathrm{1s0}}^2$ and ${\mathrm{\sigma }}_{\mathrm{10i}}^2$, respectively) were simulated for subjects and items random slopes (0 and 10,000 ms; standard deviations equal to 0 and 100 ms). An intercept-slope covariance for subjects ( ${\mathrm{\sigma }}_{U_{0s0},U_{1s0}}$) and items ( ${\mathrm{\sigma }}_{U_{00i},U_{10i}}$) was simulated equivalent to an intercept-slope correlation equal to 0 and .60. The random effects and the error term were normally distributed. Different sample sizes were simulated for subjects and items (sample size for subjects: 15, 30, 60, and 90 subjects; and sample size for items: 5, 10, 30, and 60 items). A total of 1,152 simulation conditions were considered in the present study, and 1,000 replications were generated for each simulation condition using R software. This makes a total of 1,152,000 replications.

Data Analysis

We assume that MEMs-CR naturally have random intercepts for both subjects and items. Also, we consider that the most complex MEM-CR will be the one with random intercepts and random slopes for both subjects and items. The former is called here minimal, whilst the latter is called here maximal MEMs-CR. Two intermediate random structures were also considered here: a MEM-CR with random intercepts for both subjects and items and random slopes for subjects (here, subject random slopes), and a MEM-CR with random intercepts for both subjects and items and random slopes for items (here, item random slopes). Thus, four MEMs-CR were fitted to the data using the restricted maximum likelihood (REML) approach. This makes a total of 4,608,000 analyses (1,152,000 replications x 4 MEMs-CR). All MEMs were fitted with the lme4 package (Bates et al., 2015) in R software. lme4’s AIC and BIC fit indices were used to compute model selection and model averaging. Approximately, 99.79% of the estimated models converged. Please note that the simulated data was analyzed with random intercepts even when some simulation conditions set those parameters to zero to emulate a natural analytical strategy, and to compare the recovery of random intercepts and random slopes under the same simulation conditions.

First, model selection performance of AIC and BIC fit indices was evaluated. It was computed as the proportion of true model selection: if the model used to simulate the data was selected, it was considered a correct model selection. Otherwise, it was considered as incorrect (please note that selecting the minimal MEMs-CR was also considered correct in those simulation conditions that set intercepts to zero to ease the simulation design). In the present study, these results are useful to contextualize the mean true model selection of each fit index as benchmarks.

Second, the estimations of model averaging with Akaike weights and BMA with BIC posterior probabilities were computed. For this purpose, we applied the above-mentioned procedures using the four competing models of this simulation study, namely: minimal, subject random slopes, item random slopes, and maximal MEMs-CR. Specifically, we computed the average estimations of the random effects.

Third, the bias of the estimations of the variances of the random effects was computed using different measures. Bias was computed for all the random effects using the following formula:

where θ is the population parameter and ${\widehat{\mathrm{\theta }}}_r$ is the sample estimate of each R replicate. The interpretability of bias was done considering the simulated random variance in not-null variance conditions, that is, calculating the percentage of bias regarding the 100 or 200 simulated variances. Root mean squared error (RMSE) was also computed for all the random effects using the following formula: where θ is the population parameter and ${\widehat{\mathrm{\theta }}}_r$ is the sample estimate of each R replicate. To ease the interpretability of the results, the analyses were conducted in different simulation scenarios depending on (1) symmetric vs. asymmetric simulated population variances of random intercepts and slopes; and (2) zero vs. not-zero simulated population random effects.

First Simulation Scenario: Symmetric Variances of Random Intercepts and Slopes

Second Simulation Scenario: Asymmetric Variances of Random Intercepts and Slopes

In the second simulation scenario, asymmetric variances of random intercepts and slopes (200 and 100 standard deviations, respectively) were considered for the two clusters (subjects and items). Additionally, we only considered the estimations of not-null variances because of their results did not present large differences with those of null symmetric variances of random effects. Thus, we only present the estimations of not-null asymmetric variances of random effects of subjects and items in model averaging and model selection in this section for the sake of brevity. Figure 4 presents the mean of bias and RMSE of the estimations of not-null asymmetric variances of random effects of subjects and items in model averaging and model selection with AIC and BIC.

First, it was found an absence of bias of the estimations of the random variances of the intercepts of subjects and items. In the most demanding conditions, the most severe bias was found to be around the 6% (that is, approximately, -12.5 bias was for the estimation of the 200 standard deviations). Similarly, we found more bias of the estimations of random intercepts in items than in subjects. The comparison of bias and RMSE of the random intercepts of subjects and items show that there is a relevant sample variance of the estimates that decreases as the sample sizes of both clusters increase.

Second, important bias was found for the estimations of the variances of the random slopes in the more demanding conditions. Again, the model averaging and the model selection derived from the BIC index presented a significantly worse performance than the ones of AIC index (the results of both indices are only comparable when there are larger sample sizes in both clusters). Given that the simulated random slopes had a standard deviation of 100, we can easily interpret the results as the percentage of bias. Bias was large in all the conditions where the sample sizes of both clusters were small, that is, for 5–10 items and 15–30 subjects. When one of the clusters presented a medium/large sample size (e.g., 60–90 subjects or 30–60 items), the recovery of the variances of both random effects was accurate (bias was around 8%), and the results were near (but larger than) the 10% of bias even in the more extreme conditions of small sample sizes of items. On the contrary, none of the conditions with 5 items reached an appropriate level of bias (being the bias of subjects random slopes larger than the bias of items random slopes), and the 10 item conditions presented a similar pattern of results with lower bias. Again, RMSE of the random slopes of subjects and items decrease as the sample sizes of both clusters increase. In this line, it is worth mentioning that model averaging with Akaike weights presented a significantly lower bias and RMSE than model selection with AIC index in the most demanding conditions, that is, when the sample sizes were smaller in both clusters. This means that the bias of model averaging with Akaike weights, although relevant, was lower than the one of model selection with AIC, and that the sample variance of the estimates followed the same pattern of results.

Summary of Findings

We compared the bias of the average estimations of random effects of MEMs-CR using Akaike weights and BIC posterior probabilities, using AIC and BIC model selection and fitting a full random structure (here, maximal model) as benchmarks. Additionally, two simulation scenarios were considered where the random variances of intercepts and slopes of subjects and items were symmetrical or asymmetrical. The first simulation scenario was more artificial but allows to study the recovery of the variances of random intercepts and slopes with the same metric, while the second simulation scenario was more ecological as usually the variances of random intercepts tend to be larger than the ones of random slopes in experimental psychology. In general, AIC index obtained higher true model selection performances than BIC in different demanding conditions like having small sample sizes in any of the clusters (subjects and/or items), but incorrect random structures were selected in many cases (specially, under-parametrizing the random structure, that is, setting random effects equal to zero when they are different from zero). This was found to be crucial to comprehend the advantages of model averaging in front of model selection.

In scenarios with null random effects, no relevant differences were found between model averaging and model selection. Bias was larger in more demanding conditions (e.g., small sample sizes), but it was not alarming. This means that both model averaging and model selection provide virtually unbiased estimates of population random effects that do not differ from zero, once minimal sample sizes are raise. Presumably, this simulation scenario is not very common in empirical research.

In scenarios with not-null random effects, which are, presumably, the most common empirical scenario, interesting differences were found between model averaging and model selection. First, no relevant differences were found between them when random intercepts were estimated. This means that random intercepts were estimated unbiasedly in almost all conditions of the simulation study, that is, including the presence of symmetrical and asymmetrical variances of random intercepts and slopes. Second, some differences were found between model averaging and model selection when random slopes were simulated, especially in random slopes of subjects. That is, given a minimum sample size, model averaging could estimate random slopes more accurately than model selection (the differences between the two strategies mainly appeared in random slopes of subjects, but their differences were less important for random slopes of items). But the main advantage of model averaging was using all the available information while model selection made all-or-nothing decisions that were incorrect on many occasions (in some scenarios, we even found a median of the 100% percent of bias in random slopes when a null variance was settled to a cluster). These differences were dependent of the smaller sample sizes, as they probably limited the available information for the fit indices, being model averaging capable of overcoming, at least in part, this limitation.

Some relevant differences were also found between the versions of model averaging, that is, Akaike weights and BIC posterior probabilities. In general, it was found that Akaike weights obtained a better performance than BIC posterior probabilities when there are not-null random effects (which is, in fact, the most common scenario in applied research). Again, once minimum sample sizes were raised, Akaike weights could estimate unbiased random effects while BIC posterior probabilities obtained larger bias. In any case, the estimations of both versions of model averaging were affected by lower sample sizes in any of the simulated clusters.

Additionally, interesting differences were found between model averaging with Akaike weights and fitting a full random structure (here, maximal model). Model averaging was capable of recovering the random effects under conditions with large sample sizes in the target clusters. But fitting a full random structure (without using a model selection strategy; see for example Barr et al., 2013; Martínez-Huertas et al., 2022; Matuschek et al., 2017 for different model selection strategies) was found to present large bias in the random effects that were analyzed in this simulation study if there are null population variances for random effects. Thus, fitting a full random structure (without a previous model selection) would lead to large bias in the estimation of random effects. Probably, the variability of the different clusters would be incorrectly distributed among the estimated random effects (e.g., the model could not be able to determine if the variability of a cluster is related to the intercepts or the slopes, and also the model could mix the variability between clusters). In any case, the estimates of fitting a full random structure were accurate when all the estimated random effects exist at the population level, which is not always known by applied researchers. Thus, we would like to endorse the use of model selection and model averaging instead of just fitting a full random structure to the data. In this sense, model averaging would consider all the available information of the different competing models based on their model fits, being more robust against model misspecification (see also Hinne et al., 2020) than fitting a full random structure.

Theoretical and Methodological Considerations

Applying a MEM-CR assumes decisions about their parametrization given that true population models are not known. Model averaging proposals use all the available information of different models to deal with uncertainty. Two general conclusions can be made from our analyses: model averaging and model selection show similar results for recovering null random effects, and model averaging shows less bias than model selection for recovering not-null random effects, especially for Akaike weights and random slopes of subjects. These conclusions are conditioned to minimum sample sizes in both clusters (that is, subjects and items). Thus, we would like to endorse model averaging with Akaike weights as a relevant tool to recover random effects in studies with medium to large sample sizes.

There is consensus about the necessity of including all relevant dimensions in experimental conditions to potentiate ecological validity (Hoffman, 2015). Thus, the importance of these methodological tools that allow to deal with the uncertainty of random effects is significantly emphasized. Therefore, model averaging approaches could have important theoretical implications for empirical research because they can avoid errors of statistical inference and to unbiasedly estimate random effects. Thus, model averaging allows to use the evidence in favor of each possible model without losing relevant information of their random effects by making questionable decisions about the random structure of MEMs-CR. The empirical illustration of this study aims to provide an example of the applicability of model averaging to the estimation of crossed random effects for applied researchers.

Limitations and Future Directions

The present study simulated different random structures for MEMs-CR, but only two population parameters (a null vs. not-null value) were considered to analyze the recovery of random effects. Also, a simple experimental design with two fixed effects (the intercept and the within-effect) was simulated due to our objective was to evaluate the recovery of random effects using model averaging. However, we think that model averaging would be an affordable tool to recover small random effects, which are more difficult to estimate, avoiding all-or-nothing decisions of model selection strategies. Also, given the importance of the design of the study to establish model random structures, it would be interesting to expand these results to other relevant designs like longitudinal studies (see Martínez-Huertas & Ferrer, 2022, for an illustration of the use of MEMs-CR in different longitudinal designs). Other complex designs, like models that use cross-level interaction effects or level-2 predictors, could be benefited of using model averaging approaches because these effects are largely influenced by the variability around the fixed effects. Moreover, these complex models would require considering more competing models that: (1) would generate a very complex model selection scenario, and (2) could bias all the results of fixed effects if incorrect random structures were selected (e.g., setting some random effects equal to zero). Also, future research should explore other average estimates like evidence ratios, which can be understood as the evidence in favor of a model as the ratio of the weights (Burnham & Anderson, 2002). Additionally, some readers might find some parallelism of model averaging with ensemble methods and other statistical procedures that combine the predictions of different models or estimators into a single estimate to improve the generalizability or robustness of the predictions of algorithms. Model averaging and other ensemble methods share a common underlying logic of computing a single estimate from different sources of information. In this study, we analyzed the usefulness of model averaging using AIC and BIC fit indices in MEMs-CR, but the research area of ensemble methods is more broaden. In this sense, the rationale of the present simulation study could be generalized to other interesting statistical procedures that use penalized estimation criteria to compute model estimations like penalized least squares for structural equation modeling (Huang, 2022), or regularized partial correlation networks (Epskamp & Fried, 2018).

Conclusions and Recommendations

Model averaging attempts to use all the available information of the competing models to deal with model uncertainty, beyond model selection based on AIC and BIC. Unbiased estimates of random effects were found in both model averaging with Akaike weights and BMA with BIC posterior probabilities under conditions with sufficient sample sizes in the target clusters. In general, 60 units per cluster were found to be an appropriate sample size to obtain very accurate estimations of the variances of random effects in the conditions of the present simulation study. But some differences were found in favor of model averaging with Akaike weights in demanding conditions like small sample sizes, which was also capable of estimating unbiased variances of random effects if one of the clusters presented large sample sizes. This means that larger number of subjects (e.g., 60 subjects) could lead to compensate the handicaps of smaller sample sizes of items (e.g., 10 items). Whilst the performance of model averaging is questionable under some simulation conditions, it supposes an alternative to deal with model uncertainty even in those scenarios where using model selection would require a risky all-or-nothing decision regarding to the inclusion of parameters like random slopes. Thus, we recommend using model averaging Akaike weights, and to use both model averaging approaches with small sample sizes to analyze their convergences and divergences.