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

Ordinary least squares (OLS) regression is widely employed for statistical prediction and theoretical explanation in psychology studies. However, OLS regression has a critical drawback: it becomes less accurate in the presence of outliers and non-random error distribution. Several robust regression methods have been proposed as alternatives. However, each robust regression has its own strengths and limitations. Consequently, researchers are often at a loss as to which robust regression method to use for their studies. This study uses a Monte Carlo experiment to compare different types of robust regression methods with OLS regression based on relative efficiency (RE), bias, root mean squared error (RMSE), Type 1 error, power, coverage probability of the 95% confidence intervals (CIs), and the width of the CIs. The results show that, with sufficient samples per predictor (n = 100), the robust regression methods are as efficient as OLS regression. When errors follow non-normal distributions, i.e., mixed-normal, symmetric and heavy-tailed (SH), asymmetric and relatively light-tailed (AL), asymmetric and heavy-tailed (AH), and heteroscedastic, the robust method (GM-estimation) seems to consistently outperform OLS regression.

is widely used in many areas of research in fields such as biology, business, education, computer science, psychology, and more (Anderson & Schumacker, 2003;Haupt et al., 2014;Sauvageau & Kumral, 2015;Yellowlees et al., 2016).OLS regression is used to predict dependent variables based on explanatory variables plus errors, which can be presented as where y i denotes the dependent or response variable from i = 1,…, n observations, x ij = x i1 , …, x ip (where j = 1,…, p) denote p numbers of predictors or explanatory variables, β 0 , …, β p are the p + 1 regression coefficients, and ε i represents the difference between the actual observed score and the score predicted from a statistical model.The purpose of OLS regression is to minimize the sum of squares of the difference between predicted values and actual observed values.That can be written as With the OLS regression model, researchers in psychology can examine whether predic tors can significantly predict an outcome measure, which is a highly appealing and widely-employed statistical procedure.For example, in one study (McCrone et al., 2005), OLS regression was used to predict costs based on family history of psychiatric illness.
In another study, Mossakowski (2011) used OLS regression to analyze if unfulfilled expectations predict subsequent symptoms of depression.
To use OLS regression properly, researchers need to check whether or not four critical assumptions have been met (Anderson & Schumacker, 2003;Field & Wilcox, 2017;Greene, 2003;Yellowlees et al., 2016).First, errors must follow a normal (N) distribution.Second, explanatory variables should not be correlated with error distribution.Third, homoscedasticity must be achieved; that is, residuals at each level of the predictor variable should have a common variance (Anderson & Schumacker, 2003;Greene, 2003;Yellowlees et al., 2016).Finally, the relationship between the response variable and the explanatory variable should be linear (Field & Wilcox, 2017).When all of these assumptions are met, OLS regression will be the maximal, unbiased linear estimator of the regression coefficients in the population (Field & Wilcox, 2017).However, in reality, this is not what researchers commonly face (Erceg-Hurn & Mirosevich, 2008).
In many research scenarios involving behavioral and social data, these assumptions are violated.One of the assumptions that is often violated is the assumption of homo scedasticity.Erceg-Hurn and Mirosevich (2008) claimed that the presence of heterosce dasticity (HE) is common in real data.In addition to that, normality assumptions are rarely met in practice.Micceri (1989) found that, of 440 large-sample measures related to psychology such as achievement and other common psychometric measures, none were normally distributed among the data they investigated.Approximately 15.2% of 440 distributions followed Gaussian distribution, and most of the distributions followed either a heavy or skewed distribution.
Using OLS regression under violations of assumptions gives rise to several critical problems for researchers.First, when the normality assumption is not met, OLS regres sion produces lower power and wider confidence intervals.In other words, it weakens the generalizability of data.Moreover, it may inflate the possibility of making a Type 1 error (Anderson & Schumacker, 2003;Erceg-Hurn & Mirosevich, 2008).Even when the normality assumption is not violated, the problems mentioned above can occur with the presence of heteroscedasticity (Brossart et al., 2011).As a result, OLS regression becomes inefficient and generates unstable results when underlying assumptions are not met, which is quite common in practice.
Just as violated assumptions are commonly encountered in reality, so are outliers.As can be seen in Equation 2, each observation is weighted equally by the OLS estimator, and hence, a large difference between an outlier and the mean of all other scores will have a substantial impact in the accuracy of the slope estimates.For this reason, the presence of outliers in data also markedly distorts the efficiency of OLS regression, resulting in erroneous results.Outliers can arise from different sources: man-made or random (Osborne & Overbay, 2004).That is, outliers can be caused by researchers during data entry, incorrect distribution assumptions, and sampling error, or by random chance when collecting samples from a population.
Outliers have different influences on the estimation of regression coefficients, de pending on their location.The influence of an outlier may be much more severe when it lies on the x-axis, called the leverage point, than when it lies on the y-axis (Anderson & Schumacker, 2003).The leverage point can be either good or bad (Rousseeuw & Leroy, 2003).A good leverage point, which is away from the bulk of the points but close to a regression line, reduces the standard error.However, when the location of a data point is far away from the rest of the data points and from the line of best fit, called bad leverage, it can pull the regression line towards the outlier's location.Therefore, outliers in data not only indicate the assumptions that normality and homogeneity may not hold, but also seriously impact the result based on OLS regression including inaccurate intervals, lowering statistical power, and Type 1 and Type 2 errors (Brossart et al., 2011).
One of the common approaches to handling these outliers is to transform data using logarithms or square roots (Grissom, 2000).Transformation may be a valuable option; however, it often gives rise to more problems.Transforming data often fails to restore normality and homoscedasticity and is not adequate to deal with outliers (Wilcox, 2022).Moreover, it also leads researchers to interpret data in an inaccurate way by changing the original construct (Osborne, 2003).
The other option that researchers commonly use is to discard outliers.Yellowlees et al. (2016) argued that this method, however, can only be used if an outlier can be traced back to an error that occurred during an experiment (e.g., wrong dose of a drug or product).Finney (2009, p. 51) also stated: "the danger [is] that scientists bias their conclusions by removing data that deviate markedly from current ideas of truth." They warned: "[n]ever discard an apparent outlier unless there is strong evidence that it was the product of a measurement or other form of observation that suffered a gross mistake or accident, this misfortune being unrelated to any experimental treatment under investi gation." Hence, this method may only be advisable when researchers clearly know the cause of outliers (Yellowlees et al., 2016).Unfortunately, it is hard for researchers to objectively define outliers especially when the data have many explanatory variables (Maronna et al., 2006).Moreover, the presence of outliers sometimes disguises other outliers (called a masking effect; Wilcox, 2022).On the other hand, retaining outliers could still result in a highly misleading result, if the goal of a research study is to characterize the nature of the association among the bulk of participants.
Robust regression can be an alternative method to deal with outliers and assumption violations.Despite the potential of robust regression, there is no single study that com prehensively integrates these methods and systematically evaluates their accuracy in a Monte Carlo experiment.The current research fills this research gap.The following section discusses the mathematical and computational details of these methods.

