Robust statistical methods and the credibility movement of psychological science

Linearity and additivity

Broadly speaking, assumptions made by OLS models can be divided into two categories. The first one includes assumptions that OLS shares with other statistical models regardless of the method used to estimate the parameters, like linearity and additivity. The second category concerns assumptions related to the error distribution and structure where, unlike alternative estimation methods, OLS fails to offer any modelling flexibility should violations occur.

OLS models assume that the true, or population, relationship between predictor(s) and the outcome is linear. As an equation, this is expressed as ${\displaystyle \begin{array}{c}\hfill Y_i={\beta }_0+{\beta }_1{X}_{1i}+{\varepsilon }_i\\ \hfill {\varepsilon }_i\sim N\left(0,{\sigma }^2\right)\end{array}}$ where Y_i is the outcome value for an individual i, β₀ is the unknown value of the outcome when all predictors are 0, and β₁ is the unknown parameter associated with the predictor X_1i. This parameter represents the change in the outcome variable associated with a unit change in the predictor (Fig. 1A), which for dummy-coded categorical predictors represents the difference in the mean level of the outcome between one category of the predictor and a reference category (Fig. 1B). The error terms for each i (ɛ_i), are unobservable but are assumed to be random variables that are normally distributed with a mean of 0 and constant variance of σ². This assumption is made explicit in the second line of the equation. The parameters (β₀ and β₁) are also unobservable but can be estimated from data observed in a sample. The resulting values are known as parameter estimates and are typically denoted with hats to remind us that they are estimates based on a sample. The resulting model for a sample is, therefore, expressed as ${\displaystyle \begin{array}{c}\hfill Y_i={\widehat{\beta }}_0+{\widehat{\beta }}_1{X}_{1i}+e_i\\ \hfill \end{array}}$ in which e_i is the observed residual (the difference between the observed value of the outcome and the value predicted by the model) for entity i (Fig. 2A). When there are several predictors in the model, their combined effect is assumed to be additive, i.e.: ${\displaystyle \begin{array}{c}\hfill Y_i={\beta }_0+{\beta }_1{X}_{1i}+{\beta }_2{X}_{2i}+...+{\beta }_n{X}_{ni}+{\varepsilon }_i\\ \hfill {\varepsilon }_i\sim N\left(0,{\sigma }^2\right).\end{array}}$

The parameter estimate associated with each predictor represents the change in the outcome variable associated with a unit change in the predictor when other predictors are held constant.

Linearity and additivity are crucial prerequisites for modelling linear relationships. This assumption is not unique to OLS models but generalises across estimation methods. Whether researchers apply OLS or opt for alternative methods like Maximum Likelihood Estimation or Bayesian Estimation instead, the primary consideration should be whether the linear and additive form of the model adequately reflects the relationship between the variables of interest. When the relationship being modelled is not, in reality, linear or additive the model is not fit for purpose and any estimation or inference based on it is invalid (Gelman & Hill, 2007). In this sense, linearity and additivity are the most important conditions for a useful model.

Spherical errors

The Gauss-Markov (GM) theorem states that under certain conditions the OLS estimator for an additive linear model will have lowest sampling variance compared to other unbiased estimators. These conditions are:

1. Errors are, on average, zero. More formally, the expected value of model errors is zero, E[ɛ_i] = 0.

2. Homoscedastic errors: the variance of errors is a constant, finite, value (V(ɛ_i) = σ² < ∞). We have already mentioned that model errors are assumed to be normally distributed with constant variance and this condition (referred to as homoscedasticity) reflects the constant variance.

3. Errors are uncorrelated. More formally, Cov(ɛ_i, ɛ_j) = 0,  i ≠ j.

The last two of these conditions are often referred to as the population having ‘spherical errors’, that is, errors are both independent and homoscedastic (have a constant variance). Non-spherical errors (where either independence or homoscedasticity is not met) do not bias the parameter estimate, however this estimate will not be optimal—it will produce a larger amount of error than if the assumption was met, or if an alternative more robust estimator was applied. Additionally, the standard error will be inaccurate because the formulas for computing standard errors assume independence and constant spread of errors across the levels of the predictors. This in turn affects the confidence intervals and biases the significance tests (Field, 2024).

