1693

ANALYST. DECEMBER 1989, VOL. 114

Robust Statistics- How Not to Reject OutliersPart 1. Basic Concepts
Analytical Methods Committee*
Royal Society of Chemistry, Burlington House, Piccadilly, London W1V OBN, UK

The subject of outliers has been controversial whenever analytical data have been processed. Modern
statistical theory provides an alternative t o outlier rejection, in which outlying observations are retained but
given less weight. This approach is known as robust statistics and is beginning to find favour with analyticalchemists. An introduction t o robust statistics is given and some examples are described.

Keywords: Outlier; robust statistics; mean; median; variance

The Analytical Methods Committee has received and has
approved for publication the following report from its
Statistical Sub-Committee.

The constitution of the Sub-committee responsible for the
preparation of this report was: Dr. M. Thompson (Chairman),Mr. H. M. Bee, Dr. W. H. Evans, Mr. M. J. Gardner, Dr. E. J.
Greenhow, Dr. R. Howarth, Dr. E. J. Newman, Professor
B. D. Ripley and Dr. R. Wood with Mr. J. J. Wilson as
Secretary.
The Analytical Methods Committee gratefully acknowledges the financial support given by the Ministry of Agriculture, Fisheries and Food to the work of this Sub-committee.

i

..
...
....
.....
.

0

5.0

5.0

10.0

10.0

25.0

15.0
20.0
Copper, p.p.m.

30.0

I

I

I

I

15.0

20.0

25.0

30.0

I

Copper, p.p.m.

Fig. 1. Two views of dataset (1) from the ctstistical package
MINITAB. ( a ) A dot plot; and ( b ) a box plot. * and 0,
extremeobservations;the median; and 0,
the quantiles+.

where x(l), ..., x(24) are the observations sorted into
increasing order:

Introduction
Occasionally sets of analytical data occur in which a few
observations appear discordant with the remainder. Suchobservations are known as outliers. For example, considering
the following 24 determinations of copper (vgg-1) inwholemeal flour

2.9 3.1 3.4 3.4 3.70 3.7 2.8 2.5 2.4 2.4 2.7 2.2
5.28 3.37 3.03 3.03 28.95 3.77 3.4 2.2 3.5 3.6 3.7 3.7 ( l )one value, 28.95, stands out from the remainder. In thisinstance we may be particularly suspicious of the value, as asimple explanation suggests itself. Although recording andrange errors are almost certainly the major cause of outliers,
mistakes can also occur in many other parts of the analyticalprocess and from contamination and transposition of specimens.
The almost universal practice amongst analytical chemists
has been to regard outliers as errors, and to delete them fromthe set of data. In some circumstances this is plainly wrong,
and in others there are much safer procedures. Why should we
be interested in outliers? One good reason is to catch
transcription errors while the original laboratory records areeasily accessible. In such an instance we would want to checkall the extreme results whether or not they are rejected by an
outlier test (such as the tests of Dixonl or Grubbs2). The
traditional procedure3 for dataset (1) would be to compute