Robust Regression
Robust regression is a modern method conceptualized many decades ago.In the 1950s, Siegel (1956) stated that non-parametric and robust techniques of hypothesis testing are best suited to behavioral sciences data.However, robust regression has only recently been studied due to advancements in computer technology (Anderson & Schumacker, 2003).Unlike OLS regression, which gives weights to outliers, robust regression reduces the impact of the outliers by weighing them down.This allows researchers to take outli ers into account in the statistical model rather than using other, potentially problematic methods to deal with them.
Two pivotal concepts need to be addressed to understand robust regression methods: breakdown point, and relative efficiency (Anderson & Schumacker, 2003).The break down point measures the minimum proportion of points that are needed to make a statistical estimate, such as regression slope, arbitrarily large or small.The breakdown points vary between 0, or 1/n and 50%, or n/2.In other words, an estimator with a 0% breakdown point does not efficiently prevent the regression equation from being influ enced by regression outliers or bad leverage points.OLS regression has the breakdown point of 0%; consequently, the presence of one outlier or bad leverage point can render the data inefficient.On the other hand, robust regression estimators with a 50% break down point can contain as many as 50% bad leverage points without making a statistical estimate being arbitrarily large or small.It is noteworthy that even though replacing the OLS regression with robust estimators having a high breakdown point appears to be a reasonable solution when data contains outliers or bad leverage points, some recent studies (e.g., Wilcox, 2022;Wilcox & Xu, 2023) have showed that the presence of a few outliers or bad leverage points could still have a noticeable impact on the slope estimates.Or, stated differently, robust estimators can only ensure that having a few outliers or bad leverage points will not lead to an arbitrarily large or small statistical estimate, but they do not guarantee that the slope estimates are not substantially influenced by outliers or bad leverage points, thereby leading to a misleading conclusion of the true association in practice.
Another key concept is relative efficiency, which refers to the extent to which ro bust regression performs like OLS regression when error distribution follows a normal pattern.The relative efficiency is determined by dividing OLS regression mean square error (MSE) by the robust regression MSE, which can be expressed as (Anderson & Schumacker, 2003) This is often expressed using a percentage value ranging between 0% and 100% or more.
For example, if one robust regression technique has 75% relative efficiency, this means that the method is 75% as efficient as OLS regression.While robust regression may suffer from a slightly lower efficiency than OLS regression when the normality and homoscedasticity assumptions are met, robust regression is expected to produce much more accurate (or robust) estimates and results if the assumption is violated, which is crucial for obtaining more accurate statistical results in psychological research.

