Valid Inference Corrected for Outlier Removal
Shuxiao Chen∗, Jacob Bien†

arXiv:1711.10635v2 [stat.ME] 30 Nov 2017

November 29, 2017

Abstract
Ordinary least square (OLS) estimation of a linear regression model is well-known to be highly
sensitive to outliers. It is common practice to first identify and remove outliers by looking at the data
then to fit OLS and form confidence intervals and p-values on the remaining data as if this were the
original data collected. We show in this paper that this “detect-and-forget” approach can lead to invalid
inference, and we propose a framework that properly accounts for outlier detection and removal to
provide valid confidence intervals and hypothesis tests. Our inferential procedures apply to any outlier
removal procedure that can be characterized by a set of quadratic constraints on the response vector,
and we show that several of the most commonly used outlier detection procedures are of this form.
Our methodology is built upon recent advances in selective inference (Taylor & Tibshirani 2015), which
are focused on inference corrected for variable selection. We conduct simulations to corroborate the
theoretical results, and we apply our method to two classic data sets considered in the outlier detection
literature to illustrate how our inferential results can differ from the traditional detect-and-forget strategy.
A companion R package, outference, implements these new procedures with an interface that matches
the functions commonly used for inference with lm in R.

1

Introduction

Linear regression is routinely used in just about every field of science. In introductory statistics courses,
students are shown cautionary examples of how even a single outlier can wreak havoc in ordinary least
squares (OLS). Outliers can arise for a variety of reasons, including recording errors and the occurrence of
rare phenomena, and they often go unnoticed without careful inspection (see, e.g., Belsley et al. 2005). Given
this reality, one route taken by statisticians is robust regression, in which the least-squares loss function is
replaced with other functions that are less sensitive to outliers, such as the Huber loss (Huber & Ronchetti
1981) or the least absolute deviation loss (see, e.g., Bloomfield & Steiger 1984). Yet in practice OLS remains
the predominant approach, presumably due to its accompanying inferential procedures, which are elegant
and easy to use. Thus, the most common approach is a two-step procedure, which we will refer to as
detect-and-forget:
1. detect and then remove outliers;
2. fit OLS and perform inference on the remaining data as if this were the original data set.
While this common approach is attractive for its simplicity, we will show that it can lead to confidence intervals and hypothesis tests with incorrect operating characteristics. In particular, this detect-and-forget approach is problematic for its use of the same data twice. While the term “outlier removal” might lead one to
think of Step 1 as a clear-cut, essentially deterministic step, in fact Step 1 should instead be thought of as
“potential outlier removal,” an imperfect process in which one has some probability of removing non-outliers,
a process that can alter the distribution of the data. The act of searching for and removing potential outliers must be considered as part of the data-fitting procedure and thus must be considered in Step 2 when
inference is being performed.
∗ Department
† Data

of Biological Statistics and Computational Biology, Cornell University
Sciences and Operations, Marshall School of Business, University of Southern California

1

Figure 1: Confidence intervals for the regression surface. Normal data are in black while the only outlier is
marked in red. The point with a yellow cross is the detected outlier which has the largest Cook’s distance.
The black line is the regression line fitted using the data in which the detected outlier is removed.
Similar concerns over “double dipping” are well-known in prediction problems, in which sample splitting
(into training and testing sets) is a common remedy. However, such a strategy does not translate in an obvious
way to the outlier problem: If one removes outliers on a subset of observations and performs inference on the
remaining observations, then one is of course left vulnerable to outliers in the second set that could throw
off the inference stage.
To illustrate how the detect-and-forget strategy can be problematic, consider the situation shown in Figure
1, in which there are 19 “normal” points (in black), and a single “outlier” point (in red) has been shifted
upward by different magnitudes. For this illustration, we use a well-known approach for outlier detection
called Cook’s distance (Cook 1977):
εp2
hii
Di “ i 2
,
(1)
σ p1 ´ hii q2
pp
where εpi is the i-th residual from OLS on the entire data set, σ
p2 :“ }p
ε}22 {pn ´ pq is the scaled sum of
squares, and hii is the i-th diagonal entry of the hat matrix XpX T Xq´1 X T . We declare the observation
with the largest Cook’s distance to be the outlier (indicated in the figure by a yellow cross) and then refit
the regression model with this point removed (black regression line). We then construct confidence intervals
for the regression surface in two different ways: first, using the traditional detect-and-forget strategy, which
ignores the outlier removal step, and second using corrected-exact, a method we will introduce in this paper,
which properly corrects for the removal. When the outlier is obvious (leftmost panel), our method makes no
discernible correction. With such a pronounced separation between the outlier and non-outliers, Step 1 is
unlikely to have removed a non-outlier, and thus the distribution of the data for inference is likely unaltered.
However, when the outlier is less easily distinguished from the data, our corrected confidence intervals are
noticeably different from the classical ones. In particular, the corrected-exact intervals are pulled in the
direction of the removed data point, thereby accounting for the possibility that the removed point may not
in fact have been an outlier.
While Figure 1 shows only a single realization of the two intervals in four different scenarios, Figure 2
shows the empirical coverage probability, averaged over 2000 realizations, of these two types of confidence
intervals along the regression surface for the same four scenarios. We see that when the outlier signal is strong
(leftmost panel), both detect-and-forget and corrected intervals achieve 95% coverage, as desired. However, as
the outlier signal decreases, the detect-and-forget intervals begin to break down, while our corrected intervals
remain unaffected. Indeed, we will show in this paper how all sorts of inferential statements (confidence
intervals for regression coefficients, coefficient t-tests, F -tests, etc.) can be thrown off using a detect-andforget strategy but can be corrected with a proper accounting for the outlier detection and removal step.
The machinery underlying our methodology is built on recent advances in selective inference (Fithian
et al. 2014, Taylor & Tibshirani 2015), specifically the framework introduced in Lee et al. (2016), Loftus &
Taylor (2015), and it fits within the framework of inferactive data analysis introduced by Bi et al. (2017).
2

Figure 2: Empirical coverage probability along the regression surface (across 2000 realizations). The dashed
line represents 95% coverage.
We give a brief introduction to the philosophy of selective inference in the context of outlier detection and
refer readers to Fithian et al. (2014) for more details.
We assume a standard regression setting, py, Xq P Rn ˆ Rnˆp with y „ N pµ, σ 2 In q, where µi “ xTi β ˚ for
i P M ˚ . Here M ˚ is the set of non-outliers. If M is the set of detected non-outliers, then the detect-and`
forget strategy forms the OLS estimator on the subset of observations in M , β̂ M “ XM,
¨ yM (where yM and
`
XM,¨ are formed by taking rows indexed by M , and XM, is the Moore-Penrose pseudoinverse of XM,¨ ),
¨
and then proceeds with inference assuming that
T
´1
β̂ M „ N pβ ˚ , σ 2 rXM,
¨ XM,¨ s q.

However, the above assumes that M is non-random (or at least independent of y) and that M Ď M ˚ , i.e. all
true outliers have been successfully removed. However, in practice the set of declared non-outliers is in fact
xpy, Xq, and thus to perform inference would in principle require an understanding
a function of the data, M
of the distribution of the much more complicated random variable
β̂ M py,Xq “ X `
x
x

¨ yMxpy,Xq .

M py,Xq,

xpy, Xq, such as “make plots and inspect by eye,” the above distribuFor general outlier removal procedures M
tion may be completely unobtainable. However, in this paper we define a class of outlier removal procedures
for which the conditional distribution
x
xpy, Xq
β̂ M py,Xq | M
can be precisely characterized. Access to this conditional distribution will allow us to construct confidence
intervals and p-values that are valid conditional on the set of outliers selected. For example, we will produce
a procedure for forming outlier-removal-aware confidence intervals CjM py, Xq such that
xpy, Xq “ M q ě 1 ´ α
Ppβj˚ P CjM py, Xq | M
xpy, Xq Ď M ˚ (i.e., the
for all subsets M that do not include a true outlier. If one could be certain that M
procedure is adjusted to be sufficiently conservative and outliers are known to be sufficiently large), then
such conditional coverage statements can be translated into a marginal (i.e., traditional) coverage statement:
xpy,Xq
M

Ppβj˚ P Cj

py, Xqq ě 1 ´ α.

xpy, Xq Ę M ˚ ,
However, in practice we do not know if all true outliers have been successfully removed. If M
˚
then OLS is no longer guaranteed to produce an unbiased estimate of β . OLS performed on the observations

3

in M instead estimates a parameter β M , which depends on both β ˚ and on µM zM ˚ , the mean of the true
outliers that were not detected:
`
µM .
β M :“ arg min Er}yM ´ XM,¨ β}22 s “ XM,

¨

βPRp

(2)

The goal of this paper is not to improve the performance of outlier removal procedures—certainly there is
already extensive work in the literature on outlier removal. Rather, our goal is to provide valid inferential
xpy, Xq. Thus, to
statements for someone who has chosen to use a particular outlier removal procedure, M
x
stay within the scope of this problem, we will simply acknowledge that if a procedure M py, Xq is prone to
failing to identify outliers, then one cannot hope to estimate β ˚ but must instead focus on estimating and
x
performing inference for β M py,Xq , which reflects more accurately than β ˚ the relationship between X and
xpy, Xq. For example, we will provide intervals with guaranteed
y in the data that is provided to us by M
xpy,Xq
M
coverage of β
:
xpy,Xq
M

Ppβj

xpy,Xq
M

P Cj

py, Xqq ě 1 ´ α.

We will likewise provide all the standard confidence intervals and hypothesis tests for regression but focused
x
on β M py,Xq in place of β ˚ .
The rest of the paper is organized as follows: in Section 2 we formulate the problem more precisely and
describe the class of outlier detection procedures over which our framework applies; Section 3 describes our
methodology for forming confidence intervals and extracting p-values that are properly corrected for outlier
removal; Section 4 provides empirical comparisons of the naive detect-and-forget strategy and our method,
both through comprehensive simulations and a re-analysis of two data sets previously studied in the outlier
removal literature; Section 5 gives a discussion and possible next steps.
We conclude this section by introducing some notation that will be used throughout this paper. For
n P N, we let rns :“ t1, 2, . . . , nu. For a matrix X, we let C pXq be its column space and trpXq be its trace.
We let XI,J be the submatrix formed by rows and columns indexed by I and J, respectively, and we let XI,¨
be the submatrix formed by rows indexed by I. We denote X ` as its Moore-Penrose pseudoinverse. We let
K
:“ I ´ PX . For a submatrix XI,J , we write PI,J :“ PXI,J
PX be the projection matrix onto C pXq and PX
K to denote statistical independence.
when there is no ambiguity. We use K

2
2.1

Problem Formulation
The General Setup

We elaborate on the framework described in the previous section, introducing some additional notation. We
assume y “ µ ` ε, where y P Rn , ε „ N p0, σ 2 In q and consider the mean-shift model,
µ “ Xβ ˚ ` u˚ ,

(3)

where β ˚ P Rp , and X is a non-random matrix of predictors. The set M ˚ “ rnszsupppu˚ q is the index
set of true non-outliers; equivalently, M ˚c is the index set of true outliers. By definition of M ˚ , u˚M ˚ “ 0
x : Rn Ñ 2rns , as a
and u˚i ‰ 0 for i P M ˚c . We denote a data-dependent outlier removal procedure, M
function mapping the data y to the index set of detected non-outliers (for notational ease, we suppress the
x on X since X is treated as non-random). We will assume throughout that X x
has
dependence of M
M pyq,¨
linearly independent columns.
For a fixed subset of the observations M Ď rns, the parameter β M defined in (2) represents the best linear
approximation of µM using the p predictors in XM,¨ . In what follows, we will provide hypothesis tests and
xpyq “ M u.
confidence intervals for β M conditional on the event tM
Combining (2) and (3) with the assumption that XM,¨ has linearly independent columns,
`
β M “ XM,
pXM,¨ β ˚ ` u˚M q

¨

`
“ β ˚ ` XM,
u˚M .

¨

4

Since uM ˚ “ 0, it follows that β M “ β ˚ when M Ď M ˚ . This result makes it clear that if one wishes to
x is screening out all outliers.
make statements about β ˚ , then one must ensure that the procedure M
M
xpyq “ M u. Importantly,
Our focus will be on performing inference on β conditional on the event tM
such inferential procedures in fact provide asymptotically valid inferences for β ˚ as long as one’s outlier
removal procedure asympotically detects all outliers. For example, the next proposition establishes that
xpyq “ M u do in fact achieve traditional
confidence intervals providing conditional coverage of βjM given tM
˚
(i.e., unconditional) coverage of βj asymptotically if one is using an outlier detection procedure that is
guaranteed to screen out all outliers as n Ñ 8.
Proposition 2.1. For y generated through the mean-shift model (3), consider intervals CjM , satisfying
x
x
x “ M q “ 1 ´ α. If the outlier detection procedure M
x satisfies PpM
x Ď M ˚ q Ñ 1 as n Ñ 8,
Ppβ M P C M | M
x

j

j

then we have Ppβj˚ P CjM q Ñ 1 ´ α.
x

Proof. See Appendix D.1 for a detailed proof.
This proposition is based on two simple observations: first, that conditional coverage of β M implies
x
x Ď M ˚ implies that β M
unconditional coverage of β M ; second, that M
“ β˚.
Such a screening property is reasonable to demand of an outlier detection procedure, and related results
exist in the literature (Zhao et al. 2013). For example, consider using Cook’s distance (1) to detect outliers:
xpyq “ ti : Di ă λ{nu,
M

(4)

where λ is a prespecified cutoff. In Appendix A, we provide conditions (based on a result of Zhao et al.
x “ M ˚ q Ñ 1 for an appropriate choice of λ. While λ “ 0 would trivially satisfy
2013) under which PpM
the screening property, we of course need a procedure that leaves sufficient observations for estimation and
inference.

2.2

Quadratic Outlier Detection Procedures

In this section we define a general class of outlier detection procedures for which our methodology will apply.
We then show that this class includes several of the most famous outlier detection procedures.
xpyq “ M u “ ty P EM u,
Definition 2.2. We say an outlier detection procedure is quadratic if the event tM
where
EM :“ XptEM,i uiPIM q,
n

(5)

T

EM,i :“ ty P R : y QM,i y `