[

(a)

0

Report~ ( 2 4) 422)
max 4 3 ) - ~ ( 1 )) 43)
4 2 2 ) - ~ ( 1 ) ’~ ( 2 4 -II= 0.948* Correspondence should be addressed to the Secretary, Analytical
Methods Committee, Analytical Division, Royal Society ofChemistry, Burlington House, Piccadilly, London W1V OBN, UK.

2.2 2.2 2.4 2.4 2.5 2.7 2.8 2.9 3.03 3.03 3.1 3.37
3.4 3.4 3.4 3.5 3.6 3.7 3.7 3.7 3.70 3.77 5.28 28.95
and therefore reject x(24) = 28.95. For the remaining 23
observations the Dixonl test statistic is 0.549, so ~ ( 2 3 =) 5.28
would also be rejected. Using the test yet again gives 0.133,
this being judged not significant. Grubbs’ test gives the same
results.
The traditional procedure has the merit of pointing out thesecond outiier, 5.28, but this would have been obvious fromany plot of the data (Fig. 1).The second-largest value will besignificant (at 5 % ) only if it exceeds 4.80. However, surely wewould want to check a value of 4.77 for a transcription error?
If the purpose of detecting outliers is to check the values, it istheir extremeness and the plausibility of simple explanations
that should weigh in the decision, and not the statisticalsignificance.
Outlier rejection is positively wrong when included in aprocedure to assess the variability of an analytical method.
The outlier rejection procedure used above is that of BS54974;
the illustrative dataset (1) is taken from a co-operative trial.5
The mean and variance of the whole set are 4.28 and 28.1,
respectively, whereas after outlier rejection they are 3.11 and
0.281, respectively. If the second-largest observation had been
4.77 we would have obtained 3.19 and 0.387, respectively.
From this, two conclusions can be drawn that are true
generally. The traditional procedure is (a) sensitive to theactual data values and (b) seriously underestimates thevariance that is attainable in practice. The outliers in ourdataset are only revealed because we have 24 replicates.
Duplication might throw doubt on a value of 28.95, but it isvery unlikely to do so on 5.28. On the other hand, estimating
the variance by 28.1 is also unfair, as values as large as 28195
might be spotted and are much rarer than 1 in 24.

1694

ANALYST, DECEMBER 1989, VOL. 114

Barnett and Lewis6 discussed the outlier problem inconsiderable detail and described a whole battery of outlier
rejection tests. The change of emphasis from their first edition(1978) to the second edition (1984) reflects a change instatistical practice from outlier rejection to outlier accommodation. The prevailing philosophy is known as robust statistics(or, occasionally, as resistant statistics) and is expounded in a
number of recent monographs,7-10 some of which are forbidding even to professional statisticians.
Robust statistics have been used occasionally by chemists,
especially in geochemistry.11-15 These papers concentrate on
establishing reference values, whereas robust methods can be
as useful in assessing variability as for central tendency.

Philosophy of Robust Statistics
The normal distribution pervades statistical methodology, andits very name suggests widespread applicability. Yet careful
studies show that real errors do not fit the normal distribution!
Users of statistics point to a theoretical result, the central limit
theorem, to justify the assumption, whereas theoreticiansbelieve its applicability to have been proved empirically. The
central limit theorem is, of course, a perfectly correct result
about sums of many small independent errors having (approximately) a normal distribution. The problem is that outliers
result from single large errors. It is also generally acceptedthat real error distributions have “heavier tails” than thenormal distribution, i.e., large deviations (in either direction)
are more likely than under a normal distribution. One of thebases of robust statistics is to use procedures that work well forsuch distributions.
The second basis is to protect against gross errors. We
observed from our example that recording 28.95 rather than
2.895 increased the sample mean considerably (to 4.28 from
3.19). Recording 289.5 and 2895 would give 15.1 and 123.7,
respectively. Hence the effect of a missing decimal point is
disastrous for the mean. On the other hand, the median isalmost unchanged, from 3.24 with 2.895 to 3.38 with any value
greater than 3.40. This property is shared by a trimmed mean.
Suppose we discard the smallest r and the largest r observations out of the total, n, and then take the mean of theremainder. [This is called a (lOOr/n)(%) trimmed mean.]
Discarding r = 1and r = 2 gives estimates of the mean of 3.25
and 3.21, respectively, both being insensitive to the actual sizeof 28.95. Trimmed means obey both principles of robust
statistics; they are insensitive to small numbers of gross errors,
and they work well for heavy-tailed distributions close to thenormal. The mean fails the first of these, the median thesecond.
The sample variance is even more sensitive to outliers than
the sample mean, increasing from 0.46 (when 28.95 is replacedby 2.895) to 28.1. The inter-quartile range (IQR) is thedifference between observations one quarter in from eachend, the 6th and 19th in the present example, so IQR = 1.0.
For a normal distribution the IQR would be expected to beabout 1.350, which suggests that we determine 02 by (IQW
1.35)2 = 0.55. Again this is insensitive to gross errors; it is infact unchanged if 28.95 is replaced by 2.895 or 289.5.
The trimmed mean and the IQR were developed in the days
when all data were analysed by hand. Computers have allowed
more sophisticated methods to be used. Many such procedures exist, but we will only describe some of the simplest
which are known to perform well. A whole book16 has been
devoted to comparisons of 68 procedures! The references
disagree as to the best procedure, but all accept thosedescribed here as amongst the best.

Measuring “True Values”
Both the sample mean and median are estimates of thelocation of a distribution of results. This distribution can be-C

Fig. 2.

C

A robust loss function p. Note that p(x) = x2 for 1x1 d c, but

p(x) < x2 for large 1x1

considered as a “true value”, p, plus errors, and we want tofind an estimate of p. We assume that p is the mean of the“reliable” results, but not of the whole error distribution, as inanalytical chemistry the distribution of errors will almost
always be asymmetric. Consider n data points xl, ..., x,. The
sample mean minimises the sum of squares SS = X(xi - p)2
and this is the source of its sensitivity to gross errors as large
errors inflate SS significantly. Suppose we minimise S S =Cp(xi - p) where P(E) does not weight large errors, E, as much
. good choice is the functionas ~ 2 A

illustrated in Fig. 2 where 02 is a robust variance and c is aconstant in the range 1-2. This penalises errors larger than coless severely than x2. The corresponding location estimate, $,is the mean of pseudo-values Ti{

xi

f j = p-co
p+co

if Ixi- $1 ~ c o
if x i < f i - c o..
if x j > p + c o. . (2)and also the weighted mean of x i , with weights w i
1if (xi- F(s cowi = {coilx, - fiI if Ixi - $1 > coHence extreme values can be thought of as being either
brought in or downweighted. We can compute p from either ofthese properties. To start with take any estimate fi(”), say themean or the median. At each stage compute fiG) as theweighted mean with weights min(l,co/lxi - (lU-1))) or as themean of the values fi[with p = $Q- ”1. [The function min(x,y)
denotes the smaller value of x and y,] The values PO) converge
rapidly to p.
The value of c = 1.5 has wide support. Suppose we knew oto be 0.70, then starting from the meanpG) = 4.28, 3.56, 3.27, 3.22, 3.21, ...
and starting from the median= 3.39, 3.24, 3.21, ...

It is usually unrealistic to assume that o is known, although
only a rough estimate is needed, which might be available
from past trials. One rough estimate is based on the medianabsolute deviation (MAD); MAD = median ( / x i- medianl),
6 = MAD/0.6745, which is similar to the re-scaled IQR, and in
our example gives 6 = 0.53 and fi= 3.207. (The scale factor
0.6745 is used to obtain the correct answer for normally
distributed data.) More sophisticated estimates are consideredlater. However, this simple proposal is already very reliable.
The corresponding estimator of location is sometimes known
as A15.

ANALYST, DECEMBER 1989, VOL. 114

1695

Table 1. Values of the constants 1-3 and 0 for a range of cut-off values c. Further values can be obtained from 8 =
P ( ( N (< c), f3 = 0 + cZ(1 - 0) - 2c exp(-c2/2)/V2~1where N is a standard normal deviate
C

1-3
8

1.0
0.516
0.683

1.1
0.578
0.729

1.2
0.635
0.770

1.3
0.688
0.806

1.4
0.736
0.838

1.5
0.778
0.866

1.6
0.816
0.890

1.7
0.849
0.911

1.8
0.877
0.928

1.9
0.900
0.943

2.0
0.921
0.954

Measuring “Precision”

Discussion

If it is supposed that we are in the somewhat more realistic
position of knowing p (say from a reference sample) and wish
to estimate cr, then we could use the sample variance or the
scaled IQR, or

We have observed that robust procedures can be constructed
to estimate the true value and precision of a set of data by
relatively simple iterative calculations. These are very tediousto do manually but easy to program. (The longest part of the
program will be to find the starting values, see Appendix.)
The robust estimates are completely insensitive to how
outlying the extreme data values are and obtain most of theirinformation from the values in the centre of the dataset.
The cut-off value c should in theory be chosen depending on
how frequent outliers are thought to be, although it is safer tochoose a smaller value of c if in doubt. About 1% of outliers
suggest c = 2.0 and about 5 % suggest c = 1.4. The value c =1.5 is widely used. The actual estimates obtained are not very
sensitive to c:6, = median(lxj - p.))/0.6745
A robust procedure is to solvec min(lxi - pl/a,c)2 = np..

. . (3)where again p is chosen to obtain the correct answer for
normally distributed data. Some values of p are given in Table
1.
There are a number of ways to solve (3). One of the easiest
is to compute a sequence of values 6.0’)with 6(0)= 6, andc = 1.0
p = 3.229
6 = 0.648

where

1

xi

zj= p.-c6.0’-1)
p+c6(.-1)

if (x- p( < c d - 1 )if xi<p--&0‘-1)
if x i > p + t c r ( . - I )which converges rapidly to the solution of (3).
Suppose for dataset (1) that we knew that p. = 3.68 (which isthe consensus of a much larger set of measurements), then
6i = 0.911, 0.927, 0.934, 0.938, 0.939, 0.940, 0.941, ...

so cr can be estimated by 0.941. Note that p62 is the variance ofthe pseudo-values ifi (with divisor n as p is known).
Unknown “true value”
The more sophisticated approach referred to in a previous
section involves estimating cr alongside p.. Hence at each
iteration we form pseudo-values ZLand compute their mean, X,
and variance, s2. Then @(i) = X and 6(i) = d(s2/p). This isrepeated until the values stabilise, starting from (median, 6m).
The present example gives@ = 3.385 3.255 3.213 3.206 3.205 3.205 ... 3.205
6 = 0.526 0.595 0.639 0.657 0.666 0.671 ... 0.674

The following shows the insensitivity to the outlier(s):
X(24)
P
6
28.95
3.205
0.674
2.895
3.146
0.613
0.674
289.5
3.205
and in fact (@,6)do not depend on any of the exact values
greater than @ + c6 = 4.22. Hence it is irrelevant whether thevalue 5.28 is considered an outlier or not; all that matters isthat it exceeds 4.22. This combined estimator is known as H15
or “Huber proposal 2.”
When n, the number of observations, is small, a smallsample correction should be made. The variance of xi - @ will
be about a2(n - l)/n and so the cut-off point, c, should be
reduced to c V‘(1 - l/n) in forming the pseudo-values. This
will be important in Part 2.17 (It reduces 6 to 0.662 in our
example.)

15
3.205
0.662

2.0
3.234
0.678

If we really wanted to look for outliers to check them against
the original records, a useful rule would be to check xivalues
outside @ k 26. In dataset (1) this suggests that all values
greater than 4.53 would be checked.

Use of robust estimates
Some considerable care is needed in interpreting (i and 6.
They do not estimate the mean and standard deviation of theobservations (note, not the population), and this is an asset
rather than a liability. Rare but very large outliers will affectthe theoretical mean p considerably when, as in analyticalchemistry, they will almost always occur in one direction.
Instead we should regard p as measuring the mean of the“reliable” observations, a consensus value which is the nearest
we can get to a “true value”. (This interpretation is only
possible if the “reliable” observations form the majority.
Examples do occur in which the outliers are the only valid
observations, but no statistical procedure can redeem such adisastrous trial.)
In a similar manner 6 measures the standard deviation of
the “reliable” observations. If we take m replicates then therobust measure @ obtained from these will have a variance ofabout 62/m for moderate m. [In fact 62/m X (3/CI2 where 8 =P ( [ x l - p I/a d c), where P is the probability, which can be
estimated either from the normal distribution or by the
proportion of the dataset with 2, = x,. As Table 1 shows, thecorrection factor p/W is only just larger than one.] Howeverthe (population) variance of one observation will usually muchexceed 6 2 , as outliers cannot be downweighted. Duplicatesalso do not help, as we always find that (i = (x1 x#2. At
least three replicates are needed to allow downweighting, andthis may not be sufficient unless cr is known a priori. If we hadrecorded just the three observations 2.9, 3.1 and 28.95, then+s = 14.98
Mean = 11.65
Median = 3.1
6, = 0.297
A15 = 3.222
H15 = 11.65
6 = 16.98

There are two plausible explanations for this triple ofobservations. One, favoured by (A15, a),is of two reliable

ANALYST, DECEMBER 1989, VOL. 114

1696
Program*program rob1
C

C
C

C

program from 'Robust S t a t i s t i c s
Analyst(C) B .D .Ripley- How Not t o Reject O u t l i e r s 'Cparameter (NMAX = 1 0 0 )r e a l x(NMAX), w s (NMAX)r e a l median, h15
i n t e g e r if n
r e a l xmed, xm, xs
c h a r a c t e r name*50
p r i n t *, ' F i l e nameread ( a ) , , nameopen (1, file=name, s t a t u s = , o l d f )p r i n t *, ' # d a t a p o i n t sread *, nread (1,*) ( x ( i ), i=l,n)c l o s e (1)xmed = median(x, n, w s )p r i p t *, ,median ,xxned
xm = a15(x, n, xs, w s )p r i n t *, ,a15, sigma , xm, xs
xm = h l 5 ( x , n, xs, w s )p r i n t *,h15, sigma , xm, xs
end
r e a l function median (x, n, w s )r e a l x ( n ) , ws(*), v
i n t e g e r i f j , h, n, n2
do 1 0 i = 1, n
ws(i) = x(i)10h = l
20 h = 3*h+li f (h . l e . n) goto 20
30 h = h/3
do 50 i = h + l , n
v = ws(i)j = ii f (ws(j-h) . l e . v) goto 5 040w s ( j ) = ws(j-h)
j = j-hi f (j .gt. h) goto 40
50
ws(j) = v
i f (h . g t . 1) goto 30
n2 = n/2
i f (2*n2 .eq. n) thenmedian = 0.5* ( w s (n2)+ws (n2+1) )else
median = ws (n2+1)
endif
end

smad = sm/0.6745
end
r e a l function a15(x, n, xs, w s )r e a l median, x ( n ) , w s ( * )r e a l c, xm, xm0, xs, xs0, xc, sum
integer if ndata c/1.5/
xm = median(x, n, w s )xs = $mad(=, x, n, w s )xc = c*xs

10 Xmo = xm
sum = 0 . 0
do 20 i = 1, n
20
sum = sum + min(xmO+xc, max(xm0-xc, x ( i ) ) )x m = sum/ni f (abs(xm-xmO) . g t . (l.Oe-4)*xs) go t o 1 0
a15 = xmend
r e a l f u n c t i o n h15(x, n, xsc, w s )r e a l median, x ( n ) , w s ( * ) , xsc
r e a l a , b e t a , c , c l , xm, xm0, xs, xs0, xc, sum, sum2
i n t e g e r if n
d a t a c, b e t a /l.S, 0 . 7 7 8 /cl = cwithout small-sample c o r r e c t i o n
c l = c * sqrt(l.O-l.O/n)
xm = median(x, n, w s )xs = smad(xm, x, n, w s )C

10 xmo = xm
xso = xs
sum = 0 . 0
sum2 = 0 . 0
xc = cl*xs
do 20 i = 1, n
a = min(xm+xc, max(xm-xc, x ( i ) ) )sum = sum + a
20
sum2 = sum2(a-=)*(a-xm)
xm = sum/n
xs = s q r t (sum2/(beta* (n-1) ) )i f ((abs(xm--0).gt. l.Oe-Q*xsO) .or.
h abs(xs/xsO-1.0). ~ 1 . 0.e - 4 ) go t o 1 0
h15 = xm
xsc = xsend+

r e a l function smad(mu, x, n, w s )r e a l median, mu, x ( n ) , w s ( * )r e a l sm, sum
i n t e g e r if n
do 1 0 i = l , n
10ws (i)= abs ( x ( i )-mu)s m = median(ws, n, w s )i f (sm . l e . 0 . 0 ) thensum = 0.0

do 20 i = 1, n
20
sum = sum + w s (i)sm = sum/n
endif* Copyright of Professor B. D. Ripley

observations plus one outlier. On the other hand, (H15, 6)
regards the sample as three variable observations. We know
only from other data which explanation is correct.
N o statistical method can make sense of disastrous trials.

Abbey14 quoted 31 determinations of the nickel content(pg g-l) of Canadian syenite rock:5.2 6.5 6.9 7 7 7 7.4 8 8 8 8 8.5 9 9 10
11 11 12 12 13.7 14 14 14 16 17 17 18 24 28 34 125

1697

ANALYST, DECEMBER 1989, VOL. 114

which gives
Mean = 16.01
Median = 11.00
A15 = 11.55
H15 = 11.70

s = 21.27
6m= 4.45

6 = 5.19

suggesting a standard error of the robust mean of ca. 1.0 (=6h’rn). Abbey quoted other robust estimators of p, but all
agreed to within the (considerable) uncertainty. However,this uncertainty is so large that very little has been learnt from31 determinations.
In general, outliers are not a problem when data are lookedat carefully. Increasingly data are not looked at at all. They arerecorded in machine-readable form and summarised bycomputer programs. In such circumstances robust statistics
are preferable to conventional ones, and a marked difference
between them should give a warning that the data should be
examined carefully.

AppendixComputation
The exact form of the algorithms used to calculate robust
estimates can be deduced from the FORTRAN
77 program
shown. They cover the most general case of unknown p and 0,
but are easily modified to handle other instances.
Medians are found by sorting the data by the sort algorithmof Shell.18 There are ways to find medians without sorting thatwill be faster for large values of n , but these are considerably
more awkward to program correctly. One other difficulty isthat 6, could turn out to be zero, but only if half the data are
equal to the median. In that instance we would report 2lxi medianlh, which is zero only if all the data are equal to themedian.
The programs use a common workspace (ws) which should
be as large as the data array. One trap for the unwary: medianand H15 are real qualities, despite the implicit rules of
FORTRAN.
To aid translation to other languages, all variables
used are declared. NMAX can be set as required.
Professor Peter Rousseeuw’s comments were most helpful inclarifying an earlier draft.

References
1. Dixon, W. J . , Ann. Math. Statist., 1950, 21, 488.
2. Grubbs, F. E., Technometrics, 1969, 11, 1.
3. Crow, E. L., Davis, F. A., and Maxfield, M. W., “Statistical
Manual,” Dover, New York, 1960.
4. “Precision of Test Methods,” BS 1979, No. 5497, British
Standards Institution, Milton Keynes.
5 . Analytical Methods Committee, Analyst, 1989, 114, 1489.
6. Barnett, V . , and Lewis, J., “Outliers in Statistical Data,”
Second Edition, Wiley, Chichester, 1984.
7. Huber, P. J . , “Robust Statistical Procedures,” SIAM, Pennsylvania, 1976.
8. Huber, P. J . , “Robust Statistics,” Wiley, New York, 1981.
9. Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J., andStahel, W. A., “Robust Statistics: The Approach Based on
Influence Functions,” Wiley, New York, 1986.
10. Rousseeuw, P. J., and Leroy, A. M., “Robust Regression andOutlier Detection,” Wiley, New York, 1987.
11. Lister, B., Geostand. Newsl., 1984, 8, 171.
12. Lister, B., Geostand. Newsl., 1987, 11, 11.
13. Abbey, S . , Geostand. Newsl., 1986, 10, 159.
14. Abbey, S., Geosrand. Newsl., 1988, 12, 241.
15. Rock, N. M. S . , Math. Geol., 1988, 20,243.
16. Andrews, D. F., Bickel, P. J., Hampel, F. R., Huber, P. J.,
Rogers, W. H., and Tukey, J . W., “Robust Estimates ofLocation,” Princeton University Press, 3 972.
17. Analytical Methods Committee, Analyst, 1989, 114, 1699.
18. Sedgewick, R . , “Algorithms,” Addison Wesley, Reading, MA,
1983.

Paper 9102435B
Received June 9th, 1989