Types of Robust Regressions
Least-Square-Fit-CI (lsfitci) Approach One approach to dealing with violation of the normal-error assumption is bootstrapping residuals (ε i ) instead of raw scores (x and y) to generate empirical distribution and standard error of residuals for the construction of 1 − α CI, where α is the type 1 error rate.By locating the α/2 and 1 − α/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 , where r i is the residual for i = 1,..,n, ℎ ij = x i X′X −1 x i ′, and , where x i is the ith row of X i .Consequently, the 1 − α CI surrounding β p is b j ± tS j , where t is the 1 − α/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, HC4

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 break down 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 α/2 and 1 − α/2 percentile bootstrap slopes based on the LMS-estimator for the 1 − α 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 1 2 ≤ … ≤ r ℎ 2 are the ordered squared residuals following ascending order.The breakdown point of 0.5 can be achieved with ℎ = n/2 + p + 1 /2 , which is common ly 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 α/2 and 1 − α/2 percentile bootstrap slopes based on the LTS-estimator for the 1 − α CI, which is labelled as LTS-B.

Maximum Likelihood Type Estimation (M-estimator)
The M-estimator proposed by Huber (1973) minimizes where ψ 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, ω 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, …, p. Mallow used ω i = 1 − ℎ ii based on a condition if ℎ ii > ℎ jj , ω i < ω 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, …, 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 ω 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 devel oped the S1S estimator, expressed as where ω i is determined by using the same criterion that the original Schweppe's estima tor 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 ε 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 β with high breakdown points, 0.5, using s-estimation.
Then, a robust M-estimate of scale σ of the residuals is calculated by using the initial estimation (Maronna et al., 2006).The robust scale of σ satisfies The final step is to compute the regression parameters by solving the following equation (Wilcox, 2022): where j = 0, …, p, and ψ, is Tukey's biweight which is commonly used.Tukey's biweight is expressed as where σ in Equation 10is 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 α/2 and 1 − α/2 percentile bootstrap slopes based on MM-estimator for the 1 − α 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 − …b p x ip , and the SNM approach searches the values of b i , …, b p such that the standard error estimate (S) is minimized through some measure of scale based on the values of R i …, R n .Consequently, the intercept is b 0 = M y − b 1 M 1 − …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 r i − M r , …, r n − M r .Consequent ly, for an ith point (x i , y i ) with r i − M r > 2 MAD r /.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 β 1 = σ xy /σ 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, ∑ r i .According to Koenker and Bassett (1978), researchers could estimate the qth quantile of y scores given x scores.Suppose ρ q u = u q − I u < 0 , where I is the indicator function.Hence, the regression model is estimated by minimizing ∑ ρ q r i .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 regres sion 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.

Method: A Monte Carlo Simulation
To provide a better understanding of the selection of appropriate robust regression methods, we used a Monte Carlo simulation study to compare OLS and robust regression methods under a variety of conditions that researchers commonly face.The simulation was conducted based on multiple regression models with two independent predictors that can be expressed as To represent varied research conditions, we included variations across the following research variables: sample size per predictor, slope, and different types of error distribu tion.

Sample Size Per Predictor
Several rules have been proposed for minimum required sample ratio per predictor to conduct multiple regression analysis (Miller & Kunce, 1973;Schmidt, 1971).Schmidt (1971) recommends 15 to 20 samples per predictor, whereas Miller and Kunce (1973) argue that there should be 30 samples per predictor for accurate regression analysis.In response, and to provide more representative results, we added two more sample size variables, thus including 20:1, 30:1, 50:1, and 100:1 (n; sample size per predictor) in this study.