aTM,i y

` bM,i ě 0u,

(6)

for some finite index set IM and some QM,i P Rnˆn , aM,i P Rn , bM,i P R, and X denotes a general set operator
that maps a finite family of sets to a single set.
Generally, X should be thought of as taking finite unions, intersections, and complements.
Ş The above
definition is a direct generalization of Definition 1.1 of Loftus & Taylor (2015), in which X ” . We will see
that many outlier detection procedures are quadratic in the sense of Definition 2.2. While most of the time
the definition in Loftus & Taylor (2015) will apply, there are certain cases that require our generalization
(see Appendix C.2 for a specific example). The next proposition shows that outlier detection using Cook’s
distance is quadratic in the sense of Definition 2.2.
Proposition 2.3. Outlier detection using Cook’s distance (4) is quadratic with
č
EM “
EM,i ,

(7)

iPrns

#

Cook
EM,i
if i P M
,
Cook c
pEM,i q if i P M c
"
ˆ
˙
*
λp
K
K
K
:“ y P Rn : y T
p1 ´ hi q2 PX
´ pn ´ pqhi PX
ei eTi PX
yą0 .
n

EM,i “
Cook
EM,i

5

(8)
(9)

ˆ
Proof. We may write Di “

N
K
K
y T pPX
ei eTi PX
qy

˙ ˆ
N
˙
K
y T PX
y ¨ pn ´ pqhi pp1 ´ hi q2 , where ei is i-th standard

x “ ti : Di ă λ{nu and M
xc “ ti : Di ě λ{nu
basis for Rn and hi “ pPX qii . Plugging this expression to M
gives the desired result.
As a second example, we consider the soft-IPOD method (She & Owen 2011), which identifies outliers
using a lasso penalty (Tibshirani 1996) within the context of the mean-shift model (3):
pλ q “ arg min
pβpλ , u

βPRp ,uPRn

1
}y ´ Xβ ´ u}22 ` λ}u}1 .
2n

(10)

xpyq “ ti : u
pλ , and one takes M
pλ,i “ 0u as the detected non-outliers. In
The `1 -penalty induces sparsity in u
Appendix C, this approach is shown to be a quadratic outlier detection procedure. Interestingly, βpλ can be
equivalently expressed as the minimizer of the Huber loss function (as observed in She & Owen 2011):
#
n
1 2
r
if |r| ď λ
1 ÿ
ρpyi ´ Xi,¨ β, λq
minp
where
ρprq “ 2
1
2
βPR n
λ|r| ´ 2 λ
if |r| ą λ,
i“1
xpyq corresponding to those residuals in the quadratic (rather than linear) region of
with the elements of M
Huber’s loss function. Thus, applied in this example, the methodology of this paper provides valid “inference
after robust regression.”
One more outlier detection method, called DFFITS (Welsch & Kuh 1977), is described in Appendix C,
where it is shown to be quadratic. The only outlier detection method we know of that is not quadratic is
the Θ´IPOD procedure of She & Owen (2011) when Θ is not soft-thresholding.

3

Inference Corrected for Outlier Removal

In this section, we describe how the standard inferential tools of OLS can be corrected to account for outlier
removal. The only requirement is that the outlier detection procedure be quadratic (as defined in the
xpyq “ M u and are about
previous section). The inferential statements are made conditional on the event tM
the parameter β M . As previously discussed, such statements translate to unconditional statements about
x Ď M ˚ , that is, when all true outliers are removed. Section 3.1 treats the case in which σ is
β ˚ when M
known. Section 3.2 provides procedures for the case when σ is unknown.

3.1

Confidence Intervals and Hypothesis Tests When σ Is Known

In this section, we suppose that σ is known and provide confidence intervals and hypothesis tests. In the
classical setting, inference is based on the normal and χ2 distributions and typically involves individual
regression coefficients βjM , the regression surface xT0 β M , or groups of regression coefficients βgM . We begin
by observing that both βjM and xT0 β M are of the form ν T µ for some vector ν that depends on M : βjM “
`
`
T
T M
eTj XM,
¨ IM,¨ µ and x0 β “ x0 XM,¨ IM,¨ µ. The next theorem gives a unified treatment of these two cases
that will allow us to construct confidence intervals and p-values that properly account for outlier removal.
x “ M u is quadratic as in Definition 2.2. Let
Theorem 3.1. Assume the outlier detection procedure tM
n
ν P R be a vector that may depend on M . Define
Z :“
We have

νT y
,
σ}ν}2

z :“ PνK y “ pI ´

νν T
qy.
}ν}22

ˇ
T
ˇ
x “ M, zu „ T N p ν µ , 1; EM,z q,
Z ˇˇ tM
σ}ν}2

6

(11)

where T N pξ, γ 2 ; Eq denotes a N pξ, γ 2 q random variable truncated to the set E. And we can compute EM,z
by finding the roots of a finite set of quadratic polynomials:
EM,z “ XptEM,z,k ukPIM q,

(12)
T

EM,z,k “ tZ P R :

σ 2 ν T QM,k ν 2 2σpz QM,k ν `
Z `
}ν}22
}ν}2

aTM,k νq

Z ` z T Qm,k z ` aTM,k z ` bM,k ě 0u.

E
2
Thus, letting Fξ,γ
2 be the CDF of a T N pξ, γ ; Eq random variable, we have
ˇ
ˇ
E
x “ M u „ unifp0, 1q.
ˇ tM
pZq
1 ´ F ν TM,z
µ
ˇ
σ}ν}2

,1

(13)

(14)

Proof. See Appendix D.3.
The classical analogue to the above theorem is the (much simpler!) statement that Z „ N pν T µ{rσ}ν}2 s, 1q.
This theorem is essentially a generalization of Lee et al. (2016, Theorem 5.2) and a special case of Loftus
& Taylor (2015, Theorem 3.1); however, a key difference is that these works are focused on accounting for
variable selection rather than outlier removal (which, in essence, is “observation selection”).
3.1.1

Corrected Confidence Intervals

We begin by applying Theorem 3.1 to get confidence intervals corrected for outlier removal.
Corollary 3.2. Under the conditions and notation of Theorem 3.1, if we find L and U such that
L:F

EM,z
L
σ}ν}2 ,1

pZq “ 1 ´

U :F

EM,z
U
σ}ν}2 ,1

pZq “

α
2

α
,
2

(15)
(16)

then rL, U s is a valid p1 ´ αq selective confidence interval for ν T µ. That is,
x “ M q “ 1 ´ α.
Ppν T µ P rL, U s | M

(17)

Proof. See Appendix D.3.
This result encompasses the two most common types of confidence intervals arising in regression: intervals
for the regression coefficients βjM and intervals for the regression surface xT0 β M .
`
T
T
Corollary 3.3. We write βjM “ νcoef,j
µ and xT0 β M “ νsurf
µ, where νcoef,j “ peTj XM,
IM,¨ qT and νsurf “
¨
`
pxT0 XM,
IM,¨ qT . Then Theorem 3.1 and Corollary 3.2 apply.

¨

A third type of interval common in regression is the prediction interval, intended to cover xT0 β M ` ε0 ,
x “ M, zu is a
where x0 P Rp is a new data point and ε0 „ N p0, σ 2 q is independent of ε. While ν T y | tM
T
x
truncated normal random variable, ν y ` ε0 | tM “ M, zu is not, so the strategy adopted in Theorem 3.1
does not directly apply to this case. Instead we employ a simple (but conservative) strategy.
Proposition 3.4. Let ε0 „ N p0, σ 2 q be the noise independent of y. For a given significance level α P p0, 1q,
let α
r P p0, αq. Given x0 P Rp , let rLαr , Uαr s be the p1 ´ α
rq selective confidence intervals for xT0 β M as defined
in (15) and (16). Then we have
ˇ
ˆ
˙
r
r ˇˇ x
α´α
α´α
P Lαr ´ Φ´1 p1 ´
qσ ď xT0 β M ` ε0 ď Uαr ` Φ´1 p1 ´
qσ ˇ M “ M ě 1 ´ α,
(18)
2
2
where Φ is the CDF of a standard normal distribution.
Proof. See Appendix D.4 for a detailed proof.
r so that the length of the interval is minimized.
Remark. In practice, we can optimize over α
7

3.1.2

Corrected Hypothesis Tests

Theorem 3.1 allows us to form selective hypothesis tests about the parameter ν T µ where ν may depend on
the selected index set of observations M .
Corollary 3.5. Under the conditions and notation of Theorem 3.1, the quantity
ˇ
ˇ
EM,z
x “ M u „ unifp0, 1q.
1 ´ F0,1 pZq ˇˇ tM
gives a valid selective p-value for testing H0 : ν T µ “ 0.
Proof. See Appendix D.3.
The most common application of the above would be for testing whether a specific regression coefficient
is zero, conditional on M being the selected set of non-outliers: H0 pM, jq : βjM “ 0 for j P rps.
As a generalization, we next focus on testing H0 pM, gq : βgM “ 0 for g Ď rps. We begin with an alternative
characterization of H0 pM, gq.
rM,g “ pI|M | ´ PM,gc qXM,g , where PM,gc is the projection matrix onto C pXM,gc q.
Proposition 3.6. Set X
rM,g q. Then we have
Let PrM,g be the projection matrix onto C pX
r ` µM “ 0 ô PrM,g µM “ 0.
βgM “ 0 ô X
M,g

(19)

Further, define
ˆ
PqM,g :“

PrM,g

0|M |ˆpn´|M |q

0pn´|M |qˆ|M |

0pn´|M |qˆpn´|M |q

˙
.

Then PqM,g is an orthogonal projection matrix (it is symmetric and idempotent), and we have
βgM “ 0 ô PqM,g µ “ 0.

(20)

Proof. See Appendix D.5 for a detailed proof.
Remark. This proposition characterizes H0 pM, gq as testing the projection of µ. In the non-selective case,
testing P µ “ 0 for some projection matrix P can be done based on σ ´2 y T P y „ χ2trpP q under P µ “ 0. We
would expect that in the selective case, such tests can be done based on a truncated χ2 distribution.
x “ M u is quadratic as in Definition 2.2. Define
Theorem 3.7. Assume the outlier detection procedure tM
}PqM,g y}2
,
σ
PqM,g y
PqM,g y
w :“
“
,
σX
}PqM,g y}2

X :“

(21)
(22)

K
z :“ PqM,g
y.

(23)

Under H0 pM, gq : βgM “ 0, we have
ˇ
ˇ
x “ M, w, zu „ T χ2
X ˇˇ tM
pEM,w,z q,
trpPqM,g q
2

(24)

where the R.H.S is a central χ2 random variable with df “ trpPqM,g q truncated to the set EM,w,z . And we
can compute EM,w,z by finding the roots of a finite set of quadratic polynomials:
EM,w,z “ XptEM,w,z,k ukPIM q,
2

2

(25)

T

2

T

EM,w,z,k “ tX : X ě 0, pσ w QM,k wqX ` p2σw QM,k z `
` z T QM,k z ` aTM,k z ` bM,k ě 0u.
8

