Parametric and Semi-Parametric Bootstrap-Based Confidence Intervals for Robust Linear Mixed Models

The Linear Mixed Model

For $i=1,\dots ,n$, let the model equation for the LMM be:

where ${\mathit{y}}_i$ is a vector of length $J_i$ containing the responses of participant $i$, $\mathbf{\gamma }$ is a vector of coefficients for the $p$ fixed effects, $X_i$ is the ( $J_i$ $\times$ $p$) design matrix for fixed effects, ${\mathit{b}}_i$ is the vector of random effects of length $q$, independent of the errors ${\mathbf{\epsilon }}_i$ and $Z_i$ is the ( $J_i$ $\times$ $q$) design matrix for the random effects. In the LMM language, fixed effects refer to regression coefficients that are constant between participants and thus do not need the subscript $i$ (hence $\mathbf{\gamma }$ rather than ${\mathbf{\gamma }}_i$ in Equation 1). In contrast, random effects refer to quantities that can vary between participants and thus are designated with the subscript $i$ (e.g., ${\mathit{b}}_i$ rather than $\mathit{b}$). The vector ${\mathbf{\epsilon }}_i={\left({\epsilon }_{i1},\dots ,{\epsilon }_{i{J}_i}\right)}^T$ contains the $J_i$ individual error terms for participant $i$. Generally, it is assumed that random effects and error terms are normally distributed around zero (Laird & Ware, 1982):

where $\text{Σ}=\text{Σ}\left(\mathit{\theta }\right)$ is the $q\times q$ (parametrized) covariance matrix of the random effects and ${\sigma }_{\epsilon }^{2}I$ is the $\left(J_i\times J_i\right)$ (diagonal) covariance matrix of the error term. That is, the model assumes that the random effects ${\mathit{b}}_i$ and errors ${\mathbf{\epsilon }}_i$ stem from centered normal distributions and therefore are gathered around zero.

Robust Estimation

Alternatively, robust methods follow the frequently called central model (see e.g., Koller, 2013) and consider that the data generating process is:

where $F$ is the central model assumed for the majority of the data (e.g. the normal distribution), $H$ is an arbitrary (noncentral) unknown model, different from $F$, and $\delta$ $\in$ [0 ; 1] is the amount of contamination, that is, the proportion of data not from $F$. In this approach, the goal of the analysis is to reach inferential conclusions that are valid for the central model $F$ only and are not influenced by the outliers from $H$. Commonly, this approach involves estimation weights, which are computed so that they are near 1 for the data from $F$ and smaller (with a lower bound at 0) for the data from $H$. There are several robust methods for the estimation and inference in the LMM based on the “central model”. Here, we consider the three methods compared in Agostinelli and Yohai (2016): the S-estimator (Copt & Victoria-Feser, 2006), the composite $\tau$ estimator (cTAU, Agostinelli & Yohai, 2016) and Koller's (2016) DAStau method. These methods presume that the model defined in Equation 1 and Equation 2, which assume normality for ${\mathit{b}}_i$ and ${\mathbf{\epsilon }}_i$, is true for a majority of the population (thus is the central model), but does not hold for outliers.

Conceptually, we can consider the robust methods as weighted versions of their classical counterparts, which constitute a special case where all observations’ weights are equal to 1. Broadly speaking, the weights (between 0 and 1) produced by the robust methods express how likely an observation is to belong to the bulk of the data (and, conversely, how unlikely it is to be an outlier), and the degree to which it is taken into account in the parameter estimation of the central model. The weights can be attributed at the single observation level only (as for cTAU), at the participant level only (S), or at both levels (DAStau).

Robust estimators may differ with respect to their priorities. Bounded influence estimators, such as DAStau, aim at limiting the bias in their parameter estimates, whereas high-breakdown estimators, such as S and cTAU, strive to tolerate high proportions of outliers. However, this latter class typically produces estimators that are less efficient than those from the former category (Huber, 1964).

Below we provide a brief summary of the robust estimators for the LMM (see also section ‘Technicalities of the Estimation Methods’ in the Supplementary Materials). For a detailed treatment of robustness, we refer interested readers to Maronna, Martin, and Yohai (2006).

