Approximate selective inference via maximum likelihood

TL;DR

This paper introduces an approximate maximum likelihood inference method for Gaussian data that enhances inferential power and computational efficiency in post-model selection scenarios.

Contribution

It proposes a novel convex optimization approach for approximate selective inference, reducing reliance on complex MCMC sampling and improving inference quality.

Findings

Improved inferential power demonstrated in simulations.

Method outperforms existing strategies on gene expression data.

Efficient convex optimization replaces costly sampling procedures.

Abstract

Several strategies have been developed recently to ensure valid inference after model selection; some of these are easy to compute, while others fare better in terms of inferential power. In this paper, we consider a selective inference framework for Gaussian data. We propose a new method for inference through approximate maximum likelihood estimation. Our goal is to: (i) achieve better inferential power with the aid of randomization, (ii) bypass expensive MCMC sampling from exact conditional distributions that are hard to evaluate in closed forms. We construct approximate inference, e.g., p-values, confidence intervals etc., by solving a fairly simple, convex optimization problem. We illustrate the potential of our method across wide-ranging values of signal-to-noise ratio in simulations. On a cancer gene expression data set we find that our method improves upon the inferential power of some commonly used strategies for selective inference.

Equations

$$\underset{o}{\mathrm{minimize}}\frac{1}{2}\Vert y-Xo{\Vert }_2^2+\lambda \Vert o{\Vert }_1+\frac{ϵ}{2}\Vert o{\Vert }_2^2-{\omega }^{⊺}o.$$

$$\left\{(y,\omega ):\hat{\mathrm{E}}(y,\omega )=\mathrm{E}\right\},$$

$${\hat{\beta }}_{j\cdot \mathrm{E}}^{\text{mle}}\pm z_{1-q/2}\cdot \sqrt{I_{j,j}^{-1}({\hat{\beta }}_{\mathrm{E}}^{\text{mle}})},$$

$$Y\sim N(\beta ,1),W\sim N(0,{\eta }^2),$$

$$Y+W>\tau ,\text{ where }\tau =\sqrt{1+{\eta }^2}\cdot z_{1-q}.$$

$$\left\{(y,\omega )\in {}^2:y+\omega >\tau \right\}.$$

$${\left(\bar{\Phi }\left(\frac{(\tau -\beta )}{\sqrt{(1+{\eta }^2)}}\right)\right)}^{-1}\varphi (Y;\beta ,1)\cdot \varphi (O;Y-\tau ,{\eta }^2)\cdot 1_{(0,\infty )}(O).$$

$${\left(\bar{\Phi }\left(\frac{(\tau -\beta )}{\sqrt{(1+{\eta }^2)}}\right)\right)}^{-1}\varphi (Y;\beta ,1)\cdot \bar{\Phi }\left(\frac{1}{\eta }(\tau -Y)\right).$$

$$\nabla \alpha ({\hat{\beta }}^{\text{mle}})=Y,$$

$$\alpha (\beta )=\frac{1}{2}{\beta }^2+\log \bar{\Phi }\left(\frac{(\tau -\beta )}{\sqrt{(1+{\eta }^2)}}\right).$$

$$|\text{det}({\nabla }^2\alpha ({\hat{\beta }}^{\text{mle}}))|\cdot \varphi (\nabla \alpha ({\hat{\beta }}^{\text{mle}});\beta ,1)\cdot \bar{\Phi }\left(\frac{1}{\eta }(\tau -\nabla \alpha ({\hat{\beta }}^{\text{mle}}))\right).$$

$$\mathbb{E}\left[({\hat{\beta }}^{\text{mle}}-\beta {)}^2|Y+W>\tau \right]\le (B{)}^{-1}\cdot \text{Var}(Y|Y+W>\tau ).$$

$$1-\frac{({\hat{\beta }}^{\text{mle}}-\tau )}{(1+{\eta }^2{)}^{3/2}}\cdot {\left(\bar{\Phi }\left(\frac{(\tau -{\hat{\beta }}^{\text{mle}})}{\sqrt{(1+{\eta }^2)}}\right)\right)}^{-1}\varphi \left(\frac{(\tau -{\hat{\beta }}^{\text{mle}})}{\sqrt{(1+{\eta }^2)}};0,1\right)$$

$$-\frac{1}{(1+{\eta }^2)}\cdot {\left(\bar{\Phi }\left(\frac{(\tau -{\hat{\beta }}^{\text{mle}})}{\sqrt{(1+{\eta }^2)}}\right)\right)}^{-2}{\varphi }^2\left(\frac{(\tau -{\hat{\beta }}^{\text{mle}})}{\sqrt{(1+{\eta }^2)}};0,1\right).$$

$$\bar{\Phi }\left(\sqrt{I({\hat{\beta }}^{\text{mle}})}({\hat{\beta }}^{\text{mle}}-\beta )\right).$$

$$\left\{(y,\omega )\in {}^n\times {}^p:\hat{\mathrm{S}}(y,\omega )=\mathrm{S}\right\}\subseteq \left\{(y,\omega )\in {}^n\times {}^p:\hat{\mathrm{E}}(y,\omega )=\mathrm{E}\right\}.$$

$$\mathbb{M}_{\mathrm{S}}=\Big{\{}Y\sim N_{n}(X_{\mathcal{E}}\beta_{\mathcal{E}},\sigma^{2}I),\;\;\beta_{\mathcal{E}}\in\mathbb{R}^{|\mathcal{E}|}\Big{\}}\ \text{ for a fixed }\ \sigma\in{}^{+}.$$

$${\beta }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}}={\mathcal{F}}_{\mathrm{S}}\mathbb{E}[Y]\in {}^d$$

$${\hat{\beta }}_{\mathrm{S}}\sim N({\beta }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}},{\Sigma }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}}).$$

$$\left\{(y,\omega )\in {}^n\times {}^p:\hat{\mathrm{S}}(y,\omega )=\mathrm{S}\right\}.$$

$$\Big{\{}Y\sim N_{n}(X_{\mathrm{E}}\beta_{\mathrm{E}},\sigma^{2}I),\ \beta_{\mathrm{E}}\in{}^{|\mathrm{E}|}\Big{\}},$$

$${\hat{\beta }}_{\mathrm{S}}^{\perp }=-X^{⊺}(y-X_{\mathrm{E}}{\hat{\beta }}_{\mathrm{S}}),$$

$$\omega ={\left(\begin{array}{cc}{\omega }_{\mathrm{E}}^{⊺}& {\omega }_{E^c}^{⊺}\end{array}\right)}^{⊺}=-X^{⊺}{X}_{\mathrm{E}}{\hat{\beta }}_{\mathrm{S}}+\left[\begin{array}{c}{X}_{\mathrm{E}}^{⊺}{X}_{\mathrm{E}}+ϵI\\ X_{{\mathrm{E}}^c}^{⊺}{X}_{\mathrm{E}}\end{array}\right]o_1+\left(\begin{array}{c}\lambda z_{\mathrm{E}}\\ o_2\end{array}\right)+{\hat{\beta }}_{\mathrm{S}}^{\perp },$$

$$P_{\mathrm{S}}=-X^{⊺}{X}_{\mathrm{E}},Q_{\mathrm{S}}=\left[\begin{array}{c}{X}_{\mathrm{E}}^{⊺}{X}_{\mathrm{E}}+ϵI\\ X_{{\mathrm{E}}^c}^{⊺}{X}_{\mathrm{E}}\end{array}\right],r_{\mathrm{S}}=\left(\begin{array}{c}\lambda z_{\mathrm{E}}\\ o_2\end{array}\right)+{\hat{\beta }}_{\mathrm{S}}^{\perp },$$

$$\omega ={\left(\begin{array}{cc}{\omega }_{\mathrm{E}}^{⊺}& {\omega }_{E^c}^{⊺}\end{array}\right)}^{⊺}=P_{\mathrm{S}}{\hat{\beta }}_{\mathrm{S}}+Q_{\mathrm{S}}{o}_1+r_{\mathrm{S}}.$$

$${\bar{\Sigma }}^{-1}=Q_{\mathrm{S}}^{⊺}{\Sigma }_{\mathbb{W}}^{-1}{Q}_{\mathrm{S}},A=-\bar{\Sigma }{Q}_{\mathrm{S}}^{⊺}{\Sigma }_{\mathbb{W}}^{-1}{P}_{\mathrm{S}},b=-\bar{\Sigma }{Q}_{\mathrm{S}}^{⊺}{\Sigma }_{\mathbb{W}}^{-1}{r}_{\mathrm{S}},$$

$$f({\tilde{\beta }}_{\mathrm{S}})=\int \varphi (o_1;A{\tilde{\beta }}_{\mathrm{S}}+b,\bar{\Sigma })\cdot 1_{{\mathrm{R}}_0}(o_1)d{o}_1,$$

$${\left(\int \varphi ({\tilde{\beta }}_{\mathrm{S}};J{\beta }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}}+k,\Sigma )\cdot f({\tilde{\beta }}_{\mathrm{S}})d{\tilde{\beta }}_{\mathrm{S}}\right)}^{-1}\varphi ({\hat{\beta }}_{\mathrm{S}};J{\beta }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}}+k,\Sigma )\cdot f({\hat{\beta }}_{\mathrm{S}}),$$