σaTM,k wqX
(26)

E
Further, letting Fdf
be the CDF of a T χ2df pEq random variable, we have
ˇ
ˇ
E
2 ˇ x
pX
q
1 ´ FtrpM,w,z
q
ˇ tM “ M u „ unifp0, 1q,
PM,g q

(27)

which is a valid selective p-value for testing H0 pM, gq : pβ M qg “ 0.
Proof. See Appendix D.6 for a detailed proof.
Remark. This theorem is adapted from Loftus & Taylor (2015, Theorem 3.1) to the outlier detection
context. In the special case where g “ j is a single index, direct computation can show that PqM,j “ Pνcoef,j ,
so that
T
pνcoef,j
yq2
T
, w “ signpνcoef,j
yqνcoef,j , z “ PνKcoef,j y.
X2 “ 2
σ }νcoef,j }22
Then this theorem nearly reduces to Theorem 3.1, except that in this theorem, we need to condition on the
T
sign of νcoef,j
y.

3.2

Extension to σ Unknown Case

In this section, we extend results in Section 3.1.2 to the σ unknown case. In the non-selective case, the
K
K
hypothesis H0 : βg˚ “ 0 is equivalent to H0 : µ P C pX¨,gc q. Hence under whichever H0 , pP¨K,gc ´PX
qy and PX
y
will both be centered normal random variables, and the test can be done based on F “

K
p}P¨K,gc y}22 ´}PX
y}22 q{|g|
K y}2 {pn´pq
}PX
2

„

F|g|,n´p . By analogy, we might expect H0 pM, gq :
“ 0 to be equivalent to H0 : µM P C pXM,gc q, which
would suggest that the test should be done based on a truncated F distribution; however, we will see in the
rest of the section that this is only partially true.
βgM

Proposition 3.8. We have µM P C pXM,gc q ñ βgM “ 0 but βgM “ 0 œ µM P C pXM,gc q. Moreover, if
M Ď M ˚ , then βgM “ 0 ñ µM P C pXM,gc q.
Proof. See Appendix D.7 for a detailed proof.
In order to form an F statistic, we need both the numerator and the denominator to be composed of
centered random variables. So it is necessary to assume µM P C pXM,gc q. Hence this proposition says that
testing H0 : µM P C pXM,gc q is the best we can do. Our next result adapts a truncated F significance test
from Loftus & Taylor (2015) to our purposes.
x “ M u is quadratic as in Definition 2.2. Let
Theorem 3.9. Assume the outlier detection procedure tM
K
K
R1 :“ Psub
y, R2 :“ Pfull
y, where
ˆ
˙
ˆ
˙
PM,gc
0|M |ˆpn´|M |q
0|M |ˆpn´|M |q
PM,¨
Psub :“
, Pfull :“
.
0pn´|M |qˆ|M |
Ipn´|M |q
0pn´|M |qˆ|M |
Ipn´|M |q
Define
p}R1 }22 ´ }R2 }22 q{|g|
,
}R2 }22 {p|M | ´ pq
R1 ´ R2
R2
w∆ :“
, w2 :“
, z :“ Psub y, r :“ }R1 }2 ,
}R1 ´ R2 }2
}R2 }2
d
d
|g|F{p|M | ´ pq
1
g1 pFq :“
, g2 pFq :“
.
1 ` |g|F{p|M | ´ pq
1 ` |g|F{p|M | ´ pq
F :“

Under H0 : µM P C pXM,gc q, we have
ˇ
ˇ
x “ M, w∆ , w2 , z, ru „ T F|g|,|M |´p pEM,w ,w ,z,r q,
F ˇˇ tM
∆
2
9

(28)
(29)
(30)

(31)

where the R.H.S. is a central F random variable with df1 “ |g|, df2 “ |M |´p truncated to the set EM,w∆ ,w2 ,z,r .
And we can compute EM,w∆ ,w2 ,z,r by
EM,w∆ ,w2 ,z,r “ XptEM,w∆ ,w2 ,z,r,k ukPIM q,

(32)

EM,w∆ ,w2 ,z,r,k “ tF ě 0 : y “ rpg1 pFqw∆ ` g2 pFqw2 q ` z, y T QM,k y ` aTM,k y ` bM,k ě 0u.

(33)

E
Further, letting Fdf
be the CDF of a T Fdf1 ,df2 pEq random variable, we have
1 ,df2

1´

EM,w ,w2 ,z,r
F|g|,|M∆
pFq
|´p

ˇ
ˇ
x “ M u „ unifp0, 1q,
ˇ tM
ˇ

(34)

which is a valid selective p-value for testing H0 : µM P C pXM,gc q.
Proof. See Appendix D.8 for a detailed proof.
Computing the truncation set in the σ unknown case is non-trivial since each slice is no longer a quadratic
function in F. We adopt the strategy suggested by Loftus & Taylor (2015, Section 4.1). For completeness,
we provide the details of their strategy (adapted to our notation) in Appendix E.1.
We conclude this section by noting that Theorem 3.9 does not give us a way to construct confidence
intervals for βjM . In order to form confidence intervals for βjM , one would need to be able to test for
H0 : βjM “ c0 for some non-zero constant c0 . Under this null, F does not necessarily reduce to the square of
a truncated t distribution: First, µM P C pXM,¨ q does not necessarily hold, and as a result, R2 may not even
be centered; second, the independence between F and pw∆ , w2 , z, rq may not hold. Hence the construction
of confidence intervals does not follow directly from Theorem 3.9 and is left as future work.

4

Empirical Examples

We provide simulations and real data examples in this section. We notice that our method requires evaluation
of survival functions (equivalently, the CDFs) of truncated normal, χ2 , t and F distributions. We refer readers
to Appendix E.2 for implementation details.

4.1

Simulations

In this section, we focus on the case where the outlier detection is done by Cook’s distance, and we assume
σ is unknown. We refer the readers to the supplementary materials for more detailed and comprehensive
simulations. We compare the performance of the following three inferential procedures:
• detect-and-forget: After outlier detection, refit an OLS regression model using the remaining data
pyM , XM,¨ q and do inference based on the classical (non-selective) theory (we use t and F distributions
since σ is unknown);
• corrected-est: Do selective inference as developed in Section 3.1.1 and 3.1.2, with estimated σ, and the
estimation of σ is done by
1
2
σ
pEST
“
}y ´ XAUG βpAUG }22 ,
n ´ |SAUG |
where we fit a lasso regression of y on XAUG “ pX : In q to get βpAUG P Rp`n , and SAUG is the support
of βpAUG . Reid et al. (2013) demonstrate that such a strategy gives a reasonably good estimate of σ 2
in a wide range of situations.
• corrected-exact: Do selective inference assuming unknown σ as developed in Section 3.2 (note: this
method does not give confidence intervals).
We fix n “ 100, p “ 11. Our indexing of variables starts from 0 (i.e. β0˚ corresponds to the intercept). The
first column of?X is set to be 1 and the rest of the columns are generated from i.i.d. N p0, 1q and scaled to
have `2 norm n. We fix σ “ 1.
10

To examine the coverage of confidence intervals for β1M , we let β ˚ “ p1, 2, 1, . . . , 1qT and M ˚c “
t1, 2, 3, 4, 5u. We then fix u˚M ˚c “ ps, s, s, ´s, ´sqT , and we vary s P t2, 3, 4, 5, 6u. Outliers are then detected
using Cook’s distance with different cutoffs λ P t1, 2, 3, 4u as introduced in Equation (4). For each configuration, we do the following 2000 times: we generate the response y “ Xβ ˚ ` u˚ ` ε, where ε „ N p0, σ 2 In q;
we then detect outliers and form confidence intervals. The detect-and-forget confidence intervals are set to
be
1´α{2

T
pREFIT }νcoef }2 t|M |´p s,
rνcoef
y˘σ
2
pREFIT is different from σ
pEST and, as noted in Fithian
pREFIT
where σ
“ |M1|´p }yM ´ XM,¨ βpM }22 (note that σ
et al. 2014, is generally not considered a good estimate of σ). Figure 3 shows the empirical coverage probability for β1M and β1˚ . As our theories predict, corrected-est intervals give 95% coverage of β1M , while detectand-forget intervals are off. Although without theoretical guarantees, corrected-est intervals still achieve the
desired coverage for β1˚ .

Figure 3: Empirical
coverage probability for β1M and β1˚ . The error bar of coverage probability is obtained
a
by q̂ ˘ 1.96 q̂p1 ´ q̂q{nsim, where q̂ is the empirical coverage probability and nsim “ 2000 is the number of
realizations. The dashed line represents 95% coverage.
Figure 4 shows the length of both kinds of intervals. We see that the achievement of desired coverage
comes with a price: the length of corrected-est intervals is in general wider than detect-and-forget intervals.
We next examine the power of testing H0 pM, 1q : β1M “ 0 against H1 pM, 1q : β1M ‰ 0. We let βk˚ “ 1
for k “ 0, 2, 3, . . . , 10, and we vary β1˚ smoothly. We let s “ 4 and the rest of the setup is the same as the
previous simulation. We run 2000 iterations. In each iteration, we generate the response, detect outliers,
and extract p-values. The detect-and-forget p-value is set to be
|M |´p

2Ft

p´|

T
νcoef
y
|q.
pREFIT }νcoef }2
σ

11

Figure 4: Length of confidence intervals for β1M .
where Ftdf is the CDF of a tdf distribution. For power considerations, corrected-est p-values are defined as
2 minp1´pval, pvalq, where pval is the p-value calculated by directly applying Corollary 3.5. By construction,
we are actually examining the power of testing H0 p˚, 1q : β1˚ “ 0 against H1 p˚, 1q : β1˚ ‰ 0. Figure 5 shows
the results: the two selective methods control the type I error down to 0.05 even though this correspond
to H0 p˚, 1q (recall that our theory ensures control under H0 pM, 1q), while detect-and-forgetdoes not. Both
corrected-est and corrected-exact suffer from a loss of power, although comparing to the power of detect-andforget is not meaningful since it does not control Type I error. The power for corrected-est seems acceptable,
while corrected-exact has quite a substantial loss in power, which may be the consequence of conditioning on
too much information.