Assumptions of spherical errors and normal errors (see below) are not limited to OLS models, but they can pose a greater logistical challenge for OLS compared to other estimators. Although models using e.g., maximum likelihood estimation or Bayesian estimation often require these conditions to be met, they also allow researchers to select an alternative distribution family for modelling or to define more complex error structures if spherical or normal errors cannot be assumed. For ease of expression, we refer to all the assumptions as “OLS assumptions”, however it’s worth keeping in mind that other estimators are not immune to their violations.

Normal errors

Model errors are assumed to be normally distributed with constant variance, but as we shall see this assumption has little bearing on estimating model parameters. However, when model errors are normally distributed it can be shown that the parameter estimates have a normal sampling distribution: ${\displaystyle \widehat{\beta }\sim \mathcal{N}\left(\right.\beta ,{\sigma }^2{\left({\mathbf{X}}^{\text{T}}\mathbf{X}\right)}^{-1}\left)\right.}$ in which X is an n × p matrix containing a number of rows equal to the number of observations in the sample (n), and p columns corresponding to a column of 1s (the intercept) and remaining columns represent values of each predictor variable. Aside from the observation that the sampling distribution is normal when the model errors are normal, note that the variance of the distribution of $\widehat{\beta }$ is a function of the variance of errors (σ²). The variance of errors is unobservable and must be estimated from the observed data (σ² = s²), giving us the following estimate of the variance of the parameter estimates: ${\displaystyle \widehat{\text{Var}}\left(\widehat{\beta }\right)=s^2{\left({\mathbf{X}}^{\text{T}}\mathbf{X}\right)}^{-1}.}$

The standard error of the parameter estimates is the square root of this estimate ${\displaystyle \widehat{SE}\left(\widehat{\beta }\right)=\sqrt{s^2{\left({\mathbf{X}}^{\text{T}}\mathbf{X}\right)}^{-1}}.}$

It turns out that s² follows a scaled chi-square distribution, meaning that the standard error of parameter estimates also follows this distribution. Two important statistics related to the parameter estimates depend on their standard error: (1) test statistics; and (2) confidence intervals. To test whether β is equal to a particular expected value (typically β_Expected = 0) we use a test statistic that is the ratio of the difference between the parameter estimate and the expected value to the sampling variation (signal to noise ratio): ${\displaystyle T=\frac{\widehat{\beta }-{\widehat{\beta }}_{\text{Expected}}}{\widehat{SE}\left(\widehat{\beta }\right)}.}$

As noted earlier, when errors are normally distributed parameter estimates have a normal sampling distribution. Because of the central limit theorem parameter estimates also have a normal sampling distribution when sample sizes are large (more formally, as n → ∞). In these two scenarios, the test statistic above is the result of dividing a normally distributed variable by one that has a scaled chi-square distribution (see above). Dividing a normal distribution by a scaled chi-square distribution results in a t-distribution when the null hypothesis is true, and a non-central t-distribution when it is not. Therefore, under the aforementioned conditions the test-statistic follows some form of t-distribution. This fact is used to construct confidence intervals around the parameter estimate ${\displaystyle {\widehat{\beta }}_{\text{CI boundary}}=\widehat{\beta }\pm t_{1-\left(\alpha /2\right)}\times \widehat{SE}\left(\widehat{\beta }\right).}$ where the critical t value for given alpha level is multiplied by the standard error estimate of of the parameter estimate $\widehat{\beta }$ and then subtracted or added to $\widehat{\beta }$ to create the interval.

To summarise, the assumption of “normality” refers to the requirement of OLS for a normal sampling distribution of $\widehat{\beta }$ . If the errors in the population model are normally distributed, this requirement will be met. If the errors are not normally distributed, normal sampling distribution of $\widehat{\beta }$ can only be assumed in samples that are sufficiently large to invoke the Central Limit Theorem. If the parameter estimate does not come from a normal sampling distribution, the standard errors cannot be estimated accurately and the hypothesis tests will be biased. Additionally, if the population errors are not normally distributed, the OLS estimator is not optimal. Note that under the violations of both the spherical errors and normality, the parameter estimates remain unbiased.

