The deepest regression method.
Deepest regression (DR) is a method for linear regression introduced by P. J. Rousseeuw and M. Hubert (1999, J. Amer. Statis. Assoc. 94, 388-402). The DR method is defined as the fit with largest regression depth relative to the data. In this paper we show that DR is a robust method, with breakdown value that converges almost surely to 1/3 in any dimension. We construct an approximate algorithm for fast computation of DR in more than two dimensions. From the distribution of the regression depth we derive tests for the true unknown parameters in the linear regression model. Moreover, we construct simultaneous confidence regions based on bootstrapped estimates. We also use the maximal regression depth to construct a test for linearity versus convexity/concavity. We extend regression depth and deepest regression to more general models. We apply DR to polynomial regression and show that the deepest polynomial regression has breakdown value 1/3. Finally, DR is applied to the Michaelis-Menten model of enzyme kinetics, where it resolves a long-standing ambiguity.
Introduction
Consider a dataset Z ng ae IR p . In linear regression
we want to fit a hyperplane of the form
. We denote the x-part of each data point z i by x . The
residuals of Z n relative to the fit ' are denoted as r
To measure the quality of a fit, Rousseeuw and Hubert [16] introduced the notion of
regression depth.
Research Assistant with the FWO, Belgium
Postdoctoral Fellow at the FWO, Belgium.
Definition 1. The regression depth of a candidate fit ' 2 IR p relative to a dataset Z n ae IR p
is given by
rdepth ('; Z n
where the minimum is over all unit vectors
The regression depth of a fit ' ae IR p relative to the dataset Z n ae IR p is thus the
smallest number of observations that need to be passed when tilting ' until it becomes
vertical. Therefore, we always have 0  rdepth('; Z n )  n.
In the special case of are no x-values, and Z n is a univariate dataset. For
any ' 2 IR we then have rdepth('; Z n which is the
'rank' of ' when we rank from the outside inwards. For any p  1, the regression depth of '
measures how balanced the dataset Z n is about the linear fit determined by '. It can easily
be verified that regression depth is scale invariant, regression invariant and affine invariant
according to the definitions in Rousseeuw and Leroy ([17, page 116]).
Based on the notion of regression depth, Rousseeuw and Hubert [16] introduced the
deepest regression estimator (DR) for robust linear regression. In Section 2 we give the
definition of DR and its basic properties. We show that DR is a robust method with
breakdown value that converges almost surely to 1/3 in any dimension, when the good data
come from a large semiparametric model. Section 3 proposes the fast approximate algorithm
MEDSWEEP to compute DR in higher dimensions (p  3). Based on the distribution of
the regression depth function, inference for the parameters is derived in Section 4. Tests
and confidence regions for the true unknown parameters ' are constructed. We also
propose a test for linearity versus convexity of the dataset Z n based on the maximal depth
of Z n . Applications of deepest regression to specific models are given in Section 5. First
we consider polynomial regression, for which we update the definition of regression depth
and then compute the deepest regression accordingly. We show that the deepest polynomial
regression always has breakdown value at least 1/3. We also apply the deepest regression
to the Michaelis-Menten model, where it provides a solution to the problem of ambiguous
results obtained from the two commonly used parametrizations.
2 Definition and properties of deepest regression
Definition 2. In p dimensions the deepest regression estimator DR(Z n ) is defined as the
fit ' with maximal rdepth ('; Z n ), that is
rdepth ('; Z n
(See Rousseeuw and Hubert [16].) Since the regression depth of a fit ' can only increase
if we slightly tilt the fit until it passes through p observations (while not passing any other
observations), it suffices to consider all fits through p data points in Definition 2. If several
of these fits have the same (maximal) regression depth, we take their average. Note that no
distributional assumptions are made to define the deepest regression estimator of a dataset.
The DR is a regression, scale, and affine equivariant estimator. For a univariate dataset, the
deepest regression is its median. The DR thus generalizes the univariate median to linear
regression.
In the population case, let be a random p-dimensional variable, with distribution
H on IR p . Then rdepth('; H) is defined as the smallest amount of probability mass that
needs to be passed when tilting ' in any way until it is vertical. The deepest regression
DR(H) is the fit ' with maximal depth. The natural setting of deepest regression is a
large semiparametric model H in which the functional form is parametric and the error
distribution is nonparametric. Formally, H consists of all distributions H on IR p that satisfy
the following conditions:
H has a strictly positive density and there exists a ~
med H (y j
': (H)
Note that this model allows for skewed error distributions and heteroscedasticity. Van Aelst
and Rousseeuw [22] have shown that the DR is a Fisher-consistent estimator of med(yjx)
when the good data come from the natural semiparametric model H. The asymptotic
distribution of the deepest regression was obtained by He and Portnoy [9] in simple regression,
and by Bai and He [2] in multiple regression.

Figure

1 shows the Educational Spending data, obtained from the DASL library at
http://lib.stat.cmu.edu/DASL. This dataset lists the expenditures per pupil versus the
average salary paid to teachers for regions in the US. The fits ' 1
both have regression depth 2, and the deepest regression DR(Z n
average salary
expenditures DR(Z n )2

Figure

1: Educational spending data, with observations in dimensions. The
lines
both have depth 2, and the deepest regression DR(Z n ) is the average of fits
with depth 23.
(0:17; \Gamma0:51) t is the average of fits with depth 23. Figure 1 illustrates that lines with high
regression depth fit the data better than lines with low depth. The regression depth thus
measures the quality of a fit, which motivates our interest in the deepest regression DR(Z n ).
We define the finite-sample breakdown value "
n of an estimator T n as the smallest fraction
of contamination that can be added to any dataset Z n such that T n explodes (see also Donoho
and Gasko [6]). Let us consider an actual dataset Z n . Denote by Z n+m the dataset formed
by adding m observations to Z n . Then the breakdown value is defined as
Zn+m
The breakdown value of the deepest regression is always positive, but it can be as low
as when the original data are themselves peculiar (Rouseeuw and Hubert [16]).
Fortunately, it turns out that if the original data are drawn from the model, then the
breakdown value converges almost surely to 1=3 in any dimension p.
Theorem 1. Let Z
be a sample from a distribution H on IR p
(p
\Gamma\Gamma\Gamma!
(All proofs are given in the Appendix.) Theorem 1 says that the deepest regression does
not break down when at least 67% of the data are generated from the semiparametric model
H while the remaining data (i.e., up to 33% of the points) may be anything. This result holds
in any dimension. The DR is thus robust to leverage points as well as to vertical outliers.
Moreover, Theorem 1 illustrates that the deepest regression is different from L 1 regression,
which is defined as L 1 (Z n
(')j. Note that L 1 is another generalization of
the univariate median to regression, but with zero breakdown value due to its vulnerability
to leverage points.
In simple regression, Van Aelst and Rousseeuw [22] derived the influence function of the
DR for elliptical distributions and computed the corresponding sensitivity functions. The
influence functions of the DR slope and intercept are piecewise smooth and bounded, meaning
that an outlier cannot affect DR too much, and the corresponding sensitivity functions
show that this already holds for small sample sizes.
The deepest regression also inherits a monotone equivariance property from the univariate
median, which does not hold for L 1 or other estimators such as least squares, least trimmed
squares (Rousseeuw [14]) or S-estimators (Rousseeuw and Yohai [20]). By definition, the
regression depth only depends on the x i and the signs of the residuals. This allows for
monotone transformations of the response y Assume the functional model is
with g a strictly monotone link function. Typical examples of g include the logarithmic, the
exponential, the square root, the square and the reciprocal transformation. The regression
depth of a nonlinear fit (4) is defined as in (1) but with r i
Due to the monotone equivariance, the deepest regression fit can be obtained as follows.
First we put ~
determine the deepest linear regression "
the transformed data
Then we can backtransform the deepest linear regression "
yielding the deepest nonlinear regression fit
to the original data.
Computation
In dimensions the regression depth can be computed in O(n log n) time with the
algorithm described in (Rousseeuw and Hubert [16]). To compute the regression depth of a
fit in constructed exact algorithms
with time complexity O(n p\Gamma1 log n). For datasets with large n and/or p they also give
an approximate algorithm that computes the regression depth of a fit in O(mp 3
mn log n) time. Here m is the number of (p \Gamma 1)-subsets in x-space used in the algorithm.
The algorithm is exact when all
\Delta such subsets are considered.
A naive exact algorithm for the deepest regression computes the regression depth of all
observations and keeps the one(s) with maximal depth. This yields a
total time complexity of O(n log n) which is very slow for large n and/or high p. Even
if we use the approximate algorithm of Rousseeuw and Struyf [18] to compute the depth of
each fit, the time complexity remains very high. In simple regresssion, collaborative work
with several specialists of computational geometry yielded an exact algorithm of complexity
O(n log 2 n), i.e. little more than linear time (van Kreveld et al. [24]). To speed up the
computation in higher dimensions, we will now construct the fast algorithm MEDSWEEP
to approximate the deepest regression.
The MEDSWEEP algorithm is based on regression through the origin. For regression
through the origin, Rousseeuw and Hubert [16] defined the regression depth (denoted as
rdepth 0 ) by requiring that in Definition 1. Therefore, the rdepth 0 (') of a fit ' ae IR p
relative to a dataset Z n ae IR p+1 is again the smallest number of observations that needs
to be passed when tilting ' in any way until it becomes vertical. Rousseeuw and Hubert
[16] have shown that in the special case of a regression line through the origin (p = 1), the
deepest regression (DR of the dataset Z given by the slope
med
where observations with x are not used. This estimator has minimax bias (Martin,
Yohai and Zamar [11]) and can be computed in O(n) time.
We propose a sweeping method based on the estimator (5) to approximate the deepest
regression in higher dimensions. Suppose we have a dataset Z
ng. We arrange the n observations as rows in a n \Theta p matrix
where the X i and Y are n-dimensional column vectors.
Step 1: In the first step we construct the sweeping variables X S
. We start
with
To obtain X S
(j ? 1) we successively sweep X S
out of the original
variable X j . In general, to sweep X S
k out of X l (k ! l) we compute
med
med
x il
med
ik
where J is the collection of indices for which the denominator is different from zero, and
then we replace X l by
now sweep the next variable X S
out of this new X l . If
. Thus we obtain the sweeping variables
Step 2: In the second step we successively sweep X S
out of Y . Put Y . For
med
y
med
y
med
ik
with J as before, and we replace the original Y
k . Thus we obtain
The proces (7)-(8) is iterated until convergence is reached. In each iteration step all the
coefficients "
fi are updated. The maximal number of iterations has been set to
100 because experiments have shown that even when convergence is slow, after 100 iteration
steps we are already close to the optimum. After the iteration proces, we take the median
of Y S to be the intercept "
Step 3: By backtransforming "
we obtain the regression coefficients
corresponding to the original variables . The obtained fit " '
is then slightly
adjusted until it passes through p observations, because we know that this can only improve
the depth of the fit. We start by making the smallest absolute residual zero. Then for each
of the directions we tilt the fit in that direction until it passes an observation
while not changing the sign of any other residual. This yields the fit "
'.
Step 4: In the last step we approximate the depth of the final fit "
'. Let u S
be the directions corresponding to the variables X S
, then we compute the minimum
over
instead of over all unit
vectors in the right hand side of expression (1).
Since computing the median takes O(n) time, the first step of the algorithm needs O(p 2 n)
time and the second step takes O(hpn) time where h is the number of iterations. The adjustments
in step 3 also take O(p 2 n) time, and computing the approximate depth in the last
step can be done in O(pn log n) time. The time complexity of the MEDSWEEP algorithm
thus becomes O(p log n) which is very low.
To measure the performance of our algorithm we carried out the following simulation. For
different values of p and n we generated
ng from the standard gaussian distribution. For each of these samples we computed
the deepest regression
with the MEDSWEEP algorithm and measured the
total time needed for these 10; 000 estimates. For each n and p we also computed the bias
of the intercept, which is the average of the 10; 000 intercepts, and the bias of the vector of
the slopes, which we measure by
((ave
(ave
We also give the mean squared error of the vector of the slopes, given by
where the true values and the mean squared error of the
intercept, given by 1

Table

1 lists the bias and mean squared error of the
vector of the slopes, while the bias and mean squared error of the intercept are given in

Table

2. Note that the bias and mean squared error of the slope vector and the intercept
are low, and that they decrease with increasing n. From Tables 1 and 2 we also see that the
mean squared error does not seem to increase with p.

Table

1: Bias (9) and mean squared error (10) of the DR slope vector, obtained by generating
standard gaussian samples for each n and p. The DR fits were obtained with
the MEDSWEEP algorithm. The results for the bias have to be multiplied by 10 \Gamma4 and the
results for the MSE by 10 \Gamma3 .
50 bias 18.01 18.27 28.53 22.47 21.42
MSE 54.55 52.32 52.23 52.54 58.65
100 bias 7.56 8.41 11.16 8.98 12.68
MSE 24.59 25.32 25.99 26.15 27.17
300 bias 8.11 5.98 6.29 7.39 10.89
MSE 8.11 8.27 8.26 8.41 8.42
500 bias 0.44 1.68 2.91 9.10 5.49
MSE 4.93 4.89 5.02 4.93 4.97
1000 bias 6.34 3.58 6.77 3.54 3.98
MSE 2.47 2.51 2.46 2.46 2.50

Table

3 lists the average time the MEDSWEEP algorithm needed for the computation
of the DR, on a Sun SparcStation 20/514. We see that the algorithm is very fast.
To illustrate the MEDSWEEP algorithm we generated 50 points in 5 dimensions, according
to the model e coming from
the standard gaussian distribution. The DR fit obtained with MEDSWEEP is
21. The algorithm needed 14 iterations
till convergence. In a second example, we generated 50 points according to the model
iterations
the MEDSWEEP algorithm yielded the fit
with approximate depth 21. Note that in both cases the coefficients obtained by the algorithm
approximate the true parameters in the model very well. The MEDSWEEP algorithm
is available from our website http://win-www.uia.ac.be/u/statis/ where its use is explained


Table

2: Bias and mean squared error of the DR intercepts, obtained by generating
10; 000 standard gaussian samples for each n and p. The DR fits were obtained with the
MEDSWEEP algorithm. The results for the bias have to be multiplied by 10 \Gamma4 and the
results for the MSE by 10 \Gamma3 .
50 bias -3.04 -48.70 14.38 13.23 -19.64
MSE
100 bias 5.75 -3.32 12.21 9.92 -1.70
MSE 15.78 16.01 16.92 16.77 18.14
300 bias -2.71 -7.99 -2.37 3.49 4.54
MSE 5.25 5.31 5.22 5.23 5.47
500 bias 1.35 -5.53 -4.33 -4.26 2.39
MSE 3.09 3.15 3.21 3.21 3.18
1000 bias -3.76 -3.40 -5.10 0.75 -0.01
MSE 1.56 1.56 1.55 1.58 1.62
Inference
4.1 Tests for parameters
In simple regression, the semiparametric model assumptions of condition (H) state that
and that the errors e
are independent with P
it is possible to compute F n (k) :=
n has the same fx as the actual dataset Z n . By
invariance properties,
where the e i are i.i.d. from (say) the standard gaussian. Thus we can compute F n (k) by
simulating (11). When there are no ties among the x i we can even compute F n (k) explicitly
from formula (4.4) in Daniels (1954), yielding

Table

3: Computation time (in seconds) of the MEDSWEEP algorithm for a sample of size
n with p dimensions. Each time is an average over 10; 000 samples.
1000
for k  [(n \Gamma 1)=2], and F n and each term is a
probability of the binomial distribution B(n; 1=2), which stems from the number of e i in
with a particular sign. For increasing n we can approximate B(n; 1=2) by a
gaussian distribution due to the central limit theorem, so (12) can easily be extended to
large n.
The distribution of the regression depth allows us to test one or several regression
coefficients. To test the combined null hypothesis ( ~
and the corresponding p-value equals F n (k) and can be computed from (12).
Consider the dataset in Figure 2 about species of animals (Van den Bergh [23]). The
plot shows the logarithm of the weight of a newborn versus the logarithm of the weight of
its mother. The deepest regression line has depth 19. For this dataset
yielding the p-value F 41 which is highly significant.
To test the significance of the slope
. This is easy, because we only have to compute the rdepth of all horizontal lines
passing through an observation. For the animal dataset, the maximal rdepth((0;
equals 5. Therefore the corresponding p-value is P (rdepth( ~
so H 0 is rejected. This p-value 0:00002 should be interpreted in the same way as the p-value
associated with R 2 or the F-test in LS regression. Analogously, to test ~
by considering all lines through the origin and an observation,
yielding the p-value F 41 which is also highly significant.
More generally, we can test the hypothesis H
by computing
the maximal regression depth of all lines with ' that pass through an observation.
log(weight mother)
log(weight
newborns)
DR

Figure

2: Logarithm of the weight of a newborn versus the logarithm of the weight of its
mother for species of animals, with the DR line which has depth 19.
For example, to test the hypothesis H
for the animal data (i.e., the hypothesis
that the weight of a newborn is proportional to the weight of its mother) we compute
and the corresponding p-value F 41 which is
significant at the 5% level but not at the 1% level.
These tests generalize easily to higher dimensions and situations with ties among the x i ,
but then we can no longer use the exact formula (12) which is restricted to the bivariate case
without ties in fx g. Therefore, in these cases we compute F n (k) by simulating (11).
Let us consider the stock return dataset (Benderly and Zwick [3]) shown in Figure 3 with
28 observations in dimensions. The regressors are output growth and inflation
(both as percentages), and the response is the real return on stocks. The deepest regression
obtained with the MEDSWEEP algorithm equals
depth 11. To test the null hypothesis ( ~
that both slopes are zero (this
would be done with the R 2 in LS regression) we compute the maximal rdepth((0; 0; ' 3
over all ' 3 (i.e. over all y i in the dataset). By computing the exact rdepth of these 28 horizontal
planes (by the fast algorithm of Rousseeuw and Struyf [18]) we obtain the value 6.
Simulation yields the corresponding p-value F 28 which is not significant. To test
the significance of the intercept
That is, we compute the depth of all planes through two observations and the ori-
gin. For the example this yields max corresponding p-value
F 28 (11)  1, which is not at all significant.
4.2 Confidence regions
In order to construct a confidence region for the unknown true parameter vector ~
use a bootstrap method. Starting from the dataset Z
ng 2 IR p , we construct a bootstrap sample by randomly drawing n observations, with
replacement. For each bootstrap sample Z (j)
we compute its deepest regression
. Note that there will usually be a few outlying estimates "
in the set of
which is natural since some bootstrap samples contain
disproportionally many outliers. Therefore we don't construct a confidence ellipsoid based
on the classical mean and covariance matrix of the f "
but we use the robust
minimum covariance determinant estimator (MCD) proposed by (Rousseeuw [14,15]).
The MCD looks for the h  n=2 observations of which the empirical covariance matrix
has the smallest possible determinant. Then the center "
of the dataset is
defined as the average of these h points, and the covariance matrix "
\Sigma of the dataset is a
certain multiple of their covariance matrix. To obtain a confidence ellipsoid of level ff we
compute the MCD of the set of bootstrap estimates with ff)me. The
confidence ellipsoid E 1\Gammaff is then given by
Here
is the statistic of the robust distances of the
where the robust distance (Rousseeuw and Leroy
[17]) of a bootstrap estimate "
is given by
From this confidence ellipsoid E 1\Gammaff in fit space we can also derive the corresponding
regression confidence region for "
defined as
Theorem 2. The region R 1\Gammaff equals the set
Let us consider the Educational Spending data of Figure 1. Figure 4a shows the deepest
regression estimates of 1000 bootstrap samples, drawn with replacement from the original
data. Using the fast MCD algorithm of Rousseeuw and Van Driessen [19] we find the center
in Figure 4a to be (0:19; \Gamma0:95) t which corresponds well to the DR fit
the original data. As a confidence region for ( ~
we take the 95% tolerance ellipse E 0:95
based on the MCD center and scatter matrix, which yields the corresponding confidence
region R 0:95 shown in Figure 4b. Note that the intersection of this confidence region with
a vertical line is not a 95% probability interval for an observation y at x 0 . It is the
interval spanned by the fitted values "
in a 95% confidence region
for ( ~
An example of a confidence region in higher dimensions is shown in Figure 3. It shows
the 3-dimensional stock return dataset with its deepest regression plane, obtained with
the MEDSWEEP algorithm. The 95% confidence region shown in Figure 3 was based on
samples.
4.3 Test for linearity in simple regression
If the observations of the bivariate dataset Z n lie exactly on a straight line, then
is the highest possible value. On the other hand, if Z n lies exactly
on a strictly convex or strictly concave curve, then
n=3 is at its lowest
(Rousseeuw and Hubert [16, Theorem 2]). Therefore,
can be seen as a
measure of linearity for the dataset Z n . Note that this lower bound does not depend on the
amount of curvature when the lie exactly on the curve. However, as soon as there
is noise (i.e. nearly always), the relative sizes of the error scale and the curvature come into
play.
The null hypothesis assumes that the dataset Z n follows the linear model
for some ~
and with independent errors e i each having a distribution with zero
median. To determine the corresponding p-value we generate
DR
output growth
inflation
stock
return

Figure

3: Stock return dataset with its deepest regression plane, and the upper and lower
surface of the 95% confidence region R 0:95 based on samples.

Table

4: Windmill data: Maximal rdepth k with corresponding p-value p(k) obtained by
simulation.
ng with the same x-values as in the dataset Z n and with standard gaussian
e i . For each j we compute the maximal regression depth of the dataset Z (j) . Then for each
value k we approximate the p-value P(maxrdepth  kjH 0 ) by
For example, let us consider the windmill dataset (Hand et al. [8]) which consists of 25
measures of wind velocity with corresponding direct current output, as shown in Figure 5.
The p-values in Table 4 were obtained from (17). For the actual
so we reject the linearity at the 5% level.
(a)
(b)
average salary
expenditures DR
R 0.95

Figure

4: (a) DR estimates of 1; 000 bootstrap samples from the Educational Spending data
with the 95% confidence ellipse E 0:95 ; (b) Plot of the data with the DR line and the 95%
confidence region R 0:95 for the fitted value.
5 Specific models
5.1 Polynomial regression
Consider a dataset Z ng ae IR 2 . The polynomial regression model
wants to fit the data by is called the degree of
the polynomial. The residuals of Z n relative to the fit are denoted as
We could consider this to be a multiple regression problem with regressors
and determine the depth of a fit ' as in Definition 1. But we know that the joint distribution
of degenerate (i.e. it does not have a density), so many properties of the
deepest regression, such as Theorem 1 about the breakdown value, would not hold in this
case. In fact, the set of possible x forms a so-called moment curve in IR k , so it is inherently
one-dimensional.
A better way to define the depth of a polynomial fit (denoted by rdepth k ) is to update
the definition of regression depth in simple regression, where x was also univariate. For each
polynomial fit ' we can compute its residuals r i (') and define its regression depth as
rdepth k
where . Note that this definition
only depends on the x i -values and the signs of the residuals. With this definition of depth,
we define the deepest polynomial of degree k as in (2) and denote it by DR k (Z n ).
Theorem 3. For any dataset Z ng 2 IR 2 with distinct x i the deepest
polynomial regression of degree k has breakdown value
Theorem 3 shows that with the definition of depth of a polynomial fit given in (18), the
deepest polynomial regression has a positive breakdown value of approximately 1/3, so it is
robust to vertical outliers as well as to leverage points.
In Section 4.3 we rejected the linearity of the windmill data. Therefore we now consider
a model with a quadratic component for this data, i.e. Figure 5 shows
the windmill data with the deepest quadratic fit, which has depth 12.
wind velocity
direct
current
output

Figure

5: Windmill data with the deepest quadratic fit, which has regression depth 12.
5.2 Michaelis-Menten model
In the field of enzyme kinetics, the steady-state kinetics of the great majority of the enzyme-
catalyzed reactions that have been studied are adequately described by a hyperbolic relationship
between the concentration s of a substrate and the steady-state velocity v. This
relationship is expressed by the Michaelis-Menten equation
where the constant v max is the maximum velocity and Km is the Michaelis constant. The
Michaelis-Menten equation has been linearized by rewriting it in the following three ways:
s
Km
\GammaKm
s
which are known as the Scatchard equation (Scatchard [21]), the double reciprocal equation
(Lineweaver and Burke [10]) and the Woolf equation (Haldane [7]). Each of the three relations
(21), (22), (23) can be used to estimate the constants v max and Km . In general the three
relations yield different estimates for the constants v max and Km , because the error terms
are also transformed in a nonlinear way. Cressie and Keightley [5] compared these three
linearizations of the Michaelis-Menten relation in the context of hormone-receptor assays,
and concluded that for well-behaved data the Woolf equation (23) works best, but for data
containing outliers the double reciprocal equation (22) with robust regression gives better
results.
Theorem 4 shows that applying the deepest regression to the Woolf equation (23) yields
the same estimates for v max and Km as the deepest regression applied to the double reciprocal
equation (22). This resolves the ambiguity.
Theorem 4. Let Z ng
denote the DR fit of the double reciprocal equation as DR(f( 1
Then the DR fit of the Woolf equation satisfies
In both cases we obtain the same " v
Example: In assays for hormone receptors the Michaelis-Menten equation describes
the relationship between the amount B of hormone bound to receptor and the amount F of
hormone not bound to receptor. These assays are used e.g. to determine the cancer treatment
method (see Cressie and Keightley [5]). Equation (20) now becomes
The parameters of interest are the concentration B max of binding sites and the dissociation
constant KD for the hormone-receptor interaction. Figure 6a shows the Woolf plot
of data from an estrogen receptor assay obtained by Cressie and Keightley [4]. Note that
this dataset clearly contains one outlier. To this plot we added the deepest regression line
In Figure 6b we also show the double reciprocal
plot with its deepest regression line DR(f( 1
in both cases we obtain the same estimated values "
are comparable to the least squares estimates "
obtained
from the Woolf equation based on all data except the outlier. On the other hand, least
squares applied to the full data gives "
based on the Woolf
(a)
50 100 150 200 250 300246F
DR
(b)
DR

Figure

(a) Woolf plot of the Cressie and Keightly data with the deepest regression line
and the least squares fit double reciprocal
plot of the data with the deepest regression line and the least squares
equation, and "
based on the double reciprocal equation.
These estimates are quite different. With least squares, in both cases the estimates "
KD are highly influenced by the outlying observation (both "
KD come out
too high) which may lead to wrong conclusions, e.g. when determining a cancer treatment
method.


Appendix


Proof of Theorem 1. To prove Theorem 1 we need the following lemmas.
Lemma 1. If the x i are in general position, then f'; rdepth('; Z n )  pg is bounded.
Proof: For J ae ng with we denote the fit that passes through the p
observations
Since the x i are in general position,
such a fit ' J will be non-vertical (for any J). Therefore, convf' J ; J ae
stands for the convex hull.) 2
Lemma 2. For any dataset Z n 2 IR p with the x i in general position (i.e. no more than
of the x i lie in any (p \Gamma 2)-dimensional affine subspace of IR p\Gamma1 ) the deepest regression has
breakdown value
Proof: By Lemma 1 we know that f'; rdepth('; Z n )  pg is bounded. Therefore to break
down the estimator we must add at least m observations such that rdepth(DR(Z n+m ); Z n )
it holds that rdepth(DR(Z n ); Z n )  dn=(p
Amenta et al. [1]), we obtain
from which it follows that m  n\Gammap 2 +1
. Finally, it holds that
Lemma 3. If Z n ae Z n+m then
Proof. Since Z n ae Z n+m it holds that rdepth('; Z n )  rdepth('; Z n+m ) for all '. Thus for
all ' it holds that rdepth('; Z n )  max
Lemma 4. Under the conditions of Theorem 1 we have that
a:s:
\Gamma\Gamma\Gamma!
Proof. Let us consider the dual space, i.e. the p-dimensional space of all possible fits '.
Dualization transforms a hyperplane H '
to the point '. An
observation z mapped to the set D(z i ) of all ' that pass through z i ,
so D(z i ) is the hyperplane H i given by ' In dual space, the
regression depth of a fit ' corresponds to the minimal number of hyperplanes H i intersected
by any halfline ['; '
A unit vector u in dual space thus corresponds to an affine hyperplane V in x-space
and a direction in which to tilt ' until it is vertical. For each unit vector u, we therefore
define the wedge-shaped region A ';u
in primal space, formed by tilting ' around V (in the
direction of u) until the fit becomes vertical. Further denote H n the empirical distribution
of the observed data Z n . Define the metric
)j:
If Z n is sampled from H, then it holds that
\Gamma\Gamma\Gamma!
0:
This follows from the generalization of the Glivenko-Cantelli theorem formulated in (Pollard
[13, Theorem 14]) and proved by (Pollard [13, Lemma 18]) and the fact that A ';u
can
be constructed by taking a finite number of unions and intersections of half-spaces. Now
define
its empirical version \Pi n
). It is clear that
1. Moreover we have that j\Pi n ( ~
hence
\Gamma\Gamma\Gamma!
Finally it holds that
a:s:
The latter inequality uses Theorem 7 in (Rousseeuw and Hubert [16]) which is valid since
Z n is almost surely in general position. Taking limits in (27) and using (26) then finishes
the proof. 2
Proof of Theorem 1: We will first show that lim inf
1almost surely. From Lemma 1
we know that f'; rdepth('; Z n )  pg is bounded a.s. if Z n ae IR p is sampled from H.
Therefore, to break down the estimator we must add at least m observations such that
This implies that
where the first inequality is due to Lemma 3, and thus
a:s:
Finally we apply Lemma 4 to conclude that
lim inf
a:s:
Next, we prove that lim sup
1almost surely. Since regression depth is affine
invariant, we may assume that none of the jjx i jj in the dataset are zero. We will denote a
hyperplane ' with equation by its two components fi ' and ff ' corresponding
to the slopes and the intercept. Fix two strictly positive real numbers fi 0 and ff 0 . We now
consider a point such that for any hyperplane ' passing through
data points
it holds that 1. Because we
assumed that none of the jjx i jj equals 0, we can always find such a point
x would be unbounded.
The dataset Z n+m is then obtained by enlarging the dataset Z n with
points equal to We know that the deepest fit DR(Z n+m ) must pass through at least p
different observations of Z n+m . Denote by ' J any candidate deepest fit. If ' J passes through
it is clear that rdepth(' J ; Z n+m )  On the other hand, any ' J
which passes through p data points of Z n has rdepth(' J ; Z n+m This can be seen
as follows. First consider the dataset Z n+1 which consists of the n original observations and
one copy of the point Theorem 7 in (Rousseeuw and Hubert [16]) it follows that
\Theta n+1+p
. Now there always exists a unit vector
such that x t
n. Then the number of observations passed
when tilting ' J around (u; v) as in Definition 1 plus the number of observations passed
when tilting ' J around (\Gammau; \Gammav) equals because the fit ' J passes through exactly
observations. Therefore, we can suppose that the number of observations passed when
tilting ' J around (u; v) is at most [ n+1+p]. (If not, we replace (u; v) by (\Gammau; \Gammav).) Note
that the residual r does not pass through First suppose
the residual r strictly positive. The data are in general position, therefore we
can always find a value " ? 0 such that for all it holds that x t
number of observations passed when tilting ' J around (u; is the same as when tilting
around (u; v). Finally, adding the other replications of does not change this
value. Therefore rdepth(' J ; Z n+m If the residual r strictly negative,
we replace v by in a similar way and obtain the same result.
The above reasoning shows that the deepest fit must pass through
original data points. Since we have shown that all these fits have an arbitrarily large slope
and intercept, it holds that
a:s:
\Gamma\Gamma\Gamma!
and thus
lim sup
a:s:
From (28) and (29) we finally conclude (3). 2
Proof of Theorem 2. Consider x 2 IR p\Gamma1 . We will prove that the upper and lower bounds
in expression (16) are the values y such that in the dual plot the hyperplane
is tangent to the ellipsoid
First consider the special case "
yielding the unit sphere ' t
The tangent hyperplane in an arbitrary point ~
on the sphere is given by
~
Therefore the hyperplane becomes a tangent hyperplane
' on the sphere. Together with ~
yields giving the lower bound
and the upper bound
corresponding to expression (16) for this case.
Consider the general case of an ellipsoid
is the diagonal matrix of eigenvalues of "
is the matrix of eigenvectors of "
\Sigma. We can transform this to the previous
case by the transformation ~
). The hyperplane
transforms to c(x t ; 1)P  1=2 ~
which becomes a tangent hyperplane if
This yields the lower bound y =
and the upper bound y =
of expression (16). 2
Proof of Theorem 3. From the definition of rdepth k by (18) it follows for any data set
ng with distinct x i as in Lemma 1 that f'; rdepth k ('; Z n )
is bounded. Now any bivariate linear fit
corresponds to a polynomial fit
holds for the deepest regression line DR 1 that rdepth(DR 1
(Rousseeuw and Hubert [16]), it follows that the deepest polynomial regression DR k of degree
k has rdepth k (DR
Because must add at least m observations
such that rdepth k (DR k (Z n+m ); Z n )  k to break down the estimator. We obtain
from which it follows that m  (n \Gamma 3k)=2. This yields
Proof of Theorem 4. We will show that when s holds for every
This follows from
(v i
for all switching the x-values from 1
to s i reverses their order. Therefore,
according to Definition 1 both depths are the same. 2



--R

"Regression Depth and Center Points,"
Asymptotic distributions of the maximal depth estimators for regression and multivariate location

The underlying structure of the direct linear plot with application to the analysis of hormone-receptor interactions
Analysing data from hormone-receptor assays
Breakdown properties of location estimates based on halfspace depth and projected outlyingness
Graphical methods in enzyme chemistry
A Handbook of Small Data Sets
"Applied Statistical Science III: Nonparametric Statistics and Related Topics"
The determination of enzyme dissociation constants

"On Depth and Deep Points: a Calculus,"
Convergence of Stochastic Processes
Least median of squares regression
"Mathematical Statistics and Applications, Vol. B"

New York: Wiley-Interscience
Computing location depth and regression depth in higher dimensions
A fast algorithm for the minimum covariance determinant estimator
"Robust and Nonlinear Time Series Analysis,"
The attractions of proteins for small molecules and ions
Robustness of deepest regression

"Proceedings of the 15th Symposium on Computational Geometry,"
--TR
Robust regression and outlier detection
Efficient algorithms for maximum regression depth
A fast algorithm for the minimum covariance determinant estimator
An optimal algorithm for hyperplane depth in the plane
Robustness of deepest regression
Computing location depth and regression depth in higher dimensions

--CTR
R. Wellmann , S. Katina , Ch. H. Mller, Calculation of simplicial depth estimators for polynomial regression with applications, Computational Statistics & Data Analysis, v.51 n.10, p.5025-5040, June, 2007
Christine H. Mller, Depth estimators and tests based on the likelihood principle with application to regression, Journal of Multivariate Analysis, v.95 n.1, p.153-181, July 2005
algorithm;regression depth;inference