˚
˚
Figure 5: The empirical
a power of testing H0 p˚, 1q : β1 “ 0 against H1 p˚, 1q : β1 ‰ 0. The error bar is
obtained by q̂ ˘ 1.96 q̂p1 ´ q̂q{nsim, where q̂ is the empirical power and nsim “ 2000 is the number of
realizations.

We next examine the Type I error and power of testing the global hypothesis H0 pM, gq : βgM “ 0 against
H1 pM, gq : βgM ‰ 0, where g “ t1, 2, . . . , 10u, the setup is the same as the previous simulation, except that
we let β0˚ “ 1, βk˚ “ 0 for k “ 2, 3, . . . , 10 and we vary β1˚ smoothly. The detect-and-forget p-value is set to
be
|g|,|M |´p
1 ´ FF
pFq,
where F is defined in Equation (28) and FFdf 1,df 2 is the CDF of an F distribution. We examine the power of
testing H0 p˚, gq : βg˚ “ 0 against H1 p˚, gq : βg˚ ‰ 0. Figure 6 shows the power as a funtion of β1˚ . Again we
notice the failure of detect-and-forget method to control the Type I error and the loss of power of the two
selective methods.
12

Figure 6: The empirical power of testing H0 p˚, gq : βg˚ “ 0 against H1 p˚, gq : βg˚ ‰ 0,The error bar is
a
obtained by q̂ ˘ 1.96 q̂p1 ´ q̂q{nsim, where q̂ is the empirical power and nsim “ 2000 is the number of
realizations.

4.2

Data Examples

We next apply our method on two classic data sets from the outlier detection literature. Since the number
of observations are relatively small, the estimation of variance will not be as accurate as in the previous
simulations, and we will use corrected-exact, despite the fact that we may suffer from a substantial loss of
power.
Stack Loss Data
Brownlee’s Stack Loss Plant Data (Brownlee 1965) involves measures on an industrial plant’s operation and
has 21 observations and three covariates. According to ?stackloss in R (R Core Team 2017), “Air.Flow is
the rate of operation of the plant, Water.Temp is the temperature of cooling water circulated through coils
in the absorption tower, and Acid.Conc is the concentration of the acid circulating, minus 50, times 10.”
The response, stack.loss, “is an inverse measure of the overall efficiency of the plant.” This data set is
considered by many papers in the outlier detection literature (Daniel & Wood 1999, Atkinson & Atkinson
1985, Hoeting et al. 1996). The general consensus is that observations 1, 3, 4 and 21 are outliers.
We use Cook’s distance to detect outliers, then fit the model and extract p-values, assuming σ is unknown.
The results are shown in Table 1. We see that as we detect outliers, the adjusted R2 increases, which is
an indication that the model is getting better. We also notice that corrected-exact p-values are in general
different from detect-and-forget ones, but there are several cases where the methods’ p-values coincide (e.g.
Water.Temp with cutoff 4). This is because the truncation set for the F statistic is r0, 8q in those cases. This
means that outlier removal does not have an effect on the conditional distribution of these test statistics.
Scottish Hill Races Data
This data set records the time for 35 Scottish hill races in 1984 (Atkinson 1986). There are two covariates:
“dist is the distance in miles, and climb is the total height gained during the route in feet”. The response,
time, “is the record time in hours”. This data set is also a classic one considered by many papers in the
outlier detection literature (e.g., Atkinson 1986, Hadi 1990, Hoeting et al. 1996). The consensus is that
observation 7 and 18 are obvious outliers, while observation 33 is an outlier that is masked by the other two
outliers.
Again, we use Cook’s distance to detect outliers, then fit the model and extract p-values, assuming σ is
unknown. The results are shown in Table 2. We can see the increase in adjusted R2 as outliers are detected,
and the corrected-exact p-values differ from the detect-and-forget p-values in general. Observation 33 is not
detected until the cutoff is set to 1, and observation 11 is always detected as an outlier. Atkinson (1986)
reports that observations 7, 18, 11, 33, 35 are high-leverage points but argues that only 7, 11, 33 are actual

13

Outlier Detected
Adjusted R2
Air.Flow: Estimate

Full Fit
None
0.8983
0.7156

Air.Flow: p-value

5.8 ˆ 10´5

Water.Temp: Estimate 1.2953
Water.Temp: p-value

0.00263

Acid.Conc: Estimate

´0.1521

Acid.Conc: p-value

0.34405

Cutoff = 4
21
0.9392
0.8891
0.00403
1.31 ˆ 10´6
0.8166
0.02309
0.02309
´0.1071
0.40234
0.40234

Cutoff = 3
1, 21
0.9171
0.8458
0.345
7.7 ˆ 10´6
0.8153
0.335
0.02431
´0.0881
0.376
0.49585

Cutoff = 2
1, 3, 4, 21
0.9692
0.7977
3.18 ˆ 10´4
2.48 ˆ 10´8
0.5773
0.00694
0.00408
´0.0671
0.2961
0.2961

Cutoff = 1
1, 2, 3, 4, 7, 12, 17, 21
0.9057
0.6666
0.245
1.19 ˆ 10´4
0.6357
0.792
0.01465
´0.0411
0.208
0.653

Table 1: Inference for each variable in Stack Loss data after outlier detection using Cook’s distance. The
p-values in bold font are the selective ones, while the other p-values are refitted ones.
outliers, while the others are high-leverage points that agree with the bulk of the data. But we recall that
our intent is not to concern ourselves with the accurate detection of outliers but rather with the proper
adjustment to inference based on outlier detection and removal.
Full Fit
Outlier Detected None
Adjusted R2
0.914
dist: Estimate 0.1036
dist: p-value

9.94 ˆ 10´12

climb: Estimate 1.84 ˆ 10´4
climb: p-value

6.49 ˆ 10´6

Cutoff = 4
7, 11, 18
0.9721
0.1138
1.76 ˆ 10´6
6.80 ˆ 10´15
1.28 ˆ 10´4
0.05918
9.15 ˆ 10´6

Cutoff = 3
7, 11, 18
0.9721
0.1138
1.06 ˆ 10´4
6.80 ˆ 10´15
1.28 ˆ 10´4
0.02465
9.15 ˆ 10´6

Cutoff = 2
7, 11, 18, 31
0.9723
0.1111
0.1219
6.50 ˆ 10´14
1.42 ˆ 10´4
0.06060
1.16 ˆ 10´5

Cutoff = 1
7, 11, 18, 31, 33, 35
0.9395
0.1034
6.99 ˆ 10´9
2.13 ˆ 10´13
1.17 ˆ 10´4
7.02 ˆ 10´4
1.53 ˆ 10´5

Table 2: Inference for each variable in Scottish Hill Races data after outlier detection using Cook’s distance.
The p-values in bold font are the selective ones, while the other p-values are refitted ones.

5

Discussion

In this paper, we have introduced an inferential framework for properly accounting for the removal of outliers
from a data set. The commonplace approach, detect-and-forget, makes the incorrect assumption that outlier
removal does not affect the distribution of the data. Our work is based on recent developments in the
selective inference literature, which carries out inference that properly accounts for variable selection (Lee
et al. 2016, Loftus & Taylor 2015). A key idea in that work is to characterize the event that a certain set
of variables is selected in terms of a simple to describe set of constraints on the response vector y. Doing
so makes it tractable to derive the conditional distribution of the estimator given this selection event. Our
work likewise relies on the fact that the most commonly used outlier detection procedures can be expressed
in a relatively simple form.
Our target of inference is β M , where M is the selected set of non-outliers. When M excludes all true
outliers, β M coincides with β ˚ . When the true outliers are easily detected, then (via Proposition 2.1), our
methodogy translates to inference on β ˚ . However, when the true outliers are less easily detected, this
may indicate that β M is not too far from β ˚ , meaning that our inferential statements may be translated,
approximately, to statements about β ˚ . An interesting future direction would be to characterize the regimes
(in terms of size of outlier) in which (i) all true outliers are easily detected and thus we can make inferential
statements about β ˚ and (ii) not all outliers are easily detected but β M « β ˚ so that approximate statements
about β ˚ can be made. And of course, a central question would then be whether there is a gap between
regimes (i) and (ii).

14

Another future direction would be to consider proper inference after outlier removal in the high-dimensional
setting. Our method explicitly assumes a low-dimensional setting through the assumption that XM
x, has
¨
linearly independent columns. A direct generalization of the outlier detection method (10) is to instead solve
1
}y ´ Xβ ´ u}22 ` λ1 }β}1 ` λ2 }u}1 .
βPRp ,uPRn 2n

pu
pβ,
pq “ arg min

Applying Lee et al. (2016, Theorem 4.3), one could perform inference corrected simultaneously for both
variable selection and outlier removal. Another approach would be to use a high-dimensional extension of
Cook’s distance proposed by Zhao et al. (2013) (it too can be shown to be a quadratic outlier detection
procedure). One could then do variable selection with the remaining data, for example using the lasso. In
this case our methodology would still, in principle, apply.

Acknowledgments
JB acknowledges support from an NSF CAREER grant, DMS-1653017.

Appendix A

Exact Screening Property of Cook’s Distance