Slope
For slope, we selected values of 0, .3606,and .5099(Cohen, 1988), which correspond to the zero effect, medium effect (.3606 2 = 13% of the variance of y explained by x),

Error Distribution
For error distribution, we adopted the same criteria used by Yuan and MacKinnon (2014) and Wilcox (2022).First, a normal distribution, N ~ (0, 1 2 ), was simulated.Second, follow ing Yuan and Mackinon, a mixed normal distribution (MN) with 90% random errors was generated from N ~ (0, 1 2 ) and 10% random errors generated from N ~ (0, 10 2 ).It is noteworthy that MN is a type of symmetric distributions with heavy tails on both ends, and it has been widely employed and tested in previous simulation studies (e.g., Algina et al., 2005).Third, a heteroscedastic error term, in which the variance used to generate the random error score of each simulated participant depends on his or her x ip , N ~ (0, x ip 2 ), where i = 1,..,n, and p = 1 or 2, was simulated (Yuan & MacKinnon, 2014).The remaining three distributions followed Wilcox's (2022) method for simulating asymmetric and/or heavy tailed error distribution for testing the performance of robust regressions.That is, the fourth distribution was a symmetric and heavy-tailed (SH) distribution based on a g-and-h distribution with (g, h) = (0, 0.5).The fifth distribution was an asymmetric and light-tailed (AL) distribution with (g, h) = (0.5, 0).The sixth distribution was an asymmetric and heavy-tailed (AH) distribution with (g, h) = (0.5, 0.5).
In summary, four sample sizes per predictor, nine combinations of slopes, and six types of error distributions were evaluated.This factorial design created a total of 4 × 9 × 6 = 216 different conditions.Each of the 216 conditions were replicated 1,000 times.For bootstrapping, the simulated x and y scores were resampled 1,000 times for constructing the 95% bootstrap percentile intervals.In sum, this design produced a total of 216 × 1, 000 × 1, 000 = 216,000,000 simulated data sets for evaluation.OLS-regression estimates and CI, nine robust regression estimates and analytic-based CIs (i.e., LTS-es timator, M-estimator, GM-estimator, S1S-estimator, MM-estimator, S-estimator, E-type, ROB, and QUA regressions), one robust, analytic-based CI for OLS-regression estimates (i.e., lsfitci), and five bootstrap-based CIs (i.e., HC3, HC4, LMS-B, LTS-B, and MM-B) described above were performed on these simulated data sets in order to compare the accuracy of their results.We used the statistical software, R, to conduct our simulation (R Core Team, 2023), and the code is shown in the Supplementary Materials.

Criteria
The criteria we used to compare the regression methods were relative efficacy, bias, RMSE, Type I error, power, coverage probability of the 95% CI, and width of the CI.For relative efficacy, higher percentages are desirable.For example, if the relative efficacy of a robust approach is .98,this means that it can maintain 98% of efficiency compared to the conventional OLS estimates, when the normality and homoscedasticity assumptions are met.Regarding bias, when a regression method generates a slope farther away from the true slope, this is considered bias.Consequently, bias, in this study, is defined as the difference between the mean of the 1,000 replicated slopes minus the true slope (i.e., bias = b − β, where b is the mean of1,000 replicated observed slopes in each simulated condition, and β is the true slope value manipulated in the condition).Evaluating the bias of the slope estimates is insufficient because it does not measure the variability of those values.The root mean square error (RMSE) and variance of the 1,000 replicated slope values are also included.RMSE is defined by ∑ r = 1 1, 000 b r − β 2 /1, 000, which measures the average (squared) distance between each of the 1,000 replicated slope values with the true value.Confidence width is used to evaluate the precision and sampling error of the slope estimates.A narrower confidence width indicates a more precise estimate.For Type 1 error, we set its error rate at α = .05level (two-tailed test), which is commonly used as a criterion in psychological research.When the true slope was set at 0, we examined the number of times (or probability) that a regression method would lead to an incorrect decision (i.e., rejecting the null hypothesis: β = 0) out of 1,000 replications for each of the 216 manipulated conditions.By the same token, when the true slope was set at .3606 and .5099,we evaluated the number of times (or probability) that a regression method would lead to a correct decision (i.e., rejecting the null hypothesis: β = 0) out of 1,000 replications for each of the 216 manipulated conditions.The coverage probability examines the probability a 95% CI has spanned a true slope parameter value.Theoretically speaking, of the 1,000 replications, the number of the 95% CIs that has spanned the true parameter value is expected to be 950 (or coverage probability = 95%).
In practice, sampling error exists, and hence, an observed coverage probability ranging from .925 to .975yielded by a regression method is deemed desirable (Chan & Chan, 2004).