$$\Sigma =({\Sigma }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}}^{-1}+P_{\mathrm{S}}^{⊺}{\Sigma }_{\mathbb{W}}^{-1}{P}_{\mathrm{S}}-A^{⊺}{\bar{\Sigma }}^{-1}A{)}^{-1},J=\Sigma {\Sigma }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}}^{-1},k=\Sigma (A^{⊺}{\bar{\Sigma }}^{-1}b-P_{\mathrm{S}}^{⊺}{\Sigma }_{\mathbb{W}}^{-1}{r}_{\mathrm{S}}).$$

$$\varphi ({\hat{\beta }}_{\mathrm{S}};J{\beta }_{{\mathbb{M}}_{\mathrm{S}},\mathrm{S}}+k,\Sigma )\cdot \varphi (O_1;A{\hat{\beta }}_{\mathrm{S}}+b,\bar{\Sigma }).$$

$$\displaystyle-\inf_{(\widetilde{\beta}_{\mathrm{S}},o_{1})\in\mathrm{R}}\Big{\{}\dfrac{1}{2}(\widetilde{\beta}_{\mathrm{S}}-J\beta_{\mathbb{M}_{\mathrm{\mathrm{S}}},\mathrm{\mathrm{S}}}-k)^{\intercal}$$

Full text

Keywords: Data adaptivity, Conditional inference, Maximum likelihood, Multiple queries, Post-selection inference, Randomization, Selective MLE.

1 Introduction

Querying the data has become a fairly common practice for anyone who wishes to learn a model from a range of different candidates. Naively using the same data twice, first to learn a model and then infer for the selected parameters, tends to inflate their estimated effects. As an example of a query, consider a variable selection algorithm with a shrinkage penalty (Tibshirani, 1996; Fan and Li, 2001; Yuan and Lin, 2006); the algorithm learns a set of variables (or features) into a model. Ignoring the dependence of the model (and its parameters) on the outcome of the query while calculating p values, confidence intervals, credible intervals etc. undermines inference after selection; see Benjamini and Yekutieli (2005); Leeb and Pötscher (2005, 2006); Berk et al. (2013) for a demonstration of the concerns here. The result is usually an increased chance of finding a statistically significant result when the selected variable in fact has no effect.

Various strategies for selective inference offer different solutions by characterizing the dependence between the learned models and data. Some of the strategies are easier to implement, while others fare better in inferential power. In this paper, we introduce a new method for selective inference through approximate maximum likelihood estimation. Building on the recent work by Tian and Taylor (2018), our method allows us to: (i) harness Gaussian randomization variables towards better inferential power after selection, and simultaneously (ii) bypass expensive MCMC sampling from intractable conditional distributions. Below, we provide a brief, informal overview of our method in a standard setup of linear regression.

An informal overview of our method. Consider a regression problem in which we observe a response vector $y\in{}^{n}$ and a matrix of $p$ predictors $X\in{}^{n\times p}$. Let $\omega\in{}^{p}$ be drawn from a centered Gaussian distribution with known covariance $\Sigma_{\mathbb{W}}$. For fixed values $\lambda\in{}^{+}$, $\epsilon\in{}^{+}$, we solve the following query:

[TABLE]

We call the optimization in (1) a “randomized LASSO” query; so named because of the randomization variable $\omega$ added to the objective of the canonical LASSO. Suppose the query selects a nonempty set of variables $\mathrm{E}\subseteq\{1,2,\cdots,p\}$. Following selection, we describe our response variable through the model: $y=X_{\mathrm{E}}\beta_{\mathrm{E}}+\mathrm{e},\;\mathrm{e}\sim N(0,\sigma^{2}I_{n}),$ where $I_{n}$ is the identity matrix with $n$ rows and columns.

A natural ask in the learned model is inference for the partial regression coefficients after adjusting for their dependence on data through $\mathrm{E}$. One concrete course of action is to condition on selection, specifically, base inference on the likelihood of the observed data when conditioned on the event

[TABLE]

where $\widehat{\mathrm{E}}$ represents the (data-dependent) selected set of variables. Maximizing the conditional likelihood function gives us $\widehat{\beta}^{\text{\;mle}}_{\mathrm{E}}$, the maximum likelihood estimate (MLE) for $\beta_{\mathrm{E}}$. Taking the hessian of the negative log-likelihood at the MLE yields us $I(\widehat{\beta}_{\mathrm{E}}^{\text{\;mle}})$, the observed Fisher information matrix. The (approximate) confidence intervals resulting from our method take the form

[TABLE]

where $\widehat{\beta}^{\text{\;mle}}_{j\cdot\mathrm{E}}$ is the $j^{\text{th}}$ component of $\widehat{\beta}^{\text{\;mle}}_{\mathrm{E}}$, $I^{-1}_{j,j}(\widehat{\beta}_{\mathrm{E}}^{\text{\;mle}})$ is the $(j,j)^{\text{th}}$ entry of $I^{-1}(\widehat{\beta}_{\mathrm{E}}^{\text{\;mle}})$ and $z_{1-q}$ is the $(1-q)$-th quantile of a standard normal distribution for $q\in(0,1)$.

The intervals proposed in (2) are seemingly straightforward if only we could directly calculate the two estimates in the expression, the MLE and the observed Fisher information matrix. But as it turns out, the conditional likelihood and subsequently the two estimates based on it do not admit expressions in closed forms. In the remaining development, we solve this challenge head-on in two steps. First, we construct an approximate, statistically consistent proxy for the likelihood function after conditioning on the selection event. Then, we provide a tractable system of estimating equations to obtain the MLE and the observed Fisher information matrix from our proxy likelihood. At the core of the proposed estimating equations is a fairly simple, convex optimization problem in relatively few dimensions.

Comparison with common baselines. Continuing with the regression setup, we compare our proposal with some common baselines in a simulated experiment and relate our method with existing work. In Table 1, four methods for selective inference including our method are evaluated on three criteria after selecting variables using a fixed value of tuning parameter: average coverage of interval estimates with nominal false coverage rate (FCR) level of $0.10$; average length of the interval estimates; the power of detecting true associations after applying the selective inference strategy. The data in this experiment obeys a linear model with a $300$-by-$100$ design matrix $X$ and Gaussian errors, such that the rows of $X$ are i.i.d. copies of a correlated multivariate normal vector; the model coefficient vector $\beta$ has $6$ nonzero components that are linearly varying in magnitude, and the setting corresponds to a relatively weak signal-to-noise ratio value. The simulation setting is described more precisely later in the paper and relative comparisons between all the methods are made for a wide range of values for signal-to-noise ratio.

The first baseline, “Lee et al.”, proposed by Lee et al. (2016) reduces inference to a truncated normal variable through the Polyhedral Lemma. The second baseline, “Split” uses a randomly chosen one-third of the data samples for inference after applying the LASSO to the remaining two-thirds of the data, (e.g., Cox, 1975; Hurvich and Tsai, 1990). Based on Liu et al. (2018), the third baseline “Liu et al.” conditions on (strictly) less information than “Lee et al.” by choosing to infer for the parameters associated with the selected variables in the full model: $y\sim N(X\beta,\sigma^{2}I_{n})$. All the other strategies in this example use the learned model: $y\sim N(X_{\mathrm{E}}\beta_{\mathrm{E}},\sigma^{2}I_{n})$ based on the selected set of variables.

[TABLE]

We note that FCR is (roughly) attained at the nominal level by all the strategies, except “Liu et al.” falls slightly short of the mark in this setting. Our method, namely “MLE”, delivers the shortest intervals with the highest power. In the last column of the table, we also indicate the percentage of intervals that resulted with infinite length. Among all the methods, “Lee et al.” produces some infinitely long intervals; this observation is consistent with the established fact in Kivaranovic and Leeb (2018) that the intervals based on the Polyhedral Lemma do not have a finite expected value in the Gaussian regression setting. Both “Split” and “Liu et al.” overcome the drawbacks of “Lee et al.” by setting aside more information for inference. The former strategy does so by reserving a randomly chosen subsample for inference, while the latter achieves an increase in power through a larger truncation set.

By analogy with Tian and Taylor (2018), our method uses added randomization to remedy the excessively long intervals produced by “Lee et al.”; because we do not condition on the randomization variable itself, inference does not trivially reduce to the Polyhedral Lemma. Our choice of adding a Gaussian randomization variable to the query draws motivation from data carving, a two stage situation where parameters learned on an initial data set are estimated using new samples augmented with the initial ones (Fithian et al., 2014; Panigrahi, 2019). In the analysis here, the variance of $\omega$ is chosen so that the randomized LASSO (roughly) resembles “Split” in the amount of information used towards learning the model. Our method improves upon “Split” by conditioning upon an event that implies the selection of the variables in the set $\mathrm{E}$ and therefore still consuming some information from the data used in selection. A direct relation between the power attained with our randomized method and “Liu et al.”, however, is lacking. Some gain in power reported for the above setting might be attributed to the use of the learned model by our method as opposed to the full model under which “Liu et al.” offers inference.

Other related work. Our maximum likelihood method differs from previous proposals in the tools used and the scope of inference. Existing strategies for selective inference usually require sampling from conditional distributions, either due to the intractability of their exact counterparts or due to the lack of easily available truncation regions. For example, the pivot described in Tian and Taylor (2018) lacks exact expressions and is not readily amenable for computational analyses. The truncation region in Liu et al. (2018) takes a tractable form for the LASSO; however, this form does not directly generalize to other queries. Other inferential approaches that account for the effects of selection include sampling from a selection-adjusted posterior in Panigrahi et al. (2021) and resampling-based approaches such as bootstrap in McKeague and Qian (2015); Guo and He (2020). The computing costs of these approaches are especially acute if two or more queries are applied for learning models, and the construct of valid inference must appropriately take into consideration the effect of each such query. Bypassing the requirement to sample from intractable conditional distributions after selection, our maximum likelihood method relies on the solution to a simple, convex optimization problem. Furthermore, this convex problem assumes a separable form under multiple queries which is particularly amenable for parallel computing. Much of prior work in the area of selective inference has relied on a testing-based approach for real-valued projections of parameters; see for example Yang et al. (2016); Suzumura et al. (2017); Rügamer and Greven (2018). In contrast, the scope of the present likelihood-centric approach extends to joint inference for parameter vectors within learned models.

The rest of the paper is organized as follows. In Section 2, we introduce our approximate proposal in a univariate file drawer problem. We describe in Section 3 our method of selective inference by deriving a system of estimating equations for the MLE and the observed Fisher information matrix after solving a convex query. In Section 4, we generalize our prescription to adjust for multiple, convex queries. We conduct simulations in Section 5 to study the gains with our method over existing baselines. We apply our method to gene expression data from The Cancer Genome Atlas in Section 6, corroborating some of the numerical findings in the simulated experiments. We conclude with a discussion in Section 7. In the Appendix, we include proofs for our main results and generalize our method of selective inference for multiple convex queries.

2 MLE Inference: A First Example

Before proceeding further, note, we use (i) $\bar{\Phi}(x)$ for the upper tail probability of the standard Gaussian law at $x\in\real$, and (ii) $\phi(u;\nu,\Theta)$ for the (multivariate) Gaussian density function with mean vector $\mu$ and covariance $\Theta$ at the value $u$; the special symbol $\phi(x;0,1)$ denotes the standard Gaussian density for $x\in\real$.

2.1 Univariate soft-truncated likelihood

We consider two independent random variables:

[TABLE]

where $W$ denotes a Gaussian randomization variable. We pursue inference for $\beta$ only if:

[TABLE]

We begin by describing a conditional likelihood by conditioning the Gaussian law of $Y$ upon the selection event:

[TABLE]

Define $O=Y+W-\tau$, which we call an optimization variable in our framework. Because, $O\lvert Y=y\sim N(y-\tau,\eta^{2})$ before conditioning on the event in (3) and the selection event, $y+\omega>\tau$, is equivalent to: $o>0$, the conditional likelihood for $Y$ and $O$ is given by:

[TABLE]

Marginalizing over the optimization variable $O$ yields us a likelihood function of $\beta$, that is equal to:

[TABLE]

Compared to the conditional Gaussian law in the absence of randomization, $Y\;\lvert\;Y>\tau,$ (c.f. Example 2, Fithian et al., 2014), a soft-truncating function replaces the indicator $1_{(\tau,\infty)}(Y)$. Hereafter, we refer to the resulting function as a “soft-truncated likelihood”.

2.2 Selective MLE

We are now ready to discuss the maximizer of the soft-truncated likelihood in (4) and inspect some properties of this estimate that serve to motivate our approximate pivot in the paper. Maximizing the log-likelihood gives us the “selective MLE”, $\widehat{\beta}^{\text{\;mle}}$, based on the estimating equation:

[TABLE]

where

[TABLE]

Turning to the distribution of $\widehat{\beta}^{\text{\;mle}}$, we obtain the density for the selective MLE from (4) by applying the simple variable transformation: $\widehat{\beta}^{\text{\;mle}}=\nabla\alpha^{-1}(Y)$, and note that this density is proportional to:

[TABLE]

Proposition 2.1 obtains an upper bound for the mean squared error of the selective MLE. An immediate consequence of this bound is a global (asymptotic) consistency guarantee for the selective MLE; i.e., the guarantee continues to hold for the event in (3) even when it has a vanishing probability as the sample size grows to infinity. Streamlining the main exposition to focus on finite sample results, we defer the proof for asymptotic consistency to the Appendix C.

Proposition 2.1**.**

Fix $B=(1+\eta^{2})^{-2}{\eta^{4}}$. Then, we have:

[TABLE]

Examining for now the asymptotic behavior of selective MLE and the least squares estimate, and the role of the randomization variable $W$, we undertake a simulation by letting $Y:=\sqrt{n}\bar{Y}_{n}$ with mean $\beta:=\sqrt{n}\beta_{n}$. Figures 1, 2 and 3 summarize the three primary take-aways from the simulation. We let $\tau=0$ in (3), and fix $\beta_{n}=\beta_{0}=-0.10$ which is highlighted in the figures via a dotted black line. Notice, our choice of $\beta_{n}$ results in rarer events of selection with vanishing probabilities as $n\to\infty$. For the first two figures, the randomization variance $\eta^{2}$ is equal to $1$. Based on the density in (6), Figure 1 first studies the behavior of the selective MLE. Matching our theoretical expectations, the plot demonstrates that the selective MLE is a consistent estimate for the parameter; more specifically, we observe a concentration of the estimate around $\beta_{n}$ with increasing $n$. Next, Figure 2 replaces the selective MLE in the first plot with the least squares estimate. Unlike the selective MLE, the least squares estimate fails to concentrate around the parameter of interest for the same sample sizes. In Figure 3, we reproduce Figure 1, except now we study the behavior of the selective MLE under a very low value of randomization variance, $\eta^{2}=0.04$. An empirical affirmation of the merits of randomization, this plot shows that selective MLE fails to concentrate around $\beta_{n}$ in the (almost) absence of randomization.

2.3 Approximate pivot

Prompted by a concentration of the selective MLE around the parameter of interest, we introduce an approximate pivot in the current section. We propose to approximate the distribution of the MLE by a Gaussian distribution with: (i) mean $\beta$, and (ii) variance equal to inverse of the observed Fisher information, $I(\widehat{\beta}^{\text{\;mle}})$. Taking a second derivative of the log-likelihood in (4) at the MLE gives us the value of $I(\widehat{\beta}^{\text{\;mle}})$, which is equal to:

[TABLE]

The Gaussian approximation described above gives rise to the approximate pivot:

[TABLE]

We remark that the distribution of the selective MLE, characterized exactly by the density in (6), can yield us exact maximum likelihood inference. In contrast, the pivot in (8) is only approximate in nature, but, appealingly simple in form. Inference based on the approximate pivot requires us to compute two estimates from the soft-truncated likelihood, namely, the selective MLE and the observed Fisher information.

Before turning to the general development, we explore if the proposed Gaussian approximation mimics the exact distribution of the selective MLE. In Figure 4, we represent the density of the selective MLE in (6), our benchmark, by the gray curve. The panel with $\beta=-3$ results in a rare selection event, while the panel with $\beta=1.5$ results in a highly probable selection event with little selection bias. Noteworthy, the effectiveness of our pivot is highlighted via a strong agreement of the proposed (approximate) Gaussian density with the exact (benchmark) density of the selective MLE.

3 Maximum Likelihood Inference Post Convex Queries

We develop our method of maximum likelihood inference below, focusing on the randomized LASSO in a general class of convex queries as our leading example. In the Appendix, we show that the form of our estimating equations in the primary example generalizes directly to other convex learning queries whose solutions can be similarly characterized through affine Karush-Kuhn-Tucker (K.K.T.) conditions of optimality.

3.1 Framework under linear regression

Consider solving the randomized LASSO in (1), where $y$ and $\omega$ denote the observed instances of our response variable $Y\in{}^{n}$ and randomization variable $W\sim N(0_{p},\Sigma_{\mathbb{W}})\in{}^{p}$ respectively. Let $\widehat{\mathrm{E}}(y,\omega)\subseteq\{1,2,\cdots,p\}$ denote the active set of variables selected by the randomized LASSO. At the solution of the randomized query, we record the value of the subgradient vector for the $\ell_{1}$ penalty which we represent by $\widehat{\mathrm{S}}(y,\omega)$. Notice, the collection of instances which lead us to observe $\widehat{\mathrm{S}}(y,\omega)=\mathrm{S}$ result in the active set $\widehat{\mathrm{E}}(y,\omega)=\mathrm{E}$, i.e.,

[TABLE]

We turn to a framework for selective inference, allowing our model and parameters to depend on data through the recorded output of our randomized query, $\mathrm{S}$. Consider a pre-specified mapping $\mathcal{H}:\mathrm{S}\to\mathcal{E}\subseteq\{1,2,\cdots,p\}$. Our model, after observing $\widehat{\mathrm{S}}(y,\omega)=\mathrm{S}$ and subsequently, noting $\mathcal{E}=\mathcal{H}(\mathrm{S})$, is given by:

[TABLE]

The mapping $\mathcal{H}$ grants us the flexibility to use (arbitrary) linear models informed by $\mathrm{S}$, including the selected model in the special case when $\mathcal{E}=\mathrm{E}$, the active set of variables. Suppose, we have a matrix $\mathcal{F}_{\mathrm{S}}\in\mathbb{R}^{d\times n}$, that is allowed to depend on data through $\mathrm{S}$. Then, let

[TABLE]

be our parameter vector of inferential interest.

In the next step, we form a (multivariate) likelihood function of $\beta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$ in the learned model $\mathbb{M}_{\mathrm{S}}$. To do so, for a fixed value $\mathrm{S}$, we consider the following statistic:

[TABLE]

We derive our soft-truncated likelihood by conditioning the Gaussian law of $\widehat{\beta}_{\mathrm{S}}$ upon the selection event:

[TABLE]

The event in (11) depends not just on $\widehat{\beta}_{\mathrm{S}}$, but further involves a statistic independent of $\widehat{\beta}_{\mathrm{S}}$, which we represent by $\widehat{\beta}^{\perp}_{\mathrm{S}}$. In addition to the above event, we condition on $\widehat{\beta}^{\perp}_{\mathrm{S}}$ to eliminate nuisance parameters from the likelihood.

For ease of exposition, hereafter, we specialize the above framework to a projected parameter in the selected model

[TABLE]

which we obtain by applying the specific mapping $\mathcal{H}(\mathrm{S})=\mathrm{E}$ and fixing $\mathcal{F}_{\mathrm{S}}=(X_{\mathrm{E}}^{\intercal}X_{\mathrm{E}})^{-1}X_{\mathrm{E}}^{\intercal}\in\mathbb{R}^{|\mathrm{E}|\times n}$. As noted in Berk et al. (2013); Lee et al. (2016), our parameter for inference in this model is the projection of the mean for $Y$ onto the subspace spanned by the columns in $X_{\mathrm{E}}$. Immediately, we recognize: $\widehat{\beta}_{\mathrm{S}}=(X_{\mathrm{E}}^{\intercal}X_{\mathrm{E}})^{-1}X_{\mathrm{E}}^{\intercal}y$, the least squares statistic refitted to $y$ and $X_{\mathrm{E}}$. Besides the least squares statistic, the likelihood involves

[TABLE]

which is independent of $\widehat{\beta}_{\mathrm{S}}$ under the Gaussian model in (12); we will detail this out in the following section.

3.2 Multivariate soft-truncated likelihood

Our main result in the section, Theorem 1, obtains a soft-truncated likelihood function. As seen in the file drawer example, we begin with a compact representation for our selection event in terms of optimization variables based on the randomized LASSO solution. Introducing some more notations, we let $O_{1}\in{}^{|\mathrm{E}|}$ and $O_{2}\in{}^{p-|\mathrm{E}|}$ represent the active (non-zero) components of randomized LASSO solution and the subgradient (sub-)vector for the $\ell_{1}$ penalty at the inactive indices in $\mathrm{E}^{c}$ respectively, and let $o_{1}$ and $o_{2}$ be the observed instances for these variables. Let $z_{\mathrm{E}}=\text{sign}(o_{1})\in{}^{|\mathrm{E}|}$ be the sign vector for the active components of the estimated LASSO solution. The K.K.T. conditions for the randomized LASSO are given by:

[TABLE]

where: $-\text{diag}(z_{\mathrm{E}})\ o_{1}<0$, and $\|o_{2}\|_{\infty}<\lambda$. Because, $\mathrm{S}=\begin{pmatrix}\lambda z^{\intercal}_{\mathrm{E}}&o^{\intercal}_{2}\end{pmatrix}^{\intercal},$ our selection event in (11) is equivalent to the $|\mathrm{E}|$ linear constraints: $Uo_{1}<v$, for the fixed matrices $U=-\text{diag}(z_{\mathrm{E}})$, $v=0_{|\mathrm{E}|}$. We note a resemblance with the file drawer example, wherein the selection event is equivalent to the single linear constraint: $o>0$.

Fixing the matrices:

[TABLE]

we rewrite the stationary mapping in the K.K.T. condition as

[TABLE]

Based on

[TABLE]

define:

[TABLE]

where $\mathrm{R}_{0}=\{o_{1}\in{}^{|\mathrm{E}|}:Uo_{1}<v\}$.

Theorem 1**.**

After conditioning the law of $\widehat{\beta}_{\mathrm{S}}$ upon $\widehat{\mathrm{S}}(Y,W)=\mathrm{S}$ and $\widehat{\beta}^{\perp}_{\mathrm{S}}(Y)=\widehat{\beta}^{\perp}_{\mathrm{S}}$, the soft-truncated likelihood is

[TABLE]

where $\Sigma$, $J$, $k$ are equal to:

[TABLE]

3.3 Approximate inference

The likelihood in Theorem 1, though exact, does not directly result in tractable estimating equations for the maximum likelihood estimate and the observed Fisher information matrix. This is because the normalizer for the soft-truncated likelihood lacks a closed-form expression. To circumvent the problem, we propose an approximate proxy for our soft-truncated likelihood based on an upper bound for the normalizer in Proposition 3.1. Later in the Appendix, using a large deviations principle, we prove that the approximate proxy converges to the exact likelihood with increasing sample size. Furthermore, we show the maximizer of the approximate likelihood, $\widehat{\beta}^{\;\text{mle}}_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$, is guaranteed to concentrate around the parameter in (10). Endowed with a property we expect with the exact MLE, we call the maximizer of the approximate likelihood an “approximate selective MLE”.

Proposition 3.1**.**

Let $\mathrm{R}$ be a convex and compact subset of ${}^{|\mathrm{E}|}\times{}^{|\mathrm{E}|}$. Suppose, $\widehat{\beta}_{\mathrm{S}}$ and $O_{1}$ are drawn from a Gaussian distribution with the following likelihood:

[TABLE]

Then, $\log\mathbb{P}\left[\begin{pmatrix}\widehat{\beta}_{\mathrm{S}}^{\intercal},O_{1}^{\intercal}\end{pmatrix}^{\intercal}\in\mathrm{R}\right]$ is bounded from above by

[TABLE]

Recall, $\mathrm{R}_{0}=\{o_{1}\in{}^{|\mathrm{E}|}:Uo_{1}<v\}$. We apply the bound in Proposition 3.1 to obtain the following proxy for the exact log-likelihood:

[TABLE]

after ignoring constants free of the parameter vector.

Remark 1**.**

Observe, ${}^{|\mathrm{E}|}\times\mathrm{R}_{0}$, the subset of $\;{}^{|\mathrm{E}|}\times{}^{|\mathrm{E}|}$ associated with the our selection event is clearly not compact. While compactness is a requirement to prove that the approximation in Proposition 3.1 is an upper bound for the normalizer of the likelihood, in practice, we may consider a sufficiently large compact, convex subset such that the probability of the associated event converges to the actual probability with increasing sample size.

As noted in Panigrahi et al. (2018), we can further modify the approximation in Proposition 3.1 to solve an unconstrained optimization via a barrier penalty that reflects the same constraints, but allocates a higher preference to the optimizing variables within the selection region. Letting $\mathcal{B}_{U;v}(o_{1})$ denote a barrier penalty for the constraints $Uo_{1}<v$, the final expression for our approximate log-likelihood agrees up to an additive constant with:

[TABLE]

Based on the approximate likelihood in (15), the results in Theorem 2 and Theorem 4 give us compact estimating equations for the two ingredients of approximate maximum likelihood inference.

Theorem 2**.**

Consider the optimization problem

[TABLE]

Then, maximizing the approximate log-likelihood in (15) yields us the following estimating equation for the approximate selective MLE:

[TABLE]

In line with Proposition 2.1 for the file drawer example, Theorem 3 provides a bound for the mean squared error of the approximate selective MLE. The bound in this result allows us to formalize a global consistency guarantee for our estimate in Appendix C.

Theorem 3**.**

Let the smallest eigen values for $(\Sigma^{-1}_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}+P_{\mathrm{S}}^{\intercal}\Sigma_{\mathbb{W}}^{-1}P_{\mathrm{S}})^{-1}$ and $\Sigma^{-1}_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$ be $\lambda_{0}$ and $\lambda_{1}$ respectively. Fix $B=(\lambda_{0}\cdot\lambda_{1})^{2}$. Based on the real-valued mapping:

[TABLE]

we have

[TABLE]

We provide a proxy for the observed Fisher information matrix based on the estimate in Theorem 2.

Theorem 4**.**

Let $o_{1}^{*}(\widehat{\beta}_{\mathrm{S}})$ be the solution to the optimization problem in (16). The observed Fisher information $I(\widehat{\beta}^{\;\text{mle}}_{\mathbb{M}_{\mathrm{S}},\mathrm{S}})$ for the approximate log-likelihood in (15) is:

[TABLE]

We summarize in Algorithm 1 our steps for maximum likelihood inference. Our primary computational step is the simple, |E|-dimensional, convex optimization problem (O) resulting in (S-MLE) and (FI). Emphasized earlier, the form of the estimating equations for the randomized LASSO generalizes to convex queries with affine K.K.T. conditions of optimality as in (18). In Appendix B, we illustrate how our method applies to: (i) variable screening based on marginal correlations (Lee and Taylor, 2014); (ii) variable selection via SLOPE (Bogdan et al., 2015).

4 Maximum Likelihood Inference post Multiple Queries

We now turn to the case of multiple convex queries, outlining our method after $L$ randomized LASSO queries. The randomized LASSO query yet again serves only as our primary example for ease of presentation. More generally, any of the LASSO queries can be replaced with a convex query that admits affine K.K.T. conditions of optimality at the solution.