xpyq as in (4). Of
Recalling the definition of Cook’s distance (1), we pre-specify a cutoff λ to determine M
course for sufficiently small λ, one can always assume that this outlier detection procedure screens out all
outliers with probability converging to 1. However, obviously this is ill-advised since there will not be enough
observations left for estimation and inference. The next proposition shows that, under suitable conditions,
this outlier detection procedure has the ability to screen out all outliers with probability converging to 1
with enough observations left for estimation and inference.
Proposition A.1. Consider the random design case, where each row of X comes from i.i.d. N p0, Σq,
independent of the noise ε. Let s0 :“ n ´ |M ˚ | be the number of true outliers. Assume the following
assumptions hold:
1. the covariance matrix Σ has bounded eigenvalues, and Σii “ 1 for all i P rns;
2. we are in a low-dimensional regime: p6 {n Ñ 0;
3. the outliers are sparse: s0 p{n1{3 Ñ 0;
4. the outlying signals do not vanish asymptotically: n2{3 {|u˚i | Ñ 0 for any i P M ˚c ;
5. the outlying signals are of the same order: u˚max {u˚min “ Op1q, where u˚max and u˚min are max u˚i and
min u˚i respectively for i P M ˚c .
Then for any i P M ˚c and any constant K ą 0, we have
Ppns0 Di ´ max˚ ns0 Dj ą Kq Ñ 1,

(35)

jPM

as n Ñ 8, where Di is the Cook’s distance of i-th observation.
Proof. See Appendix D.2 for a detailed proof.
Remark. This proposition says that under suitable conditions, outliers and non-outliers are well-separated
asymptotically in terms of their Cook’s distances, hence there exists a suitable λ that can screen out all the
outliers with probability converging to 1 with enough observations left at the same time. But notice that
this proposition does not give a explicit choice of λ. We also notice that this result can be generalized to the
form PpminiPM ˚c ns0 Di ´ maxjPM ˚ ns0 Dj ą Kq Ñ 1, using a simple union bound and assuming stronger
sparsity on s0 , but we do not pursue the details here as it is not our major focus.

15

Appendix B

Details of the Simulation in Section 1

We describe the details of the simulation given in Section 1. Fix n “ 20, p “ 2. The first column of X is
set to be 1 and the second column is set to be twenty equally spaced values between ´3 and 3. We assume
σ “ 1 to be known. We let β ˚ “ p1, 2qT and M ˚c “ t2u. Then we let u˚ “ p0, s, 0, . . . , 0qT and we vary
s P t5, 2, 1.5, 0.5u. The outlier detection procedure is defined as
xc pyq “ arg max Di ,
M

(36)

iPrns

where Di is the Cook’s distance (1) of i-th data point. Ignoring ties, we see that the outlier detection
procedure defined above is quadratic as in Definition 2.2 with
č
EM “
EM,i ,
iPrnszM

EM,i “ ty P Rn :

K
K
K
K
y pn ´ pqhi
y T PX
ei eTi PX
y T PX
eM eTM PX
y pn ´ pqhM
¨
ě
¨
u.
K
K
2
T
T
pp1 ´ hM q
pp1 ´ hi q2
y PX y
y PX y

We can simplify the expression of EM,i :
K
EM,i “ ty P Rn : y T PX
p

hM eM eTM
hi ei eTi
´
qP K y ě 0u.
p1 ´ hM q2
p1 ´ hi q2 X

Then we can use the theory developed in Section 3.1.1 and 3.1.2 to do selective inference.
We do the following 2000 times: generate y “ Xβ ˚ ` u˚ ` ε, where ε „ N p0, σ 2 q; then we detect outliers
using (36); then detect-and-forget and corrected intervals for β1M , xT0 β M are constructed respectively. Since
we assume σ is known, detect-and-forget confidence intervals are set to be
T
rνcoef
y ˘ Φ´1 p1 ´ α{2qσ}νcoef }2 s,
T
rνsurf
y ˘ Φ´1 p1 ´ α{2qσ}νsurf }2 s,

where Φ is the CDF of the standard normal distribution, and νcoef , νsurf are defined in Corollary 3.3. corrected confidence intervals are produced based on Corollary 3.3.

Appendix C
C.1

More Examples of Quadratic Outlier Detection Procedures

DFFITS

DFFITS is a measure of influential observations similar to Cook’s distance. Following Belsley et al. (2005,
Chapter 2.1), DFFITS for observation i can be written as
c
hi
DF F IT Si “
ri,ext ,
1 ´ hi
where hi “ pPX qii and
ri,ext :“

K
eTi PX
y
?
.
ppiq 1 ´ hi
σ

This can be equivalently expressed as
a
K
hi pn ´ p ´ 1q
eTi PX
y
b
DF F IT Si “
.
1 ´ h1
K y ´ y T P K e eT P K y{p1 ´ h q
y T PX
i
X i i X
We specify a cutoff λ (a typical choice being λ “ 4), and set
xpyq “ ti : DF F IT S 2 ă λp{pn ´ pqu.
M
i
16

Using similar arguments as in the proof of Proposition 2.3, we can show that outlier detection using DFFITS
is quadratic with
č
EM “
EM,i ,
iPrns

#

DFFITS
EM,i
if i P M
,
DFFITS c
pEM,i
q if i P M c
"
ˆ
˙
*
λp
λp
hi pn ´ p ´ 1q
K
K
T K
`
:“ y P Rn : y T
PX
´p
qP
e
e
P
y
ą
0
.
i
n´p
p1 ´ hi q2
pn ´ pqp1 ´ hi q X i X

EM,i “
DFFITS
EM,i

C.2

Soft-IPOD

xpyq “ ti : u
Consider solving the soft-IPOD problem (10) for outlier detection with M
pλ,i “ 0u. Following the
transformation suggested by Tibshirani et al. (2011, Section 3), we can first optimize for β with u fixed:
p
βpuq
“ X ` py ´ uq.
p
Plugging βpuq
back into (10) we are left with a lasso program,
u
p “ arg min
uPRn

1
q 2 ` λ}u}1 ,
}q
y ´ Xu}
2
2n

K
q “ P K . By Lee et al. (2016, Theorem 4.3), the event tM
x “ M, signpp
uq “ su can be
where yq “ PX
y and X
X
characterized by
"
ˆ
˙
ˆ
˙*
A0 pM, sq
b0 pM, sq
LASSO
n
yq ě
EM,s
:“ y P R :
,
A1 pM, sq
b1 pM, sq

where ě is understood to apply entry-wise and
˛
¨
˜
¸
K
qT
q T pX
q T c q` s
´1 ` X
| ,M c
1 ˝´X¨,M PX
,M
,M
¨
¨
¨
‚
, b0 pM, sq :“
,
A0 pM, sq :“
T
`
qT P K
q T pX
nλ
X
´1 ´ X
¨,M X|¨,M c
¨,M q¨,M c q s
´1
qT cX
q ` c , b1 pM, sq :“ nλdiagpsqpX
A1 pM, sq :“ diagpsqX
,M
¨,M q¨,M c q s.

¨

Since yq is linear in the original response y, and taking the union over all possible sign patterns, we see that
outlier detection using soft-IPOD is quadratic (it is actually affine) with
ď
LASSO
EM,s
EM “
.
sPt´1,1u|M |

x “ M, signpp
We only in fact need to consider tM
uq “ su for valid inference, i.e. we only need to compute one
LASSO
x “ M, signpp
slice EM,s . The reason is that any probabilistic statement conditioning on tM
uq “ su is also
x
valid conditioning on tM “ M u by marginalizing over all possible sign patterns (see Lee et al. 2016, Section
5 for a detailed discussion).

17

Appendix D
D.1

Proofs

Proof of Proposition 2.1

Consider the following decomposition:
ÿ
x
x “ M qPpM
x “ Mq
Ppβj˚ P CjM q “
Ppβj˚ P CjM | M
M :M ĎM ˚

ÿ
`

x “ M qPpM
x “ Mq
Ppβj˚ P CjM | M

M :M ĘM ˚

ÿ

x “ Mq
p1 ´ αqPpM

“
M :M ĎM ˚

ÿ
`

x “ M qPpM
x “ Mq
Ppβj˚ P CjM | M

M :M ĘM ˚

ÿ

x Ď M ˚q `
“ p1 ´ αqPpM

x “ M qPpM
x “ M q.
Ppβj˚ P CjM | M

M :M ĘM ˚

x “ M q “ Ppβ M P C M | M
x“
The second equality uses that for M Ď M ˚ , βj˚ “ βjM and so Ppβj˚ P CjM | M
j
j
M q “ 1 ´ α. We can get lower and upper bounds for the desired probability:
x Ď M ˚ q ď Ppβ ˚ P C M q ď p1 ´ αqPpM
x Ď M ˚ q ` PpM
x Ę M ˚ q.
p1 ´ αqPpM
j
j
x

x Ď M ˚ q Ñ 1 and PpM
x Ę M ˚ q Ñ 0, hence Ppβ ˚ P C M q Ñ 1 ´ α.
By assumption, we have PpM
j
j
x

D.2

Proof of Proposition A.1

This proof is modified from the proof of Zhao et al. (2013, Theorem 2), where they essentially show the same
results assuming there is only one outlier. For notational convenience, we let S ˚ :“ M ˚c , and we write each
row of X to be Xi and each column to be xj . We let X̃i “ Xi {}Xi }2 be the unit vector. We denote the
maximal and minimal eigenvalue of matrix A as λmax pAq and λmin pAq, respectively. We also write u˚ ” u
when there is no ambiguity. We begin by proving two useful lemmas.
Lemma D.1. Under the same assumptions as in Proposition A.1, the following statements hold:
1. n´1 maxiPrns }Xi }22 “ op p1q;
2. for any i P S ˚ , we have |ui |´1 maxjPrns |εj |}Xj }2 “ op p1q;
3. for any i P S ˚ , we have |ui |´1 }X T ε}2 “ op p1q;
4. there exists some 0 ă Cmin ď Cmax ă 8 such that PpCmin ď λmin pX T X{nq ď λmax pX T X{nq ď
Cmax q Ñ 1;
ř
ř
5. iPS ˚ jPS ˚ ui uj XiT pX T Xq´1 Xj “ op ps0 u2max q;
ř
u
r T pX T X{nq´1 X
rj ;
6. for any i P S ˚ , we have n´1 }Xi }2 jPS ˚ uji }Xj }2 Hij “ op p1q, where Hij “ X
i
ř
u
7. for any k P S ˚ , we have n´1 jPS ˚ maxiPM ˚ ukj }Xi }22 }Xj }2 “ op p1q;
r T pX T X{nq´1 X T ε{}X T ε}2 and ηri :“ X T pX T Xq´1 X T ε, then maxiPrns ηi “
8. for any i P rns, define ηi :“ X
i
i
Op p1q and ηri “ op p|ui |q.
Proof. Parts 1–4 are proved in the proof of Zhao et al. (2013, Theorem 2) (Section 3.1, step 1 in the
supplements). Proofs for part 5–8 are modified from the proof of Zhao et al. (2013, Theorem 2).

18

ri pX T X{nq´1 X
rj . Define the
Part 5: We can write XiT pX T Xq´1 Xj “ n´1 }Xi }2 }Xj }2 Hij , where Hij :“ X
´1
event B :“ tmaxi,jPrns Hij ą Cmin u. By part 4 of this lemma, PpBq Ñ 0. Then for any δ ą 0, we have
ˇ
˙
˙
ˆ
ˆˇ ÿ ÿ
ÿ ÿ
ˇ
ˇ
}Xi }2 }Xj }2 Hij ą s0 δ
ui uj XiT pX T Xq´1 Xj ˇˇ ą s0 u2max δ ď P n´1
P ˇˇ
iPS ˚ jPS ˚

iPS ˚ jPS ˚

ď

Pps20 n´1

max }Xi }2 }Xj }2 Hij ą s0 δq

i,jPS ˚

“ Pp max˚ }Xi }2 }Xj }2 Hij ą ns´1
0 δq
i,jPS
ˆ"
*č ˙
ďP
max˚ }Xi }2 }Xj }2 Hij ą ns´1
δ
Bc
0
i,jPS

` PpBq.
On the event B c , we have
´1
max }Xi }2 }Xj }2 Hij ď max˚ }Xi }2 }Xj }2 max Hij ď max˚ }Xi }2 }Xj }2 Cmin
.

i,jPS ˚

i,j

i,jPS

i,jPS

Hence we have
ˆ"
*č ˙
P
max˚ }Xi }2 }Xj }2 Hij ą ns´1
δ
B c ď Pp max˚ }Xi }2 }Xj }2 ą ns´1
0
0 Cmin δq
i,jPS
i,jPS
ÿ ÿ
ď
Pp}Xi }2 }Xj }2 ą ns´1
0 Cmin δq.
iPS ˚ jPS ˚

For i “ j, we have
´1
2
Pp}Xi }2 }Xj }2 ą ns´1
0 Cmin δq “ Pp}Xi }2 ą ns0 Cmin δq
2
2
ď Pp}Xi }42 ą n2 s´2
0 Cmin δ q
´2 ´2 ´2 2
ď Cmin
δ n s0 Er}X1 }42 s.
4
Recall that X1j „ N p0, 1q for any j, so ErX1j
s “ 3. We then have

Er}X1 }42 s “ Erp

p
ÿ

2 2
4
X1j
q s ď p2 max ErX1j
s “ 3p2 .
jPrps

j“1

´2 ´2 ´2 2 2
Hence for i “ j, we get Pp}Xi }2 }Xj }2 ą ns´1
n s0 p . Now for i ‰ j, we have
0 Cmin δq ď 3Cmin δ
2
2
2 ´2 2
2
Pp}Xi }2 }Xj }2 ą ns´1
0 Cmin δq ď Pp}Xi }2 }Xj }2 ą n s0 Cmin δ q
´2 ´2 ´2 2
ď Cmin
δ n s0 Er}Xi }22 sEr}Xj }22 s
´2 ´2 ´2 2 2
“ Cmin
δ n s0 p .

Putting those two cases together, we have
ˆ"
*č ˙
´2 ´2 ´2 2 2
P
max˚ }Xi }2 }Xj }2 Hij ą ns´1
δ
B c ď ps20 ` 2s0 qpCmin
δ n s0 p q
0
i,jPS

´2 ´2 ´2 4 2
ď 2Cmin
δ n s0 p

for n large

Ñ0
by Assumption 3 in Proposition A.1.
Part 6: By Assumption 5 in Proposition A.1, DK ą 0 such that |uj {ui | ď K for n large. Then for any
i P S ˚ , δ ą 0 and n large, we have
ˇ
ˆˇ
˙
ÿ uj
ÿ
ˇ
ˇ
ˇ
ˇ
P ˇ}Xi }2
}Xj }2 Hij ˇ ą nδ ď Pp}Xi }2
}Xj }2 Hij ą nK ´1 δq
ui
jPS ˚

jPS ˚

ˆ"
*č ˙
ÿ
´1
P }Xi }2
}Xj }2 Hij ą nK δ
B c ` PpBq.
jPS ˚

19

Then we have
ˆ"
*č ˙
ÿ
ÿ
P }Xi }2
}Xj }2 Hij ą nK ´1 δ
B c ď Pp}Xi }2
}Xj }2 ą nK ´1 Cmin δq
jPS ˚

jPS ˚

ď Pps0 }Xi }2 max
}Xj }2 ą nK ´1 Cmin δq
jPS ˚
ÿ
´1
ď
Pp}Xi }2 }Xj }2 ą s´1
Cmin δq
0 nK
jPS ˚

ÿ

2 ´2 2
Pp}Xi }22 }Xj }22 ą s´2
Cmin δ 2 q
0 n K

ď
jPS ˚

ÿ

´2 ´2 ´2 2
ď K 2 Cmin
δ n s0

Er}Xi }22 }Xj }22 s.

jPS ˚

