MAGIC: a general, powerful and tractable
method for selective inference

arXiv:1607.02630v1 [math.ST] 9 Jul 2016

Xiaoying Tian, Nan Bi and Jonathan Taylor∗
Department of Statistics
Stanford University
Sequoia Hall
Stanford, CA 94305, USA
e-mail: xtian@stanford.edu
nbi@stanford.edu
jonathan.taylor@stanford.edu
Abstract: Selective inference is a recent research topic that tries to perform valid inference after using the data to select a reasonable statistical model. We propose MAGIC,
a new method for selective inference that is general, powerful and tractable. MAGIC
is a method for selective inference after solving a convex optimization problem with
smooth loss and `1 penalty. Randomization is incorporated into the optimization problem to boost statistical power. Through reparametrization, MAGIC reduces the problem
into a sampling problem with simple constraints. MAGIC applies to many `1 penalized
optimization problem including the Lasso, logistic Lasso and neighborhood selection
in graphical models, all of which we consider in this paper.
Keywords and phrases: selective inference, statistical power, sampling.

1. Introduction
There are a great deal of sophisticated statistical learning methods that allow us to
search through a large number of models and look for meaningful patterns. Having
done this search, we naturally want to judge the apparent associations that have been
found. This has spawned a new area of research called selective inference [5, 6, 14,
15]. Loosely speaking selective inference recognizes the inherent selection biases in
reporting the most “significant” results from various statistical models and attempts to
adjust for the bias.
At a high level, selective inference involves two stages: First, solve a convex optimization problem, usually some penalized loss function. Second, perform inference in
the statistical model suggested by the result of the optimization problem. For example,
we first use the data to solve the Lasso problem, and then want to form confidence
intervals for the variables that are nonzero in the Lasso solution. Adjustment for selection results in some constraints on the underlying distribution. Although various
such problems have been studied, most of the papers only focus on one specific optimization problem. This is necessary as different loss functions in the optimization
problems result in different geometry of the constraints. In this paper, we introduce
a method, called ”MAGIC”, Monte-carlo Algorithm for General Inference with Constraints, which provides valid selective inference for optimization problems with any
∗ Supported

in part by NSF grant DMS 1208857 and AFOSR grant 113039.
1

Tian et al./MAGIC

2

smooth loss functions. The advantage of MAGIC compared to previous selective inference methods are generality, statistical power and tractability. We elaborate each in the
following passage:
Generality: The generality of MAGIC lies in two aspects: arbitrary smooth loss
function in the penalized optimization problem and the data distribution from any exponential family. In comparison, the authors in [6] considered only inference after solving Lasso; The work [5] considered some exponential families with simple selection
rules, but also noticed the difficulty for inference after solving more complex optimization problems. Finally, the work [13] is the closest in generality to this work, but shows
substantially weaker statistical power, which we discuss below.
Statistical power: Earlier work [6, 13] has provided valid inference after selection,
but sometimes lacks power. Other work [15, 5] suggested introducing randomness in
the optimization algorithm which produces much improved power. This is the approach
we take in this work. In simulation, we show that MAGIC produces much improved
power over [6, 13].
Tractability: The earlier work [6, 13] computes valid p-values in closed form, thus
involving the least computation cost. The framework in [5, 15] involves sampling in
a constrained subset in the sample space. Both work used hit-and-run algorithm proposed in [1], which is a method to generate distributions on a subset of the space. The
constrained subsets as described in [5, 15, 6, 13] can be quite complicated and depend
on the particular loss function. Algorithms that do not use MCMC, such as [6, 13], do
not suffer from this problem too much as they only need to compute the boundary once,
but the methods in [5, 15] need to compute the boundary at each step of simulation, resulting in much more computation cost. MAGIC, however, transforms the constrained
subset to a canonical set through reparametrization, removing the computational cost
involved in computing the boundary at each step of sampling. Thus it is more tractable
than [5, 15].
In Section 2, we introduce the general form of our randomized optimization problem, and describe the inference method as well as theory for MAGIC. Section 3 gives
applications of MAGIC to different statistical learning problems. To demonstrate the
applicability of MAGIC, we give three distinct examples: Lasso, `1 penalized logistic
regression and neighbourhood selection [8], which are applied in regression, classification and Gaussian graphical models respectively. Section 4 includes the comparisons of
MAGIC with existing selective inference methods both in terms of statistical power and
confidence intervals. All proofs are collected in Section 5 and the sampling methods
are covered in Section 6. We conclude with discussions about future work in Section 7
1.1. Related works
Most of the theoretical work on high-dimensional data focuses on consistency, either
the consistency of solutions [11, 17] or the consistency of the models [19, 21].
In the post selection literature, the authors in [2] proposed the PoSI approach, which
reduce the problem to a simultaneous inference problem. Because of the simultaneity,
it prevents data snooping from any selection procedure, but also results in more conservative inference. In addition, the PoSI method has extremely high computational

Tian et al./MAGIC

3

cost, and is only applicable when the dimension p < 30 or for very sparse models.
The authors [10] proposed a method for computing p-values that controls false discovery rate (FDR) among all variables. This is quite different from the hypothesis testing
framework of this work, as the hypotheses tested in selective inference are chosen as
a function of the data. Hence, the hypotheses tested are not directly comparable. Furthermore, compared with [10], MAGIC has the advantage of being able to construct
confidence intervals for the selected variables.
2. Randomized selective inference
2.1. A randomized selection algorithm
Many statistical learning problems can be cast as convex optimization problems. Specifically, we solve the following randomized convex optimization.
`(β; S) + λkβk1 − ω T β,
β̂ = minimize
p
β∈R