Parametric Bootstrap

First, the procedure estimates fixed effects ( $\mathbf{\gamma }$) and variance components $\left({\sigma }_{\epsilon }^2,\mathit{\theta }\right)$ parameters from model (1) on the original data $\left({\mathit{y}}_i,X_i,Z_i;i=1,\dots ,n\right)$ to produce $\widehat{\mathbf{\gamma }}$, ${\widehat{\sigma }}_{\epsilon }^2$, and $\text{Σ}\left(\widehat{\mathit{\theta }}\right)$. Then, for each of the $B$ bootstrap samples:

The estimator used at step 2 can either be the same as, or differ from, that used on the original sample. Because we generated $B=5000$, and robust LMM estimation can be computationally extremely intensive¹, we choose to apply ML at step 2 also when a robust estimator was used at step 1. We believe this procedure legitimate because the bootstrap samples are unlikely to contain outliers.

Wild Bootstrap

The wild bootstrap method replaces step 1 of the parametric bootstrap computing “residual” ${\tilde{\mathit{\upsilon }}}_i$ as follows (Modugno & Giannerini, 2013):

with $H_i=X_i{\left(X^{T}X\right)}^{-1}{X}_i^T$, where the operator “ $\circ$” denotes the (element-wise) Hadamard product.

For $i=1,\dots ,n$, sample independently $w_i^{*}$ from the following standardized random variable (Mammen, 1993):

According to Modugno and Giannerini (2013), these weights produce better results than the simple weighting scheme usually used with classical linear models, consisting in sampling either 1 or -1 with equal probability of .50. Then construct individual responses

Uncontaminated Data ( $C_0$)

We first generated data without contamination from the following model:

with $\mathbf{1}$ a vector of 1’s of dimension $J_i$, $\mathbf{\gamma }={\left(250,10\right)}^T$ and with the elements of ${\mathit{x}}_i$ representing the day of measurement (hence, from 0 to $J$-1, with $J_i=J$). We additionally assumed

and

We then generated contaminated data under two different conditions: with outlying errors ${\mathbf{\epsilon }}_i$ ( $C_{{\mathbf{\epsilon }}_i}$) and with outlying random slope effects $b_{1i}$ ( $C_{b_{1i}}$), as described below. In the Supplementary Materials (see Figures S7 and S10), we also provide results with outlying random intercept effects $b_{0i}$ ( $C_{b_{0i}}$), but do not discuss them here, because they were analogous to those for $C_{b_{1i}}$.

Contamination of the Errors ( $C_{{\mathbf{\epsilon }}_i}$)

This condition simulated within-participant outliers, with outlying single observations, randomly distributed across participants. While we assumed that a proportion (1- $\delta$) of the generated data points follow Equation 9, we also assumed that the remaining proportion $\delta$ had outlying individual error terms, so that the full data generating mechanism was defined by Equation 7 and Equation 8 with ${\mathbf{\epsilon }}_i$, such that

We centered the contaminated distribution at −80, thus four SDs ( $=4\sqrt{400}$) below zero, because in the sleepstudy application the ML estimate of ${\sigma }_{\epsilon }$ was $25.59$ and the smallest residual was $-101.18$ (hence roughly 4 SDs below 0). We arbitrarily chose the variance of the contaminated errors at $0.25$ to limit their variability. We expect degraded performance of CIs associated with ML for ${\gamma }_0$ and ${\sigma }_{\epsilon }^2$. For all methods, effects of within-participant outliers should be more deleterious when $J=4$ than $J=8$. This should especially be the case for ${\sigma }_1^2$, because of the greater influence that a single point out of four, compared to eight, can exert on the estimation of this parameter. In general, we expect these effects to be greater in classical than in robust estimators. Finally, because S downweights data only at the participant-level, we expect CIs to have lower coverage than when associated to cTAU, which is also a high-breakdown estimator but weights data at the observation level and thus should be more efficient.

Contamination of Random Slope Effects ( $C_{{\mathit{b}}_{1i}}$)

This condition implemented between-participant outliers, where all values of a participant were particularly deviant. To produce a $\delta$-proportion of outliers at the random slope level, we modified the mean of their distribution by shifting it four SDs to the left, generating data from Equation 7 and Equation 9, using the following distributional assumptions for the random effects:

We expect contamination effects only on the slope parameters. In particular, ${\gamma }_1$ should be underestimated, because of the negative contamination, whereas ${\sigma }_1^2$ should be overestimated, because of an overestimation of the between-individual variability on the slope reducing the coverage rates. Again, across all methods, the estimation of the individual slopes should improve with $J=8$ compared to $J=4$. Indeed, with only four individual observations there is less signal to determine the atypical nature of a trajectory.

Study Design

Our contamination method is inspired by Richardson and Welsh (1995), and reproduces the generating process of the central model: most of the data follow a normal distribution, whereas outliers are generated from an alternative normal distribution. Ultimately, our design contained 2 ( $\delta$ = 0.05, 0.1) $\times$ 2 ( $J$ = 4, 8) $\times$ 3 ( $C_0$, $C_{{\mathbf{\epsilon }}_i}$, $C_{b_{1i}}$) conditions, for a total of 12 simulation cell conditions.

For each condition we generated $250$ data sets (which, according to Equation 1 of Morris, White, and Crowther (2019), is generally enough to obtain a sufficient degree of precision), for a total of $3000$ data sets, each to which we applied the possible combinations of the four estimators and the 2 CI estimation methods, producing 8 types of CIs. For the percentile CIs based on bootstrap, we generated $B=5000$ bootstrap replications (cf. Deen & de Rooij, 2019) and all CIs were calculated for a .95 confidence level.

Simulation Analysis

For each of the six parameters ( ${\gamma }_0$, ${\gamma }_1$, ${\sigma }_{\epsilon }^2$, ${\sigma }_0^2$, ${\sigma }_1^2$, ${\sigma }_{10}$), we comparatively assessed the coverage of each type of CI. For each parameter, in each cell condition, we defined the actual CI coverage as the proportion of the 250 CIs including the true population value. It is important that the observed coverage approaches the nominal coverage (e.g., .95) to allow for sound inferential conclusions.

Uncontaminated Data ( $C_0$)

For the fixed effects, when associated to ML, no major differences in coverage emerged between the parametric and the wild bootstrap as we can see on Figure 1. However, when associated to robust estimators (mostly S and cTAU), the wild bootstrap obtained better coverage than the parametric bootstrap for ${\gamma }_1$, whereas coverage was similar for ${\gamma }_0$. Indeed, as we can see on the upper major panel of Figure 2, for ${\gamma }_1$, CIs were much larger with the wild than with the parametric bootstrap, resulting in a lower coverage for the latter when associated to S and cTAU. Also, ML (with parametric and wild bootstrap), cTAU and DAStau (both with wild bootstrap only) had coverage close to the nominal threshold for both fixed-effect parameters in the uncontaminated condition.

For the variance components (see Figure 1), again, the two types of bootstrap when associated to ML obtained coverage close to the nominal threshold, with a slight difference for ${\sigma }_{\epsilon }^2$, where the parametric bootstrap was closest to .95. The wild bootstrap obtained better coverage when associated to the robust estimators, particularly for ${\sigma }_1^2$. The lower panel of Figure 2 reveals that the bias estimation of ${\sigma }_1^2$ with S, cTAU, and also, when $J$ = 4, with DAStau, affected the parametric bootstrap method and explained the coverage differences noted above.

Contamination of ${\mathbf{\epsilon }}_i$ ( $C_{{\epsilon }_i}$)

This contamination condition affected mainly ${\gamma }_0$ and ${\sigma }_{\epsilon }^2$, whereas results for the remaining four parameters were similar to those observed in $C_0$ (mainly for ${\gamma }_1$, ML, cTAU and DAStau obtained closer to nominal coverage than S, especially when associated to the wild bootstrap; for the remaining variance components, all estimators performed well when associated to the wild bootstrap; cf. Figures S5 and S8). Figure 3 thus displays coverage CIs for ${\gamma }_0$ and ${\sigma }_{\epsilon }^2$ only. ML obtained the worst coverage for ${\gamma }_0$ in all conditions (cf. upper panel of Figure 3). When $\delta =0.05$, the two bootstrap methods obtained the highest coverage for ${\gamma }_0$ with S and cTAU, with slightly better coverage for S with the parametric bootstrap, followed by cTAU with the wild bootstrap. When $\delta =0.1$, all coverages for ${\gamma }_0$ fell drastically, except for S when $J=4$. The upper panel of Figure 4 reveals that differences in coverage for ${\gamma }_0$ are principally explained by differences in estimation, with more bias for ML than for other estimators, and when $\delta =0.1$ compared to $\delta =0.05$. As in $C_0$, the two bootstrap methods behaved similarly for ${\gamma }_0$.