Absence of outliers and influential cases

Outliers are typically defined as extreme cases with a large residual (Darlington & Hayes, 2017). The term ‘residual’ is equivalent to model error when we’re talking about a sample model rather than a population model. Simply put, it is the distance between the outcome value predicted by the model and the actual outcome value observed in a sample. Because the goal of the OLS estimator is to always minimise the amount of error in the model, the estimates produced by a model containing outliers will still be optimal, however they will be biased in the direction of the outlier. Influential cases are cases that have substantial influence on the parameter estimates. Influential cases can often be less obvious than outliers, especially if they have enough influence to bias the estimate to the point where their residual is minimised (Fig. 2A) and would not be detected as an outlier using conventional methods.

Are violations of assumptions common in practice?

This lackluster image of statistical practice could be justifiable if violations of assumptions were not something occurring in practice, or if the OLS estimation method remained robust despite violated assumptions. Unfortunately, the evidence suggests to the contrary on both fronts. Initial studies of distributional characteristics focused on a set of educational achievement scores, demonstrating that high skewness and kurtosis are a common feature of these variables (Lord, 1955; Cook, 1959; Burt, 1963). In a landmark early meta-research, Micceri (1989) followed-up this work by studying the distributions of 440 educational and psychometric measures. Micceri (1989) found that almost all measures were significantly non-normal due to heavy tails, exponential levels of skewness, or multi-modality.

In the decades that followed, non-normal distributions have been recorded across different areas of psychological research. Skewness is now a well-documented characteristic in various time-dependent outcomes like reaction times (Kranzler, 1992; Mewhort, Braun & Heathcote, 1992; Juhel, 1993; Reber, Alvarez & Squire, 1997; Leclaire, Osmon & Driscoll, 2020), visual search and recognition tasks (Hockley, 1984; Miyaoka, Iwamori & Miyaoka, 2018), or eye-tracking responses (Unsworth et al., 2011), with exacerbated distribution tail in clinical populations (Leth-Steensen, Elbaz & Douglas, 2000; Waschbusch, Sparkes & Northern Partners in Action for Child and Youth Services, 2003; Reckess et al., 2014). A mixture of exponential and normal distribution is often reported as the most fitting distributional family for these measures (Hockley, 1984; Juhel, 1993; Reber, Alvarez & Squire, 1997; Leclaire, Osmon & Driscoll, 2020). Additionally, Haslbeck, Ryan & Dablander (2023) found that multi-modality and high levels of skewness are common in emotion time-series, while Ho & Yu (2015) confirmed the presence of non-normality in educational scale scores in a conceptual replication of Micceri’s work. Focusing on broader psychological measures, Blanca et al. (2013) found that skewness and kurtosis are common, although they are not as extreme as originally reported by Micceri (1989). Cain, Zhang & Yuan (2017) extends this work to multivariate designs. These studies have a key limitation, in that they focus on the distribution of raw variables. As noted above, the assumption of normality refers to the distribution of the model errors, not the variables included in the model. However, when studying the residuals extracted from OLS models, Sladekova & Field (2024c) found similar levels of skewness and kurtosis as those reported by Blanca et al. (2013). In scenarios where distributions are non-normal, the distributional families researchers should expect to find most often are gamma, negative binomial, multinomial, binomial, log-normal, exponential, and Poisson, respectively (Bono et al., 2017). The only case where the normal distribution is a reasonable expectation seem to be the within-person T-scores of cognitive measures in non-clinical populations, which tend to be symmetrical with only some measures showing high kurtosis (Buchholz et al., 2024). Outside of this context, assuming normal distribution is, at best, wishful thinking. It could be argued that the non-normal distributions are not really a cause for concern—as we outlined above, we can often assume normal sampling distributions due to the central limit theorem, as long as the sample size is large enough. While some psychological fields report average sample sizes in studies to be over 100 (Fraley & Vazire, 2014; Reardon et al., 2019; Sassenberg & Ditrich, 2019; Fraley et al., 2022), others lag behind (Holmes, 1983; Marszalek et al., 2011; Hussey, 2023), with samples sizes in some areas being as low as 10 (Schrimp et al., 2022). We elaborate this discussion below with reference to statistical power.