Let $\widehat{\mathrm{S}}^{(l)}$ be the subgradient vector for the $\ell_{1}$ penalty at the solution of the $l^{\text{th}}$ randomized LASSO query:

[TABLE]

for $l\in\{1,2,\cdots,L\}$, and let

[TABLE]

Consider the event in Equation 11. Theorem 5 and 6 give the estimating equations for the approximate selective MLE and the observed Fisher information matrix after accounting for the effect of each query on inference. Notably, the system of estimating equations relies on $L$ separable optimization problems that can be solved in parallel.

Invoking the framework for selective inference in the preceding section, suppose, we apply the mapping $\mathcal{H}(\mathrm{S})=\mathrm{E}=\cup_{l=1}^{L}\mathrm{E}^{(l)}$ and fix $\mathcal{F}_{\mathrm{S}}=(X_{\mathrm{E}}^{\intercal}X_{\mathrm{E}})^{-1}X_{\mathrm{E}}^{\intercal}$. We let $o_{1}^{(l)}$ and $o_{2}^{(l)}$ be the realized instances for the active components of the randomized LASSO solution and the subgradient (sub-)vector of the $\ell_{1}$ penalty at the inactive indices in $\mathrm{E}^{c}$, respectively. Let $z_{\mathrm{E}}^{(l)}=\text{sign}(o_{1}^{(l)})$. The K.K.T. conditions of optimality for each randomized LASSO query are given by:

[TABLE]

together with the constraints: $-\text{diag}(z_{\mathrm{E}})\ o^{(l)}_{1}<0$, and $\|o^{(l)}_{2}\|_{\infty}<\lambda$, where $\widehat{\beta}_{\mathrm{S}}=(X_{\mathrm{E}}^{\intercal}X_{\mathrm{E}})^{-1}X_{\mathrm{E}}^{\intercal}y$, $\widehat{\beta}^{\perp}_{\mathrm{S}}=-X^{\intercal}(y-X_{\mathrm{E}}\widehat{\beta}_{\mathrm{S}}),$ and

[TABLE]

Our selection event in Equation 11 is equivalent to $L\cdot|\mathrm{E}|$ linear constraints: $U^{(l)}o^{(l)}_{1}<v^{(l)}$ for $l\in\{1,2,\cdots,L\}$, where $U^{(l)}=-\text{diag}(z_{\mathrm{E}^{(l)}})$ and $v^{(l)}=0_{|\mathrm{E}|}$.

Using the matrices

[TABLE]

[TABLE]

and ignoring an additive constant, an approximate proxy for the exact log-likelihood agrees with

[TABLE]

In line with Section 3 of the paper, the approximate proxy is motivated by an upper bound for the normalizer of the soft-truncated likelihood, given in the following Proposition.

Proposition 4.1**.**

Suppose, $\widehat{\beta}_{\mathrm{S}}$ and $O^{(l)}_{1}$ for $l=1,2,\cdots,L$ are drawn from a Gaussian distribution with the following likelihood:

[TABLE]

For a convex and compact set $\mathrm{R}$, $\log\mathbb{P}\left[\begin{pmatrix}\widehat{\beta}_{\mathrm{S}}^{\intercal}&(O^{(1)}_{1})^{\intercal}&\cdots&(O^{(L)}_{1})^{\intercal}\end{pmatrix}^{\intercal}\in\mathrm{R}\right]$ is bounded from above by

[TABLE]

The proof of the proposition closely follows the steps in the proof for Proposition 3.1; we thus omit the proof here.

Theorem 5**.**

Consider solving

[TABLE]

for $l\in\{1,2,\cdots,L\}$. Then, the maximizer of the approximate likelihood in (19) is

[TABLE]

Theorem 6**.**

Consider solving the $L$ optimization problems in Theorem 5. Then, the observed Fisher information matrix for the approximate likelihood (19) is

[TABLE]

We provide our steps for maximum likelihood inference in Algorithm 2. Our approximate selective MLE is a linear combination of the solutions from $L$ separable optimization problems, which we state in Step ($\mathbf{O^{(l)}}$) of the algorithm. In Section 5 of the paper, we apply the algorithm to the analysis of a cancer gene expression data for eliminating selection bias from multiple queries.

5 Simulation Experiments

We explore the potential of our method for a wide range of signal-to-noise ratio (SNR) values in a linear regression setting. Following closely the setup in Hastie et al. (2017), in every round of experiment, we simulate each $\mathbb{R}^{p}$-valued row of the design matrix $X$ from $N(0,\Sigma(\rho))$ where the $(i,j)$-th entry of $\Sigma(\rho)=\rho^{|i-j|}$. We then draw the response as: $Y\lvert X\sim N(X\beta,\sigma^{2}I)$. Fixing $n=300$, $p=100$, $\rho=0.35$, we consider a linearly-varying coefficient vector with $s=6$ non-zero equally spaced components that have magnitudes: $-10$,$-6$,$-2$,$2$,$6$,$10$. We vary the noise level $\sigma^{2}$ to match the SNR value: $\text{SNR}=\sigma^{-2}\cdot(\beta^{\intercal}\Sigma\beta)$, and vary SNR in the set $\{0.15,0.21,0.26,0.31,0.42,0.71,1.22,2.07,3.52\}$.

Dividing our experiments into two regimes, namely randomized and non-randomized, we run a canonical LASSO query (without randomization)

[TABLE]

and a randomized LASSO query in (1) with $\omega\sim N(0,\eta^{2}I_{p})$ and $\epsilon=n^{-1/2}$. The randomization variance $\eta^{2}$ is chosen so that $(\hat{\sigma})^{-2}\eta^{2}=0.50$, using the estimated noise level in the data: $\hat{\sigma}^{2}=(n-p)^{-1}\|(I-X(X^{T}X)^{-1}X^{T})y\|^{2}$. Based on an asymptotic equivalence between data splitting and a Gaussian randomization scheme (see Proposition 4.1 Panigrahi et al., 2021), our choice of randomization variance roughly matches the amount of information used up in selection when two-thirds of the samples are allocated for the LASSO. For each query, we carry out three common schemes to choose $\lambda$. Our first choice is a theoretical value proposed in Negahban et al. (2009) and is given by $\lambda_{\text{theory}}=\mathbb{E}[\|X^{T}\Psi\|_{\infty}]$ where $\Psi\sim\mathbb{N}(0,\hat{\sigma}^{2}I)$. Our second and third choices are obtained from cross-validation. Denoted by $\lambda_{\text{cv.min}}$ and $\lambda_{\text{cv.1se}}$, the tuning parameters are associated with the lowest cross-validated error and error within 1 standard error of the best model, respectively.

5.1 Methods and metrics

In our experiments, we illustrate maximum likelihood inference for two sets of parameter vectors after selection: (i) the partial regression coefficients in the selected model, obtained by letting $\mathcal{H}(\mathrm{S})=\mathrm{E}$ and $\mathcal{F}_{\mathrm{S}}=(X_{\mathrm{E}}^{\intercal}X_{\mathrm{E}})^{-1}X_{\mathrm{E}}^{\intercal}$; (ii) the selected set of parameters in the full model, obtained by letting $\mathcal{H}(\mathrm{S})=\{1,2,\cdots,p\}$ and $\mathcal{F}_{\mathrm{S}}=\mathcal{L}_{\mathrm{E}}(X^{\intercal}X)^{-1}X^{\intercal}$, where $\mathcal{L}_{\mathrm{E}}\in{}^{|\mathrm{E}|\times p}$ is a matrix of all zeros except for the indices

[TABLE]

We call the former parameter vector “Partial” and the latter “Full” in our depictions.

In the linear regression setting described above, we compare the relative performance of our proposed method with “Lee et al.”, “Liu et al.”, “Split”, and “Naive”. Our method follows Algorithm 1 with the barrier penalty

[TABLE]

where $U_{j}$ is the $j^{\text{th}}$ row of $U$ and $v_{j}$ is the $j^{\text{th}}$ component of $v$. First, we report the average coverage of the interval estimates produced by each strategy across all simulations. Each simulation records the proportion of intervals that cover our target parameters in a single round of experiment. The nominal level of FCR aimed by the interval estimates in $10\%$. Note, selective inference produced by “Liu et al.” is tied only to the “Full” parameter vector. The coverage for “Naive” intervals, not adjusted for any selection, underscores the extent of selection bias for a specific value of SNR. Next, we provide a breakdown of all methods in terms of their inferential power. To this end, we record the average lengths of the interval estimates and their power which we define to be the proportion of signals detected by a strategy from the ones successfully screened by the query. Detected variables here count the variables for which the corresponding interval estimates do not cover zero. Faced with two different regimes, we depict the power comparisons for the randomized and non-randomized estimates under certain best-case scenarios for each estimate. The best-case scenarios in our experiment are led by an assessment of predictive risks for a point estimate associated with each strategy, measured through the relative risk metric:

[TABLE]

where $\widehat{o}^{\lambda}$ is the estimate and $\beta$ is the parameter vector. We consider the LASSO solution as a natural point estimate for the parameter vector $\beta$ in the non-randomized regime and associate the LASSO estimate with both strategies “Lee et al.” and “Liu et al.”. For “Split”, we consider the least squares estimate obtained after refitting the selected model to the remaining one-third of the data samples and append it with zeros for the indices not selected by the LASSO. Note, the selective MLE, an immediate byproduct of Algorithm 1, appended with zeros for the inactive indices serves as a point estimate for our proposal. We discuss the detailed findings of our experiments next.

5.2 Findings and interpretation

