Robust statistics

The plot below shows a density plot of the speed-of-light data, together with a rug plot (panel (a)). Also shown is a normal Q–Q plot (panel (b)). The outliers are clearly visible in these plots.

Panels (c) and (d) of the plot show the bootstrap distribution of the mean (c) and the 10% trimmed mean (d). The trimmed mean is a simple robust estimator of location that deletes a certain percentage of observations (10% here) from each end of the data, then computes the mean in the usual way. The analysis was performed in R and 10,000 bootstrap samples were used for each of the raw and trimmed means.

The distribution of the mean is clearly much wider than that of the 10% trimmed mean (the plots are on the same scale). Also whereas the distribution of the trimmed mean appears to be close to normal, the distribution of the raw mean is quite skewed to the left. So, in this sample of 66 observations, only 2 outliers cause the central limit theorem to be inapplicable.

Robust statistical methods, of which the trimmed mean is a simple example, seek to outperform classical statistical methods in the presence of outliers, or, more generally, when underlying parametric assumptions are not quite correct.

Whilst the trimmed mean performs well relative to the mean in this example, better robust estimates are available. In fact, the mean, median and trimmed mean are all special cases of M-estimators. Details appear in the sections below.

The outliers in the speed-of-light data have more than just an adverse effect on the mean; the usual estimate of scale is the standard deviation, and this quantity is even more badly affected by outliers because the squares of the deviations from the mean go into the calculation, so the outliers' effects are exacerbated.

The plots below show the bootstrap distributions of the standard deviation, the median absolute deviation (MAD) and the Rousseeuw–Croux (Qn) estimator of scale.[2] The plots are based on 10,000 bootstrap samples for each estimator, with some Gaussian noise added to the resampled data (smoothed bootstrap). Panel (a) shows the distribution of the standard deviation, (b) of the MAD and (c) of Qn.

The distribution of standard deviation is erratic and wide, a result of the outliers. The MAD is better behaved, and Qn is a little bit more efficient than MAD. This simple example demonstrates that when outliers are present, the standard deviation cannot be recommended as an estimate of scale.

Traditionally, statisticians would manually screen data for outliers, and remove them, usually checking the source of the data to see whether the outliers were erroneously recorded. Indeed, in the speed-of-light example above, it is easy to see and remove the two outliers prior to proceeding with any further analysis. However, in modern times, data sets often consist of large numbers of variables being measured on large numbers of experimental units. Therefore, manual screening for outliers is often impractical.

Outliers can often interact in such a way that they mask each other. As a simple example, consider a small univariate data set containing one modest and one large outlier. The estimated standard deviation will be grossly inflated by the large outlier. The result is that the modest outlier looks relatively normal. As soon as the large outlier is removed, the estimated standard deviation shrinks, and the modest outlier now looks unusual.

This problem of masking gets worse as the complexity of the data increases. For example, in regression problems, diagnostic plots are used to identify outliers. However, it is common that once a few outliers have been removed, others become visible. The problem is even worse in higher dimensions.

Robust methods provide automatic ways of detecting, downweighting (or removing), and flagging outliers, largely removing the need for manual screening. Care must be taken; initial data showing the ozone hole first appearing over Antarctica were rejected as outliers by non-human screening.[3]

Although this article deals with general principles for univariate statistical methods, robust methods also exist for regression problems, generalized linear models, and parameter estimation of various distributions.

Intuitively, the breakdown point of an estimator is the proportion of incorrect observations (e.g. arbitrarily large observations) an estimator can handle before giving an incorrect (e.g., arbitrarily large) result. For example, given ${\displaystyle n}$ independent random variables ${\displaystyle (X_{1},\dots ,X_{n})}$ and the corresponding realizations ${\displaystyle x_{1},\dots ,x_{n}}$, we can use ${\displaystyle {\overline {X_{n}}}:={\frac {X_{1}+\cdots +X_{n}}{n}}}$ to estimate the mean. Such an estimator has a breakdown point of 0 because we can make ${\displaystyle {\overline {x}}}$ arbitrarily large just by changing any of ${\displaystyle x_{1},\dots ,x_{n}}$.

The higher the breakdown point of an estimator, the more robust it is. Intuitively, we can understand that a breakdown point cannot exceed 50% because if more than half of the observations are contaminated, it is not possible to distinguish between the underlying distribution and the contaminating distribution Rousseeuw & Leroy (1986) harvtxt error: no target: CITEREFRousseeuwLeroy1986 (help). Therefore, the maximum breakdown point is 0.5 and there are estimators which achieve such a breakdown point. For example, the median has a breakdown point of 0.5. The X% trimmed mean has breakdown point of X%, for the chosen level of X. Huber (1981) and Maronna, Martin & Yohai (2006) contain more details. The level and the power breakdown points of tests are investigated in He, Simpson & Portnoy (1990).