Likewise, non-constant variance is common in published research. Grissom (2000) reviews a range of clinical outcomes and illustrates how heteroscedastic variance can be a direct by-product of study design in psychological research. In a review of published reports, Wilcox (1987) identified group variance ratios as large as 16—that is, the largest variance among the groups being compared was 16 times that of the smallest variance. In a more recent research investigating heterogeneity of variance in a range of factorial designs, Ruscio & Roche (2012) found that variance ratio at the 50th percentile of the samples was 2.74, going up to 5.10 in the 90th percentile for all samples, and up to 9.43 for designs with at least four groups. Reviewing both published and unpublished work, Sladekova & Field (2024c) reports similar findings for residual variances while also noting heteroscedasticity in designs with continuous predictors.

Are OLS models robust to violated assumptions?

Researchers therefore routinely neglect assumption checks even though they should realistically expect their violations in practice. This discrepancy is not entirely surprising—if our goal is to answer the question “Are OLS-based models robust to violations of assumptions?” with a simple yes or no, the messaging scattered in the past five decades of the robustness literature can appear mixed. If, however, we consider the effects violations of different assumptions on specific metrics of robustness, the general message is more consistent.

Perhaps the most influential work on the effects of violated assumptions on the ANOVA F-test is a review by Glass, Peckham & Sanders (1972). In the rare instances where researchers report violated assumptions, Glass et al.’s work is often cited in support of researchers’ decision not to take any remedial action. One of the conclusions presented in the review was the following: the rate of false positives of the F-test is unaffected by skewness. This conclusion has been consistently supported in further reviews and simulation studies—as long as skewness is the only concern, false positives will remain close to nominal 5% error rate (Harwell et al., 1992; Schminder et al., 2010; Liu, 2015; Blanca et al., 2017; Yang, Tu & Chen, 2019; Delacre et al., 2019). False positives will, however, be inflated in the presence of heteroscedasticity (Glass, Peckham & Sanders, 1972; Rogan & Keselman, 1977; Harwell et al., 1992; Hsiung & Olejnik, 1996; Moder, 2010; Blanca et al., 2018; Nguyen et al., 2019; Yang, Tu & Chen, 2019) and the effects on the analysis can compound when heteroscedasticity is combined with non-normal errors (Delacre et al., 2019; Sladekova & Field, 2024a). In factorial designs, the detrimental effects on false positives are further exacerbated when the group sample sizes are unequal (Lantz, 2013; Delacre et al., 2019)—a situation that often occurs in practice (Sladekova & Field, 2024c).

In a frequentist framework, an estimator’s ability to control the rate of false positives is instrumental and should be the starting point in evaluating robustness if the goal is to test the null hypothesis. Other metrics of robustness which have come to the forefront of considerations for statistical analyses only in the past couple of decades, include statistical power, estimation accuracy, and the coverage of confidence intervals (Cumming, 2014). Power is the probability of correctly rejecting a false null hypothesis. Informally, it is the probability of detecting an effect of a specified magnitude as statistically significant at a given alpha level (Cohen, 2013), and it is a crucial aspect of the robustness of OLS models. Both skewness and high kurtosis can reduce statistical power of OLS models below the recommended threshold of 80% (Glass, Peckham & Sanders, 1972; Delacre et al., 2019; Nwobi & Akanno, 2021; Kim & Li, 2023; Sladekova & Field, 2024a)—this fact is also highlighted in Glass, Peckham & Sanders (1972) but is often overlooked by researchers. Kurtosis produces distributions with heavier than normal tails—sampling from these distributions also reduces the coverage of confidence intervals (Kim & Li, 2023), and increases the probability of encountering outliers and influential cases, which can bias the parameter estimates (Kim & Li, 2023; Sladekova & Field, 2024a). Power plays an important role in the recent efforts to improve the credibility of psychology as a science (Open Science Collaboration, 2015; Chambers, 2017; Vazire, 2018), because inadequate power combined with poor control over false positives can alter the landscape of the available research findings. Studies with low power are less likely to reach statistical significance and are therefore less likely to be published (Rosenthal, 1979; Sterling, Rosenbaum & Weinkam, 1995; Fanelli, 2010a). Conversely, underpowered research that does reach statistical significance is more likely to reflect a ‘lucky sample’ and is less likely to replicate, but will stand a greater chance of becoming part of the published record (Sterne, Gavaghan & Egger, 2000). Attempts to synthesise or reproduce the effects found in published literature are seriously undermined by this kind of publication bias.