Coverage rates for ${\sigma }_{\epsilon }^2$ were always exactly null except for S and cTAU with parametric bootstrap when $\delta =0.05$ and $J=4$ (but rates were close to .10 only; cf. lower panel of Figure 3). We can see on the lower panel of Figure 4 that all estimates were biased. Overall, bias was larger with $\delta =0.1$ than with $\delta =0.05$, parametric bootstrap CIs were shorter than with wild bootstrap, and for all estimation methods, variability in estimation was higher when $J=4$ than when $J=8$.

Contamination of $b_{1i}$ ( $C_{b_{1i}}$)

This contamination resulted mainly in effects on CIs for ${\gamma }_1$ and ${\sigma }_1^2$ (see Figure 5 and Figure 6). Again, for the other parameters, results were similar to those obtained in the other simulation conditions and are available in the Supplementary Materials (see Figures S6 and S9), which also contains results when contamination occurred on $b_{0i}$, affecting coverage rates on ${\gamma }_0$ and ${\sigma }_0^2$ (see Figures S7 and S10 for $C_{b_{0i}}$).

Generally, the wild bootstrap outperformed the parametric bootstrap with robust estimators and, to a lesser extent, also with ML. Again, ML obtained the lowest coverage for ${\gamma }_1$, especially when $J=8$. When $\delta =0.05$, S obtained coverage close to the nominal threshold, but only when $J=8$. cTAU and then DAStau followed next in coverage. Again, coverage dropped drastically when $\delta =0.1$, but S retained the highest rates. Interestingly, only robust estimators’ CIs had better coverage when $J=8$ than $J=4$. The upper panel of Figure 6 shows that all CIs were shorter when $J=8$ than $J=4$, especially with the parametric bootstrap, which may explain the reduced coverage for ML when $J=8$. On the contrary, with robust estimators associated to the wild bootstrap, especially when $\delta =0.05$, CIs appear slightly shifted to the right, towards the true value, thereby resulting in higher coverage.

For ${\sigma }_1^2$, all coverage rates were extremely low (see lower panel of Figure 5). The lower panel of Figure 6 shows that parametric CIs with S and cTAU were very close to zero and extremely small, resulting in nearly null coverage. Again, CIs were larger when $J=4$ than when $J=8$ and closer to the expected value when $\delta =0.05$ than $\delta =0.1$, explaining the best coverage rate with the wild bootstrap when $\delta =0.05$ and $J=4$.

Major Statistical Considerations

First, we found that the contamination effects were different across within- and between-participant outliers, but in both cases CI coverage was mainly reduced for specific estimates. For the former, CIs of ${\gamma }_0$ and ${\sigma }_{\epsilon }^2$ were strongly affected, whereas for the latter, CIs of ${\gamma }_1$ and ${\sigma }_1^2$ were impacted under $C_{b_{1i}}$ (and analogously, for ${\gamma }_0$ and ${\sigma }_0^2$ under $C_{b_{0i}}$). But across all contamination conditions, coverage of CIs of the remaining parameters were similar to those obtained without outliers. And, unsurprisingly, the more the outliers, the stronger the effects of the contamination on CI coverage rates.

Second, with ML estimation and uncontaminated data, both parametric and wild CIs obtained excellent coverage across all parameters. Moreover, wild CIs were larger for variance components, especially for ${\sigma }_{\epsilon }^2$, with rates close to 1, as in Modugno and Giannerini (2013). On the contrary, with robust estimators, the wild bootstrap nearly always outperformed the parametric method in terms of coverage, particularly for ${\gamma }_1$ and the variance components.