Supporting the validity of inference after selection, Figure 5 highlights the performance of the different interval estimates for the “Full” and “Partial” parameters. The three columns in the plot are associated with the three different choices of tuning parameter: $\lambda_{\text{theory}}$, $\lambda_{\text{cv.1se}}$, and $\lambda_{\text{cv.min}}$. We remark that none of the methods adjust for the adaptivity involved in the choice of the cross-validated tuning parameters. Coherent with expectations, the interval estimates for all the methods approximately attain the nominal FCR level $10\%$ at $\lambda_{\text{theory}}$; “Lee et al.” fails to yield valid inference at $\lambda_{\text{cv.1se}}$ and $\lambda_{\text{cv.min}}$ since it does not account for the fact that the tuning parameters for the LASSO were chosen based on the specific data through cross-validation. The effectiveness of the normal approximation for the proposed method “MLE” is largely ascribed to the soft-truncated likelihood, due to the use of randomization in the LASSO query. Besides, the accuracy of the large deviations-type approximation for the exact likelihood is maintained under moderate dimensions. Interestingly, “MLE” and “Liu et al.” recover the nominal levels even at the cross validated tuning parameter; this is seen in the second and third columns of the plot. Although there lacks a formal justification for this observation, heuristically, the use of randomization in the LASSO query limits the role of the data-dependent regularization. A similar justification possibly holds for “Liu et al.”, which achieves the same goal by constructing inference for parameters that are less affected by selection.

Figure 6 highlights the averaged lengths of the interval estimates produced by different methods. The bars in red depict the percentage of intervals with infinite length for “Lee et al.”, confirming the conclusions in Kivaranovic and Leeb (2018). The interval estimates produced by all other methods are bounded in length. Consistent with the example presented in the introduction, the interval estimates based on “Liu et al.” and “MLE” are comparable for the “Full” parameters. The estimates for the “Partial” parameters produced by “MLE” are shorter than those for the “Full” parameters; this gain in power is in part due to inference in the learned model as opposed to the full model. Assuredly, the new proposal dominates the simple “Split” estimates which roughly use the same information in selection as the randomized query in our method.

Emphasized earlier, for a fair comparison of power between the randomized and non-randomized estimates, we use their relative risks to guide us to a best-case scenario within each regime. Purely from a predictive lens, the risk assessment across the different values of SNR suggest running the randomized LASSO at $\lambda_{\text{cv.1se}}$ and the canonical LASSO at $\lambda_{\text{cv.min}}$. Taking on direct comparisons for “MLE” after the randomized LASSO at $\lambda_{\text{cv.1se}}$, and “Liu et al.”, “Lee et al.”, “Split” after the usual LASSO at $\lambda_{\text{cv.min}}$, Figure 7 depicts their relative risks and power under the respective best-case situations. For the moderately high SNR values, the selective MLE proves to be a competing estimate when compared against the LASSO estimate. We note a far superior predictive performance of the LASSO at $\lambda_{\text{cv.min}}$ than the selective MLE in the lower range of SNR values. Our proposal, however, turns out to be a better choice for inference across the range of SNR values, outperforming the non-randomized alternatives in terms of power.

6 Real Data Example

We apply our method to investigate associations between gene expressions and patient survival times for Gliomas (a type of brain tumor) in the TCGA data. With survival times ranging between $1$ to $15$ years, and some of these tumors quickly progressing to Glioblastoma, genetic associations are increasingly utilized for prognostic decisions (e.g., Zhang et al., 2019; Panigrahi et al., 2020). In our analysis of $441$ samples, we use log-transformed survival times as our response. As potential predictors, we choose the top $2500$ predictors with the largest sample variation from a candidate pool of $17500$ molecular measurements of gene expression values (mRNAseq). Before running a meaningful LASSO query, we account for the presence of strongly correlated predictors by further pruning the $2500$ predictors to a subset of $140$ predictors. We do so by applying the hierarchical clustering scheme in Bien and Tibshirani (2011), followed by collecting the prototype representatives for each resulting cluster of predictors. We consider the following algorithms: (i) the LASSO; (ii) two runs of the LASSO; (iii) a marginal screening of the predictors followed by the LASSO, adding a Gaussian randomization variable $\omega\sim N(0,\eta^{2}I_{p})$ to the queries for our method. Consistent with the simulations, we fix $(\hat{\sigma})^{-2}\eta^{2}=1$. We use $\lambda_{\text{theory}}$ in Section 5 to tune the LASSO penalty. We conduct a marginal screening of variables at the nominal level $q=0.20$ and let the screening threshold be $\zeta=z_{1-q/2}\cdot\sqrt{\hat{\sigma}^{2}\text{diag}(X^{T}X)+\eta^{2}1_{p}}$ for the randomized version of this query.

Figure 8 showcases the distribution of the lengths of interval estimates produced by “MLE”, “Lee et al.” and “Liu et al.” after solving (i). For inference post (ii) and (iii), we compare “MLE $+$ double-LASSO” and “MLE $+$ MS-LASSO” against “Split $+$ LASSO” where half of the samples are reserved for inference. The conditional prescriptions in “Lee et al.” and “Liu et al.” do not directly apply to accommodate multiple queries at the time of selection. Corroborating our findings in the simulations, the lengths of the estimates using our proposal are way shorter than those based on “Liu et al.” and “Lee et al.”. Observe, “MLE $+$ double-LASSO” and “MLE $+$ MS-LASSO” outperform “Split $+$ LASSO” with shorter intervals, despite querying the data twice before inference.

7 Discussion

We investigate in the current paper a method for selective inference via maximum likelihood estimation. Amenable to a large class of convex queries at the time of selection, we rely on an optimization problem whose solution yields us estimating equations for the MLE and the observed Fisher information matrix, the two main ingredients for the proposed method. The estimating equations easily generalize to multiple convex queries at the time of selection and assume a separable form across the queries. The appeal of our method is two fold: (1) the computing costs for selective inference are reduced by orders of magnitude in comparison with MCMC sampling-based alternatives, and (2) statistical power for inference is preserved despite querying the data multiple times through randomized queries.

Future extensions of our proposal include a development of theory to use the method beyond Gaussian data. Along this direction, uniform guarantees for coverage (e.g., Leeb and Pötscher, 2005, 2006) require closer investigation. The framework for selective inference in the paper adjusts for queries with affine K.K.T. conditions at the solution. The ability of the proposal to account for learning queries that present non-affine representations for the K.K.T. conditions, e.g., the Group LASSO, remain to be explored in the future.

8 Acknowledgements

S.P. would like to sincerely thank and acknowledge Veera Baladandayuthapani and Yujia Pan for their inputs in the analysis of the TCGA dataset. S.P. is immensely thankful to Xuming He and Liza Levina for offering valuable comments on an initial draft of the paper. The authors thank the anonymous reviewers for their many insightful suggestions on earlier drafts of the paper.

Appendix A Proofs

Proof.

Proposition 2.1.  Recall,

[TABLE]

Then, we observe $\nabla^{2}\alpha(\beta)$ is equal to:

[TABLE]

From the above display, we note $\nabla^{2}\alpha(\beta)$ bounded from below by

[TABLE]

which implies that the function $\alpha$ is strongly convex. Let $\alpha^{*}$ denote the conjugate function for $\alpha:\real\to\real$. Based on the estimating equation for the selective MLE in Equation 5, we now have:

[TABLE]

Noting that the conjugate of a strongly convex function with index $K$ is $L$-Lipschitz smooth, where $L=K^{-1}$, we have:

[TABLE]

The finite sample bound in (21) now leads to our claim after taking conditional expectations on both sides of the inequality. ∎

Proof.

Theorem 1.   For a fixed value $\mathrm{S}$, denote the density of $\widehat{\beta}^{\perp}_{\mathrm{S}}$ by $\ell^{\perp}(\cdot)$ and the density of $\omega$ conditional on $\widehat{\beta}_{\mathrm{S}}$ and $\widehat{\beta}^{\perp}_{\mathrm{S}}$ by $\ell_{\mathbb{W}}(\cdot\;\lvert\;\widehat{\beta}_{\mathrm{S}},\widehat{\beta}^{\perp}_{\mathrm{S}})$. We note that the unconditional multivariate likelihood based on the joint density of $\widehat{\beta}_{\mathrm{S}}$, $\widehat{\beta}^{\perp}_{\mathrm{S}}$ and $\omega$ is proportional to:

[TABLE]

In obtaining (22), we have used the mutual independence between $\widehat{\beta}_{\mathrm{S}}$, $\widehat{\beta}^{\perp}_{\mathrm{S}}$ and $\omega$. Consider a change of variables mapping: $\omega\stackrel{{\scriptstyle\pi^{-1}}}{{\longrightarrow}}o,$ based on the invertible mapping $\pi$:

[TABLE]

Ignoring the constant Jacobian, the likelihood function based on the density of $\widehat{\beta}_{\mathrm{S}}$, $\widehat{\beta}^{\perp}_{\mathrm{S}}$, $o_{1}$ and $o_{2}$ agrees (up to proportionality constants) with:

[TABLE]

To complete the proof, let:

[TABLE]

and let

[TABLE]

Conditioning upon $\widehat{\mathrm{S}}(Y,W)=\mathrm{S}$ and $\widehat{\beta}^{\perp}_{\mathrm{S}}(Y)=\widehat{\beta}^{\perp}_{\mathrm{S}}$, the likelihood derived from the truncated (joint) density of $\widehat{\beta}_{\mathrm{S}}$ and $o_{1}$ is equal to:

[TABLE]

In the previous display, we use the fact that our conditioning event is equivalent to the $|\mathrm{E}|$-dimensional linear constraints: $Uo_{1}<v$. In the final step, we marginalize over the optimization variables $o_{1}$, which gives us the expression of our soft-truncated likelihood. ∎