Thus far, discussions around power analyses have focused on the process of specifying the right sample size for a given effect size in order to reach sufficient statistical power for a hypothesis test (Cohen, 2013). Historically, sample sizes collected in psychological research have been low. In a series of investigations, Holmes (1979), Holmes, Holmes & Fanning (1981) and Holmes (1983) reported no changes in sample sizes between 1950 and 1970, with median samples of 55 across four areas of psychology (Abnormal Psychology, Applied Psychology, Developmental Psychology, and Experimental Psychology), while Marszalek et al. (2011) found the median sample size across these areas to be even lower (n = 40) when replicating Holmes’s (1983) study three decades later. Since then, sample sizes have improved in some areas of psychology. Sassenberg & Ditrich (2019) reports doubling of sample sizes between 2009 and 2018, with the lowest increase being from 122 to 185. Clinical research has seen the largest increase, with average sample size around 180 and statistical power to detect a moderate effect just below 90% in 2019 (Reardon et al., 2019). Doubling of sample sizes has also been observed in social and personality psychology between 2011 and 2019 (Fraley et al., 2022), with average sample size of 104 reported in 2014 (Fraley & Vazire, 2014). Other areas are still lagging on improvements—in applied and experimental behavioral research, most studies report samples up to 10 (Schrimp et al., 2022), while Hussey (2023) reports an average sample of n = 64 in Implicit Relational Assessment Procedure research, translating into statistical power of about 34%.

Sample sizes have therefore been at the forefront of focus in response to wide-spread replicability failures (Open Science Collaboration, 2015). Issues that can further affect statistical power, like violated assumptions, have not been given much spotlight because when assumptions are considered, it is usually as an afterthought rather than something that is accounted for in advance.

Violated assumptions as a source of analytic flexibility

The possibility of violated assumptions or the presence of outliers and influential cases should be accounted for during the preparation of an analytic plan. This can prove to be challenging in practice. Methods that rely on pre-defined cut-off points for decision making (like significance tests) are easy to plan for in advance, but, as highlighted above, they are often not an appropriate tool for the job. Other methods—like for example diagnostic plots of residuals and fitted values—might provide a more nuanced picture of the problems, but they also rely on the subjective judgement of the researcher and often highlight only the most blatant violations (Hayes & Cai, 2007a; Darlington & Hayes, 2017). Decisions based on these methods can only be made after the data have been observed, and are inevitably subject to bias (Fanelli, 2009; Steegen et al., 2016).

While data-driven decision making is not problematic in its own right and is often utilised, for example, as part of Bayesian estimation and hypothesis testing (Kruschke & Liddell, 2018; McElreath, 2020), it is at odds with null hypothesis significance testing and the use of p-values as long run probabilities. Any analytic decisions made after the data have been observed can invalidate the p-value associated with test statistic for the hypothesis. The p-value is conditional on the decisions made about the study design and the analytic strategy. If these decisions are made before the data are observed, the p-value remains valid as a long run probability and the finding is potentially replicable, as long as the steps of the original study are being reproduced. If the analytic decisions are made after observing the data, the p-value becomes conditional on those decisions. In such situations, the meaning of the p-value changes to the probability of detecting the effect of the observed or greater magnitude if the null hypothesis is true given a specific set of decisions following a specific line of reasoning based on specific characteristics of the sample (Greenland et al., 2016). Replicating a finding based on a p-value conditional on one researcher’s reasoning might be an impossible undertaking. On the one hand, there is great amount of flexibility when it comes to performing statistical analyses even in the simplest scenarios. Studies of the researchers’ analytic degrees of freedom (Steegen et al., 2016) show that multiple researchers will come to divergent conclusions when faced with the same data and the same hypothesis (Silberzahn et al., 2018; Botvinik-Nezer, 2020; Bastiaansen et al., 2020; Breznau et al., 2022). The present paper alone has so far outlined only a small number of methods that could be combined and applied in a variety of ways, leading to divergent analytic paths. On the other hand, there is no guarantee that the same type of assumption violation or the same degree of violation will be detected in the sample of the replication study, in which case a new set of sample-based decisions needs to be made.