(1)

where data S ∼ F , ` can be the negative log-likelihood for F , but generally just needs
to be some convex loss function in β, the randomization variable ω ∼ G, a distribution
on Rp independent of F , λ, is fixed. Without randomization, that is to set G = δ0 , the
point mass at 0, (1) includes many classical statistical learning problems, e.g. lasso [16],
elastic net [22], `1 penalized logistic regression, neighbourhood selection [8]. Although
it might seem strange to add noise to data for model selection, it is seen in other forms
in literature and applications. Common use of data splitting is an example [4, 20], as
a random subset of data is used for model selection. The form of our randomization is
also related to [9]. We can control the amount of randomization through the variance
of G, usually just a little randomization will produce much improved power.
We define the variable selection map as
Ê(S, ω) = supp(β̂(S, ω)).
For the observed data Sobs and an instance of ωobs both considered fixed, we define
E = Ê(Sobs , ωobs ) which is the active set of (1) and consider it fixed hereafter.
After having solved the above problem, we now consider inference for parameters
chosen on the basis of this set of non-zero coefficients E. Suppose the data S ∼ F
is a member of an exponential family with parameters b ∈ Rp and sufficient statistics
T (S) ∈ Rp . In particular, its density fb (s) has the following form,
dfb
(s) = exp(bT T (s) − A(b))
dµ
where µ is the reference measure on the sample space of S and A is the normalizing
constant with µ, A known. Having observed a set of selected variables E, we can and
often do then consider a submodel of the above model with b−E = 0. If E ⊇ supp(b),
then our model is correctly specified. This is the scenario we always consider hereafter. For treatment of misspecified models, see [6, 5]. Under this submodel, the joint
distribution of (S, ω) is fully specified. Our target of inference is now bE .

Tian et al./MAGIC

4

Since E is not given a priori, but selected by the data, it seems to be only fair to
consider (S, ω) such that Ê(S, ω) = E. This is equivalent to condition on the event
{(S, ω) : Ê(S, ω) = E}. This is the general approach taken in [6, 5, 15] to provide
valid (selective) inference in the above model.
Let A be the region where {(S, ω) ∈ A} ⇐⇒ {Ê(S, ω) = E}, then this general
approach to selective inference requires us to describe the conditional distribution
S | (S, ω) ∈ A,

(S, ω) ∼ F × G.

(2)

We first state the following result,
Theorem 1. Suppose (S, ω) ∼ F × G, F is the exponential family specified above,
with the parameters b satisfying supp(b) ⊆ E. G is a distribution on Rp and A is
defined as above. Then for any variable j ∈ E, there exists a p-value function Pj :
supp(F ) → [0, 1], such that
PF ×G [Pj (S; A) ≤ α | (S, ω) ∈ A] ≤ α,

(3)

under the null hypothesis H0j : bj = 0. The function Pj only depends on data S and
A.
In some cases, equality holds in (3), we will discuss the conditions in the proof. In
this case, the test proposed above is the Uniformly Most Powerful Unbiased test [5],
providing theoretical ground for the power of MAGIC. Theorem 1 gives a construction of the p-value, which we can use to reject the null hypothesis at level α. We will
give the exact construction of Pj in the proof of Theorem 1, which is an multivariate
integral and is hard to compute in general. We instead try to acquire samples from (2)
and approximate the multivariate integral. The constrained region A is the bottleneck
for the sampling, which is complicated and specific to the loss function `. However,
through a reparametrization, we can form the constrained region as a simple set that is
independent of `.
2.2. Augmented parameter space
Once we solve the optimization (1), we get β̂ the solution and ẑ the subgradient of
kβ̂k1 . β̂, ẑ are functions of (S, ω). We formally define the optimization map as follows:
θ̂

(s, ω) 7→ (s, β̂(s, ω), ẑ(s, ω)) ∈ S F (`),

(4)

where
def

S F (`) =