Relative Efficiency (RE)
The results of RE of each robust regression method compared to OLS regression, with normal error distribution, are presented in Table 1 1 (for all tables, see Supplementary Ma terials).When errors were normally distributed, three types of RE results were observed.The first type was observed by LMS and LTS, which produced the least efficient RE with a range from .761 to .778 when the sample size was small (20).The second type was observed by S, E, and QUA, and they were regarded as moderate RE, which ranged from .935 to .960 with a small sample size of 20.The third type was found in the M, GM, S1S, MM, and ROB approaches leading to the highest RE, which ranged from .958 to .996 with a small sample size of 20.Comparing the effects of different manipulated factors, the sample size was found to be the most influential.As the sample size got larger, the discrepancy between the efficiency of robust regression and OLS regression got smaller.When the sample size was larger (e.g., n = 100), most of the robust regression methods (except LMS and LTS) were nearly the same in efficiency, ranging from 0.982 to 0.999.Therefore, this result indicates that there is no obvious or substantial loss of efficiency based on the robust regression methods compared to OLS regression, especially when n is greater than 100.
Second, HC3, HC4, and W-B appear to offer good adjustment when errors are HE.In particular, the means (or medians) of the Type I error rates were .056,.051,.055(or .060,.054,and .053)for HC3, HC4, and W-B respectively.The coverage probabilities were also improved: the means (or medians) were .944, .949, and .947 (or .943, .948, and .949)for HC3, HC4, and W-B, respectively.On the other hand, HC3, HC4, and W-B still have the same limitation as in the conventional OLS-based CIs, meaning that the power rates were small, and the widths of the CIs were wide for MN, SH, AL, and AH.Hence, HC3, HC4, and W-B could potentially solve the issue of the Type I error and coverage of the true parameter for HE, but they may not be the most appropriate approach to use in practice when errors are non-normal (MN, SH, AL, AH, and HE).
Third, of the remaining robust methods, only M, GM, MM, MM-B, and ROB could be considered for MN, SH, AL, and AH, and only GM, E, and QUA could be considered for HE because of their superiority of the coverage probabilities that span the true parameter slope values.Comparatively, MM produced Type I error rates slightly larger than the criterion of .05,and it also consistently led to "no solution" for HE.Hence, it is not the most appropriate approach in practice.MM-B, a modified approach based on MM, seemed to overcome the no-solution issue with MM for HE, but its Type I error rates seem to fall in the conservative side of .05(e.g., mean = .032,median = .031for HE).ROB's performance is also similar to MM-B's performance, meaning that ROB behaved appropriately for N, MN, SH, AL, and AH, except for more conservative Type 1 error rates (mean = median = .028)for HE.Conversely, QUA seems to have good, protected Type 1 error rates for HE (mean = .047,median = .046),but they became noticeably smaller (or more conservative) for N, MN, SH, AL, and AH.E consistently produced smaller Type I error rates for all the six error conditions.
In sum, GM seems to have the best all-round performance.The Type 1 error rates remain slightly conservative for N, MN, SH, AL, and AH (means = .043,.042,.042,.040,and .041;medians = .043,.043,.041,.039,and .040,respectively) The power rates are reasonable and comparable to those obtained by M, MM, MM-B, and ROB for MN, SH, AL, and AH.When errors are HE, the Type 1 error rates are still appropriate (mean = .048,median = .047),and the power rates are desirable (mean = .625,median = .586)and comparable to those obtained by other robust methods.

Discussion
OLS regression is a widely employed statistical method in psychology (Anderson & Schumacker, 2003;Erceg-Hurn & Mirosevich, 2008;Haupt et al., 2014).However, its efficiency is often distorted in practice due to violated assumptions and the presence of outliers (Erceg-Hurn & Mirosevich, 2008;Micceri, 1989;Wilcox, 1998).When researchers have outliers in their data, robust regression may be a valuable alternative to handle out liers without causing other potential problems caused by other methods (e.g., changing the original construct or distribution).However, due to the existence of several types of robust regression estimators, which have their own strengths and weaknesses, it may be confusing and challenging for researchers to choose which robust regression method is appropriate for their research-without a clear guideline based on empirical evidence from a simulation study.