Flexible analytic decision making based on sample characteristics is therefore problematic in its own right, and it has recently received increased attention as a factor potentially contributing to the publication of spurious unreplicable effects (Simmons, Nelson & Simonsohn, 2011; Wagenmakers et al., 2011). The existence of publication bias favouring statistically significant results (Ferguson & Heene, 2012) coupled with the incentive structures in academia and the ‘publish or perish’ culture (Fanelli, 2010b) puts researchers under a great amount of pressure to produce ‘publishable’ results. As such their decision making during the analytic process after observing the data may become biased towards obtaining a statistically significant result, even in the absence of deliberate attempts to manipulate the data. The efforts to improve the credibility of psychology have resulted in a number of proposed methodological reforms to address this. One such reform was the introduction of preregistration (Nosek et al., 2018), where researchers submit their analytic plans into a time-stamped online repository prior to collecting or accessing the data, and follow through with this plan when performing the analysis. Registered reports as publication format take this idea further—researchers submit a detailed study protocol, including the analysis plan, for peer-review and acceptance in principle is granted before the results are known, as long as the protocol is followed (Chambers et al., 2015).

There’s growing evidence that registered reports are successful at combating publication bias (Scheel, Schijen & Lakens, 2021), however, preregistering a realistic plan can prove challenging when we consider the complications arising with violations of OLS conditions. As discussed, current methods for detecting violations require data inspection and post hoc decision making, while the most commonly applied ‘remedies’ (rank-based methods and data transformations) transform the data into a form that no longer maps onto the original hypotheses, introducing unnecessary degrees of freedom into the analytic process.

An ideal solution would be to preregister an alternative estimator with consistent performance across a variety of error distributions which also retains interpretation of parameter estimates comparable to OLS estimates. While no single method that meets these conditions outperforms all possible alternatives in all situations, a class of methods known as robust statistical methods offers a range of estimators that outperform OLS in a vast majority of scenarios in terms of statistical power and accuracy of estimates, and can therefore complement the efforts to improve the credibility of findings in psychology.

What are robust methods?

A statistical method can be considered robust if (1) it has adequate control over the rate of false positive findings, (2) it retains adequate power to detect a true effect under a variety of scenarios, (3) the point estimates produced by the method are not overly sensitive to outliers—i.e., they remain unbiased (4) the estimates accurately describe the typical individual in the sample, and (5) the previous points hold true regardless of whether or not the assumptions of OLS are violated.

“Robust methods” is an umbrella term encompassing a wide variety of methods that meet these conditions. This section will focus on some examples of robust methods that can be split into two general categories—methods improving the estimation of standard errors and confidence intervals constructed around the original OLS point estimates, and methods improving the estimation of point estimates as well as the estimates of standard errors and confidence intervals. The first category includes bootstrapping (Efron & Tibshirani, 1993) and heteroscedasticity-consistent standard errors (Hayes & Cai, 2007a), while the second category includes M-estimators (Huber, 1964), and robust trimming (Yuen, 1974; Wilcox, 1998b). These four methods represent only small selection of all available robust methods. In an ideal world, the researcher would select a method that is known for optimal performance given the problem at hand. In the real world, it is unreasonable to expect the researchers to have an in depth knowledge about the performance and application of every single method there is, and statisticians or methodologists who would be able to advise on the issue while also understanding the research context may not be readily available in most psychology departments (Golinski & Cribbie, 2009). We focus on the methods outlined above because they are fairly intuitive extensions of the standard methods, easily applied to common research designs in psychology, and guidance already exists for their application (Hayes & Cai, 2007a; Wilcox, 2017; Field & Wilcox, 2017). The following section provides an introduction to these methods, followed by a discussion of the benefits of their application.