Proof.

Proposition 3.1.   Based on the distribution of the Gaussian variables in the Proposition, we have the following bound:

[TABLE]

The bound above holds for all $\lambda_{1},\lambda_{2}$. Optimizing over $\lambda_{1},\lambda_{2}$ in (23) leads us to observe

[TABLE]

The penultimate step holds due to a minimax equality, because of the convexity and compactness of $\mathrm{R}$. The last claim follows by plugging in the log-moment generating function for the multivariate Gaussian distribution of the variables $\widehat{\beta}_{\mathrm{S}}$ and $O_{1}$, followed by computing its conjugate function at $\begin{pmatrix}{\widetilde{\beta}_{\mathrm{S}}}^{\;\intercal}&o_{1}^{\intercal}\end{pmatrix}^{\intercal}$. ∎

Proof.

Theorem 2.   Let $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}=\Sigma^{-1}(J\beta_{\mathbb{M}_{\mathrm{\mathrm{S}}},\mathrm{\mathrm{S}}}+k)$, and let $h^{*}(\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}},0)$ be the convex conjugate of

[TABLE]

when evaluated at $\begin{pmatrix}(\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}})^{\intercal}&0^{\intercal}_{|\mathrm{E}|}\end{pmatrix}^{\intercal}$. Based on the expression in Equation 15, note that our approximate log-selection probability, in terms of $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$ is equal to:

[TABLE]

Then, the approximate selective MLE for $\beta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$ is given by

[TABLE]

where

[TABLE]

In the remaining proof, we derive an expression for $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}$. Introducing duplicate variables $\eta^{\prime}_{\mathbb{M}_{\mathrm{S}},\mathrm{S}},u^{\prime}$, we rewrite the problem in (24) as

[TABLE]

with the linear constraints $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}=\eta^{\prime}_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$ and $u=u^{\prime}$. The dual of the optimization in the preceding display is equivalent to solving

[TABLE]

where the Lagrangian is as follows:

[TABLE]

Because,

[TABLE]

[TABLE]

[TABLE]

the dual problem further simplifies as:

[TABLE]

Denoting the optimal variables for $u$, $\zeta$ and $o_{1}$ by $u^{*}$, $\zeta^{*}$ and $o^{*}_{1}(\widehat{\beta}_{\mathrm{S}})$ respectively, the K.K.T. conditions of optimality yield us:

[TABLE]

where

[TABLE]

We conclude that $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}$ is obtained by looking at the first $|\mathrm{E}|$ coordinates of $\nabla h(\widehat{\beta}_{\mathrm{S}},o^{*}_{1}(\widehat{\beta}_{\mathrm{S}}))$. This gives us:

[TABLE]

which then results in the estimating equation:

[TABLE]

∎

Proof.

Theorem 3.    Let $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}=\Sigma^{-1}(J\beta_{\mathbb{M}_{\mathrm{\mathrm{S}}},\mathrm{\mathrm{S}}}+k)$. We simplify our approximate expression for the log-selection probability:

[TABLE]

by writing it equivalently as the following $|\mathrm{E}|$-dimensional optimization:

[TABLE]

where

[TABLE]

We obtain the above equivalence by optimizing over the unconstrained variables, $\widetilde{\beta}_{\mathrm{S}}$ in the former objective. The equivalent formulation then leads us to observe

[TABLE]

where $\Lambda^{*}$ is the conjugate function of

[TABLE]

Further, we note

[TABLE]

Clearly,

[TABLE]

is positive semi-definite, and $\alpha$ is strongly convex, satisfying specially,

[TABLE]

Lastly, using the fact that

[TABLE]

where $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}=\nabla\alpha^{-1}(\widehat{\beta}_{\mathrm{S}})$, we have:

[TABLE]

The ultimate inequality relies on the $L$-Lipschitz nature of $\alpha^{*}$ where $L\leq(\lambda_{0})^{-1}$; this fact is a direct consequence of the strong convexity noted for $\alpha$. Taking conditional expectations on both sides gives us the finite-sample bound in our claim. ∎

Proof.

Theorem 4.   Consider the notations in the proof of Theorem 2. Recall, the expression for our approximate log-likelihood is given by:

[TABLE]

The second derivative of the negative (log-) likelihood with respect to $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$ is equal to

[TABLE]

where $\widehat{\beta}^{*}_{\mathrm{S}}$, dependent on $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$, satisfies the following set of equations

[TABLE]

[TABLE]

Based on (28), taking a derivative with respect to $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$, we have:

[TABLE]

and therefore:

[TABLE]

We obtain an expression for $\dfrac{\partial o_{1}^{*}}{\partial\widehat{\beta}^{*}_{\mathrm{S}}}$ from (29) after taking a derivative with respect to $\widehat{\beta}^{*}_{\mathrm{S}}$, which gives us:

[TABLE]

Lastly, we note $\widehat{\beta}^{*}_{\mathrm{S}}=\widehat{\beta}_{\mathrm{S}}$ when $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}=\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}$, the optimizer defined in (24). Thus, \dfrac{\partial\widehat{\beta}^{*}_{\mathrm{S}}}{\partial\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}}(\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}})\Big{\lvert}_{\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}=\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}} is equal to

[TABLE]

where $o_{1}^{*}(\widehat{\beta}_{\mathrm{S}})$ is the solution of (29) when $\widehat{\beta}^{*}_{\mathrm{S}}=\widehat{\beta}_{\mathrm{S}}$. Using the reparameterization

[TABLE]

we deduce that observed Fisher information matrix, $I(\widehat{\beta}^{\;\text{mle}}_{\mathbb{M}_{\mathrm{S}},\mathrm{S}})$, equals

[TABLE]

∎

Proof.

Theorem 5.   Let $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}=\Sigma^{-1}(J\beta_{\mathbb{M}_{\mathrm{\mathrm{S}}},\mathrm{\mathrm{S}}}+k)$. Based upon the approximation in Equation 19, the MLE problem solves the optimization

[TABLE]

where $h^{*}(\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}},0)$ is the convex conjugate of $h$ defined as

[TABLE]

Our selective MLE is then given by:

[TABLE]

such that $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}$ is the minimizer of the optimization in (30). Our proof is complete by noting that the dual of the optimization in (30) is equal to:

[TABLE]

Observe, the objective in the dual problem is separable in the optimizing variables $o_{1}^{(1)},\cdots,o_{1}^{(L)}$. The K.K.T. conditions of optimality indicate that $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}$ is obtained by looking at the first $|\mathrm{E}|$ coordinates of

[TABLE]

where $o^{*(l)}_{1}(\widehat{\beta}_{\mathrm{S}})$ for $l\in\{1,2,\cdots,L\}$ are the optimal variables for the dual problem. This gives us the estimating equation for the approximate selective MLE. ∎

Proof.

Theorem 6.   Borrowing notations from the proof of Theorem 5, the Hessian of the negative log-likelihood

[TABLE]

with respect to $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}$ is equal to

[TABLE]

In the above display, $\widehat{\beta}^{*}_{\mathrm{S}}$ satisfies the following system of equations

[TABLE]

and

[TABLE]

Taking derivatives of equations (32) and (33), we have

[TABLE]

[TABLE]

Observe, $\widehat{\beta}^{*}_{\mathrm{S}}=\widehat{\beta}_{\mathrm{S}}$ when $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}=\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S}}^{*}$, the minimizer for the problem in (30). Therefore, we have

[TABLE]

where $o^{*(l)}_{1}(\widehat{\beta}_{\mathrm{S}})$ satisfies (33) when $\widehat{\beta}^{*}_{\mathrm{S}}=\widehat{\beta}_{\mathrm{S}}$. Finally, using the reparameterization

[TABLE]

we obtain our expression for the observed Fisher information matrix. ∎

Appendix B Soft-truncated likelihood: examples

In this section, we illustrate how our method applies to: (i) variable screening based on marginal correlations (Lee and Taylor, 2014); (ii) variable selection via SLOPE (Bogdan et al., 2015).

Example B.1**.**

In this example, we consider a randomized variable screening algorithm based on marginal correlations. Fix a level $q\in(0,1)$ and consider drawing a randomization variable $\omega\sim N(0_{p},\eta^{2}I_{p})$. Let

[TABLE]

Suppose, we solve:

[TABLE]

where $\mathcal{P}^{\zeta}(o)=\displaystyle\sum_{j}\chi_{\{|o_{j}|\leq\zeta_{j}\}}(o_{j})$, and $\chi_{R}(o)$ takes the value [math] if $o\in R$ and is infinity otherwise. Let the solution of (34) be given by: $o^{(\zeta)}$. Then, let $\mathrm{E}$ represent the indices of the selected variables, i.e., $\mathrm{E}=\{j:|o^{(\zeta)}_{j}|=\zeta_{j}\}$, and let $z_{E}$ be the signs of the vector: $X_{E}^{\intercal}y+\omega_{E}$. Let $o_{1}$ be the subgradient (subvector) of the penalty $\mathcal{P}^{\zeta}$ at the indices in $\mathrm{E}$, and let $o_{2}=X_{\mathrm{E}^{c}}^{\intercal}y+\omega_{\mathrm{E}^{c}}$.