Statistics with high breakdown points are sometimes called resistant statistics.[4]

In the speed-of-light example, removing the two lowest observations causes the mean to change from 26.2 to 27.75, a change of 1.55. The estimate of scale produced by the Qn method is 6.3. We can divide this by the square root of the sample size to get a robust standard error, and we find this quantity to be 0.78. Thus, the change in the mean resulting from removing two outliers is approximately twice the robust standard error.

The 10% trimmed mean for the speed-of-light data is 27.43. Removing the two lowest observations and recomputing gives 27.67. Clearly, the trimmed mean is less affected by the outliers and has a higher breakdown point.

If we replace the lowest observation, −44, by −1000, the mean becomes 11.73, whereas the 10% trimmed mean is still 27.43. In many areas of applied statistics, it is common for data to be log-transformed to make them near symmetrical. Very small values become large negative when log-transformed, and zeroes become negatively infinite. Therefore, this example is of practical interest.

Tukey's biweight function

The empirical influence function is a measure of the dependence of the estimator on the value of one of the points in the sample. It is a model-free measure in the sense that it simply relies on calculating the estimator again with a different sample. On the right is Tukey's biweight function, which, as we will later see, is an example of what a "good" (in a sense defined later on) empirical influence function should look like.

In mathematical terms, an influence function is defined as a vector in the space of the estimator, which is in turn defined for a sample which is a subset of the population:

1. ${\displaystyle (\Omega ,{\mathcal {A}},P)}$ is a probability space,
2. ${\displaystyle ({\mathcal {X}},\Sigma )}$ is a measurable space (state space),
3. ${\displaystyle \Theta }$ is a parameter space of dimension ${\displaystyle p\in \mathbb {N} ^{*}}$,
4. ${\displaystyle (\Gamma ,S)}$ is a measurable space,

For example,

1. ${\displaystyle (\Omega ,{\mathcal {A}},P)}$ is any probability space,
2. ${\displaystyle ({\mathcal {X}},\Sigma )=(\mathbb {R} ,{\mathcal {B}})}$,
3. ${\displaystyle \Theta =\mathbb {R} \times \mathbb {R} ^{+}}$
4. ${\displaystyle (\Gamma ,S)=(\mathbb {R} ,{\mathcal {B}})}$,

The definition of an empirical influence function is: Let ${\displaystyle n\in \mathbb {N} ^{*}}$ and ${\displaystyle X_{1},\dots ,X_{n}:(\Omega ,{\mathcal {A}})\rightarrow ({\mathcal {X}},\Sigma )}$ are i.i.d. and ${\displaystyle (x_{1},\dots ,x_{n})}$ is a sample from these variables. ${\displaystyle T_{n}:({\mathcal {X}}^{n},\Sigma ^{n})\rightarrow (\Gamma ,S)}$ is an estimator. Let ${\displaystyle i\in \{1,\dots ,n\}}$. The empirical influence function ${\displaystyle EIF_{i}}$ at observation ${\displaystyle i}$ is defined by:

${\displaystyle EIF_{i}:x\in {\mathcal {X}}\mapsto n\cdot (T_{n}(x_{1},\dots ,x_{i-1},x,x_{i+1},\dots ,x_{n})-T_{n}(x_{1},\dots ,x_{i-1},x_{i},x_{i+1},\dots ,x_{n}))}$

What this actually means is that we are replacing the i-th value in the sample by an arbitrary value and looking at the output of the estimator. Alternatively, the EIF is defined as the (scaled by n+1 instead of n) effect on the estimator of adding the point ${\displaystyle x}$ to the sample.

Instead of relying solely on the data, we could use the distribution of the random variables. The approach is quite different from that of the previous paragraph. What we are now trying to do is to see what happens to an estimator when we change the distribution of the data slightly: it assumes a distribution, and measures sensitivity to change in this distribution. By contrast, the empirical influence assumes a sample set, and measures sensitivity to change in the samples.[5]

Let ${\displaystyle A}$ be a convex subset of the set of all finite signed measures on ${\displaystyle \Sigma }$. We want to estimate the parameter ${\displaystyle \theta \in \Theta }$ of a distribution ${\displaystyle F}$ in ${\displaystyle A}$. Let the functional ${\displaystyle T:A\rightarrow \Gamma }$ be the asymptotic value of some estimator sequence ${\displaystyle (T_{n})_{n\in \mathbb {N} }}$. We will suppose that this functional is Fisher consistent, i.e. ${\displaystyle \forall \theta \in \Theta ,T(F_{\theta })=\theta }$. This means that at the model ${\displaystyle F}$, the estimator sequence asymptotically measures the correct quantity.

Let ${\displaystyle G}$ be some distribution in ${\displaystyle A}$. What happens when the data doesn't follow the model ${\displaystyle F}$ exactly but another, slightly different, "going towards" ${\displaystyle G}$?