Recall that in Er}Xi }42 s ď 3p2 and Er}Xi }22 }Xj }22 s “ p2 for i ‰ j. Hence we have
ˆ"
*č ˙
ÿ
´2 ´2 ´2 2
P }Xi }2
}Xj }2 Hij ą nK ´1 δ
B c ď K 2 Cmin
δ n s0 ps0 ` 2qp2
jPS ˚
´2 ´2 ´2 3 2
ď 3K 2 Cmin
δ n s0 p

Ñ0
by Assumption 3 in Proposition A.1.
Part 7: Again by Assumption 5 in Proposition A.1, DK ą 0 such that |uj {uk | ď K for any j, k P S ˚ when
n large. So for any k P S ˚ , δ ą 0, when n large, we have
ˇ
˙
ˆˇ ÿ
ÿ
ˇ
ˇ
uj
2
ˇ
ˇ
max }Xi }22 }Xj }2 ą nK ´1 δq
P ˇ
max˚ | |}Xi }2 }Xj }2 ˇ ą nδ ď Pp
uk
iPM ˚
iPM
˚
˚
jPS

jPS

´1
ď Pp max
}Xi }22 }Xj }2 ą s´1
δq
0 nK
jPS ˚ ,iPM ˚
ÿ ÿ
´1
ď
Pp}Xi }22 }Xj }2 ą s´1
δq
0 nK
iPM ˚ jPS ˚

ÿ

ÿ

ď

2 ´2 2
Pp}Xi }42 }Xj }22 ą s´2
δ q
0 n K

iPM ˚ jPS ˚

ď K 2 δ ´2 s20 n´2

ÿ

ÿ

iPM ˚

jPS ˚

Er}Xi }42 sEr}Xj }22 s

ď K 2 δ ´2 s20 n´2 pns0 qp3p2 pq
“ 3K 2 δ ´2 n´1 s30 p3
Ñ0
by Assumption 3 in Proposition A.1.
ˆ
Part 8: maxiPrns ηi “ Op p1q directly follows from part 4 of this lemma. We can write ηri “

˙ˆ
n

´1

}Xi }2

˙
T

}X ε}2 ηi “

op p1q ¨ op p|ui |q ¨ Op p1q “ op p|ui |q by part 1 and 3 of this lemma.
Lemma D.2. Under the same assumptions as in Proposition A.1, we have
p2 “ pn ´ pq´1 s0 u2min pOp1q ` op p1qq,
σ
K
where σ
p2 :“ }p
ε}22 {pn ´ pq and εp :“ PX
y.

Proof. This proof is adapted from the proof of Zhao et al. (2013, Theorem 2). Consider the decomposition

20

K
εp “ V1 ´ V2 ` V3 , where V1 “ PX
ε, V2 “

ř
iPS ˚

XpX T Xq´1 Xi ui , V3 “ u. Direct computation shows that

V1T V2 “ 0,
ÿ ÿ

V2T V3 “ }V2 }22 “

ui uj XiT pX T Xq´1 Xj ,

iPS ˚ jPS ˚

V1T V3

ÿ
“

pεi ui ´ ui ηri q,

iPS ˚

}V3 }22 “

ÿ

u2i ,

iPS ˚
K
where we define ηri :“ XiT pX T Xq´1 X T e. Since PX
only has 0 and 1 as its eigenvalues, we have }V1 }22 “
2 T K
}ε}2 εr PX εr “ Op pnqOp p1q “ Op pnq, where εr :“ ε{}ε}2 . By Assumption 4 in Proposition A.1, we have
n{c2max Ñ 0, and hence }V1 }22 “ op pc2max q. By Part 5 of Lemma D.1, we have V2T V3 “ op ps0 c2max q. We also
have
˙
ÿ ˆ
V1T V3 “
Op p1qui ´ ui op p|ui |q ,
iPS ˚

where we call ηri :“ XiT pX T Xq´1 X T ε “ op p|ui |q by part 8 of Lemma D.1. Hence
˙
ÿ ˆ
T
2
V1 V3 “
Op pui q ´ op pui q
iPS ˚

ÿ ˆ
“

˙
op pu2i q ´ op pu2i q

iPS ˚

“ op ps0 u2i q
“ op ps0 u2max q.
As for }V3 }22 , we have }V3 }22 “ Ops0 u2max q. Hence we have
p2 “ pn ´ pq´1 p}V1 }22 ´ }V2 }22 ` }V3 }22 ` 2V1T V3 q
σ
ˆ
˙
“ pn ´ pq´1 op pu2max q ´ op ps0 u2max q ` op ps0 u2max q ` Ops0 u2max q
ˆ
˙
“ pn ´ pq´1 s0 u2max Op1q ` op p1q
ˆ
˙
“ pn ´ pq´1 s0 u2min Op1q ` op p1q ,
where the last step follows from Assumption 5 in Proposition A.1.
Equipped with the above lemmas, the rest of the proof is very similar to that of Zhao et al. (2013,
Theorem 2). Recall that we can write Cook’s distance as
Di “

εp2i
hii
,
2
pp
σ p1 ´ hii q2

where
K
εp “ PX
y “ pu ` εq ´ XpX T Xq´1 X T pu ` εq,

σ
p2 “ }p
ε}22 {pn ´ pq,
hij “ XiT pX T Xq´1 Xj .
For i P S ˚ , we can write
εpi “ εi ´ XiT pX T Xq´1 X T ε ` ui ´ XiT pX T Xq´1

ÿ
jPS ˚

21

Xj uj .

For i P M ˚ , we can write
εpi “ εi ´ XiT pX T Xq´1 X T ε ´ XiT pX T Xq´1

ÿ

Xj uj .

jPS ˚

ri :“ Xi {}Xi }2 . So we can write
Recall that we define X
T
T
r T p X X q´1 X e
XiT pX T Xq´1 X T e “ n´1 }Xi }2 }X T e}2 X
i
T
n
}X e}2
looooooooooooomooooooooooooon
:“ηi

:“ ηri .
r T pX T X{nq´1 X
rj . Thus, we can rewrite εp as:
Also recall that we define Hij :“ X
i
ÿ
for i P S ˚ , εpi “ εi ´ ηri ` ui ´
uj n´1 }Xi }2 }Xj }2 Hij
jPS ˚

ÿ uj
εi ´ ηri
´
n´1 }Xi }2 }Xj }2 Hij q;
ui
ui
˚
jPS
ÿ
˚
for i P M , εpi “ εi ´ ηri ´
uj n´1 }Xi }2 }Xj }2 Hij
“ ui p1 `

jPS ˚

“ uk p

ÿ uj
εi ´ ηri
´
n´1 }Xi }2 }Xj }2 Hij q,
uk
uk
˚
jPS

where k is any index in S ˚ . Expressing hii in terms of Hii , we have
εp2i n´1 }Xi }22 Hii
.
2
pp
σ p1 ´ n´1 }Xi }22 Hii q2

Di “
Define
An,ui :“ p1 `

1{2
ÿ uj
}Xi }2 Hii
εi ´ ηri
,
´
n´1 }Xi }2 }Xj }2 Hij q
ui
ui
|1 ´ n´1 }Xi }22 Hii |
˚
jPS

Gn,uk :“ max˚ p
iPM

1{2
ÿ uj
εi ´ ηri
}Xi }2 Hii
.
´
n´1 }Xi }2 }Xj }2 Hij q
uk
uk
|1 ´ n´1 }Xi }22 Hii |
˚
jPS

It can be verified that
Di “ u2i ppnp
σ 2 q´1 A2n,ui ,
max Di “ u2k ppnp
σ 2 q´1 G2n,uk ,

iPM ˚

i P S˚
k is any index in S ˚ .

We then start to show that
p

1{2

An,ci Ñ }Xi }2 Hii ,
p

Gn,ck Ñ 0,
p

where we use Ñ to denote convergence in probability. We first consider Gn,uk . For any δ ą 0, we define the
event
"
*č"
*
XT X
XT X
´1
2
A0 :“ n max }Xi }2 ă δCmin
Cmin ď λmin p
q ď λmax p
q ď Cmax .
n
n
iPrns
By part 4 of Lemma D.1, we know PpAc0 q Ñ 0. On the event A0 , we have
max Hij ď λmax pp

i,jPrns

X T X ´1
´1
q q ď Cmin
,
n

22

and for small enough δ, we have
1{2

1{2

´1{2

}Xi }2 Hii
}Xi }2 Cmin
}Xi }2 Hii
“
ď
.
|1 ´ n´1 }Xi }22 Hii |
1 ´ n´1 }Xi }22 Hii
1´δ
Hence on A0 , we have
ˆ
|Gn,uk | ď

ηi |}Xi }2 `
|uk |´1 max˚ |εi |}Xi }2 ` |uk |´1 max˚ |r
iPM
iPM
loooooooooooomoooooooooooon
loooooooooooomoooooooooooon
p1q

p2q

:“Gn,uk

:“Gn,uk

uj
n
max˚ | |}Xi }22 }Xj }2
u
iPM
k
jPS ˚
loooooooooooooooooooooomoooooooooooooooooooooon
´1

´1
Cmin

´1{2

˙

ÿ

¨

Cmin
.
1´δ

p3q

:“Gn,uk
p1q

By part 2 of Lemma D.1, we have Gn,uk “ op p1q. Recalling the definition of ηri , we have
´1 ´1
Gp2q
n }Xi }22 }X T ε}2 |ηi |
n,uk “ max |uk |
iPM
ˆ
˙ˆ
˙ˆ
˙
ď max˚ |uk |´1 }X T e}2
max˚ n´1 }Xi }22
max˚ |ηi |
iPM

iPM

iPM

“ op p1q ¨ op p1q ¨ Op p1q,
p2q

where the last step is true by part 3, 1, 8 of Lemma D.1 respectively. Hence Gn,uk “ op p1q. Meanwhile,
p3q
Gn,uk “ op p1q by part 7 of Lemma D.1. So Gn,uk “ op p1q on A0 . We then notice that
ˆ"
*č ˙
Pp|Gn,uk | ą δq ď P |Gn,uk | ą δ
A0 ` PpAc0 q Ñ 0.
p

So Gn,uk Ñ 0. We then consider An,ui . It is easy to see that εi “ Op p1q. Since |ui | Ñ 8 by Assumption 4
in Proposition A.1, we have εi {ui “ op p1q. Also we know ηri “ op p|ui |q by part 8 of Lemma D.1. Now notice
that in the definition of An,ci ,
ÿ uj
n´1 }Xi }2 }Xj }2 Hij “ op p1q
u
i
˚
jPS

by part 6 of Lemma D.1, and n }Xi }22 Hii “ op p1q ¨ Op p1q “ op p1q by part 1 and 4 of Lemma D.1. Putting
p
1{2
the above pieces together, we conclude that An,ui Ñ }Xi }2 Hii , so that
´1

p

A2n,ui ´ G2n,uk Ñ }Xi }22 Hii .
In the proof of Zhao et al. (2013, Theorem 2) (step 4.2 in the supplements), they show that p´1 }Xi }22 “
´1
1 ` op p1q, so Ppp´1 }Xi }22 ą 1 ´ δq Ñ 1 for small δ ą 0. Combining with PpHii ą Cmax
q Ñ 1, we have
´1
PpA2n,ui ´ G2n,uk ą pp1 ´ δqCmax
q Ñ 1.
´2
p2 pn ´ pqs´1
Lemma D.2 tells that σ
0 umin “ Op1q ` op p1q. Hence for any δ ą 0 we have
´2
σ 2 pn ´ pqs´1
Ppp
0 umin ´ θn ą δq Ñ 0,

where we use θn to denote the Op1q term. Suppose that θn is asymptotically bounded above by Kmax ą 0,
then for n large we have
´2
´2
σ 2 pn ´ pqs´1
σ 2 pn ´ pqs´1
Ppp
0 umin ą Kmax ` δq ď Ppp
0 umin ´ θn ą δq Ñ 0.

Take δ “ 1, we get
´2
Ppp
σ 2 pn ´ pqs´1
0 umin ą Kmax ` 1q Ñ 0.

23

Equivalently, we have
ˆ
˙
2 ´1 ´2
´1
p
P σ
s
u
p
ă
pK
`
1qppn
´
pq
Ñ 1.
max
0
min
loooooooooooooooooooooooomoooooooooooooooooooooooon
:“D1

PpA2n,ui

G2n,uk

´1
δqCmax
q

´1
Recall that
´
ą pp1 ´
Ñ 1. Since pp1 ´ δqCmax
Ñ 8 and ppn ´ pq´1 Ñ 0 by
Assumption 2 in Proposition A.1, we have that, for any K ą 0,
ˆ
˙
2
2
´1
P looooooooooooooooooooooooooomooooooooooooooooooooooooooon
An,ui ´ Gn,uk ą pKmax ` 1qppn ´ pq K Ñ 1.
:“D2

For i P S ˚ , define
have

piq
D3

piq

:“ tA2n,ui ´ G2n,ui ą 0u. It is easy to see that PpD3 q Ñ 1. For any i P S ˚ , K ą 0, we
ˆ

˙
s0 c2i ppp
σ 2 q´1 A2n,ui

Ppns0 Di ´ max˚ ns0 Dj ą Kq “ P
jPM

´

s0 c2k ppp
σ 2 q´1 G2n,uk

ˆ

ąK

˙

σ2
“ P s0 u2i pA2n,ui ´ G2n,ui q ą Kpp
take k “ i
ˆ"
*č
˙
piq
2
2
2
2
ě P s0 ui pAn,ui ´ Gn,ui q ą Kpp
σ
D3
˙
ˆ"
*č
piq
2
2
2
2
ě P s0 umin pAn,ui ´ Gn,ui q ą Kpp
σ
D3
ˆ"
*č
˙
piq
2
2
2 ´1 ´2
σ s0 umin p
“ P An,ui ´ Gn,ui ą Kp
D3
piq

ě PpD1 X D2 X D3 q
Ñ 1,
which completes the proof.

D.3

Proof of Theorem 3.1 and Corollaries 3.2 and 3.5