Third, without contamination, CI coverage rates with cTAU and DAStau were similar to those of ML for fixed effects and variance components parameters. With S, coverage was lower for the fixed effects parameters, but similar for the variance components. Thus, with uncontaminated data, ML with both bootstrap methods, and cTAU and DAStau with the wild bootstrap produced the best CIs.

Fourth, our results confirm that in the presence of relatively many (i.e., 5% or 10%) outliers, CIs for ML estimates can obtain very low coverage rates. As expected, with 5 $\%$ of outlying data, robust estimators obtain higher CI coverage than ML. Consistently with Agostinelli and Yohai (2016), who observed that DAStau fixed effect estimates were less accurate than those of S and cTAU under within- and between-cluster outliers, we found that the wild bootstrap obtained slightly lower coverage with DAStau than with S and cTAU for the parameters affected by the contamination condition. With 10 $\%$ of outliers in the data, only S maintained a high coverage for ${\gamma }_0$ under $C_{{\mathbf{\epsilon }}_i}$ and $C_{b_{0i}}$ when $J=4$. But under $C_{b_{1i}}$, the coverage for ${\gamma }_1$ dropped across all estimators, even if it still obtained the best performance. For the variance components affected by contamination, coverage were very low for all methods, independently of the proportion of outliers.

In conclusion, in the realm of our simulation, when there are no or a few outliers, the wild bootstrap associated to the cTAU estimator provides excellent coverage generally for all parameters, except for the variance component parameters particularly affected by the contamination condition. However, when the proportion of outliers becomes larger, CIs obtained with the S estimator and the wild bootstrap method appear superior.

Limitations, Extensions, and Further Applications

Motivated by the empirical sleepstudy example that is often used as a didactic example for R (Bates et al., 2015; Koller, 2016), we based the population values of our simulation on LMM results from this data set. A characteristic of the empirical LMM results was the nearly null intercept-slope covariance estimate ( ${\sigma }_{10}=-8.5$, implying $r_{10}=-.05$). Of the results we obtained, we are most wary of those about this parameter. As usual, the extension of the conclusions beyond the specific simulation conditions has to be made cautiously.

Also, our example enjoyed balanced data, in that each participant was assessed at the same occasions and had complete data. Two of the robust estimation methods we examined cannot be applied to unbalanced data (namely S and cTAU), so that analysts with such data cannot apply all robust estimators to the LMM. Extending S and cTAU to handle unbalanced data would greatly increase their applicability. Nevertheless, based on the current results, it appears that the use of DAStau with the wild bootstrap when analyzing unbalanced data that possibly contain outliers might be useful. Indeed, this method obtains CIs for most parameters with larger coverage than those based on classical estimators, with and without outliers.

The example did not include explanatory covariates, other than time. Typically, researchers are interested not only in testing for sample heterogeneity in intercept and slope values, but also in trying to explain such interindividual differences by means of person-specific covariates. For instance, in LMM applications to experimental data a main covariate is group membership, which distinguishes participants in the treatment from those in the control group. It seems plausible that the presence of outliers in the response variable may differ between the groups, to the point of possibly clouding the treatment effect. Knowledge about robust estimation of the LMM in such situations should further increase the attractiveness of applying the LMM to experimental data.

Conclusions

Nowadays, several methods for the robust estimation of the LMM are available. In this work, which we believe to be the most comprehensive on CIs with robust estimation in the LMM, we have shown that cTAU with the wild bootstrap can produce sound inferential conclusions, both with and without data contamination.

We are not claiming that all repeated-measure (or, more generally, nested or crossed) data sets necessarily contain outliers that bias results when classical estimation methods are used. Nevertheless, we hope to have raised awareness about the possible detrimental effects of outliers in the context of the LMM and to have provided useful suggestions about alternative estimation methods for such data. We provide R code (see Supplementary Materials) to implement these methods, accompanied by detailed comments, so that most users familiar with the R language and environment ought to be able to use our code to estimate their LMMs when suspicion arises about the presence of outliers. Our wish is that research in this domain continues to provide tools aimed at reducing the potentially detrimental effects of outliers, and that researchers may find the applications of such tools meaningful to their work.