Implications of the Findings
The primary purpose of this study was to provide applied researchers with empirical evi dence and guidelines to select more appropriate regression methods by comparing OLS and robust regression under different conditions.This simulation study suggests that when the normality assumption is met, OLS regression outperforms robust regression methods in terms of bias, RMSE, Type 1 error, power, coverage probabilities and confi dence width of the CIs.This is because OLS regression achieves maximum efficiency when the normality assumption is met (Andersen, 2008;Field & Wilcox, 2017).Consis tent with previous research (Andersen, 2008;Anderson & Schumacker, 2003;Mercer et al., 2015), the current results show that when the sample size per independent variable is large (n = 100), robust regression methods are quite comparable to OLS regression.That is, as the sample size increases, the efficiency of robust regression methods also increases, thereby paralleling the OLS regression method.Therefore, with sufficiently large samples, robust regression methods may be an appropriate alternative.
These research findings concur with those of other researchers (Andersen, 2008;Anderson & Schumacker, 2003;Brossart et al., 2011;Finger, 2010;Mercer et al., 2015;Sauvageau & Kumral, 2015;Yellowlees et al., 2016), and indicate that robust regression methods, in general, are better options than OLS regression when the normality and homoscedastic assumptions are violated.More specifically, when errors were non-normal (MN, SH, AL, AH, HE), HC3 and HC4 provide good adjustment for HE, but the associated power rates are still noticeably small and the widths of the CIs are wide, leading to a subpar approach in practice.Comparatively, GM is the most all-around and appropriate method in terms of the Type I error, power, coverage probability, width of the CI, as well as the precision of the point estimates.There are alternative robust methods which researchers could consider if they know that the errors are either with long tails (i.e., MM-B and ROB for MN, SH, AL, and AH) or heteroscedastic (QUA), if they would like to have slightly higher power with reasonably protected Type 1 error rates.
In addition to reporting the slope estimates, researchers often report and interpret the associated standard error and p-value to evaluate the sampling error or precision of the estimates.According to Wilcox and Keselman (2004), heteroscedasticity does not impact the slope but the standard error, which, in turn, influences the p-value.These findings provide empirical evidence that concurs with Wilcox and Keselman's explanation of het eroscedasticity.Therefore, although all regression methods perform well in terms of bias, only GM handles the Type 1 error rate effectively with desirable coverage probability of the CI.Therefore, GM seems to be an appropriate option for researchers to deal with heteroscedasticity.

Limitation and Future Directions
Although the results of this study provide a valuable guideline for future research regarding the use of appropriate robust regression methods, researchers need to consider the limitations of this study.First, although the Monte Carlo simulation tool allowed us to compare OLS and different types of robust regression methods under a variety of conditions, the findings were based on the simulated data.Therefore, like all other Monte Carlo simulation studies, the degree to which the findings may generalize to real data is uncertain and needs further examination.On the other hand, by including a variety of sample sizes, slopes, and error distribution conditions, based on previous studies, the research findings may serve as a guideline for future researchers when they select an appropriate robust regression method for their research.Future research may explore additional manipulated factors, such as larger sample size per predictor (n > 100) and other types of non-normal distributions; it could also compare OLS and robust regression methods based on real-world data in published studies.Second, this study shows that, broadly speaking, the GM estimation seems to perform appropriately across all the simulated conditions.Indeed, other studies (e.g., Wilcox, 2022) illustrated that it is overly simplistic for researchers to suggest one single estimator that is robust to all different levels or patterns of heteroscedasticity.Additional research is needed to further examine the details regarding the performance of various robust estimators under different levels and patterns of heteroscedasticity.

Conclusion
Due to the accessibility and advancement of computing power and the existence of free statistical software R and packages, researchers are readily able to conduct research using robust regression.When researchers suspect outliers in their data, but don't know the exact source, robust regression methods may be a valuable option to consider when addressing the specific outliers in their analysis.Based on the current research findings, researchers may be able to deal with outliers in their data efficiently using robust regression methods, especially if they use the GM-estimator.With sufficiently large sample sizes (n = 100), robust regression methods could be used by default instead of OLS regression without worrying about whether the error scores meet or violate the normality assumption.On the other hand, when the normality and homoscedasticity assumptions are met, OLS is found to offer small advantage in terms of slightly improved power, more precise width of the CIs, and protected Type I error as compared with other robust methods.