This proof is modified from the proof of Loftus & Taylor (2015, Theorem 3.1). We begin by notice that
y “ σZν{}ν}2 `z. Plugging this expression to the quadratic constraints in Definition 2.2 gives (12)(13). Now
K z by orthogonality and properties of normal distribution.
notice that for the fixed index set M , we have Z K
Then it follows that
ˇ
ˇ
ˇ
ˇ
ˇ
ˇ
d
d
x “ M, z “ z0 u “
Z ˇˇ tZ P EM,z0 , z “ z0 u “ Z ˇˇ tZ P EM,z0 u,
Z ˇˇ tM
d

where “ denotes equivalence in distribution, and we drop
ˇ tz “ z0 u in the last equation by independence.
ˇ
x “ M, z “ z0 u is a truncated normal random
Notice that there is no randomness in EM,z0 , hence Z ˇˇ tM
variable as in (11). Applying inverse survival function transform (equivalently, the inverse CDF transform)
gives
ˇ
ˇ
E
0
x “ M, z “ z0 u „ unifp0, 1q.
ˇ tM
1 ´ F ν TM,z
pZq
µ
ˇ
σ}ν}2

,1

Then we can marginalize over all possible values of z0 and conclude (14). To show (17), we first notice
E
that Fξ,γ
2 is a decreasing function in ξ (this follows from the fact that the truncated normal distribution has
monotone likelihood ratio in the mean parameter, see Lee et al. (2016, Lemma A.1) for a detailed discussion).
Hence we have
ˇ
˙
ˆ
ˇ
EM,z
EM,z
x
x “ M q “ P F EM,z
ˇ
M
“
M
PpL ď ν T µ ď U | M
pZq
ď
F
pZq
ď
F
pZq
L
U
νT µ
ˇ
,1
σ}ν}2 ,1
σ}ν}2 ,1
σ}ν}2

α
α x
“ Mq
“ Pp ď unifp0, 1q ď 1 ´ | M
2
2
“ 1 ´ α.
24

D.4

Proof of Proposition 3.4

x “ M u, B “ t|ε0 | ď Φ´1 p1 ´
Denote the event A “ tLαr ď xT0 β M ď Uαr | M
PpAq “ 1 ´ α
r and PpBq “ 1 ´ pα ´ α
rq. Then

r
α´α
2 qσ

x “ M u. We have
|M

PpA X Bq “ PpAq ` PpBq ´ PpA Y Bq
ě PpAq ` PpBq ´ 1
“ 1 ´ α.
α
r
α
r
T M
x
Since A X B ñ tLαr ´ Φ´1 p1 ´ α´
` ε0 ď Uαr ` Φ´1 p1 ´ α´
2 qσ ď x0 β
2 qσ | M “ M u, we have
ˇ
ˆ
˙
α´α
r
α´α
r ˇˇ x
P Lαr ´ Φ´1 p1 ´
qσ ď xT0 β M ` ε0 ď Uαr ` Φ´1 p1 ´
qσ ˇ M “ M ě PpA X Bq ě 1 ´ α.
2
2

D.5

Proof of Proposition 3.6

The desired result follows from direct computations as follows. By definition we have
ˆ T
˙´1 ˆ T ˙
T
XM,g XM,g XM,g
XM,gc
XM,g
βM “
µM .
T
T
T
XM,g
XM,g
XM,g
c XM,g
c XM,g c
c
ˆ
˙
A B
Write K “
to be the block matrix above. The inversion formula gives
C D
pK ´1 q11 “ pA ´ BD´1 Cq´1 ,

pK ´1 q12 “ ´pA ´ BD´1 Cq´1 BD´1 .

We are only interested in
ˆ
M

pβ qg “

˙
pK

´1

T
q11 XM,g

` pK

´1

T
q12 XM,g
c

µM .

We have
T
T
T
´1 T
pA ´ BD´1 Cq´1 “ pXM,g
XM,g ´ XM,g
X
XM,gc XM,g q´1
M,g c pXM,g c XM,g c q
loooooooooooooooooomoooooooooooooooooon
PM,gc

“
“

T
pXM,g
pI ´ PM,gc qXM,g q´1
T
K
´1
pXM,g
PM,g
.
c XM,g q

So we have
T
T
K
´1 T
pK ´1 q11 XM,g
“ pXM,g
PM,g
XM,g
c XM,g q
T
T
K
´1 T
T
´1 T
pK ´1 q12 XM,g
XM,g X
XM,gc
c “ ´pXM,g PM,g c XM,g q
M,g c pXM,g c XM,g c q
loooooooooooooooooomoooooooooooooooooon
PM,gc

Hence
T
K
´1 T
pβ M qg “ pXM,g
PM,g
XM,g pI ´ PM,gc qµM
c XM,g q
T
K
´1 T
K
“ pXM,g
PM,g
XM,g PM,g
c XM,g q
c µM

r ` µM ,
“X
M,g
which shows the first iff in Equation (19). For the second iff in Equation (19), we have
r ` µM “ 0 ñ X
rM,g X
r ` µM “ 0 ñ PrM,g µM ;
• X
M,g
M,g
rM,g X
r ` µM “ 0 ñ X
r` X
r
r`
r`
• PrM,g µM “ 0 ñ X
M,g
M,g M,g XM,g µM “ 0 ñ XM,g µM “ 0.
Then we consider P9M,g . Its symmetry and idempotency are easy to check. Now notice that
P9M,g µ “ PrM,g µM ,
which establishes the iff in Equation (20).
25

D.6

Proof of Theorem 3.7

This proof is modified from the proof of Loftus & Taylor (2015, Theorem 3.1). We begin by proving two
useful lemmas.
Lemma D.3. Consider y „ N p0, σ 2 Iq. Then the length }y}2 is a complete sufficient statistic for σ ´2 and
the direction y{}y}2 is ancillary for σ. Moreover, we have }y}2 K
K y{}y}2 .
Proof. Notice that y is in the exponential family with densities fY pyq9 exptσ ´2 }y}22 u. By Shao (2003,
Proposition 2.1), }y}22 is a complete sufficient statistic for σ ´2 . Hence by the factorization theorem and the
definition of completeness, the length }y}2 is complete and sufficient for σ. To show the direction y{}y}2 is
ancillary, we notice that y{}y}2 “ σ ´1 y{}σ ´1 y}2 , and σ ´1 y „ N p0, Iq, so the distribution of the direction
does not depend on σ. Then the independence between the length and the direction follows from Basu’s
theorem.
Lemma D.4. Consider y „ N pµ, σ 2 Iq. Let L Ď Rn be a linear subspace and let PL be the projection
matrix onto L. If PL µ “ 0, then the length }PL y}2 is a complete sufficient statistic for σ and the direction
PL y{}PL y}2 is ancillary for σ. Moreover, we have }PL y}2 K
K PL y{}PL y}2 .
Proof. This proof is modified from the proof of Yang et al. (2016, Lemma A.1). Take V to be the matrix whose
columns form an orthonormal basis of L. We then have V T y „ N pV T µ, σ 2 Iq. Notice that PL “ V V T and we
have V V T µ “ 0 by our assumption. Multiplying both sides by V T gives V T µ “ 0. Hence V T y „ N p0, σ 2 Iq
and by Lemma D.3, we know }V T y}2 is complete and sufficient for σ and V T y{}V T y}2 is ancillary. Now
notice that }V T y}2 “ }PL y}2 , which gives the completeness and sufficiency of }PL y}2 . On the other hand, we
have PL y{}PL y}2 “ V V T y{}V T y}2 , which establishes the ancillarity of PL y{}PL y}2 . Then the independence
follows from Basu’s theorem.

K pw, zq for
By proposition 3.6, PqM,g µ “ 0 under H0 pM, gq : βgM “ 0. Then by Lemma D.4, we have Z K
the fixed index set M . Now notice that y “ σX w ` z. Plugging this expression to the quadratic constraints
in Definition 2.2 gives (25)(26). We then have
ˇ
ˇ
ˇ
ˇ
d
x “ M, w “ w0 , z “ z0 u “
X 2 ˇˇ tX 2 P EM,w0 ,z0 , w “ w0 , z “ z0 u
X 2 ˇˇ tM
ˇ
ˇ
d
“ X 2 ˇˇ tX 2 P EM,w0 ,z0 u,
where we drop tw “ w0 , z “ z0 u in the last equation by independence. Since EM,w0 ,z0 has no randomness, we
can conclude that X 2 is a truncated χ2 random variable as in (24). Then applying inverse survival function
transform and marginalizing over all possible values of pw0 , z0 q gives (27).

D.7

Proof of Proposition 3.8

r ` µM “ pX T P K c XM,g q´1 X T P K c µM “ 0, which is equivalent to
By proposition 3.6, βgM “ 0 ô X
M,g M,g
M,g M,g
M,g
T
K
K
M
XM,g
PM,g
“ 0.
c µM “ 0. So µM P C pXM,g c q ñ PM,g c µM “ 0 ñ βg
T
K
To show the opposite direction does not hold, we need to find µM such that XM,g
PM,g
c µM “ 0 and µM R
K
T
K
C pXM,gc q. Notice that µM P C pXM,gc q ñ PX,gc µM “ 0, so it suffices to find µM such that XM,g
PM,g
c µM “
K
T
K
0 and PX,gc µM ‰ 0. This reduces to showing nullpXM,g q X C pPM,gc q ‰ t0u, since if this is true, DµM P
T
K
T
K
R|M | , 0 ‰ γ P R|M | with XM,g
γ “ 0, PM,g
c µM “ γ. Notice that nullitypXM,g q “ |M | ´ |g| and rankpPM,g c q “
T
K
|M | ´ |g c | “ |M | ´ p ` |g|. Hence the dimension of nullpXM,g
q and C pPM,g
c q adds up to 2|M | ´ p ą |M |,
hence the intersection of the two spaces must be non-trivial.
Now if instead we assume M Ď M ˚ and βgM “ 0, then βg˚ “ 0 and u˚M “ 0 and the proof is completed
by noticing that µM “ XM,g βg˚ ` XM,gc βg˚c ` u˚M “ XM,gc βg˚c “ XM,gc βgMc .

26

D.8

Proof of Theorem 3.9

This proof is modified from the proof of Loftus & Taylor (2015, Theorem 4.1). We begin by noticing that
K
K
Psub , Psub
, Pfull , Pfull
are all orthogonal projection matrices. We also have
˙
ˆ
PM,gc PM,¨ 0|M |ˆpn´|M |q
“ Psub ,
Psub Pfull “
0pn´|M |qˆ|M |
Ipn´|M |q
K
K 2
and similarly Pfull Psub “ Psub . It follows that pPfull ´ Psub q2 “ Pfull ´ Psub and pPsub
´ Pfull
q “ pPfull ´
2
K
K
K
K
Psub q “ Pfull ´ Psub “ Psub ´ Pfull , which means Pfull ´ Psub and Psub ´ Pfull are both orthogonal projection
matrices. Then we have
K
K
K
K
K
K
}R1 }22 ´ }R2 }22 “ }Psub
y}22 ´ }Pfull
y}22 “ y T pPsub
´ Pfull
qy “ }pPsub
´ Pfull
qy}22 “ }R1 ´ R2 }22 .

Under the null, we have
ˆ
K
K
pPsub
´ Pfull
qµ “

K
PM,g
c µM
0

˙

ˆ
´

K
PM,
¨ µM
0

˙
“ 0,

which means the numerator of F is a centered random variable. Similarly, under the null, we have
˙
ˆ K
PM,¨ µM
K
“ 0,
Pfull
µ“
0
K
K
K
which means the denominator of F is also a centered random variable. Moreover, since pPsub
´Pfull
qPfull
“ 0,
the numerator and the denominator of F are orthogonal to each other. In the non-selective setup, F will be
a centered F random variable. Now let us notice that we can write

y “ rpg1 pFqw∆ ` g2 pFqw2 q ` z.
Plugging the above expression to Definition 2.2 gives the expression for EM,w∆ ,w2 ,z,r as in (32)(33). Then
we have
ˇ
ˇ
x “ M, w∆ “ w0 , w2 “ w0 , z “ z0 , r “ r0 u
F ˇˇ tM
∆
2
ˇ
ˇ
d
0
0
0 ,w 0 ,z ,r , w∆ “ w , w2 “ w , z “ z0 , r “ r0 u.
“ F ˇˇ tF P EM,w∆
∆
2
2 0 0
0
We claim that F K
K pw∆ , w2 , z, rq for the fixed M , so that we can drop tw∆ “ w∆
, w2 “ w20 , z “ z0 , r “ r0 u
in the R.H.S., and a similar argument as in the proof of Theorem 3.1 and Theorem 3.7 will give the desired
result. Next we establish the independence.
Independence between F and w∆ : We have shown R1 ´ R2 is orthogonal to, hence independent of R2
for fixed M . So w∆ K
K }R2 }22 . By Lemma D.4, w∆ K
K }R1 ´ R2 }22 . Hence F K
K w∆
Independence between F and w2 : The independence between pR1 ´ R2 , R2 q implies w2 K
K }R1 ´ R2 }22
2
and by Lemma D.4, w2 K
K }R2 }2 . Hence F K
K w2 .
K
K
Independence between F and z: Since Psub Psub
“ 0 and Psub Pfull
“ 0, we have z K
K pR1 , R2 q. Hence
K z.
FK
Independence between F and r: By Lemma D.4, r is complete and sufficient for σ. Under the null, F
is ancillary. Hence F K
K r by Basu’s theorem.

Appendix E
E.1