We're looking at: ${\displaystyle dT_{G-F}(F)=\lim _{t\rightarrow 0^{+}}{\frac {T(tG+(1-t)F)-T(F)}{t}}}$,

which is the one-sided Gateaux derivative of ${\displaystyle T}$ at ${\displaystyle F}$, in the direction of ${\displaystyle G-F}$.

Let ${\displaystyle x\in {\mathcal {X}}}$. ${\displaystyle \Delta _{x}}$ is the probability measure which gives mass 1 to ${\displaystyle \{x\}}$. We choose ${\displaystyle G=\Delta _{x}}$. The influence function is then defined by:

${\displaystyle IF(x;T;F):=\lim _{t\rightarrow 0^{+}}{\frac {T(t\Delta _{x}+(1-t)F)-T(F)}{t}}.}$

It describes the effect of an infinitesimal contamination at the point ${\displaystyle x}$ on the estimate we are seeking, standardized by the mass ${\displaystyle t}$ of the contamination (the asymptotic bias caused by contamination in the observations). For a robust estimator, we want a bounded influence function, that is, one which does not go to infinity as x becomes arbitrarily large.

Properties of an influence function which bestow it with desirable performance are:

1. Finite rejection point ${\displaystyle \rho ^{*}}$,
2. Small gross-error sensitivity ${\displaystyle \gamma ^{*}}$,
3. Small local-shift sensitivity ${\displaystyle \lambda ^{*}}$.

${\displaystyle \rho ^{*}:=\inf _{r>0}\{r:IF(x;T;F)=0,|x|>r\}}$

${\displaystyle \gamma ^{*}(T;F):=\sup _{x\in {\mathcal {X}}}|IF(x;T;F)|}$

${\displaystyle \lambda ^{*}(T;F):=\sup _{(x,y)\in {\mathcal {X}}^{2} \atop x\neq y}\left\|{\frac {IF(y;T;F)-IF(x;T;F)}{y-x}}\right\|}$

This value, which looks a lot like a Lipschitz constant, represents the effect of shifting an observation slightly from ${\displaystyle x}$ to a neighbouring point ${\displaystyle y}$, i.e., add an observation at ${\displaystyle y}$ and remove one at ${\displaystyle x}$.

M-estimators do not necessarily relate to a probability density function. Therefore, off-the-shelf approaches to inference that arise from likelihood theory can not, in general, be used.

It can be shown that M-estimators are asymptotically normally distributed, so that as long as their standard errors can be computed, an approximate approach to inference is available.

Since M-estimators are normal only asymptotically, for small sample sizes it might be appropriate to use an alternative approach to inference, such as the bootstrap. However, M-estimates are not necessarily unique (i.e., there might be more than one solution that satisfies the equations). Also, it is possible that any particular bootstrap sample can contain more outliers than the estimator's breakdown point. Therefore, some care is needed when designing bootstrap schemes.

Of course, as we saw with the speed-of-light example, the mean is only normally distributed asymptotically and when outliers are present the approximation can be very poor even for quite large samples. However, classical statistical tests, including those based on the mean, are typically bounded above by the nominal size of the test. The same is not true of M-estimators and the type I error rate can be substantially above the nominal level.

These considerations do not "invalidate" M-estimation in any way. They merely make clear that some care is needed in their use, as is true of any other method of estimation.

It can be shown that the influence function of an M-estimator ${\displaystyle T}$ is proportional to ${\displaystyle \psi }$,[6] which means we can derive the properties of such an estimator (such as its rejection point, gross-error sensitivity or local-shift sensitivity) when we know its ${\displaystyle \psi }$ function.

${\displaystyle IF(x;T,F)=M^{-1}\psi (x,T(F))}$

with the ${\displaystyle p\times p}$ given by:

${\displaystyle M=-\int _{\mathcal {X}}\left({\frac {\partial \psi (x,\theta )}{\partial \theta }}\right)_{T(F)}\,dF(x).}$

In many practical situations, the choice of the ${\displaystyle \psi }$ function is not critical to gaining a good robust estimate, and many choices will give similar results that offer great improvements, in terms of efficiency and bias, over classical estimates in the presence of outliers.[7]

Theoretically, ${\displaystyle \psi }$ functions are to be preferred, and Tukey's biweight (also known as bisquare) function is a popular choice. Maronna, Martin & Yohai (2006) recommend the biweight function with efficiency at the normal set to 85%.

For the speed-of-light data, allowing the kurtosis parameter to vary and maximizing the likelihood, we get

${\displaystyle {\hat {\mu }}=27.40,{\hat {\sigma }}=3.81,{\hat {\nu }}=2.13.}$

Fixing ${\displaystyle \nu =4}$ and maximizing the likelihood gives

${\displaystyle {\hat {\mu }}=27.49,{\hat {\sigma }}=4.51.}$