Robust standard errors and confidence intervals

Robust point estimates

A point estimate for the relationship between variables that we are studying should represent a typical individual in a given population as accurately as possible and produce standard error of size that is not detrimental to statistical power. No estimator beats OLS when sampling from a normal distribution. However as we have seen, such distributions are an exception rather than the rule in applied settings. OLS can be especially volatile in heavy-tailed distributions that are likely to contain outliers. The variance computed using OLS will be too large and therefore the power to detect a real statistically significant effect will be low. Additionally, OLS becomes biased when outliers are present. This is related to its low finite sample breakdown point, which is the proportion of extreme values in the sample that can cause the estimate to be biased in the direction of these values. Each estimator has its own breakdown point, and for OLS this value is 1/n, n being the sample size. In proportional terms, the value gets smaller as the sample size increases, therefore large sample sizes are not protected from the effects of outliers.

M-estimators

Among the robust methods that show substantially better performance than the mean (which is an OLS estimator) are the trimmed mean and M-estimators. M-estimators are among the robust methods that show substantially better performance in situations where OLS becomes biased or produces non-optimal estimates. During M-estimation, extreme observations are weighted down to lessen their impact on the point estimate. Which observations are weighted down and to what extent is determined algorithmically by the computer based on the properties of the model in question. The researcher can however select a function according to which the down-weighting should be applied.

One example of a weight function is Huber’s ψ (Huber, 1964), which applies smaller weights to observations with large residuals, where the residuals that are beyond a specific point are given a weight of zero—this the equivalent of the scores being trimmed off and not weighing in on the computation of the point estimate. Figure 3A illustrates one form Huber’s ψ can take. Observations with residuals between −1.2 and +1.2 are assigned the same weights as they would have in and OLS estimation, whereas the observations with the residuals outside of this boundary are given a weight of zero—this is represented by the horizontal lines at the tails of the weight function. As such, M-estimation is always completed in two or more steps. First, an OLS model is fitted and the weights are determined based on the residuals in this model, then an M-estimate is computed. The process is iteratively repeated until the model converges (Susanti et al., 2014). The breakdown point of M-estimators is 0.5, meaning that 50% of scores can be at the extreme ends of the curve and the point estimate will remain at the same location. This is a maximum possible breakdown point an estimator can have. As a result of the weighting procedure, the standard errors in the model that uses M-estimation are smaller than they would be in an OLS model when sampling from heavy tailed distribution, thus increasing the power of the associated statistical significance tests (Wilcox, 2010).

M-estimators using Huber’s ψ for weighting can however be inefficient under heteroscedasticity and when there are influential cases in the sample (Croux, Dhaene & Hoorelbeke, 2003; Toka & Cetın, 2011). Extensions of this estimator have been derived over the years to deal with these issues. We focus on the S-estimator, MM-estimator, and DAS-estimator, although other versions of the M-estimator exist (Wilcox, 2017). Each of these estimators builds on the previous version by iterating through the original steps and then applying a unique adjustment. The S-estimator affects the variance and aims to minimise the dispersion of the residuals around the M-estimate (Rousseuw & Yohai, 1984; Toka & Cetın, 2011). The MM-estimator (Yohai, 1987) applies an alternative redescending function (often Tukey’s biweight function, Fig. 3B) to weigh the residuals and re-estimate the point estimate originally produced by M-estimation and the variance estimate produced by S-estimation in the previous steps. The weakness of the MM-estimator is that it becomes inefficient in designs with small samples (below 20), but also in larger samples if the ratio of the number of predictors in the model to the number of observations (p/n ratio) is greater than 1/10. In other words, if there are less than 10 observations for each predictor in the model, the estimate suffers. High p/n ratio also causes bias in the estimation of variance (Croux, Dhaene & Hoorelbeke, 2003). The Design Adaptive Scale estimator (DAS; Koller & Stahel, 2011) was developed to correct this bias and lack of efficiency, and can do so successfully for ratios of up to p/n = 1/3. The function applied by this estimator is also less steep than the typical the biweight function typically applied in MM-estimation, which accounts for overly aggressive down-weighting of observations as outliers in heteroscedastic models.