Implementation Details

Computing the Truncation Set for F Statistic

Recall that we want to solve the set
tF ě 0 : y “ rpg1 pFqw∆ ` g2 pFqw2 q ` z, y T QM,k y ` aTM,k y ` bM,k ě 0u,
27

where all the quantities are defined in Theorem 3.9. Plugging the expression for y to the quadratic constrataints, we can show the desired set are those F such that
2 T
2 T
2 T
rlooooooomooooooon
w∆ QM,k w∆ g12 pFq ` 2r
w∆ QM,k w2 g1 pFqg2 pFq ` rlooooooomooooooon
w2 QM,k w2 g22 pFq
loooooooomoooooooon
x11

x12

x22

T
T
T
T
` rp2w
∆ QM,k z ` aM,k w∆ q g1 pFq ` rp2w
2 QM,k z ` aM,k w2 q g2 pFq
looooooooooooooomooooooooooooooon
loooooooooooooomoooooooooooooon
x1

`

T

zlooooooooooooooomooooooooooooooon
QM,k z ` aTM,k z ` bM,k

x2

ě 0,

x0

where we define px11 , x12 , x22 , x1 , x2 , x0 q with a little abuse of notation. We then let
π
|g|F{p|M | ´ pq “ tan2 pθq, 0 ď θ ă .
2
With this change of parameterization and using eiθ “ cos θ ` i sin θ, we have
1 ´ cosp2θq
1
“ p2 ´ ei2θ ´ e´i2θ q,
2
4
1
`
cosp2θq
1
g22 pFq “
“ p2 ` ei2θ ` e´i2θ q,
2
4
1
i i2θ
g1 pFqg2 pFq “ sinp2θq “ ´ pe ´ e´i2θ q,
2
4
i iθ
g1 pFq “ ´ pe ´ e´iθ q,
2
1 iθ
g2 pFq “ pe ` e´iθ q.
2
g12 pFq “

Let t “ eiθ . Then the original inequality in F becomes
r2
r2
r2
r2
r2
r2
x11 ´ x12 i ` x22 qt2 ` p´ x11 ` x12 i ` x22 qt´2
2
4
4
4
4
4
r
r
r
r
´1
` p x2 ´ x1 iqt ` p x2 ` x1 iqt
2
2
2
2
r2
` px11 ` x22 q ` x0 ě 0.
2
We numerically solve the solution of the above polynomial in t with complex coefficients. Direct transformation from t to the original domain in F can be numerically unstable, so we transform the solution in t to
a small interval that contains the roots (in F), and then use the uniroot function in R to numerically solve
for the solution in F. This strategy was adopted by Loftus & Taylor (2015) and is numerically stable.
p´

E.2
E.2.1

Evaluating Truncated Survival Functions
The Naive Implementation

Let X be a one-dimensional random variable. Let F be its (non-truncated) CDF and let f be its (nontruncated) PDF. In our case, we will have a truncation set E Ď R, and we are interested in the truncated
random variable X|tX P Eu, whose PDF is proportional to f pxq ¨ 1tx P Eu. What we want is the survival
function (or CDF) of this truncated random variable. Explicitly, we want
PpX ě q | X P Eq “

PpX P E X rq, 8qq
,
PpX P Eq

where q is some quantile. The naive implementation to compute this quantity is to compute the numerator
and denominator separately and then take the fraction. Empirically, we find that under the null, this
strategy works quite well. However, when the truth is far away from the null, we find empirically that both
the numerator and the denominator will be extremely close to zero, so that taking the fraction causes severe
numerical issues, and sometimes R will give NaN.
28

E.2.2

The Approximation

Let us restate the problem explicitly. We want to compute
PpX P E1 q
PpX P E2 q
in a numerically stable way, where E1 , E2 are unions of intervals.
Truncated Normal Distribution
Once we solve this problem for the normal distribution, we will also have solved the problem for χ2 , t,
and F distributions since these can be efficiently and accurately approximated by normal distributions.
The strategy for the truncated normal distribution was first adopted in the selectiveInference package
(Tibshirani et al. 2017), where in their case, E1 , E2 are both a single interval. We extend their strategy to
the more general case where E1 , E2 are both unions of intervals.
Let Z „ N p0, 1q. WLOG we can assume E1 and E2 are unions of disjoint intervals with positive endpoints.
By Equation p14q in Bryc (2002), for z P r0, 8q, we have
PpZ ą zq « hpzq expt´z 2 {2u,
hpzq “ ?

2πz 3

z 2 ` 5.575192695z ` 12.77436324
.
` 14.38718147z 2 ` 31.53531977z ` 2 ˆ 12.77436324

Recall that we want
ř
Ppli ď Z ď ui q
PpZ P E1 q
“ ř iPI
PpZ P E2 q
jPJ Pplj ď Z ď uj q
ˆ
˙
ř
2
2
hpu
q
expt´u
{2u
´
hpl
q
expt´l
{2u
i
i
i
i
iPI
ˆ
˙,
«
ř
2
2
jPJ hpuj q expt´uj {2u ´ hplj q expt´lj {2u
where I, J are finite index sets. In order to avoid the situation when expt´u2i {2u is too close to zero, we
introduce a constant C, and write
ˆ
˙
ř
2
2
hpu
q
expt´u
{2
`
Cu
´
hpl
q
expt´l
{2
`
Cu
i
i
i
i
iPI
PpZ P E1 q
ˆ
˙.
«
ř
PpZ P E2 q
2
2
jPJ hpuj q expt´uj {2 ` Cu ´ hplj q expt´lj {2 ` Cu
To ensure the numerical stability, we loop C over the 0%, 20%, 40%, 60%, 80%, 100% quantiles of the set
tu2k {2, lk2 {2 : k P I Y Ju, and empirically this choice of C gives numerically stable results.
Truncated χ2 Distribution
Canal (2005) gives a normal approximation for the χ2 distribution. Consider X „ χ2r . Then the random
variable
1 X
1 X
X
LpXq :“ p q1{6 ´ p q1{3 ` p q1{2
r
2 r
3 r
is approximated normal distributed with mean
5
1
7
25
´
´
`
2
6 9r 648r
2187r3
and variance

1
1
37
`
´
.
18r 162r2
11664r3
Taking derivatives, one can easily show that L is a monotone increasing function of X. So PpX ě qq “
PpLpXq ě Lpqqq. Then we can use the approximation for the normal distribution to compute the desired
quantity.
29

Truncated t Distribution
Ling (1978) studies the accuracy of approximation formulas for the t, χ2 , and F CDFs. According to this
work, the normal approximation of the t distribution given by Peizer & Pratt (1968) performs well in most
situations. Consider X „ tr . By Ling (1978, Equation (2.6)), we have
d
˙
ˆ
1
logp1 ` q 2 {rq
2
q¨
,
PpX ě qq « P Z ě pr ´ `
3 10r
r ´ 5{6
where Z „ N p0, 1q. Then we can again use the approximation for truncated normal distribution to compute
the desired quantity.
Truncated F Distribution
Again by Ling (1978), the normal approximation for F distribution given by Peizer & Pratt (1968) performs
well in most situations. Consider X „ Fr1 ,r2 . By Ling (1978, Equation (4.3)), we have
g
f
ˆ
f 1 ` Q ¨ gp NSP q ` P ¨ gp NRQ q ˙
,
PpX ě qq « P Z ě D ¨ e
pN ` 16 qP Q
where Z „ N p0, 1q, and
r1 ´ 1
r2 ´ 1
r1 ` r2 ´ 2
,S “
,N “
,
2
2
2
r2
P “
, Q “ 1 ´ P,
r1 q ` r2
1
Q
P
Q ´ 0.5
1
`
q,
D “ S ` ´ pN ` qP ` 0.04p ´
6
3
r2
r1
r1 ` r2
$ 1´x2 `2x log x
,
0ăxă1
’
’ 2 1´x2
’
& ´x `1`2x log x
, xą1
x2 ´1
gpxq “
.
’1,
x“0
’
’
%
0,
x“1
R“

Using the approximation for the truncated normal distribution gives the desired quantity.

References
Atkinson, A. (1986), ‘[influential observations, high leverage points, and outliers in linear regression]: Comment: Aspects of diagnostic regression analysis’, Statistical Science 1(3), 397–402.
Atkinson, A. C. & Atkinson, A. C. (1985), Plots, transformations and regression; an introduction to graphical
methods of diagnostic regression analysis, Technical report.
Belsley, D. A., Kuh, E. & Welsch, R. E. (2005), Regression diagnostics: Identifying influential data and
sources of collinearity, Vol. 571, John Wiley & Sons.
Bi, N., Markovic, J., Xia, L. & Taylor, J. (2017), ‘Inferactive data analysis’, arXiv preprint arXiv:1707.06692
.
Bloomfield, P. & Steiger, W. L. (1984), ‘Least absolute deviations: theory, applications and algorithms’.
Brownlee, K. A. (1965), Statistical theory and methodology in science and engineering, Vol. 150, Wiley New
York.
Bryc, W. (2002), ‘A uniform approximation to the right normal tail integral’, Applied mathematics and
computation 127(2), 365–374.
30

Canal, L. (2005), ‘A normal approximation for the chi-square distribution’, Computational statistics & data
analysis 48(4), 803–808.
Cook, R. D. (1977), ‘Detection of influential observation in linear regression’, Technometrics 19(1), 15–18.
Daniel, C. & Wood, F. S. (1999), Fitting equations to data: computer analysis of multifactor data, John
Wiley & Sons, Inc.
Fithian, W., Sun, D. & Taylor, J. (2014), ‘Optimal inference after model selection’, arXiv preprint
arXiv:1410.2597 .
Hadi, A. (1990), A stepwise procedure for identifying multiple outliers in linear regression, in ‘American
Statistical Association Proceedings of the Statistical Computing Section’, Vol. 137, p. 142.
Hoeting, J., Raftery, A. E. & Madigan, D. (1996), ‘A method for simultaneous variable selection and outlier
identification in linear regression’, Computational Statistics & Data Analysis 22(3), 251–270.
Huber, P. & Ronchetti, E. (1981), ‘Robust statistics, ser’, Wiley Series in Probability and Mathematical
Statistics. New York, NY, USA, Wiley-IEEE 52, 54.
Lee, J. D., Sun, D. L., Sun, Y., Taylor, J. E. et al. (2016), ‘Exact post-selection inference, with application
to the lasso’, The Annals of Statistics 44(3), 907–927.
Ling, R. F. (1978), ‘A study of the accuracy of some approximations for t, χ2, and f tail probabilities’,
Journal of the American Statistical Association 73(362), 274–283.
Loftus, J. R. & Taylor, J. E. (2015), ‘Selective inference in regression models with groups of variables’, arXiv
preprint arXiv:1511.01478 .
Peizer, D. B. & Pratt, J. W. (1968), ‘A normal approximation for binomial, f, beta, and other common,
related tail probabilities, i’, Journal of the American Statistical Association 63(324), 1416–1456.
R Core Team (2017), R: A Language and Environment for Statistical Computing, R Foundation for Statistical
Computing, Vienna, Austria.
URL: https://www.R-project.org/
Reid, S., Tibshirani, R. & Friedman, J. (2013), ‘A study of error variance estimation in lasso regression’,
arXiv preprint arXiv:1311.5274 .
Shao, J. (2003), ‘Springer texts in statistics’.
She, Y. & Owen, A. B. (2011), ‘Outlier detection using nonconvex penalized regression’, Journal of the
American Statistical Association 106(494), 626–639.
Taylor, J. & Tibshirani, R. J. (2015), ‘Statistical learning and selective inference’, Proceedings of the National
Academy of Sciences 112(25), 7629–7634.
Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical
Society. Series B (Methodological) pp. 267–288.
Tibshirani, R. J., Taylor, J. et al. (2011), ‘The solution path of the generalized lasso’, The Annals of Statistics
39(3), 1335–1371.
Tibshirani, R., Tibshirani, R., Taylor, J., Loftus, J. & Reid, S. (2017), selectiveInference: Tools for PostSelection Inference. R package version 1.2.2.
URL: https://CRAN.R-project.org/package=selectiveInference
Welsch, R. E. & Kuh, E. (1977), ‘Linear regression diagnostics’.
Yang, F., Barber, R. F., Jain, P. & Lafferty, J. (2016), Selective inference for group-sparse linear models, in
‘Advances in Neural Information Processing Systems’, pp. 2469–2477.
Zhao, J., Leng, C., Li, L., Wang, H. et al. (2013), ‘High-dimensional influence measure’, The Annals of
Statistics 41(5), 2639–2667.
31