Suppose, we condition on the observed value of $\mathrm{E}$ along with the additional conditioning on $z_{E}$ and $o_{2}$. That is, $\mathrm{S}$ is the observed active set, alongside the sign vector $\text{sign}(X_{E}^{\intercal}y+\omega_{E})$ and the observed vector $X_{\mathrm{E}^{c}}^{\intercal}y+\omega_{\mathrm{E}^{c}}$. Fixing $\mathcal{H}(\mathrm{S})=\mathrm{E}$, $\mathcal{F}_{\mathrm{S}}=(X_{\mathrm{E}}^{\intercal}X_{\mathrm{E}})^{-1}X_{\mathrm{E}}^{\intercal}\in\mathbb{R}^{|\mathrm{E}|\times n}$, we work with $\widehat{\beta}_{\mathrm{S}}$ and $\widehat{\beta}^{\perp}_{\mathrm{S}}$ defined within Section 3 for the randomized LASSO query. Fix the following matrices

[TABLE]

Then, the K.K.T mapping associated with the solution of (34) is

[TABLE]

Our event of selection is clearly equivalent to: $Uo_{1}<v$, for the fixed matrices $U=-\text{diag}(z_{\mathrm{E}})$, $v=0_{|\mathrm{E}|}$. Now, we can follow the exact same analysis as prescribed for the randomized LASSO query.

Example B.2**.**

We consider a randomized version of the SLOPE (Bogdan et al., 2015), given by:

[TABLE]

The randomization instance $\omega$ is drawn from a Gaussian distribution as our previous examples. The SLOPE penalty is an ordered version of the $\ell_{1}$ penalty. Assuming the $p$ tuning parameters are distinct, we let $\mathrm{E}$ be the indices of the selected variables with the sign vector $z_{\mathrm{E}}$ based on the solution of (35) . Because, the SLOPE penalty generates ties in the absolute values of (the nonzero components of) the solution vector, we let $|\breve{\mathrm{E}}|$ be the number of distinct nonzero components for the solution to (35). Consider the vector containing the magnitudes of these distinct nonzero components, $O_{1}\in{}^{|\breve{\mathrm{E}}|}$; that is,

[TABLE]

Suppose, each distinct solution $O_{1}^{(k)}$ is associated with the following cluster of variables: $\mathcal{A}_{k}$. Denote $\bar{X}_{k}=\sum_{k\in\mathcal{A}_{k}}z_{k}X_{k}$, where $z_{k}$ is the sign of the estimated solution for the $k^{\text{th}}$ variable. Let $O_{2}$ stand for the inactive components of the subgradient for the SLOPE penalty. We denote the observed instances for $O_{1}$ and $O_{2}$ by $o_{1}$ and $o_{2}$, respectively.

In this example, we define our conditioning event to be the set of realizations that lead us to observing the active set of variables $\mathrm{E}$, along with the signs $z_{\mathrm{E}}$ and the inactive subgradient $o_{2}$. Working with the specific choices: $\mathcal{H}(\mathrm{S})=\mathrm{E}$ and $\mathcal{F}_{\mathrm{S}}=(X_{\mathrm{E}}^{\intercal}X_{\mathrm{E}})^{-1}X_{\mathrm{E}}^{\intercal}$, define $\widehat{\beta}_{\mathrm{S}}$ and $\widehat{\beta}^{\perp}_{\mathrm{S}}$ as we have done in the previous example. We recognize that the K.K.T. mapping at the solution is given by:

[TABLE]

Let $V\in{}^{|\breve{\mathrm{E}}|-1\times|\breve{\mathrm{E}}|}$ be a matrix of all zeros, except the indices

[TABLE]

In particular, our selection event is equivalent to the following linear constraints on $o_{1}$: $Uo_{1}<v$ for the fixed matrices $U=\text{diag}(-I_{|\breve{\mathrm{E}}|},V)$, $v=0_{2|\breve{\mathrm{E}}|-1}$. We are now ready to apply Theorem 2 and 4 in order to obtain the estimating equations for the approximate selective MLE and observed Fisher information matrix, respectively.

Appendix C Asymptotic Properties of Selective MLE

We first turn to the univariate example in Section 2. An application of Proposition 2.1 results in a global consistency guarantee for the selective MLE that we formalize below. We consider the selection rule

[TABLE]

such that $\sqrt{n}\bar{Y}_{n}\sim N(\sqrt{n}\beta_{n},1)$ and $W\sim N(0,\eta^{2})$ is our randomization variable, independent of $\bar{Y}_{n}$. We denote the corresponding selective MLE by $\widehat{\beta}_{n}^{\text{\;mle}}$.

Proposition C.1**.**

Fix $\delta>0$. Then,

[TABLE]

as $n\to\infty$.

Proof.

Using the bound for the mean squared error of the selective MLE in Proposition 2.1, we have

[TABLE]

The inequality on the last display uses a reduction in variance when the variable is truncated to a convex region (Kanter and Proppe, 1977). The proof of consistency is now immediate by applying the Chebyshev’s inequality. ∎

We devote the remaining section to study the asymptotic guarantees for our approximate soft-truncated likelihood in Section 3. Establishing a convergence of our approximation in Equation 15 to the exact counterpart in Theorem 7, our main result in Theorem 8 proves a global consistency guarantee for the approximate selective MLE.

Suppose, we observe $\widehat{\mathrm{S}}_{n}=\mathrm{S}$; the subscript $n$ specifically indicates the dependence of $\widehat{\mathrm{S}}_{n}$ on the sample size. In the asymptotic setting, let

[TABLE]

for any fixed value $\mathrm{S}$. We denote our selective MLE by $\widehat{\beta}_{\mathbb{M}_{\mathrm{S}},\mathrm{S},n}^{\text{\;mle}}$. Consider the following sequence of parameters

[TABLE]

such that $n^{-1/2}\bar{b}_{n}=O(1)$, $\bar{b}_{n}\to\infty$ as $n\to\infty$. Based on the K.K.T. conditions of optimality, our optimization variables $O_{1,n}$ are now defined as

[TABLE]

where $\sqrt{n}W_{n}\sim N(0,\Sigma_{\mathbb{W}})$, and the selection event is equivalent to the linear constraints

[TABLE]

In a fixed $p$ and $n\to\infty$ regime, Theorem 7 verifies that the upper bound in Proposition 3.1 consistently approximates the exact selection probability in the following sense.

Theorem 7**.**

Assume

[TABLE]

converges to [math] as $n\to\infty$, whenever $v_{n}=O(1)$. Then, we have

[TABLE]

Before we provide the proof, note, Theorem 7 suggests approximating

[TABLE]

by

[TABLE]

Observe, the display on the right-hand side is equivalent to

[TABLE]

The approximate value in (40) coincides with the bound in Proposition 3.1, as we highlight the dependence of our variables and the parameter vector on $n$ in our asymptotic setting.

Proof.

We write the variables $\widehat{\beta}_{\mathrm{S},n}$ and $O_{1,n}$ as a sum of i.i.d. variables

[TABLE]

where $\begin{pmatrix}(\widehat{\beta}^{j}_{\mathrm{S},n})^{\intercal}&(O^{j}_{1,n})^{\intercal}\end{pmatrix}^{\intercal}$ are distributed as Gaussian variables with the covariance matrix:

[TABLE]

and the mean vector:

[TABLE]

Marginally, the covariance and mean for $O^{j}_{1,n}$ are equal to

[TABLE]

and

[TABLE]

respectively. For the remaining proof, we let $O^{j;(c)}_{1,n}$ be the centered version of $O^{j}_{1,n}$.

We first examine the limit of

[TABLE]

Applying a large deviations limit, we have

[TABLE]

which can be rewritten as

[TABLE]

Completing our proof, observe, the probability of selection simplifies as

[TABLE]

Using our assumption in (39), we have

[TABLE]

which proves the claim. ∎

Let $\mathcal{B}_{U;v}(\cdot)$ be a barrier penalty such that $(\bar{b}_{n})^{-2}\mathcal{B}_{U;v}(\bar{b}_{n}o_{1})$ converges to [math] in a pointwise sense only if $Uo_{1}\leq\frac{1}{\bar{b}_{n}}v$, and assumes the value infinity otherwise. Based on Theorem 7, the second term in the limit trivially agrees with

[TABLE]

as $n\to\infty$.

Theorem 8**.**

Fix $\delta>0$. Consider the assumptions in Theorem 7. Then

[TABLE]

as $n\to\infty$.

Proof.

Let $\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S},n}=\Sigma^{-1}\left(J\beta_{\mathbb{M}_{\mathrm{\mathrm{S}}},\mathrm{\mathrm{S}}}+\frac{1}{\bar{b}_{n}}k\right)$, let $\alpha_{n}(\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S},n})$ be given by

[TABLE]

and let $\bar{\alpha}_{n}(\eta_{\mathbb{M}_{\mathrm{S}},\mathrm{S},n})$ assume the value

[TABLE]

To simplify notations in the proof, let us denote the conditional expectation

[TABLE]

by $\overline{\mathbb{\mathrm{E}}}[\cdot]$. Using the finite sample bound in Theorem 3, we have

[TABLE]

The last display is obtained after using a triangle inequality. The first term on the right hand side is

[TABLE]

Applying the variance reduction result in Kanter and Proppe (1977), this term is further bounded by

[TABLE]

and converges to [math] as $n\to\infty$. The convergence of the second term uses the convexity and differentiability of the functions $\alpha_{n}$ and $\bar{\alpha}_{n}(\cdot)$ in addition to the fact

[TABLE]

as $n\to\infty$. These facts imply

[TABLE]

Our claim follows immediately by observing

[TABLE]

∎