Which Robust Regression Technique Is Appropriate Under Violated Assumptions? A Simulation Study

Least-Square-Fit-CI (lsfitci) Approach

One approach to dealing with violation of the normal-error assumption is bootstrapping residuals ( ${\epsilon }_i$) instead of raw scores (x and y) to generate empirical distribution and standard error of residuals for the construction of $1-\alpha$ CI, where $\alpha$ is the type 1 error rate. By locating the $\alpha /2$ and $1-\alpha /2$ percentiles of the bootstrap residuals, it is expected that the CI could be asymmetrical, thereby adjusting for any non-normal residuals. This method is known as “least square fit CI” (lsfitci) in Wilcox (2022), and its performance is evaluated in this study.

Heteroscedastic-Consistent (HC) Standard-Error Approaches

HC standard-error approaches (Huber, 1967; Long & Ervin, 2000; White, 1980) can be used to fit a regression model that contains heteroscedastic errors. Among them, HC3 method was developed for studies with large sample sizes (Wilcox, 2022). This method can effectively replace the bootstrap standard error estimate by $\text{HC}3=S={\left(X\text{'}X\right)}^{-1}{X}^{\prime }diag\left[r_i^2/{\left(1-h_{ij}\right)}^2\right]X{\left(X^{\prime }X\right)}^{-1}$, where $r_i$is the residual for i = 1,..,n, $h_{ij}={\mathbf{x}}_{\mathbf{i}}{\left(\mathbf{X}\mathrm{\text{'}}\mathbf{X}\right)}^{-1}{\mathbf{x}}_{\mathbf{i}}^{\mathrm{\text{'}}}$, and

where ${\mathbf{x}}_{\mathbf{i}}$ is the ith row of X_i. Consequently, the $1-\alpha$ CI surrounding ${\beta }_p$ is $b_j\pm t{S}_j$, where t is the $1-\alpha /2$ quantile of t distribution with n – p – 1 degrees of freedom. Another approach is called HC4, which is a modified version of HC3 that was found to be better for more general use. That is, $\text{HC}4=S={\left(X\text{'}X\right)}^{-1}{X}^{\prime }diag\left[r_i^2/{\left(1-h_{ij}\right)}^{d_{ii}}\right]X{\left(X^{\prime }X\right)}^{-1}$.

HC-Robust Wild Bootstrap (W-B) Approach

The W-B approach was originally developed by Wu (1986), and it can be used to produce unbiased estimates of regression models with heteroscedastic errors. The W-B approach resamples the multiway, clustered heteroscedastic error terms to estimate the bootstrap dependent variable scores for constructing the CI surrounding the slope parameters. Roodman et al. (2019) have advanced and modified Wu’s (1986) approach by developing a fast W-B approach that can efficiently calculate bootstrap test statistics and implements a HC-robust W-B procedure for constructing the CI via their developed R package (fwildclusterboot; Fischer et al., 2023), and the performance of this approach is evaluated in this study.

Least Median of Squares (LMS) Estimator

The LMS estimator was developed by Rousseeuw (1984). Unlike OLS regression, using the sum of the squared residuals, the LMS-estimator uses the median of the squared residuals. This can be expressed as

where M is the median. The LMS-estimator is the first method that achieves the breakdown point 0.5; therefore, it is resistant to outliers. The LMS function in R does not provide an analytic method for the standard error, but it can be estimated through bootstrapping that locates the $\alpha /2$ and $1-\alpha /2$ percentile bootstrap slopes based on the LMS-estimator for the $1-\alpha$ CI (called LMS-B in this study). However, it has a critical limitation: the relative efficiency of the LMS-estimator to OLS regression is 0 due to $n^{-1/3}$ convergence. For this reason, the LMS-estimator is not practically useful, but it plays a significant role in other robust methods such as the MM-estimator, which is described below (Andersen, 2008).

Least Trimmed Squares (LTS) Estimator

Another robust estimator developed by Rousseeuw (1984) is called LTS, which is defined as

where $r_{\left(1\right)}^2\le \dots \le r_{\left(h\right)}^2$ are the ordered squared residuals following ascending order. The breakdown point of 0.5 can be achieved with $h=\left[n/2\right]+\left[(p+1)/2\right]$, which is commonly used (Wilcox, 2022). Although the LTS-estimator may have a high breakdown point depending on h, its efficiency is very low, about 8% (Stromberg et al., 2000). Nevertheless, this method has some value insofar as it is used as an initial estimate for other robust methods (Andersen, 2008). Another related approach is the use of bootstrapping that locates the $\alpha /2$ and $1-\alpha /2$ percentile bootstrap slopes based on the LTS-estimator for the $1-\alpha$ CI, which is labelled as LTS-B.

Maximum Likelihood Type Estimation (M-estimator)

The M-estimator proposed by Huber (1973) minimizes

where $\psi$ is a robust loss function with a unique minimum at zero. The robustness of the M-estimator depends on a robust loss function that researchers choose. One commonly used function is Huber’s p function, expressed as

where c is a tuning constant which can be adjusted to control asymptotic efficiency. When outliers lie in the y-axis, the Huber M-estimator, in general, is more efficient than OLS regression against outliers. However, it does not consider a leverage point. If there is an outlier in the x-axis, the Huber’s M-estimator is not better than OLS regression; therefore, the breakdown point is $1/n$(Wilcox, 2022).

Generalized M-Estimator (GM-Estimator)

Due to the limitation of the M-estimator which does not consider leverage points, a generalized M-estimator was developed to guard against leverage points by adding some weight, ${\omega }_i$, to $x_i$ values. Mallow (1973, as cited in Krasker & Welsch, 1982, p. 596) proposed the GM-estimator, which can be expressed as

where $x_{io}=1$, and $j=0,\dots ,\mathrm{ }p$. Mallow used ${\omega }_i= \sqrt{1-h_{ii}}$ based on a condition if $h_{ii}>h_{jj}, {\omega }_i<{\omega }_j$. That is, high leverage points to $x_i$ receive less weight than low leverage points to $x_i$. As a result, this method gives less weight to good leverage points resulting in a loss of efficiency (Andersen, 2008).

To solve the low efficiency issue by using Mallow’s weight, Schweppe proposed a different solution (Handschin et al., 1975), expressed as

where $j=0,\dots ,p$. The idea of Equation 8 is to use different weight values according to the size of the residuals. In other words, Schweppe tried to solve the limitation of Mallow’s weight by dividing $r_i$ by ${\omega }_i$. Even though Schweppe’s estimator may provide a better option for dealing with leverage points than the regular M-estimator, its break point is less than 0.5 (Maronna et al., 2006), and especially low with a large number of predictors (Andersen, 2008).

Schweppe One-Step (S1S) Estimator

Coakley and Hettmansperger (1993) expanded from the Schweppe’s estimator and developed the S1S estimator, expressed as

where ${\omega }_i$ is determined by using the same criterion that the original Schweppe’s estimator used. The S1S-estimator is different from the two GM-estimators mentioned above in that it can achieve 95% efficiency when the error term is normally distributed. Moreover, it achieves a breakdown point of 0.5 by using the LTS-estimator as an initial estimator (Andersen, 2008). Wilcox (2022) states that the S1S-estimator can be effective when the sample size is large, and $\mathrm{\epsilon }$ follows a normal distribution. However, it becomes inefficient when samples sizes get smaller.

MM-Estimator

Another popular robust technique derived by Yohai (1987) is MM-estimator, so called because it calculates the final estimates by employing more than one M-estimation. Three steps are involved to find the MM-estimator. The first step is to compute initial estimates of the coefficients $\widehat{\beta }$ with high breakdown points, 0.5, using s-estimation. Then, a robust M-estimate of scale $\widehat{\sigma }$ of the residuals is calculated by using the initial estimation (Maronna et al., 2006). The robust scale of $\widehat{\sigma }$ satisfies

The final step is to compute the regression parameters by solving the following equation (Wilcox, 2022):

where $\mathit{j}=0,\dots ,\mathit{p}$, and $\psi$, is Tukey’s biweight which is commonly used. Tukey’s biweight is expressed as

where $\widehat{\sigma }$ in Equation 10 is a robust M-estimate of scale, and Tukey’s biweight is used as a redescending function in Equation 10.

When c is 4.685, the relative efficiency of the MM-estimator is 95% to OLS regression. Under normality, it has a high breakdown point, 0.5, and relative efficiency, 95%, to OLS regression (Wilcox, 2022). Another related approach is the use of bootstrapping that locates the $\alpha /2$ and $1-\alpha /2$ percentile bootstrap slopes based on MM-estimator for the $1-\alpha$ CI; this method is called MM-B in this study.

S-Estimators

Rousseeuw and Yohai (1984) proposed the S-estimators that estimate slope and intercept values with the goal to minimize some measure of scale corresponding to the residuals. Indeed, the conventional, least squares approach is regarded as one type of S-estimator that minimizes the variance of the residuals. In this case, replacing the variance with some measure of location that is robust to outliners is the goal of using S-estimators for estimating robust slopes and intercepts. Wilcox (2022) mentioned that S-estimators may have some practical value, but no study has examined their empirical performance via simulation studies. One approach is the Nelder-Mead method (SNM; e.g., Olsson & Nelson, 1975). According to Wilcox (2022); let $R_i=y_i-b_i{x}_{1i}-\dots b_p{x}_{ip}$, and the SNM approach searches the values of $b_i,\dots ,b_p$ such that the standard error estimate (S) is minimized through some measure of scale based on the values of $R_i\dots ,R_n$. Consequently, the intercept is $b_0=M_y-b_1{M}_1-\dots b_p{M}_p$, where $M_y$ and $M_j$ are the medians of the y values and of the $x_{ij}$ scores, where i = 1,…,n.

E-Type Skipped Estimator

The purpose of using an E-Type (or error-type) skipped estimator is to remove or decrease the influence of any outliers existing in a dataset on fitting a regression model. This estimator often begins by running preliminary fit to search for any outlier residuals, removing or downweighing those outliers, and fitting a regression model based on the remaining data. Wilcox (2022) supposes that $M_r$ is the median of the residuals, and MAD_r (median absolute deviation) is the median of the values $\left|r_i-M_r\right|,\dots ,\left|r_n-M_r\right|$. Consequently, for an ith point ( $x_i,y_i$) with $\left|r_i-M_r\right|>2\left(MA{D}_r\right)/.6745$, this point is deemed an outlier. The slope and intercept estimates will then be based on the data points that are not declared as outliers.

Methods Based on Robust Covariances

Replacing conventional covariances with robust covariances in fitting a regression model is a general approach that leads to robust intercept and slope parameter estimators of the model (ROB; Huber, 1981). In its simplest form with one predictor, the slope of the OLS regression line is ${\beta }_1={\sigma }_{xy}/{\sigma }_x^2$, in which the numerator can be replaced by a robust covariance estimate between x and y, and the denominator can be replaced by a robust variance estimate of x. One common approach for estimating robust variance and covariance is the biweight midcovariance that estimates the variability and co-variability of the scores based on the robust location (medians) and robust distance (MAD) that exist in a dataset.

Quantile (QUA) Regression

Another robust approach is estimating regression parameters based on minimizing the summation of the absolute of the residual scores, $\sum \left|r_i\right|$. According to Koenker and Bassett (1978), researchers could estimate the qth quantile of y scores given x scores. Suppose ${\rho }_q\left(u\right)=u\left(q-I_{u<0}\right)$, where I is the indicator function. Hence, the regression model is estimated by minimizing $\sum {\rho }_q\left(r_i\right)$. When q = 0.5 (or 50th quantile), it refers to the least absolute value of the estimator which, in turn, leads to an estimate of the median of y scores, a robust measure of location, given x scores.

In sum, it has been suggested that robust regression methods outperform OLS regression when outliers exist in the data (Andersen, 2008; Anderson & Schumacker, 2003; Brossart et al., 2011; Finger, 2010; Maronna et al., 2006; Mercer et al., 2015; Sauvageau & Kumral, 2015; Wilcox, 2022; Wilcox & Keselman, 2004; Yellowlees et al., 2016). However, no robust regression technique is universally superior, because each regression method has strengths and limitations. Depending on the situation, one may be more appropriate than another. In terms of handling leverage points, the S1S-estimator may be the best option. However, the S1S-estimator becomes less efficient with a small sample size, in which case the MM-estimator may be more appropriate (Andersen, 2008). In general, MM-type robust estimation may be the preferred choice in terms of relative efficiency, bias, and testing the null hypothesis (Anderson & Schumacker, 2003). When outliers are located in the y-axis, the Huber M-estimator would be appropriate as well, because it has a .5 breakdown point regarding outliers in the y-axis (Yellowlees et al., 2016). Therefore, the purpose of this research study is to compare robust regression methods using different settings to provide other researchers with information that will help them select the most appropriate robust regression methods when errors violate normality and homoscedasticity assumptions in psychological research.

Type I Error

When errors were normal (Table 20), the Type 1 error rates were well protected for the OLS method, range = .044, .062, mean = median = .051. Some robust methods (e.g., HC3, HC4, W-B, MM-B, and ROB) produced Type 1 error rates similar to those obtained by OLS. Of the remaining robust methods, some resulted in slightly smaller (or more conservative) Type 1 error rates, e.g., LMS-B: range (.000, .007), mean = median = .003; LTS-B: range (.004, .036), mean = .015, median = .012; GM: range (.031, .060), mean = median = .043; S1S: range (.013, .064), mean = .032, median = .003; S: range (.010, .033), mean = median = .021; E: range (.007, .039), mean = .021, median = .020; QUA: range (.018, .042), mean = .028, median = .026, while others produced higher Type 1 error rates, e.g., lsfitci range (.056, .095), mean = .074, median = .073; LTS: range (.115, .163), mean = .146, median = .148; M: range (.036, .082), mean = .059, median = .061; MM: range (.,035, .080), mean = .058, median = .059.

When errors were mixed-normal, SH, AL, and AH (Tables 21–24), the Type I error rates produced by OLS were surprisingly comparable with those observed in the normal distribution: range (.038, .065), mean = .051, median = .050. Four robust methods, W-B, range (.031, .061), mean = .047, median = .047; GM, range (.029, .059), mean = .041, median = .040; MM-B, range (.032, .065), mean = .047, median = .046; and ROB, range (.033, .067), mean = .046, median = .045, yielded reasonable Type 1 error rates close to the nominal level of .05. Nine robust methods, HC3, HC4, LMS-B, LTS, LTS-B, S1S, S, E, and QUA, produced Type I error rates slightly smaller than the nominal level: range (.000, .066), mean = median = .023. Three robust methods, LTS, N, and MM, yielded Type I error rates higher than the nominal level of .05, range (.031, .135), mean = .072, median = .062.

When errors were heteroscedastic (Table 25), OLS performed poorly with less protected (or higher) Type I error rates: range (.361, .409), mean = median = .378. Indeed, only four robust methods, HC3, range (.043, .065), mean = .056, median = .060; HC4, range (.041, .057), mean = .051, median = .054; W-B, range (.047, .068), mean = .055, median = .053; and GM, range (.037, .058), mean = .048, median = .047, produced reasonable Type I error rates. Eight robust methods, LMS-B, LTS-B, S1S, MM-B, S, E, ROB and QUA, yielded smaller Type I error rates: range (.004, .059), mean = .027, median = .028. The remaining three robust methods, lsfitci, LTB, M, and MM resulted in larger Type I error rates: range (.072, .606), mean = .390, median = .505.

Power

When errors were normal (Table 26), the power rates produced by OLS were the largest: range (.566, 1), mean = .882, median = .949. Of the remaining robust methods, some of them, including lsfitci, HC3, HC4, W-B, LTS, M, and MM, yielded power rates similar to those obtained by OLS: range = (.512, 1), mean = .871, median = .935. GM, MM-B, and ROB produced power rates slightly smaller than OLS: range (.403, 1), mean = .816, median = .878. The rest of the robust methods, LMS-B, LST-B, S1S, S, E, and QUA, produced much smaller power rates: range (.019, 1), mean = .585, median = .584. The factors that influenced the power rates were the same for all the approaches. First, when n increases, power rates increase. Second, when ${\beta }_1$ increases, power rates increase. Third, there is no obvious relationship between the size of ${\beta }_2$ and power rates.

When errors were mixed normal, SH, AL, and AH (Table 27–30), the OLS produced the smallest power rates, range (.125, 1), mean = .47, median = .329. This pattern of results is reasonable because the standard error should be more biased with the longer tails, and hence, this affects the precision of the estimates and likelihood of detecting significant results. Of the remaining methods, most of them produced power rates larger than those obtained by OLS. The LTS and MM led to the relatively larger power rates in the mean range of .80, i.e., range (.402, 1), mean = .829, median = .874. The M, GM, MM-B, and ROB produced power rates in the mean range of .70, i.e., range (.262, 1), mean = .745, median = .789. The S1S, S, E, and QUA fell in the mean range of .60, i.e., range (.153, 1), mean = .658, median = .654. The remaining robust methods, lsfitci, HC3, HC4, W-B, LMS-B, LTS-B, resulted in power rates that fell in the mean range of .4, range (.024, 1), mean = .461, median = .399.

When errors were heteroscedastic (Table 31), the OLS produced power rates ranging from .458 to .865, with a mean (or median) of .634 (or .614), although one should be cautious about the use of the OLS because of the highly unprotected Type I error rates. Of the remaining robust methods, LTS, M, S1S resulted in the largest power rates, range (.301, 1), mean = .859, median = .928. The LTS-B, GM, S, E, and QUA resulted in power rates like those observed in OLS, range (.240, .995), mean = .629, median = .607. The other approaches, lsfitci, HC3, HC4, W-B, LMS-B, MM-B, and ROB produced power rates smaller than those obtained by OLS, range (.134, .973), mean = .402, median = .344.

Coverage Probability and Width of CI

When errors were normal (Table 32), the coverage probabilities yielded by the OLS method were desirable: range (.928, .962), mean = median = .946. Of the 36 conditions, 36 (or 100%) produced coverage probabilities within the criteria of (.925, .975). The widths of the CI ranged from .278 to .654 with mean = .460 and median = .457 (Table 33). Of the remaining robust methods, eight of them, HC3, HC4, W-B, M, GM, MM, MM-B, and ROB, led to coverage probabilities comparable to the OLS approach, range (.918, .972), mean = .950, median = .951. Of the 36 conditions, 36 (or 100%) yielded by HC3, HC4, W-B, GM, MM, MM-B, and ROB fell within (.925, .975), whereas 35 (or 97.2%) produced by M fell within (.925, .975). The widths were slightly wider than those obtained by OLS, range (.282, .837), mean = .513, median = .516, and this pattern is reasonable as the OLS-based CI is expected to be the most precise with normal errors. The remaining seven robust methods, lsfitci, LMS-B, LTS-B, S1S, S, E, and QUA, were subpar. LMS-B, LTS-B, S1S, S, E, and QUA produced coverage probabilities larger than expected, range (.936, 1), mean = .981, median = .982. Of the 36 conditions, only around 11.5 conditions (or 31.9%), on average, fell within the criteria of (.925, .975). The widths were wider than those obtained by OLS, range (.322, 2.294), mean = .963, median = .883. On the contrary, lsfitci and LTS resulted in smaller coverage probabilities, range (.814, .932), mean = .879, median = .886. Of the 36 conditions, none of the LTS’s (or one of the lsfitci’s) conditions fell within the criteria of (.925, .975). The widths were narrower than those observed in OLS, range (.253, .612), mean = .425, median = .420.

When errors were mixed-normal, SH, AL, and AH (Table 34 to 37), the OLS method still produced good coverage probabilities: range (.925, .965), mean = median = .946. Of the 144 conditions, 144 (or 100%) yielded coverage probabilities within (.925, .975). As shown in Table 38–41, the widths ranged from .333 to 3.752 with mean = 1.514 and median = 1.550. Five robust methods, M, GM, MM, MM-B, and ROB, yielded desirable coverage probabilities: range (.919, .981), mean = .951, median = .952. Of the 144 conditions, 142.4 (or 98.9%), on average, led to coverage probabilities within (.925, .975). The widths were noticeably narrower than those found in OLS, range (.281, 1.225), mean = .625, median = .581. The W-B method yielded desirable coverage probabilities, range (.936, .969), mean = .952, median = .953, but its widths were noticeably wider, range = (.337, 3.610), mean = 1.401, 1.314. Eight robust methods, HC3, HC4, LMS-B, LTS-B, S1S, S, E and QUA, produced larger coverage probabilities, range (.934, 1), mean = median = .977. Of the 144 conditions 64.1 (or 44.5%) resulted in coverage probabilities within (.925, .975). The widths ranged from .289 to 3.585 with mean = .625, median = .581. Two robust methods (lsfitci and LTS) produced smaller coverage probabilities, range (.838, .943), mean = median = .895. Of the 144 conditions, only 6 (or 4.2%) yielded coverage probabilities within (.925, .975). The widths ranged from .304 to 2.960 with mean = 1.277 and median = 1.298.

When errors were heteroscedastic (Table 42), the OLS produced undesirable coverage probabilities: range = (.584, .654), mean = .615, median = .612. Of the 36 conditions, 0 fell within (.925, .975). As shown in Table 43, the widths ranged from .464 to .989 with mean = .728 median = .735. Six robust methods, HC3, HC4, W-B, GM, E, and QUA, led to desirable coverage probabilities, range (.930, .979), mean = .951, median = .950. Of the 36 conditions, 35.7 (or 99.1%), on average, fell within (.925, .975). It is noteworthy that GM, range (.937, .967), mean = median = .950; HC3, range (.930, .960), mean = .944, median = .943; HC4, range (.937, .965), mean = .949, median = .948; and W-B, range (.932, .960), mean = .947, median = .949, resulted in coverage probabilities much closer to the true value of .950. The widths ranged from .421 to 2.366 with mean = 1.219 and median = 1.126. Six robust methods, LMS-B, LTS-B, S1S, MM-B, S, and ROB, led to over-coverage probabilities, range (.954, .997), mean = median = .977. Of the 36 conditions, 16 (or 44.4%), on average, fell within (.925, .975). The widths ranged from .399 to 1.449 with mean = .864 and median = .815. Three robust methods (lsfitci, LTS, and M) resulted in under-coverage probabilities, range (.394, .905), mean = .602, median = .489. Of the 36 conditions, 0 fell within (.925, .975). The widths ranged from .142 to 1.586 with mean = .576 and median = .299.