(s, β, z) : s ∈ supp(F ), `(β; s) < ∞, kβk1 < ∞, z ∈ ∂kβk1 .

(5)

S F (`) is the set of possible values (s, β, z) where there will be a solution to (1). We call
S F (`) the augmented parameter space. Note β̂ and ẑ are random variables (through the
randomness in (S, ω)). One way to describe their distribution, is to find the inverse of
the map θ̂ to reconstruct ω from (S, β̂, ẑ). In the following passage, we denote β̂, ẑ as
the random variables and β, z as the corresponding integration variables when writing
out the density.

Tian et al./MAGIC

5

2.3. Reconstruction and description of the constrained set
Let β̂E , ẑE be β̂ and ẑ restricted to E, and ẑ−E the subgradients restricted to E c . To
make the notation easier, we define the gradient map γ
γ(s, β) = ∂β `(β; s).
Lemma 2. Through the reparametrization in the optimization map (4), the selection
event {Ê(S, ω) = E} is equivalent to
(
γ(S, β̂) + λ · ẑ − ω = 0,
(6)
sign(β̂E ) = ẑE , kẑ−E k∞ ≤ 1.
Lemma 2 provides a reconstruction of ω using S, β̂ and ẑ. The reconstruction map
is defined as,
def
ψ(s, β, z) = (s, γ(s, β) + λ · z) = (s, ω).
It is thus easy to see that the distribution of (S, β̂, ẑ) follows satisfies the following
distributional condition,
(S, γ(S, β̂) + λ · ẑ) ∼ F × G.
Moreover, using Lemma 2, the distribution for inference (2) can be rewritten as
S | (β̂(S, ω), ẑ(S, ω)) ∈ B,
(S, ω) ∼ F × G,


B = β̂−E = 0, sign(β̂E ) = ẑE , kẑ−E k∞ < 1 .
Note that B is a much nicer set than A in that it only requires β̂E to be in a certain
quadrant and kẑ−E k∞ < 1.
Combining these two observations above, we have the following theorem. We denote by TE\j ∈ R|E|−1 the sufficient statistics T ∈ Rp restricted to the set E − {j}.
Theorem 3 (Sampling for MAGIC). Through change of variables (4), the law for
selective inference (2) is equivalent to
S|(β̂, ẑ) ∈ B,

(S, γ(S, β̂) + λ · ẑ) ∼ F × G.

(7)

Moreover, suppose F , G has densities f and g respectively, the joint distribution of
(S, β̂, ẑ) has density proportional to
f (s) · g(γ(s, β) + λ · z) · Jψ(s, β, z) · 1B (β, z)

(8)

with the Jacobian denoting the determinant of the derivative of the map ψ with respect
to (β, z) on the fiber over s.
Furthermore, assuming the assumptions in Theorem 1, Pj (S) can be computed (with
approximation) with samples from (7) and further conditional on the sufficient statistics
TE\j (S).

Tian et al./MAGIC

6

Theorem 3 gives the explicit density of the law (7) up to a constant. In the proof
we specify how to use the samples from (7) to approximate the p-value function Pj .
A natural choice of sampling would be the Metropolis-Hastings method or perhaps
the projected Langevin method [3]. To condition on TE\j (S), we just need to make
sure the proposal does not move TE\j (S) in each step. Such choice of the proposal is
usually natural, for examples see Section 6. The boundary constraint is B, which needs
small adjustment from the original Metropolis-Hasting method. Detailed description is
in Section 6. After acquiring such samples, we can use them to approximate the p-value
function Pj in Theorem 1.
All the previous work on selective inference also conditions on the observed signs
ẑE .


S | (S, ω) ∈ A, sign β̂E (S, ω) = zE,obs
where zE,obs = sign(β̂(Sobs , ωobs )) is considered fixed. The work [6] explains that any
inference valid under this law, would be valid under (2). Note the additional constraint
simply requires β̂E to be in the quadrant specified by zE,obs . In what follows, we also
condition on ẑE .
3. Examples
3.1. Randomized Lasso
Consider linear regression setting where data y ∼ N (Xb, σ 2 I), X ∈ Rn×p is fixed,
σ 2 is known. Instead of solving the original Lasso proposed by [16], we solve the
following randomized version of it,
1
β̂ = argmin ky − Xβk22 + λkβk1 − ω T β.
p
2
β∈R

(9)

The gradient of the loss γ(y, β̂) = −X T (y − X β̂) and ẑ is the subgradient for kβ̂k1 .
The reconstruction map ψ(y, β, z) = y, λ·z −X T (y −Xβ) . Suppose E is the active
set of (9), then we model the data by F = N (XE bE , σ 2 I), S = y and bE is the target
for inference.
Corollary 4 (Randomized Lasso sampler). If E ⊇ supp(b), then conditioning on
(E, zE,obs ), the joint distribution of (y, β̂E , ẑ−E ) can be used for inference (for bE ).
Further, it has density proportional to

  


ky − XE bE k22
zE,obs
T
exp −
·g λ
− X (y − XE βE )
(10)
z−E
2σ 2
supported on sign(βE ) = zE,obs and kz−E k∞ < 1.
We thus can obtain samples (y, β̂E , ẑ−E ) for any bE in the null hypothesis and use
Theorem 1 and Theorem 3 to construct valid p-values. Detailed algorithm is specified
in Section 6.

Tian et al./MAGIC

7

3.2. L1-penalized logistic regression
In practice, many statistical learning problems are classification problems, e.g. spam
iid
classification, tumor analysis, etc. Suppose xi ∼ FX , xi ∈ Rp , yi |xi ∼ Bernoulli(xTi b),
FX is unknown and p fixed, S = (X, y). The logistic loss is
#
" n
1 X
`(β) = − √
yi log π(xi β) + (1 − yi ) log(1 − π(xi β)) ,
n i=1
where π(x) = exp(x)/(1 + exp(x)). The randomized logistic regression solves the
following problem

β̂ = argminβ∈Rp `(β) + λkβk1 − ω T β + kβk22
(11)
2
with  > 0 small and fixed. The addition of the term with  is to ensure the existence of
the solution to (11). We explicitly express the  term, but in general it can be absorbed
into the loss function.
Suppose E is the active set of (11), then bE is the target of inference. With slight
abuse of notation, we allow π : Rn → Rn , x 7→ π(x) to be the function applied on
each coordinate of x ∈ Rn . With some algebra, we have the reconstruction map for ω
h
i
1
ω = λ · ẑ − √ X T y − π(X β̂) + β̂
n
To sample (X, y) jointly is not feasible when FX is unknown. Two observations help
us circumvent it and even make the sampling more efficient. First, the reconstruction
map for ω only involve the random vector
1
∇`(β̂E ) = − √ X T (y − π(XE β̂E ))
n

1 T
≈ − √ X y − π(XE β̄E ) − W (XE β̄E )XE (β̂E − β̄E )
n
where β̄E is the MLE for the unpenalized logistic regression with only the variables in
E and W (Xβ) = diag(π(Xβ)(1 − π(Xβ))) is the weight matrix. Alternatively, we
might take β̄E to be the one-step estimator in the selected model starting from β̂E [13].
The gradient ∇`(β̂E ) can be reconstructed, up to a Taylor remainder, from β̂ and the
random vector


β̄E
T =
.
T
X−E
(y − π(XE β̄E ))
Moreover, when p is fixed, pre-selection, our random vector T properly scaled is
asymptotically normal and when the selected model is correct (E ⊇ supp(b)):

 
1
b
D
√ T− E
→ N (0, Σ)
(12)
0
n
where Σ is estimable from the data [13]. Since asymptotically T is from an exponential
family with parameters bE , Theorem 1 states the p-value is a function of T only. Thus
instead of sampling (X, y), we only need to sample the distribution T .

Tian et al./MAGIC

8

Theorem 5. Suppose E ⊇ supp(b) and conditioning on (E, zE,obs ), the joint distribution of (T, β̂E , ẑ−E ) can be used for inference. Then the distribution of (T, β̂E , ẑ−E )
asymptotically (with p fixed, n → ∞) has density



T
1
XE
W (XE β̄E )XE (β̂E − β̄E )
φ(T ) · g √
T
T
n X−E
W (XE β̄E )XE (β̂E − β̄E ) − X−E
(y − π(XE β̄E ))
(13)


 
zE,obs
β̂E
+λ
+
z−E
0
where φ is the density for N ((bE , 0), Σ).
3.3. Neighborhood Selection
Gaussian graphical models have recently become a very popular way to study network
structures. In particular, it has often been used on many types of genome data (e.g. gene
expression, metabolite concentrations etc.) Suppose the data we observe is X ∈ Rn×p ,
where each row of X is independently distributed as N (µ, Σ), µ ∈ Rp , Σ ∈ Rp×p .
It is of interest to study the conditional independence structure of the variables
{1, 2, . . . , p}. The conditional independence structure is conveniently represented by
an undirectional graph (Γ, E), where the nodes Γ = {1, 2, . . . , p}, and there is an edge
between (i, j) if and only if xi 6⊥ xj conditional on all the other variables Γ\{i, j}.
Moreover, assuming the covariance matrix Σ is not singular, we denote the inverse
covariance matrix Θ = Σ−1 , then
xi ⊥ xj |XΓ\{i,j} ⇐⇒ Θij = 0.
In many applications of Gaussian graphical models, we assume the sparse edge structure, where we can hope to recover the edgeset E even when n < p2 . The authors in
[8] proposed neighborhood selection with the Lasso to achieve this goal. The algorithm
can be formulated as the following optimization problem, for any node i

(14)
β̂ i,λ = argmin n−1 kxi − Xβk22 + λkβk1 ,
β:βi =0

where xi is the i-th column of X, λ is chosen according to Chapter 3 of [8] and considered fixed. Denote B̂ = (β̂ 1 , β̂ 2 , . . . , β̂ p ), we propose the randomized version of
(14),
B̂ = argmin kX − XBk2F + λkBk1 − ΩB,
(15)
B:Bii =0

i.i.d

where Ω = (ω 1 , . . . , ω p ), ω i ∼ G. Let E i = supp(β̂ i ), and E = (E 1 , . . . , E p ).
Since E is usually not symmetric, we instead look at the set
E ∨ = {(i, j)|Eij = 1 or Eji = 1}.
Our target for inference is {Θij , (i, j) ∈ E ∨ }. Note (15) is the matrix form of (1), and
the reconstruction maps are decomposable across the p nodes; Therefore, we have the
following corollary,

9

Tian et al./MAGIC

Corollary 6. Suppose E is the active set for (15), and ẑE is the corresponding signs
of B̂E , then conditioning on (E, zE,obs ), the distribution of (X, B̂E , ẑ−E ) can be used
for inference. Furthermore, if we assume Θij = 0, i 6= j and (i, j) 6∈ E ∨ , then the
joint distribution of (X, B̂E , ẑ−E ) has the following density,


p
X
X
1
Θii kxi k2 +
Θij xTi xj 
exp −
2 i=1
(i,j)∈E ∨
(16)

Y  z i 
T
i
T
·
g λ Ei ,obs − X−i
· det(XE
(xi − XE i βE
i)
i XE i ).
z−E i
i∈Γ

4. Simulation
Theorem 1 states that our p-values should be valid at level α, for any α ∈ [0, 1], see (3).
In fact, all the three examples above satisfy the condition such that the Type-I error for
any level-α test would be equal to (or asymptotically equal to) α. That is equivalent as
saying the p-values follow Unif(0, 1) distribution. To validate Theorem 1 and Theorem
3, we ran the following simulations for each of the examples in Section 3. Our data is
generated as follows, for Lasso,
y ∼ N (Xb, σ 2 I), X ∈ Rn×p , fixed, kbk0 = s,
where s  p. The framework works for arbitrary n and p. To demonstrate the applicability of our framework in high dimensions, we set n = 50, p = 100, s = 7. For
logistic Lasso problem,
xi ∼ N (0, I),

yi |xi ∼ Bernoulli(π(xi b)),

π=

exp(x)
,
1 + exp(x)

kbk0 = s.

The framework for logistic regression is fixed p and n → ∞. Thus we take n =
500, p = 50, s = 5. For both of the examples above, the signal to noise ratio (snr) is
7. Finally, for neighborhood selection, the data matrix is X ∈ Rn×p , each row of X is
i.i.d from N (0, Θ−1 ). We take n = 100, p = 30, note this is a high-dimensional setting
since we have 30 × 30 unknown parameters. But only 1% of off-diagonal elements of
Θ is non-zero, and the non-zero off-diagonal entries of Θ are taken to be ρ = 0.245
and the diagonal elements are 1. ρ = 0.245 is chosen because any value less than 0.25
would ensure Θ is positive definite [8].
For each j ∈ E, we test the hypothesis H0j : bj = 0, against a two-sided alternative hypothesis. We call the p-values the null p-values when the null hypothesis is
true and alternative p-values otherwise. When the active set E (or E ∨ ) from the problem covers supp(b) (or supp(Θ)), the null p-values should follow Unif(0, 1). Figure
1 is the plot for the empirical cdf for the null p-values computed from Lasso, logistic
Lasso and neighborhood selection. We see that all the null-pvalues follow the uniform
distribution, verifying our Theorem 1 and Theorem 3.

10

Tian et al./MAGIC

lasso

1.0

logistic lasso

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0
0.0

neighborhood selection

0.2

0.4

0.6

0.8

1.0

Fig 1: Empirical cdf of null p-values, generated from 100 instances of each problem.
We use Laplace noise for randomization.
4.1. Comparisons of statistical powers
As we mentioned in Section 1, randomization significantly boosts power. This is shown
in both hypothesis testing and confidence intervals. We describe what it means in both
aspects. For a valid selective level-α test, Type-I error is controlled at α conditional
on selection. We hope to achieve valid tests with high power. In the selective inference
framework, statistical power is simply defined as the power in the selected model [5,
15]. If E ⊇ supp(b), then for any j ∈ E,


power = P reject H0j | H1j is true, E is selected .
The selective inference framework also offers confidence intervals by inverting a valid
test, for examples, see [6]. We want short confidence intervals which have the desired
coverage guarantees. MAGIC enjoys higher statistical power (shorter intervals), the
tradeoff is slightly worse selected models as we added randomization for model selection. However, the tradeoff is highly in favor of MAGIC. Usually just a small amount
of randomization will dramatically increase statistical power. In the linear regression
case, this has been shown in [15] with simulated data. In the following passage, we
give numerical comparisons on both a real dataset and simulated data.
4.1.1. In vitro HIV drug resistance
In [12], the authors study the genetic basis of drug resistance in HIV, using markers
of inhibitor mutations to predict a quantitative measurement of susceptibility to several antiretroviral drugs. The hope is to find the mutations highly correlated with the
susceptibility to drugs. We apply Lasso to the protease inhibitor subset of their data
and select the potential mutations set for one of the drugs, Lamivudine (3TC). We then
compute the OLS estimator in the selected set of gene mutations, and form confidence
intervals for the coefficients (Figure 2). The grey bars are the OLS estimates with only
the selected mutations. The confidence intervals are adjusted for selection and should
have the desired coverage 90%. We report the estimators together with the confidence
intervals. The procedure in left panel 2a is the same as [6] without randomization in

11

Tian et al./MAGIC

Fig 2: Confidence intervals for selected genes in 3TC DATA
3

MAGIC

Parameter values

2

1

0

1

(a) Selective intervals without randomization

P219R

P215F

P215Y

P190A

P210W

P184V

P181C

P118I

P151M

P90I

P115F

P77L

P83K

P69i

P75I

P65R

P67N

3

P62V

2

(b) Selective intervals with randomization

selecting the mutations. The right panel 2b in contrast uses the MAGIC framework
2
for LASSO with randomization ω ∼ N (0, 0.1σcv
), where σcv is the noise level estimated by cross-validation. Note the mutations selected by the two methods only differ
by 3 mutations, with small effects, and the OLS estimator for the common mutations
are very close. But the randomized selection procedure gives much shorter confidence
interval across all mutations, demonstrating the advantage and practicality of our methods.
4.1.2. Statistical power comparison with simulated data
In this section, we compare more specifically the tradeoff between power and model
selection using simulated data. The authors in [13] offered explicit calculations of
p-values after the model is selected by `1 penalized logistic regression or graphical
Lasso. Both examples can be considered in the MAGIC framework. Simulations in
[13] showed that graphical Lasso has worse power than `1 penalized logistic regression. Therefore, we compare our framework to the latter. We assume the same setup as
before, our randomization noise is ω ∼ N (0, 0.1σ 2 ) and  = 0.02. The proportions of
selecting the “true” models (E ⊇ supp(b)) is 0.91 without randomization and 0.852
in MAGIC. Much more different is the power of the two procedures; for a level-0.05
test, the statistical powers defined above is 0.176 without randomization and 0.887 in
MAGIC. Figure 3 is the histograms for the alternative p-values with or without randomization.
5. Proofs
5.1. Proof for Theorem 1
Proof. Let S be the space for S, then (S, ω) ∈ S × Rp . The joint distribution of
(S, ω) conditional on (S, ω) ∈ A has the following density with respect to the measure
µ(ds)G(dω)


exp bT T (s) 1{(s, ω) ∈ A}
h(s, ω) = R
.
(17)
exp [bT T (s)] 1{(s, ω) ∈ A}µ(ds)G(dω)
S×Rp

12

Tian et al./MAGIC

Fig 3: The alternative p-values computed from the MAGIC framework highly concentrated around 0, while without randomization the p-values are more evenly distributed between [0, 1], The statistical powers are 0.887 for MAGIC v.s. 0.176 for nonrandomized procedure with a level-0.05 test.
0.8

non-randomized
MAGIC

0.7

Frequencies

0.6
0.5
0.4
0.3
0.2
0.1
0.0
0.0

0.2

0.4
0.6
Alternative p-values

0.8

1.0

Since the denominator is merely a normalizing constant, (17) is also an exponential
distribution with parameters b, sufficient statistics T (S) and a slightly different reference measure 1{(s, ω) ∈ A}µ(ds)G(dω). Since E ⊇ supp(b), bT T (s) = bTE TE (s),
where bE , TE ∈ R|E| are b and T restricted to set E. Thus (17) can be seen as an exponential family with sufficient statistics TE and parameters bE . To test any hypothesis
H0j : bj = 0, j ∈ E, Chapter 4 of [7] states that Uniformly Most Powerful Unbiased
tests can be constructed using the statistic Tj and conditioning on all the other sufficient statistics TE\j ∈ R|E|−1 . Thus the conditional density of the one dimensional
distribution for Tj is
h
i R
exp bj tj + bTE\j tE\j · A 1{TE\j (s) = tE\j }µ(ds)G(dω)
h
i R
hj (tj ; tE\j ) =
exp bTE\j tE\j · A exp [bj Tj (s)] 1{TE\j (s) = tE\j }µ(ds)G(dω)
R
exp(bj tj ) A 1{TE\j (s) = tE\j }µ(ds)G(dω)
=R
exp [bj Tj (s)] 1{TE\j (s) = tE\j }µ(ds)G(dω)
A
(18)
Thus (18) is the density for the distribution
Tj (S)|TE\j (S), (S, ω) ∈ A, (S, ω) ∼ F × G.

(19)

Note (18) involves only the parameter bj , thus it can be used to test the composite
hypothesis H0j : bj = 0, with bE\j taking arbitrary values.
R tj
Let Hj denote the c.d.f of the above law: Hj (tj ; TE\j ) = −∞
hj (r; TE\j ) dr.
Then we can construct our function P̃j : Rp → R as
R
exp(bj Tj (s))1{Tj (s) > tj }1{TE\j (s) = tE\j }µ(ds)G(dω)
R
P̃j (t) = A
exp [bj Tj (s)] 1{TE\j (s) = tE\j }µ(ds)G(dω)
(20)
A
= 1 − Hj (tj ; tE\j ).

Tian et al./MAGIC

13

Under the null hypothesis, we take bj = 0, thus P̃j depends only on T (s) and A. We
define Pj (s) = P̃j (T (s)). Now we prove the level-α control (3). Note


PF ×G [Pj (S) ≤ α | (S, ω) ∈ A] = E PF ×G [Pj (S) ≤ α | TE\j , (S, ω) ∈ A] ,
it suffices to prove the quantity inside the expectation has the level-α control for any
TE\j . Since Hj is the c.d.f of the conditional law (19),
PF ×G [Pj (S) ≤ α | TE\j , (S, ω) ∈ A]
=PF ×G [1 − Hj (Tj ; TE\j ) ≤ α | TE\j , (S, ω) ∈ A]
=PF ×G [Tj ≥ Hj−1 (1 − α) | TE\j , (S, ω) ∈ A]
=1 − Hj [Hj−1 (1 − α)] ≤ α,
where Hj−1 generalized inverse for Hj , the equality holds when Hj is strictly increasing in tj .
5.2. Proof for Lemma 2
Proof. Equation (1) is a convex optimization problem. The solution β̂ and subgradient
of the `1 norm ẑ satisfy the KarushKuhnTucker conditions (KKT), which are sufficient
and necessary.
(
∂β `(β̂; S) + ẑ − ω = 0,
ẑ ∈ ∂kβ̂k1 .
The equations are simply the differentiation of the optimization objective function.
This gives the equation part in (6) of the lemma. Note the penalty term kβ̂k1 is differentiable except at 0, its subgradient at 0 is [−1, 1]. Thus, conditioning on the active set
Ê(S, ω) = E it is equivalent to:
(
ẑj = sign(β̂j ), ∀ j ∈ E,
|ẑj | ≤ 1, ∀ j 6∈ E.
Combining the above two, we have the conclusion of the lemma.
5.3. Proof for Theorem 3
Proof. Per the discussion above Theorem 3, it is not hard to see the distributional constraint on (S, β̂, ẑ) is that γ(S, β̂) + λẑ ∼ G and is independent of S. Moreover, (β̂, ẑ)
are constrained to be in the region B. Thus the law (7) is the marginal law of S conditional on selection.
Now we investigate the joint density of (s, β, z). Through the reconstruction map,


ψ(s, β, z) = s, γ(s, β) + λz ,

Tian et al./MAGIC

14


the density of (s, β, z) is simply the product f (s)g γ(s, β)+λz times the determinant
of the Jacobian matrix. Standard multivariate calculus yields the form of the Jacobian
matrix of ψ as follows,


I
0
Ds ψ(s, β, z) D(β,z) ψ(s, β, z)
with determinant det D(β,z) ψ(s, β, z). Thus we have (8).
Notice the construction (20), Pj (s) = 1−Hj (Tj (s), TE\j (s)) and Hj is the c.d.f for
the conditional distribution (19). It is equivalent to sampling (7) while further conditional on TE\j (S). After we acquire m samples {S (1) , . . . , S (m) }, we can approximate
the integral in (20) as the percentile of Tj (Sobs ) among {Tj (S (1) ), . . . , Tj (S (m) )}.
5.4. Proof for Corollary 4
Proof. Notice that once we condition on the active set E, and the signs zE,obs ,
 


βE
zE,obs
β=
,
z=
.
0
z−E

(21)

Therefore, the density of (y, β, z) is equivalent to that of (y, βE , z−E ), through the
construction in (21). Note the Jacobian matrix is

 T
XE XE
0
DβE ,z−E ψ(y, βE , z−E ) =
,
T
XE λI
X−E
where I is the identity matrix of dimension p − |E|. Then the Jacobian Jψ(y, β, z) =
T
XE ). Since the Jacobian is a constant only depending on X. Thus, if
λp−|E| det(XE
we plug in β, z in Theorem 3, the density of (y, βE , z−E ) has the form in Corollary
4.
5.5. Proof for Theorem 5
Proof. We first reconstruct the gradient ∇`(β̂E ) from β̂E and


β̄E
T =
.
T
X−E
(y − π(XE β̄E ))
The Taylor expansion of ∇`(β̂E ) at β̄E is
1
∇`(β̂E ) = − √ X T (y − π(XE β̂E ))
n

1 T
= − √ X y − π(XE β̄E ) − W (XE β̄E )XE (β̂E − β̄E ) + R,
n
where R = op (1). Since β̄E is the minimizer of the logistic regression with E variables,
the gradient at β̄E is zero,
T
XE
(y − π(XE β̄E )) = 0.

15

Tian et al./MAGIC

Thus we can rewrite ∇`(β̂E ) in terms of β̂E and T via the following map,

T
XE
W (XE β̄E )XE (β̂E − β̄E )
+ R.
T
T
X−E
W (XE β̄E )XE (β̂E − β̄E ) − X−E
(y − π(XE β̄E ))
(22)
Notice that β̄E is the MLE for the negative logistic likelihood, and thus satisfy the
asymptotic normality, with asymptotic mean bE , when E ⊇ supp(b). Moreover, the
following part has asymptotically mean 0,


1 T
E √ X−E
(y − π(XE β̄E ))
n




1 T
1 T
=E √ X−E (y − π(XE bE )) − E √ XE
W (XE bE )XE (β̄E − bE ) + op (1) → 0.
n
n
1
∇`(β̂E ) = √
n



Thus we have the asymptotic normality as in (12). Moreover, since ω = ∇`(β̂E ) +
λẑ + β̂E , then we have asymptotically,
d

(T, ∇`(β̂E ) + λẑ + β̂E ) → F × G.
T
The Jacobian
is det(XE
W (XE β̄E )X
 E +I) which by law of large numbers converges


T
W (XE bE )XE ) + I , a constant. Therefore, we have the density (13)
to det E(XE

if we plug in the map (22) for ∇`(β̂).
5.6. Proof for Corollary 6
Proof. For every node i, the i-th coordinate of β i is held to be zero, and (14) is in fact
T
a regression of dimension p − 1, thus γ(X, β̂ i ) = −X−i
(xi − X β̂ i ) ∈ Rp−1 , and the
reconstruction map,
ψ : (X, B̂, ẑ) 7→ (X, γ(X, B̂) + λẑ),
where

γ(X, B̂) = (γ(X, β̂ 1 ), . . . , γ(X, β̂ p )) ∈ R(p−1)×p ,

ẑ = (ẑ 1 , ẑ 2 , . . . , ẑ p ), B̂ = (β̂ 1 , . . . , β̂ p ),
 i

zE,obs
i
and ẑ =
is the subgradient of the optimization problem (14). Since ω i ’s are
i
ẑ−E
independent, and the Jacobian
Y
T
Jψ(X, z, B) =
det(XE
i XE i ),
i∈Γ

density (16) follows.

16

Tian et al./MAGIC

Algorithm 1 Metropolis Hastings sampler for randomized Lasso
Set: b = 0 for distribution fb , compute the explicit expression h.
†
Compute: P = XE\j XE\j
, R = I − P,
0 , ẑ 0 ) ← (y, β̂ , ẑ
Initialize: (y 0 , β̂E
E −E ),
−E
(n+1)
Step data: y
← P y (n) + an · Rτ,
(n)

(n)

τ ∼ N (0, I), compute the acceptance ratio r =

(n)

h(y (n+1) ,β̂E ,ẑ−E )
(n)

h(y (n) ,β̂E ,ẑ−E )
y (n+1) .

, accept

(n+1)

Step coefficient: β̂E

(n+1) (n)
g(y
,β̂E
,ẑ−E )
(n) (n)
g(y (n+1) ,β̂E ,ẑ−E )

y (n+1)

with probability r, otherwise y (n+1) ← y (n) . If r > 1, accept
(n)

← sE |β̂E

+ cn · ν|, ν ∼ G, compute the acceptance ratio r =

(n+1)

, and accept/reject accordingly.

Step subgradient: compute the upper and lower limits,
(n+1)

T
∆+ = −X−E
(y (n+1) − XE β̂E

∆

−

=

T
−X−E
(y (n+1)

−

) + λ1,

(n+1)
XE β̂E
)

− λ1,

(n+1) ind

sample λẑ−E

∼ G∆− ,∆+ .

6. Monte-Carlo sampler
Theorem 3 gives an explicit way of computing the density for the law of selective inference. We can use a Gibbs sampler to rotate through sampling (S, β̂, ẑ). For sampling S
and β̂, we can take a Metropolis-Hastings step and use the density to compute the acceptance probability. For sampling ẑ, it is even simpler as we recognize the conditional
distribution of ẑ|S, β̂ is simply a truncated G distribution. To illustrate our sampler, we
describe the sampling scheme of some of our examples in more details.
6.1. Randomized Lasso sampler
Without loss of generality, we assume the density of added noise G is symmetric and
each coordinate of ω is independent. This is in fact what we use a lot in practice.
Also denote G∆− ,∆+ as truncated distribution G with ∆− , ∆+ as the lower and upper
truncation points, and h(y, β̂E , ẑ−E ) to be the density in (10). Then to test the null
hypothesis H0j : bj = 0, we propose Algorithm 1. Note the step sizes an and cn in
Algorithm 1 is chosen through [18] to achieve the desired acceptance rate.
6.2. Neighborhood selection
Similar to the scheme in Section 6.1, we use a Gibbs sampler to sample X, B̂ and ẑ
respectively. The sampling for B̂ and ẑ are analogous to that of Section 6.1, and we only
need a proposal distribution for X. As mentioned in Section 3.3, to test the hypothesis
H0,ij : Θij = 0, we condition on {xTi0 xj 0 , (i0 , j 0 ) 6= (i, j)}. To sample the data matrix
X, we rotate through its columns, sampling one column at a time, keeping all the others
as constant. More specifically, for column i, we sample from the distribution,
xi |X−i , kxi k2 , xTi0 xj 0 , (i0 , j 0 ) ∈ E ∨ , (i0 , j 0 ) 6= (i, j).

Tian et al./MAGIC

17

Note the graph structure gives a natural partition of the nodes into different connected
components, let ne(i) be the nodes in the connected component of i, then xi ⊥ xk , ∀k 6∈
ne(i), conditioning on all the other xj ’ in ne(i). Thus the above law is equivalent to,
xi |kxi k2 , xj , xTi xj , j ∈ ne(i).

(23)

We can sample the above law (23) by sampling uniformly from a sphere with radius
kxi k, holding the projections onto the xj ’s constant. After sampling a new column of
X, we compute the accept ratio, accept/reject accordingly and move to the next column.
As for the sampling of B̂ and ẑ, we can develop an algorithm similar to Algorithm 1.
7. Discussion
MAGIC has the following limitations that we hope to remove in future work. First,
the penalty in our convex program have to be `1 penalty. Second, we assume parametric models, more specifically in the exponential family setting. Third, in the setting
for Section 3.2, we require the dimension p to be fixed, leaving the high-dimensional
problem p > n as an interesting problem.
References
[1] Claude JP Bélisle, H Edwin Romeijn, and Robert L Smith. Hit-and-run algorithms for generating multivariate distributions. Mathematics of Operations Research, 18(2):255–266, 1993.
[2] Richard Berk, Lawrence Brown, Andreas Buja, Kai Zhang, and Linda Zhao. Valid
post-selection inference. The Annals of Statistics, 41(2):802–837, April 2013.
[3] S. Bubeck, R. Eldan, and J. Lehec. Sampling from a log-concave distribution
with Projected Langevin Monte Carlo. ArXiv e-prints, July 2015.
[4] DR Cox. A note on data-splitting for the evaluation of significance levels.
Biometrika, 62(2):441–444, 1975.
[5] William Fithian, Dennis Sun, and Jonathan Taylor. Optimal inference after model
selection. arXiv:1410.2597 [math, stat], October 2014. arXiv: 1410.2597.
[6] Jason D. Lee, Dennis L. Sun, Yuekai Sun, and Jonathan E. Taylor. Exact postselection inference with the lasso. arXiv:1311.6238 [math, stat], November 2013.
[7] Erich L. Lehman. Testing Statistical Hypotheses. Springer-Verlag, New York,
2nd edition, 1997.
[8] Nicolai Meinshausen and Peter Bühlmann. High-dimensional graphs and variable
selection with the lasso. The Annals of Statistics, pages 1436–1462, 2006.
[9] Nicolai Meinshausen and Peter Bühlmann. Stability selection. Journal of the
Royal Statistical Society: Series B (Statistical Methodology), 72(4):417–473,
2010.
[10] Nicolai Meinshausen, Lukas Meier, and Peter Bühlmann. P-values for highdimensional regression. Journal of the American Statistical Association, 2012.
[11] Sahand Negahban, Pradeep Ravikumar, Martin J Wainwright, and Bin Yu. A
unified framework for high-dimensional analysis of $m$-estimators with decomposable regularizers. arXiv:1010.2731, October 2010.

Tian et al./MAGIC

18

[12] Soo-Yon Rhee, Jonathan Taylor, Gauhar Wadhera, Asa Ben-Hur, Douglas L Brutlag, and Robert W Shafer. Genotypic predictors of human immunodeficiency
virus type 1 drug resistance. Proceedings of the National Academy of Sciences,
103(46):17355–17360, 2006.
[13] Jonathan Taylor and Robert Tibshirani. Post-selection inference for l1-penalized
likelihood models. arXiv preprint arXiv:1602.07358, 2016.
[14] Xiaoying Tian, Joshua R. Loftus, and Jonathan E. Taylor. Selective inference with
unknown variance via the square-root LASSO. arXiv:1504.08031 [math, stat],
April 2015. arXiv: 1504.08031.
[15] Xiaoying Tian and Jonathan E. Taylor. Selective inference with a randomized
response. arXiv:1507.06739 [math, stat], July 2015. arXiv: 1507.06739.
[16] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of
the Royal Statistical Society: Series B, 58(1):267–288, 1996.
[17] Sara A Van de Geer. High-dimensional generalized linear models and the lasso.
The Annals of Statistics, pages 614–645, 2008.
[18] Matti Vihola. Robust adaptive metropolis algorithm with coerced acceptance rate.
Statistics and Computing, 22(5):997–1008, 2012.
[19] Martin J Wainwright. Sharp thresholds for high-dimensional and noisy sparsity
recovery using-constrained quadratic programming (lasso). Information Theory,
IEEE Transactions on, 55(5):2183–2202, 2009.
[20] Larry Wasserman and Kathryn Roeder. High dimensional variable selection. Annals of statistics, 37(5A):2178, 2009.
[21] Peng Zhao and Bin Yu. On model selection consistency of lasso. The Journal of
Machine Learning Research, 7:2541–2563, 2006.
[22] Hui Zou and Trevor Hastie. Regularization and variable selection via the elastic
net. Journal of the Royal Statistical Society: Series B (Statistical Methodology),
67(2):301–320, 2005.

