2.1 Standard practice: variable selection with one covariate

For this section, we assume an outcome $Y_{i}$ and single covariate, $X_{i}$ , both scaled as described above. We are going to consider variable selection in this simplified setting, with model

(2) $$\begin{eqnarray}\displaystyle Y_{i}=X_{i}\unicode[STIX]{x1D6FD}^{o}+\unicode[STIX]{x1D716}_{i} & & \displaystyle\end{eqnarray}$$

where the goal is how to decide whether or not to zero out $\widehat{\unicode[STIX]{x1D6FD}}$ . We focus on this simplified setting because it provides analytic results not available in the multivariate setting, and we use these results to convey the basic intuitions of variable selection. Later we return to the multivariate setting.

With a single covariate, variable selection is normally done in two stages: first, the effect is estimated and then some $p$ -value or $t$ -statistic threshold for statistical significance is used to determine whether the effect can be differentiated from zero. A standard estimate for $\unicode[STIX]{x1D6FD}^{o}$ is the least squares estimate,

(3) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}^{LS} & = & \displaystyle \mathop{\text{argmin}}_{\widetilde{\unicode[STIX]{x1D6FD}}}\frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}\widetilde{\unicode[STIX]{x1D6FD}})^{2}\end{eqnarray}$$

(4) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{\mathop{\sum }_{i=1}^{N}Y_{i}X_{i}}{\mathop{\sum }_{i=1}^{N}X_{i}^{2}}=\frac{\mathop{\sum }_{i=1}^{N}Y_{i}X_{i}}{N-1}\end{eqnarray}$$

which gives a point estimate. The point estimate is then compared to its standard error, $\widehat{\unicode[STIX]{x1D70E}}_{\widehat{\unicode[STIX]{x1D6FD}}}:$

(5) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D70E}}_{\widehat{\unicode[STIX]{x1D6FD}}}=\frac{\widehat{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D716}}}{\sqrt{N-1}};\quad \widehat{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D716}}=\sqrt{\frac{\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}\widehat{\unicode[STIX]{x1D6FD}}^{LS})^{2}}{N-2}}. & & \displaystyle\end{eqnarray}$$

If the $t$ -statistic is larger in magnitude than some critical value, normally 1.96, the effect is considered statistically significant. In this framework, estimation and selection occur in two separate steps. Estimation and selection cannot, in fact, be simultaneous: the least squares estimate is never zero, outside of pathological cases.

2.2 LASSO with a single covariate

Next we introduce the LASSO of Tibshirani (Reference Tibshirani1996) in the case with a single covariate. LASSO is an acronym for Least Absolute Shrinkage and Selection Operator. In the one-parameter case, the LASSO estimator is the solution to

(6) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}^{L}=\mathop{\text{argmin}}_{\widetilde{\unicode[STIX]{x1D6FD}}}\frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}\widetilde{\unicode[STIX]{x1D6FD}})^{2}+\unicode[STIX]{x1D706}|\widetilde{\unicode[STIX]{x1D6FD}}|. & & \displaystyle\end{eqnarray}$$

The first term in this objective function is the residual sum of a squares. The second term has two components: the tuning parameter $\unicode[STIX]{x1D706}$ , indexed by sample size, and the penalty term $|\widetilde{\unicode[STIX]{x1D6FD}}|$ . The least absolute component references the structure of the penalty, $\unicode[STIX]{x1D706}|\widetilde{\unicode[STIX]{x1D6FD}}|$ .

Denote the sign of the least squares estimate as $\widehat{s}^{L}=\text{sgn}(\widehat{\unicode[STIX]{x1D6FD}}^{LS})\in \{-1,1\}$ . With one parameter, the LASSO estimate is (Tibshirani Reference Tibshirani1996, Section 2.2)

(7) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}^{L}=\left(\widehat{\unicode[STIX]{x1D6FD}}^{LS}-\widehat{s}^{L}\frac{\unicode[STIX]{x1D706}}{N-1}\right)\mathbf{1}\left(|\widehat{\unicode[STIX]{x1D6FD}}^{LS}|>\frac{\unicode[STIX]{x1D706}}{N-1}\right). & & \displaystyle\end{eqnarray}$$

The first part of Equation (7) illustrates the shrinkage: the estimate is the least squares estimate biased toward zero by $\unicode[STIX]{x1D706}/(N-1)$ . The second part of Equation (7) illustrates the selection component. If $\widehat{\unicode[STIX]{x1D6FD}}^{LS}$ is less than $\unicode[STIX]{x1D706}/(N-1)$ in magnitude, then the LASSO estimate is simply zero.

Equation (7) reveals three shortcomings of LASSO estimation. First, $\unicode[STIX]{x1D706}_{}$ must be estimated from the data. Researchers commonly turn to an AIC or BIC criterion or to cross-validation to select this tuning parameter. Unfortunately, these three methods may give quite different results, with little theoretical guidance to choose among them. Second, both large effects and small effects are biased toward zero by the same amount, $\unicode[STIX]{x1D706}/(N-1)$ . A more accurate estimator would shrink smaller effects but leave larger effects untouched. The third problem arises due to competing demands on $\unicode[STIX]{x1D706}/(N-1)$ . In the left side of the equation, it is a bias term. As such, we want it to go to zero quickly. In the right side of the equation, though, we see that if $\unicode[STIX]{x1D706}/(N-1)$ goes to zero too quickly, then we will end up not zeroing out any effects. As we show in Appendix A,Footnote ³ no LASSO model indexed by a single $\unicode[STIX]{x1D706}$ can produce estimates that are sparse, consistent, and asymptotically normal. The LASSOplus method proposed in this paper helps to resolve all three issues.

2.3 Single covariate adaptive LASSO

Zou (Reference Zou2006) introduced the adaptive LASSO, which extends the LASSO by including parameter-specific weights that cause differential shrinkage (see also Kenkel and Signorino Reference Kenkel and Signorino2012). These weights serve to shrink larger effects less than smaller effects. The adaptive LASSO is a two-stage estimator. In the first stage, weights that are inversely related to $\unicode[STIX]{x1D6FD}$ are constructed. In the second stage, these weights are used to weigh a LASSO problem. The weights are constructed as

(8) $$\begin{eqnarray}\displaystyle \widehat{w}=1/|\widehat{\unicode[STIX]{x1D6FD}}^{1}|^{\unicode[STIX]{x1D6FE}} & & \displaystyle\end{eqnarray}$$

where $\widehat{\unicode[STIX]{x1D6FD}}^{1}$ is a first-stage, root-N consistent estimate of $\unicode[STIX]{x1D6FD}^{o}$ and $\unicode[STIX]{x1D6FE}>0$ . With these weights, the estimator becomes

(9) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}^{L}(\unicode[STIX]{x1D706})=\mathop{\text{argmin}}_{\widetilde{\unicode[STIX]{x1D6FD}}}\frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}\widetilde{\unicode[STIX]{x1D6FD}})^{2}+\unicode[STIX]{x1D706}\widehat{w}|\widetilde{\unicode[STIX]{x1D6FD}}|. & & \displaystyle\end{eqnarray}$$

In this setting, the univariate adaptive LASSO estimator is

(10) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}^{aL}(\unicode[STIX]{x1D706},w)=\left(\widehat{\unicode[STIX]{x1D6FD}}^{LS}-s\frac{w\unicode[STIX]{x1D706}}{N-1}\right)\mathbf{1}\left(|\widehat{\unicode[STIX]{x1D6FD}}^{LS}|>\frac{w\unicode[STIX]{x1D706}}{N-1}\right). & & \displaystyle\end{eqnarray}$$

For a fixed value of $\unicode[STIX]{x1D706}$ , the adaptive LASSO estimator has a lower bias for larger effects than smaller ones. With one effect, the adjustment is trivial, but with hundreds of possible effects, the gains from differential shrinkage can be substantial.

Several problems emerge with adaptive LASSO estimation. The first is that any number of methods can return a root-N consistent estimate of $\unicode[STIX]{x1D6FD}$ : least squares, ridge regression, or Bayesian regression models (e.g., Gelman et al. Reference Gelman, Jakulin, Pittau and Su2008). Second, the decay parameter $\unicode[STIX]{x1D6FE}$ must be either assumed or estimated from the data, which can grow computationally costly. Third, the adaptive LASSO inherits the same uncertainty over tuning parameter selection as the LASSO.

2.4 An empirical process approach with a single variable

We turn next to LASSO estimation through an empirical process approach. Early seminal works include Donoho and Johnstone (Reference Donoho and Johnstone1994), Candes and Tao (Reference Candes and Tao2007), and Bickel et al. (Reference Bickel, Ritov and Tsybakov2009); the approach has been recently popularized in economics in work by Victor Chernozhukov and colleagues (Belloni and Chernozhukov Reference Belloni and Chernozhukov2013; Chernozhukov et al. Reference Chernozhukov, Fernández-Val and Melly2013; Belloni et al. Reference Belloni, Chen, Chernozhukov and Hansen2012, Reference Belloni, Chernozhukov and Hansen2011). Harding and Lamarche (Reference Harding and Lamarche2016) have extended this work to estimating individual-level heterogeneous effect in quantile regression panel models. We recommend Buhlmann and van de Geer (Reference Buhlmann and van de Geer2013) as a particularly accessible introduction to this line of work.

A central result in this literature is deriving an “Oracle Inequality,” a bound showing that for a particular value of $\unicode[STIX]{x1D706}$ , the excess risk goes to zero at rate $1/N$ , up to a penalty incurred for not knowing the true covariates. In the one-parameter case, this inequality is given below:

Lemma 1 (Oracle Inequality for the LASSO in the Single-Parameter Case).

For $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D70E}\times t\times \sqrt{(N-1)}$ , the single-parameter LASSO estimator satisfies the Oracle InequalityFootnote ⁴

(11) $$\begin{eqnarray}\displaystyle \frac{1}{N}\left\{\mathop{\sum }_{i=1}\left(X_{i}(\widehat{\unicode[STIX]{x1D6FD}}^{L}-\unicode[STIX]{x1D6FD}^{o})\right)^{2}+\unicode[STIX]{x1D706}|\widehat{\unicode[STIX]{x1D6FD}}^{L}-\unicode[STIX]{x1D6FD}^{o}|\right\}\leqslant \frac{C_{L1}\unicode[STIX]{x1D70E}^{2}t^{2}}{N} & & \displaystyle\end{eqnarray}$$

with probability at least $1-2\exp \{-t^{2}/2\}$ .

We denote as $t$ the variable by which we control the probability of the bound holding, i.e., the variable $t$ enters into both the bound $C_{L1}\unicode[STIX]{x1D70E}^{2}t^{2}/(N-1)$ and the probability of it holding $1-2\exp \{-t^{2}/2\}$ . We will use $C_{\cdot }$ to denote constants that do not change in $N$ or $K$ .

Though we state the result with $K=1$ , this approach is most useful in the $K>N$ setting, a point to which we return below.

2.5 The multivariate LASSO and its variants

We now turn to the multivariate setting, where $X_{i}^{\top }$ is a vector of $K$ observation-level covariates, $i\in \{1,2,\ldots ,N\}$ . These covariates may include the values of treatment variables, pretreatment moderators, and interactions within and between the two. We assume $X_{i}$ is of length $K$ , where $K$ is a finite but possibly large number, say in the hundreds or thousands. The $k$ th element of $X_{i}$ and $\unicode[STIX]{x1D6FD}$ are $X_{ik}$ and $\unicode[STIX]{x1D6FD}_{k}$ , and we also assume all fourth moments of $[Y_{i}X_{i}^{\top }]$ exist.Footnote ⁵ We assume the data are generated as

(12) $$\begin{eqnarray}\displaystyle Y_{i}=X_{i}^{\top }\unicode[STIX]{x1D6FD}^{o}+\unicode[STIX]{x1D716}_{i} & & \displaystyle\end{eqnarray}$$

where we desire a sparse representation of the $K$ -dimensional vector $\unicode[STIX]{x1D6FD}$ . We are not assuming that elements of $\unicode[STIX]{x1D6FD}^{o}$ are zero, but that some effects are of negligible magnitude.Footnote ⁶ Instead, we are seeking the best representation of the model in which most effects are estimated as zero, so as to allow the researcher to focus on relevant effects. We return to this point more fully below.

In the multivariate setting, the LASSO estimate is the solution to

(13) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}^{L}(\unicode[STIX]{x1D706})=\mathop{\text{argmin}}_{\widetilde{\unicode[STIX]{x1D6FD}}}\frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\widetilde{\unicode[STIX]{x1D6FD}})^{2}+\unicode[STIX]{x1D706}\mathop{\sum }_{k=1}^{K}|\widetilde{\unicode[STIX]{x1D6FD}}_{k}|. & & \displaystyle\end{eqnarray}$$

The LASSO tends to over-select small effects. Previous research has addressed this problem in two ways. First, first-stage adaptive weights can be incorporated into the estimation. This leads to the adaptive LASSO, given above. Alternatively, small coefficients can be trimmed ex post and OLS fit to the surviving covariates. This approach has been developed in the empirical process framework discussed above and we focus on the popular LASSO $+$ OLS method of Belloni and Chernozhukov (Reference Belloni and Chernozhukov2013) below and in our simulations.

We start with the first. In the multivariate setting, the adaptive LASSO gives each parameter its own weight in the penalty $\widehat{w}_{k}=1/|\widehat{\unicode[STIX]{x1D6FD}}_{k}^{0}|^{\unicode[STIX]{x1D6FE}}$ , for some $\unicode[STIX]{x1D6FE}>0$ and $\unicode[STIX]{x1D6FD}_{k}^{0}$ a root-N consistent estimate. The model is now

(14) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}^{aL}(\unicode[STIX]{x1D706})=\mathop{\text{argmin}}_{\widetilde{\unicode[STIX]{x1D6FD}}}\frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\widetilde{\unicode[STIX]{x1D6FD}})^{2}+\unicode[STIX]{x1D706}\mathop{\sum }_{k=1}^{K}\widehat{w}_{k}|\widetilde{\unicode[STIX]{x1D6FD}}_{k}|. & & \displaystyle\end{eqnarray}$$

We move now to the second approach. LASSO $+$ OLS proceeds in two steps. First, a LASSO fit is used to select an initial set of possible effects. Second, OLS is conducted on the all subsets of the initial fit, and the OLS fit with residual variance close to the residual variance of the LASSO fit is selected (Belloni and Chernozhukov Reference Belloni and Chernozhukov2013). From our experience, LASSO $+$ OLS serves as a helpful tool for sparse modeling and variable selection, so we include it both in our software package and in our simulations below. We provide a more complete discussion of LASSO $+$ OLS in the online supplemental materials.

These methods, pre-estimation weighting and post-estimation selection, raise several concerns. First, each have tuning parameters that must be selected and, ideally, estimated. As we describe in the online supplemental materials, LASSO $+$ OLS has three tuning parameters that must be selected. The authors provide reasonable defaults and algorithms on how to estimate the tuning parameters, but provide no theoretical guidance as to how to choose them. Similarly, with the adaptive LASSO, the user must choose a method for estimating the first-stage weights and for the exponent in the weights. A second area of concern is in generating confidence intervals. The probabilistic bounds in LASSO $+$ OLS generate confidence intervals (Belloni et al. Reference Belloni, Chernozhukov and Hansen2011), but require user-tuning of several parameters. For the adaptive LASSO methods, generating confidence intervals is still a field of active research.Footnote ⁷

2.6 Two statistical properties of sparse models

We now discuss two statistical properties of sparse models that have been described in the literature. We present these properties and their relationship to existing methods here in order to set up our discussion in Section 4 of how and when LASSOplus has these properties. Though the two properties have similar names, the Oracle Property and Oracle Inequality, they are actually distinct statistical concepts. The Oracle Property requires that, as $N$ grows, the estimator select the correct effects and converge in distribution to OLS using only variables from the true subset. This property was initially cast in the fixed- $K$ , growing- $N$ framework (Fan and Li Reference Fan and Li2001; Zou Reference Zou2006), but we focus on the Oracle Property in when $N$ and $K$ can both be thought of as growing (see also Fan and Peng Reference Fan and Peng2004). The Oracle Inequality is a predictive property, such that the fitted values are close to the fitted values from OLS estimated only on the true subset variables. Satisfying the Oracle Inequality is weaker than satisfying the Oracle Property, as estimating fitted values well is easier than estimating $\unicode[STIX]{x1D6FD}^{o}$ correctly. The key advantage to the Oracle Inequality is that it allows for analysis when $K>N$ .

2.6.1 The Oracle Property

The adaptive LASSO satisfies the Oracle Property of Fan and Li (Reference Fan and Li2001, p. 1353, Thm 2). An Oracle estimator must satisfy two conditions. First, it must be consistent for variable selection, so in-truth nonzero effects are selected and in-truth zero effects are zeroed out. Second, the estimate has to be consistent and asymptotically normal, with variance equal to a model fit only to covariates associated with the nonzero effects. In other words, the estimator achieves the same asymptotic performance as an identical estimator told ex ante by an “oracle” which effects are zero and which are not. We give a formal definition in Definition 1.

Definition 1 (The Oracle Property).

Assume the model $Y_{i}=X_{i}^{\top }\unicode[STIX]{x1D6FD}^{o}+\unicode[STIX]{x1D716}_{i}$ with $\unicode[STIX]{x1D716}_{i}$ mean zero with four finite moments. Let the set $S$ denote the set of indices of in-truth nonzero elements of $\unicode[STIX]{x1D6FD}$ , i.e. $S=\{k:k\in \unicode[STIX]{x1D6FD}_{k}^{o}\neq 0\}$ .

An Oracle estimator $\widehat{\unicode[STIX]{x1D6FD}}^{oracle}$ has the following two properties (Zou Reference Zou2006, p. 1418):

1. 1. Consistent Variable Selection: $\lim _{N\rightarrow \infty }\{k:\widehat{\unicode[STIX]{x1D6FD}}_{k}^{oracle}\neq 0\}=S$ .

2. 2. Optimal Estimation Rate: $\sqrt{N}(\widehat{\unicode[STIX]{x1D6FD}}_{S}^{oracle}-\unicode[STIX]{x1D6FD}_{S}^{o})\stackrel{\text{d}}{\rightarrow }{\mathcal{N}}(\mathbf{0}_{|S|},\unicode[STIX]{x1D6F4}_{S}^{\ast })$ where $\unicode[STIX]{x1D6F4}_{S}^{\ast }$ is the asymptotic variance matrix from the true subset model.

Satisfying the Oracle Property is desirable because it offers asymptotic efficiency gains over the normal least squares estimator.Footnote ⁸ Denote the asymptotic relative efficiency of two estimators of vector $\unicode[STIX]{x1D703}$ , $\widehat{\unicode[STIX]{x1D703}}_{1}$ and $\widehat{\unicode[STIX]{x1D703}}_{2}$ , as

(15) $$\begin{eqnarray}\displaystyle ARE_{\unicode[STIX]{x1D703}}(\widehat{\unicode[STIX]{x1D703}}_{1},\widehat{\unicode[STIX]{x1D703}}_{2})=\lim _{N\rightarrow \infty }\frac{\mathbb{E}\left\{\Vert \widehat{\unicode[STIX]{x1D703}}_{1}-\unicode[STIX]{x1D703}\Vert _{2}^{2}\right\}}{\mathbb{E}\left\{\Vert \widehat{\unicode[STIX]{x1D703}}_{2}-\unicode[STIX]{x1D703}\Vert _{2}^{2}\right\}}. & & \displaystyle\end{eqnarray}$$

The asymptotic relative efficiency of an estimator with the Oracle Property will never perform worse than the least squares estimate, a result we state below:

Proposition 1. Assume the least squares estimator exists and is unique. Then, an estimator with the Oracle Property is at least as efficient as the non-Oracle least squares estimator, asymptotically:

(16) $$\begin{eqnarray}\displaystyle ARE_{\unicode[STIX]{x1D6FD}^{o}}(\widehat{\unicode[STIX]{x1D6FD}}^{LS},\widehat{\unicode[STIX]{x1D6FD}}^{Oracle})\geqslant 1 & & \displaystyle\end{eqnarray}$$

with equality if and only if none of the elements of $\unicode[STIX]{x1D6FD}^{0}$ are 0.

Proof. See the online supplemental materials.

We note that the Oracle Property has been extended by Fan and Peng (Reference Fan and Peng2004). The authors consider the general case of a penalized likelihood and give conditions for the estimates to satisfy the Oracle Property. Among these are that the bias induced by shrinkage disappear and that the likelihood eventually swamp the prior. The key assumption, from an applied standpoint, is that $K^{5}/N\rightarrow 0$ . Of course this holds when $K$ is fixed. When $K$ is allowed to grow, the result illustrates that the Oracle Property—basically, getting the model right pointwise and in distribution—requires quite a bit of data.Footnote ⁹

2.6.2 The Oracle Inequality

Again, we take as our benchmark the least squares estimator fit to only the in-truth nonzero effects, with the goal of producing an estimator that has similar properties. The least squares estimator, fit only to true covariates, achieves a predictive risk of order $\mathbb{E}\big(\frac{1}{N}\Vert X_{S}(\unicode[STIX]{x1D6FD}_{S}^{o}-\widehat{\unicode[STIX]{x1D6FD}}_{S}^{LS})\Vert _{2}^{2}\big)=\unicode[STIX]{x1D70E}^{2}|S|/N$ , where $|S|$ denotes the number of in-truth nonzero effects. An Oracle Inequality bounds the predictive order at a rate going to zero as $1/N$ , so it performs comparably to OLS on the true subset. Often in these inequalities, though, the bound will include a penalty that grows in $K$ , as the true model is not known in advance.

The Oracle Inequality in the multivariate case requires two additional constructions. Denote as $X_{S}$ the subset of $X$ corresponding with in-truth nonzero effects and denote as $\unicode[STIX]{x1D719}_{o}$ the smallest eigenvalue of $\frac{1}{N}\sum _{i=1}^{N}X_{S}^{\top }X_{S}$ . The Compatibility Condition holds when $\unicode[STIX]{x1D719}_{o}>0$ .Footnote ¹⁰

Proposition 2 (Oracle Inequality for the LASSO in the Multivariate Case).

For

(17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}=C_{\unicode[STIX]{x1D716}}\widehat{\unicode[STIX]{x1D70E}}\sqrt{\left(t+\log (K)\right)\times (N-1)} & & \displaystyle\end{eqnarray}$$

the LASSO estimator satisfies the Oracle Inequality

(18) $$\begin{eqnarray}\displaystyle \frac{1}{N}\left\{\Vert X_{i}(\widehat{\unicode[STIX]{x1D6FD}}^{L}-\unicode[STIX]{x1D6FD}^{o})\Vert _{2}^{2}+\unicode[STIX]{x1D706}\Vert \widehat{\unicode[STIX]{x1D6FD}}^{L}-\unicode[STIX]{x1D6FD}^{o}\Vert _{1}\right\}\leqslant \frac{C_{L}\widehat{\unicode[STIX]{x1D70E}}^{2}\left(t+\log (K)\right)}{\unicode[STIX]{x1D719}_{0}^{2}N} & & \displaystyle\end{eqnarray}$$

with probability at least $1-2\exp \{-t^{2}/2\}-\Pr (C_{\unicode[STIX]{x1D716}}\widehat{\unicode[STIX]{x1D70E}}\leqslant \unicode[STIX]{x1D70E})$ .

Corollary. The LASSO risk is consistent for the population risk when $\log (K)/N\rightarrow 0$ .

Proof. See Buhlmann and van de Geer (Reference Buhlmann and van de Geer2013, ch. 6.2).

We focus on two insights from the Oracle Inequality. First, it is achieved when $\unicode[STIX]{x1D706}$ is of order $\sqrt{N\log (K)}$ . Doing so actually allows for closed form calculation of the tuning parameters. The results also highlights that we are paying a “cost” of $\log (K)$ for not knowing the true model in advance. The second is that the requirements are quite mild for consistency. For $K$ covariates, we only need the sample size $N$ to be of order $\log (K)$ . For example, then, going from $K=100$ to $K=200$ would only require a $15\%$ increase in sample size $(=\log (200)/\log (100))$ to maintain the Oracle Inequality bound. Relative to the Oracle Property, we see that the Oracle Inequality requires a good bit less data, but it also guarantees less than the Oracle Property.

2.7 The Bayesian LASSO

LASSOplus is a Bayesian estimator and thus with the above review of recent frequentist based LASSO methods in mind, we turn to the existing Bayesian LASSO literature. In a Bayesian framework, the LASSO can be interpreted as the maximum a posteriori (MAP) estimate of a model with a double-exponential prior $\Pr (\unicode[STIX]{x1D6FD}_{j}|\unicode[STIX]{x1D706})=\frac{1}{2\unicode[STIX]{x1D706}}\exp (-\unicode[STIX]{x1D706}|\unicode[STIX]{x1D6FD}_{j}|)=DE(\unicode[STIX]{x1D706})$ . The Bayesian LASSO model of Park and Casella (Reference Park and Casella2008, hereafter PC),Footnote ¹¹ can be written as

(19) $$\begin{eqnarray}\displaystyle Y_{i}|X_{i},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D70E}^{2} & {\sim} & \displaystyle {\mathcal{N}}(X_{i}^{\top }\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D70E}^{2})\end{eqnarray}$$

(20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}|\unicode[STIX]{x1D706},\unicode[STIX]{x1D70E}^{2} & {\sim} & \displaystyle DE(\unicode[STIX]{x1D706}/\unicode[STIX]{x1D70E})\end{eqnarray}$$

(21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{2} & {\sim} & \displaystyle \unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D6FF},\unicode[STIX]{x1D70C})\end{eqnarray}$$

where we denote as $DE(a)$ the double-exponential density $f(x;a)=\frac{a}{2}\exp (-a|x|)$ . PC show that parameterizing the prior on $\unicode[STIX]{x1D6FD}_{k}$ with $\unicode[STIX]{x1D706}/\unicode[STIX]{x1D70E}$ instead of $\unicode[STIX]{x1D706}$ ensures a unimodal posterior. The prior on the tuning parameter is over $\unicode[STIX]{x1D706}^{2}$ rather than $\unicode[STIX]{x1D706}$ in order to maintain conjugacy in the augmented model, given below. Any positive value for the shape ( $\unicode[STIX]{x1D6FF}$ ) and rate parameters ( $\unicode[STIX]{x1D70C}$ ) will give a proper prior; PC take $(\unicode[STIX]{x1D6FF},\unicode[STIX]{x1D70C})=(1,1.78)$ . PC complete the hierarchy by assuming $\Pr (\unicode[STIX]{x1D70E}^{2})\propto 1/\unicode[STIX]{x1D70E}^{2}$ , though any gamma prior on $1/\unicode[STIX]{x1D70E}^{2}$ will maintain conjugacy.

The posterior mode of the Bayesian LASSO is a LASSO estimator. The negative log-posterior of $\unicode[STIX]{x1D6FD}$ under this model is, up to an additive constant that does not depend on $\unicode[STIX]{x1D6FD}$ ,

(22) $$\begin{eqnarray}\displaystyle -\!\log \left(\Pr \left(\unicode[STIX]{x1D6FD}|\unicode[STIX]{x1D706},\unicode[STIX]{x1D70E}^{2},{\mathcal{D}}_{N}\right)\right) & = & \displaystyle \frac{1}{\unicode[STIX]{x1D70E}^{2}}\left\{\frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\unicode[STIX]{x1D6FD})^{2}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D70E}\mathop{\sum }_{k=1}^{K}|\unicode[STIX]{x1D6FD}_{k}|\right\}\end{eqnarray}$$

(23) $$\begin{eqnarray}\displaystyle & \propto & \displaystyle \frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\unicode[STIX]{x1D6FD})^{2}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D70E}\mathop{\sum }_{k=1}^{K}|\unicode[STIX]{x1D6FD}_{k}|.\end{eqnarray}$$

Factoring out $\frac{1}{\unicode[STIX]{x1D70E}^{2}}$ reveals the posterior mode of the PC model is exactly a LASSO estimate with tuning parameter $\widetilde{\unicode[STIX]{x1D706}}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D70E}$ .

2.8 Shrinkage priors and scale mixtures

The normal likelihood above is not conjugate with the double-exponential prior on $\unicode[STIX]{x1D6FD}$ . In order to restore conjugacy, PC augment the parameter space by representing the double-exponential distribution as a scale mixture of normals with exponential mixing density (see also West Reference West1987):

(24) $$\begin{eqnarray}\displaystyle \frac{\widetilde{\unicode[STIX]{x1D706}}}{2}e^{-\widetilde{\unicode[STIX]{x1D706}}|\unicode[STIX]{x1D6FD}_{k}|}=\int _{0}^{\infty }\frac{1}{\sqrt{2\unicode[STIX]{x1D70B}\unicode[STIX]{x1D70F}_{k}^{2}}}e^{-\unicode[STIX]{x1D6FD}_{k}^{2}/(2\unicode[STIX]{x1D70F}_{k}^{2})}\frac{\widetilde{\unicode[STIX]{x1D706}}^{2}}{2}e^{-\widetilde{\unicode[STIX]{x1D706}}^{2}\unicode[STIX]{x1D70F}_{k}^{2}/2}d\left(\unicode[STIX]{x1D70F}_{k}^{2}\right). & & \displaystyle\end{eqnarray}$$

This suggests the following augmented representation of the double-exponential prior:

(25) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}\sim DE(\unicode[STIX]{x1D706}/\unicode[STIX]{x1D70E})\Rightarrow \unicode[STIX]{x1D6FD}_{k}|\unicode[STIX]{x1D70F}_{k}^{2},\unicode[STIX]{x1D70E}^{2}\sim {\mathcal{N}}(0,\unicode[STIX]{x1D70F}_{k}^{2}\unicode[STIX]{x1D70E}^{2});\;\;\unicode[STIX]{x1D70F}_{k}^{2}\sim \exp (\unicode[STIX]{x1D706}^{2}/2). & & \displaystyle\end{eqnarray}$$

Under the augmented parameterization, the likelihood and prior for $\unicode[STIX]{x1D6FD}$ are both normal and hence conjugate. Let $D_{\unicode[STIX]{x1D70F}}=\text{diag}(\unicode[STIX]{x1D70F}^{2})$ and $A=(X^{\top }X+D_{\unicode[STIX]{x1D70F}}^{-1})^{-1}$ . The Gibbs updates are:

(26) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}|\cdot & {\sim} & \displaystyle {\mathcal{N}}(AXY,\unicode[STIX]{x1D70E}^{2}A)\end{eqnarray}$$

(27) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}^{2}|\cdot & {\sim} & \displaystyle \text{InvGamma}\left((N-1)/2+K/2,\frac{1}{2}\left\{\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\unicode[STIX]{x1D6FD})^{2}+\mathop{\sum }_{k=1}^{K}\unicode[STIX]{x1D6FD}_{k}^{2}/\unicode[STIX]{x1D70F}_{k}^{2}\right\}\right)\end{eqnarray}$$

(28) $$\begin{eqnarray}\displaystyle 1/\unicode[STIX]{x1D70F}_{k}^{2}|\cdot & {\sim} & \displaystyle \text{InvGaussian}\left(\unicode[STIX]{x1D706}\unicode[STIX]{x1D70E}/|\unicode[STIX]{x1D6FD}_{k}|,\unicode[STIX]{x1D706}^{2}\right)\end{eqnarray}$$

(29) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{2}|\cdot & {\sim} & \displaystyle \unicode[STIX]{x1D6E4}(K+\unicode[STIX]{x1D6FF},\mathop{\sum }_{k=1}^{K}\unicode[STIX]{x1D70F}_{k}^{2}/2+\unicode[STIX]{x1D70C})\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}(a,b)$ denotes a Gamma distribution with shape parameter $a$ and rate parameter $b$ .

Additional methods have implemented different mixing densities within the scale-mixture representation of shrinkage priors (Polson and Scott Reference Polson and Scott2012). The “horseshoe prior” of Carvalho et al. (Reference Carvalho, Polson and Scott2010) is

(30) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}|\unicode[STIX]{x1D706},\unicode[STIX]{x1D706}_{k} & {\sim} & \displaystyle {\mathcal{N}}(0,\unicode[STIX]{x1D706}_{k}^{2}\unicode[STIX]{x1D706}^{2})\end{eqnarray}$$

(31) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{k} & {\sim} & \displaystyle C^{+}(0,s)\end{eqnarray}$$

where $s$ is taken as either $1$ or $\unicode[STIX]{x1D70E}^{2}$ and $C^{+}(a,b)$ denotes the half-Cauchy distribution (Gelman Reference Gelman2006). The model is so-named because the posterior density places most of its mass at either no shrinkage or full shrinkage, giving the posterior density a horseshoe shape. The horseshoe prior has proven to be an excellent default choice in sparse modeling, so we include it in our simulation study below.Footnote ¹²

2.9 Shrinkage without selection

The Bayesian estimators have shown better performance, in terms of mean-squared error and prediction, than their frequentist counterparts; see Kyung et al. (Reference Kyung, Gill, Ghosh and Casella2010) as well as our simulations below. These estimators, however, are not sparse. By a sparse Bayesian model, we mean one where either the mean, median, or mode of the conditional posterior density of $\Pr (\unicode[STIX]{x1D6FD}_{k}|\cdot )$ takes on a value of zero with nonzero probability. By the Bernstein–Von Mises theorem, we know that the data will swamp the prior and the posterior density will converge to the same density as the maximum likelihood estimate. The maximum likelihood estimate for a continuous parameter is never sparse, outside of pathological cases, and therefore neither are these Bayesian methods. In order to move toward a sparse estimate, the LASSOplus estimator “slows down” the rate at which the prior disappears asymptotically.

Under existing Bayesian methods, variable selection occurs in one of two ways. One, the variable may be selected by examining the posterior density of the estimated effect size; see Figure 8 or Hahn and Carvalho (Reference Hahn and Carvalho2015, Section 3). Two, effects may be selected off a summary statistic of the posterior density. Kyung et al. (Reference Kyung, Gill, Ghosh and Casella2010) propose fitting a frequentist LASSO such that the sums of absolute values of the frequentist method agree with the posterior mean of the sum of absolute values of the parameters. Carvalho et al. (Reference Carvalho, Polson and Scott2010) suggest a threshold for selecting off the parameter-specific weight parameters, a process shown to satisfy the Oracle Property (Datta and Ghosh Reference Datta and Ghosh2013). Hahn and Carvalho (Reference Hahn and Carvalho2015) suggest fitting a nonsparse Bayesian model and then selecting a frequentist model closest to these fitted values.

LASSOplus selects effects off a statistic of the posterior density. It is a sparse approximation to an underlying posterior density, constructed to achieve the Oracle Property.

3.1 The consistent model

We constructed the prior structure for the LASSOplus with two goals in mind. First, the log-posterior takes the same functional form as an adaptive LASSO problem. We show this property below. Second, the posterior mean of $\unicode[STIX]{x1D706}$ grows as $N^{1/4}K^{1/2}$ . We show in the next section how this growth rate helps the LASSOplus estimator achieve the Oracle Property and satisfy and Oracle Inequality.

The consistent model for LASSOplus can be written as

(32) $$\begin{eqnarray}\displaystyle Y_{i}|X_{i},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D70E}^{2} & {\sim} & \displaystyle {\mathcal{N}}(X_{i}^{\top }\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D70E}^{2})\end{eqnarray}$$

(33) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}|\unicode[STIX]{x1D706},w_{k},\unicode[STIX]{x1D70E} & {\sim} & \displaystyle DE\left(\unicode[STIX]{x1D706}w_{k}/\unicode[STIX]{x1D70E}\right)\end{eqnarray}$$

(34) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{2}|N,K & {\sim} & \displaystyle \unicode[STIX]{x1D6E4}\left(K\left(\sqrt{N}-1\right),\unicode[STIX]{x1D70C}\right)\end{eqnarray}$$

(35) $$\begin{eqnarray}\displaystyle w_{k}|\unicode[STIX]{x1D6FE} & {\sim} & \displaystyle \text{generalized}\;\text{Gamma}(1,1,\unicode[STIX]{x1D6FE})\end{eqnarray}$$

(36) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE} & {\sim} & \displaystyle \exp (1)\end{eqnarray}$$

with the generalized Gamma density $f(x;a,d,p)=\frac{p/a^{d}}{\unicode[STIX]{x1D6E4}(d/p)}x^{d-1}\exp \{-(x/a)^{p}\}$ .

The prior on $\unicode[STIX]{x1D6FD}$ is a reweighted version of that in the PC model. The tuning parameter, $\unicode[STIX]{x1D706}$ , was constructed to grow in $N$ . This growth is evident in the prior on $\unicode[STIX]{x1D706}^{2}$ : we replace the $\unicode[STIX]{x1D6FF}$ parameter in the PC model with $K(\sqrt{N}-1)$ . Any value $\unicode[STIX]{x1D70C}>0$ returns a proper prior; we take $\unicode[STIX]{x1D70C}=1$ . The Gamma prior on $\unicode[STIX]{x1D706}^{2}$ returns the Gibbs update in Formula (39). Lastly, the priors on the weights were derived so that joint posterior of $(\unicode[STIX]{x1D6FD},\{w_{k}\}_{k=1}^{K})$ would resemble the adaptive LASSO model. To see this, note that up to an additive constant that does not depend on $\unicode[STIX]{x1D6FD}$ or the weights,

(37) $$\begin{eqnarray}\displaystyle -\!\log \left(\Pr \left(\unicode[STIX]{x1D6FD},\{w_{k}\}_{k=1}^{K}|\unicode[STIX]{x1D706},\unicode[STIX]{x1D70E}^{2},\unicode[STIX]{x1D6FE}\right)\right)=\frac{1}{\unicode[STIX]{x1D70E}^{2}}\left\{\frac{1}{2}\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\unicode[STIX]{x1D6FD})^{2}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D70E}\mathop{\sum }_{k=1}^{K}w_{k}|\unicode[STIX]{x1D6FD}_{k}|\right\}+\mathop{\sum }_{k=1}^{K}w_{k}^{\unicode[STIX]{x1D6FE}} & & \displaystyle\end{eqnarray}$$

which combines the elements of equations (8) and (9).

Figure 1.

Comparison of Priors. This figure compares our unconditional LASSOplus prior to other priors (left) and shows the behavior of our conditional LASSOplus prior under different values of $\unicode[STIX]{x1D6FE}$ (right). The $x$ -axis contains values of $\unicode[STIX]{x1D6FD}$ and the $y$ -axis contains the prior density, $\Pr (\unicode[STIX]{x1D6FD})$ , on a log scale. We also include the two limiting distributions of the $t$ density, the normal and Cauchy, as well as the LASSO and horseshoe priors. The sparse priors (LASSO, horseshoe, LASSOplus) place greater mass at zero than the normal, an effect particularly pronounced for the horseshoe and LASSOplus priors. The normal prior places the least mass, and enacts the most shrinkage, on large effects. LASSOplus has the fattest tails, and thereby will have the least impact on posterior inference for large effects. The right-hand plot shows how the decay parameter $\unicode[STIX]{x1D6FE}$ affects the conditional prior, $\text{Pr}(\unicode[STIX]{x1D6FD}|\unicode[STIX]{x1D6FE})$ . The larger $\unicode[STIX]{x1D6FE}$ , the more aggressive the shrinkage. For smaller $\unicode[STIX]{x1D6FE}$ , the more the model is adjusting to large effects. As $\unicode[STIX]{x1D6FE}$ increases, less mass is placed in the tails and more toward zero. The parameter $\unicode[STIX]{x1D6FE}$ is estimated endogenously within the model.

Our model differs from earlier implementations of adaptive weights from Leng et al. (Reference Leng, Tran and Nott2014), Alhamzawi et al. (Reference Alhamzawi, Yu and Benoit2012), Griffin and Brown (Reference Griffin and Brown2012, Reference Griffin and Brown2010), Kang and Guo (Reference Kang and Guo2009) by placing a prior over $w_{k}^{\unicode[STIX]{x1D6FE}}$ rather than $w_{k}$ and $\unicode[STIX]{x1D6FE}$ separately. Like existing methods, we estimate the decay parameter $\unicode[STIX]{x1D6FE}$ from the data. As we show in Figure 1, adjusting the parameter allows the model to adapt to the global level of sparsity in the data. Taking $w_{k}=1$ for all $k$ returns the Bayesian LASSO model prior structure on $\unicode[STIX]{x1D6FD}_{k}$ ;Footnote ¹³ this model is the implementation of LASSOplus that does not have the Oracle Property.

3.2 Estimation

Estimation is nearly identical to the augmented PC model. We augment the model as

(38) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}|\unicode[STIX]{x1D706},w_{k},\unicode[STIX]{x1D70E}\sim DE(\unicode[STIX]{x1D706}w_{k}/\unicode[STIX]{x1D70E})\Rightarrow \unicode[STIX]{x1D6FD}_{k}|\unicode[STIX]{x1D70F}_{k}^{2},\unicode[STIX]{x1D70E}^{2},w_{k}^{2}\sim {\mathcal{N}}(0,\unicode[STIX]{x1D70F}_{k}^{2}\unicode[STIX]{x1D70E}^{2}/w_{k}^{2});\;\;\unicode[STIX]{x1D70F}_{k}^{2}\sim \exp (\unicode[STIX]{x1D706}^{2}/2). & & \displaystyle\end{eqnarray}$$

There are only two adjustments to the PC Gibbs sampler:

(39) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{2}|\cdot & {\sim} & \displaystyle \unicode[STIX]{x1D6E4}\left(K\sqrt{N},\mathop{\sum }_{k=1}^{K}\unicode[STIX]{x1D70F}_{k}^{2}/2+\unicode[STIX]{x1D70C}\right)\end{eqnarray}$$

(40) $$\begin{eqnarray}\displaystyle 1/\unicode[STIX]{x1D70F}_{k}^{2}|\cdot & {\sim} & \displaystyle \text{InvGaussian}\left(\unicode[STIX]{x1D706}w_{k}\unicode[STIX]{x1D70E}/(|\unicode[STIX]{x1D6FD}_{k}|),\unicode[STIX]{x1D706}^{2}w_{k}^{2}\right).\end{eqnarray}$$

We update $w_{k}$ and $\unicode[STIX]{x1D6FE}$ using a Griddy Gibbs sampler (Tierney Reference Tierney1994).

3.3 The LASSOplus estimator

As with existing methods, the data-generating process above does not return a sparse mode. The LASSOplus estimate is constructed from the estimate $\unicode[STIX]{x1D6FD}_{k}$ and a thresholding function that zeroes out sufficiently small values of $|\unicode[STIX]{x1D6FD}_{k}|$ . The threshold was constructed such that the final estimate achieves the Oracle Property, a point we return to after defining the estimator itself.

In order to guarantee that we zero out effects in the limit of $N$ , we sample an inflated variance component, $\unicode[STIX]{x1D70E}_{sp}^{2}$

(41) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{sp}^{2}|\cdot \sim \text{InvGamma}\left(\left(N^{1-2\unicode[STIX]{x1D6FC}}-1\right)/2+K/2,\frac{1}{2}\left\{\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\unicode[STIX]{x1D6FD})^{2}+\mathop{\sum }_{k=1}^{K}\unicode[STIX]{x1D6FD}_{k}^{2}/\unicode[STIX]{x1D70F}_{k}^{2}\right\}\right) & & \displaystyle\end{eqnarray}$$

that will enter into the threshold function. The parameter $\unicode[STIX]{x1D70E}_{sp}$ , which is central to the theoretical properties of LASSOplus, grows approximately as $N^{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D70E}_{}$ . We implement the model at $\unicode[STIX]{x1D6FC}=1/4$ , as this value achieves several desirable theoretical properties as we explain below.

The LASSOplus estimate is constructed from the consistent model and the inflated variance term $\unicode[STIX]{x1D70E}_{sp}$ . Define as

(42) $$\begin{eqnarray}\displaystyle V_{i}^{k}=Y_{i}-X_{i,-k}^{\top }\unicode[STIX]{x1D6FD}_{-k} & & \displaystyle\end{eqnarray}$$

the outcome less the estimated values from all effects except the $k$ th. Next, denote the conditional least squares estimate $\widehat{\unicode[STIX]{x1D6FD}}_{k}^{ols}$ as

(43) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D6FD}}_{k}^{ols}=\frac{\mathop{\sum }_{i=1}^{N}X_{ik}V_{i}^{k}}{\mathop{\sum }_{i=1}^{N}X_{ik}^{2}}. & & \displaystyle\end{eqnarray}$$

Conditional on all other parameters in both models, the LASSOplus estimate for the $k$ th element is then defined as

(44) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}^{plus}|\cdot =\unicode[STIX]{x1D6FD}_{k}\mathbf{1}\left(|\widehat{\unicode[STIX]{x1D6FD}}_{k}^{ols}|\geqslant \frac{\unicode[STIX]{x1D706}_{}\unicode[STIX]{x1D70E}_{sp}w_{k}}{N-1}\right). & & \displaystyle\end{eqnarray}$$

As LASSOplus is a Bayesian model, it returns estimates of the full posterior for all parameters. In the examples below, we select and focus on effects for which the median of the posterior density is nonzero, $\text{med}(\unicode[STIX]{x1D6FD}_{k}^{plus}|\cdot )\neq 0$ . We show that selecting off the median LASSOplus estimates serves as a conservative and powerful rule when trying to identify nonzero effects.

3.4 Approximate confidence intervals

A crucial contribution of LASSOplus is uncertainty estimates for model parameters. While the LASSOplus method returns posterior uncertainty estimates, the following discussion shows how to calculate approximate confidence intervals. We focus on confidence intervals because credible intervals are not calibrated to achieve nominal coverage, except in the limit. And as with earlier work, we found coverage to be subnominal with credible intervals (Kyung et al. Reference Kyung, Gill, Ghosh and Casella2010).

In returning approximate confidence intervals, we sample from the approximate sampling distribution of the LASSOplus estimator. To do so, we approximate Equation (44) as

(45) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}^{plus}|\cdot \approx \unicode[STIX]{x1D6FD}_{k}\boldsymbol{\unicode[STIX]{x1D6F7}}\left(\left|\frac{\widehat{\unicode[STIX]{x1D6FD}}_{k}^{ols}}{\widehat{\unicode[STIX]{x1D70E}}_{k}}\right|\geqslant \frac{\unicode[STIX]{x1D706}_{}\unicode[STIX]{x1D70E}_{sp}w_{k}}{\widehat{\unicode[STIX]{x1D70E}}_{k}(N-1)}\right) & & \displaystyle\end{eqnarray}$$

with $\widehat{\unicode[STIX]{x1D70E}}_{k}$ the variance of $\unicode[STIX]{x1D6FD}_{k}^{ols}$ . We apply the delta method to estimate $\unicode[STIX]{x1D70E}_{ci}^{2}$ (see the online supplemental materials for details and Efron (Reference Efron2015) for a more general argument).

The asymptotic normal approximation of the delta method may not hold with small samples. To correct for this, we estimate the error degrees of freedom using Satterthwaite’s approximation,

(46) $$\begin{eqnarray}\displaystyle df^{ci}=\frac{\left(\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\unicode[STIX]{x1D6FD})^{2}\right)^{2}}{\mathop{\sum }_{i=1}^{N}(Y_{i}-X_{i}^{\top }\unicode[STIX]{x1D6FD})^{4}}. & & \displaystyle\end{eqnarray}$$

With an estimate of error degrees of freedom in hand, we exploit the representation of a $t$ -density as a scale mixture of normals with inverse-gamma mixing density. We draw $v_{ci}\sim \text{InvGamma}(df^{ci}/2,df^{ci}/2)$ and use this value to inflate $\unicode[STIX]{x1D70E}_{ci}^{2}$ .Footnote ¹⁴

The sampling density of a nonzero univariate LASSO estimate is truncated-normal, conditional on the sign of the mode (Potscher and Leeb Reference Potscher and Leeb2009). Our approximation of the sampling density is then

(47) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{k}^{ci}|\cdot \sim \left\{\begin{array}{@{}l@{}}{\mathcal{N}}\left(\unicode[STIX]{x1D6FD}_{k}^{plus},\unicode[STIX]{x1D70E}_{ci}^{2}v_{ci}\right);\quad \unicode[STIX]{x1D6FD}_{k}^{plus}=0\quad \\ {\mathcal{T}}{\mathcal{N}}\left(\unicode[STIX]{x1D6FD}_{k}^{plus},\unicode[STIX]{x1D70E}_{ci}^{2}v_{ci},0,\infty \right);\quad \unicode[STIX]{x1D6FD}_{k}^{plus}>0\quad \\ {\mathcal{T}}{\mathcal{N}}\left(\unicode[STIX]{x1D6FD}_{k}^{plus},\unicode[STIX]{x1D70E}_{ci}^{2}v_{ci},-\infty ,0\right);\quad \unicode[STIX]{x1D6FD}_{k}^{plus}<0\quad \end{array}\right. & & \displaystyle\end{eqnarray}$$

where ${\mathcal{T}}{\mathcal{N}}(\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E}^{2},l,u)$ denotes the truncated-normal density with mean $\unicode[STIX]{x1D707}$ , variance $\unicode[STIX]{x1D70E}^{2}$ , and support on $(l,u)$ .

The approximate confidence interval is taken from the posterior distribution of $\unicode[STIX]{x1D6FD}_{k}^{ci}$ . For $K$ discovered effects, we take $\widetilde{K}=\max (K,1)$ and approximate the $1-\unicode[STIX]{x1D6FC}_{0}\%$ confidence interval as

(48) $$\begin{eqnarray}\displaystyle CI_{\unicode[STIX]{x1D6FC}_{0},\widetilde{K}}=\left(q_{\unicode[STIX]{x1D6FC}_{0}/(2\widetilde{K})},q_{1-\unicode[STIX]{x1D6FC}_{0}/(2\widetilde{K})}\right) & & \displaystyle\end{eqnarray}$$

where $q_{\widetilde{\unicode[STIX]{x1D6FC}}}$ is the estimated $\widetilde{\unicode[STIX]{x1D6FC}}$ quantile of $\unicode[STIX]{x1D6FD}_{k}^{ci}$ , with a Bonferroni correction for the discovered effects. Benjamini and Yekutieli (Reference Benjamini and Yekutieli2005, esp. 74–75) show that implementing a Bonferroni correction off discovered effects will maintain at least nominal coverage across all discovered effects.

3.5 Random effects and uncertainty estimates for designs using repeated unit-level observations

When researchers have repeated units in their sample, ignoring within-unit correlation can produce incorrect uncertainty estimates. Our substantive application, which uses a conjoint experiment by design, features repeated observations at the unit level. Often researchers will utilize some form of clustered standard errors. LASSOplus implements unit-level random effects for the same purpose.Footnote ¹⁵

Specifically, assume observation $i\in \{1,2,\ldots ,N\}$ as above. Now, assume each observation was generated by experimental unit $j$ , $j\in \{1,2,\ldots ,J\}.$ The function $j[i]$ maps each observation back to one of the experimental units (Gelman and Hill Reference Gelman and Hill2007).

We include random effects $u_{i}$ as

(49) $$\begin{eqnarray}\displaystyle u_{i} & = & \displaystyle a_{j[i]}\end{eqnarray}$$

(50) $$\begin{eqnarray}\displaystyle a_{j[i]} & {\sim} & \displaystyle {\mathcal{N}}(0,\unicode[STIX]{x1D70E}_{a}^{2}).\end{eqnarray}$$

We take the Jeffreys’ prior $1/\unicode[STIX]{x1D70E}^{2}$ on $\unicode[STIX]{x1D70E}^{2}$ , though a folded- $t$ density may be used (Gelman Reference Gelman2006).

3.6 Parametric extensions

Most experimental studies implement a linear (mean) model, given its connection to causal estimands. Researchers may prefer alternative models however, such as a probit model for a binary outcome. In this case, LASSOplus models this alternative data-generating process in a straightforward way by using the latent variable representation of the probit model (Albert and Chib Reference Albert and Chib1993). Briefly, the probit regression models the probability of a positive outcome as

(51) $$\begin{eqnarray}\displaystyle \Pr (Y_{i}=1|X_{i},\unicode[STIX]{x1D6FD})=\unicode[STIX]{x1D6F7}(X_{i}^{\top }\unicode[STIX]{x1D6FD}) & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6F7}(a)$ representing the cumulative distribution for a standard normal random variable. An observation-specific random variable is introduced, $z_{i}^{\ast }$ , and the problem transforms to

(52) $$\begin{eqnarray}\displaystyle z_{i}^{\ast }=X_{i}^{\top }\unicode[STIX]{x1D6FD}+e_{i} & & \displaystyle\end{eqnarray}$$

with $e_{i}\stackrel{\text{i.i.d.}}{{\sim}}{\mathcal{N}}(0,1)$ . Conditional on $z_{i}^{\ast }$ , the probit regression is a least squares problem with known variance, and LASSOplus can be applied as above. The latent variable representation also allows straightforward extension to ordinal, categorical, and censored data (Jackman Reference Jackman2009).

3.7 Expectation–maximization implementation

While we have certainly worked to produce an efficient implementation of LASSOplus, the MCMC method given above may be time-consuming and computationally expensive for large datasets with a large number of possible effects. For practical research and initial model fitting, we also implement an EM version of the LASSOplus model. We take $\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D70E}^{2}$ as the parameters to be estimated and the remaining parameters, $\{\unicode[STIX]{x1D706},\{1/\unicode[STIX]{x1D70F}_{k}^{2}\}_{k=1}^{K},\{w_{k}\}_{k=1}^{K},\unicode[STIX]{x1D6FE}\}$ , as “missing.”Footnote ¹⁶ As we have calculated all of the conditional posterior densities or kernels, the EM implementation is straightforward. We defer the details to the online supplemental materials.

4.1 Preliminary results

We present here two preliminary results. First, we consider the role of the weights, $\widehat{w}_{k}$ . Each weight is associated with a parameter in the penalty, where we want to downweight the penalty associated with large effects and upweight the penalty associated with small effects. We first show that the weights and magnitude of the effect estimates are indeed inversely related.

Lemma 2. The partial derivative of $\widehat{w}_{k}$ with respect to $|\widehat{\unicode[STIX]{x1D6FD}}_{k}|$ is

(53) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\widehat{w}_{k}}{\unicode[STIX]{x2202}|\widehat{\unicode[STIX]{x1D6FD}}_{k}|}=-\widehat{\unicode[STIX]{x1D706}}\sqrt{\widehat{\frac{1}{\unicode[STIX]{x1D70E}^{2}}}}\text{Var}(w_{k}|\cdot ) & & \displaystyle\end{eqnarray}$$

which is strictly less than zero.

This will allow us to associate the largest weights with the smallest effect estimates, and vice versa.

Second, we bound $\widehat{\unicode[STIX]{x1D706}}$ :

The bound on $\widehat{\unicode[STIX]{x1D706}}$ is a crucial component of both proofs below, as the growth rate of $\widehat{\unicode[STIX]{x1D706}}$ and $\widehat{w}_{k}$ determine the LASSOplus-EM estimator and, hence, its statistical properties. We turn next to the first result.

4.2 Oracle Property and Oracle Inequality for LASSOplus

We now give conditions on $N,K$ for which LASSOplus achieves the Oracle Property. We then show that both the estimator from the consistent model and LASSOplus each satisfy their own Oracle Inequality.

4.3 Comparison to existing priors

The Oracle Property and Oracle Inequality are both theoretical results. As we designed LASSOplus for use on real data, we next move on to finite-sample consideration by examining the behavior of the prior structure over the parameter weights $w_{k}$ . This enables us to compare the prior used in LASSOplus to other priors used in the literature.

There is, of course, no prior structure that performs well in all situations and for all data sets. We have generated a prior structure with four properties. First, the prior is concentrated at zero. This is appropriate for a setting where the researcher confronts hundreds or thousands of effects and wants to winnow these down to a small subset of relevant ones. Second, the prior places a large probability on the existence of large effects. For example, a standard normal prior places a $5.0\%$ prior probability on observing a value larger than 1.96 in magnitude. For a Cauchy prior, this value is $30.0\%$ ; for the horseshoe, $20.0\%$ , and for LASSOplus, $31.39\%$ . The more mass in the tails, the less posterior inference on large effects will be impacted by the prior. Third, the decay parameter $\unicode[STIX]{x1D6FE}$ allows the prior to adjust to the level of sparsity implied by the data. Fourth, the prior is less informative than several existing sparse priors. Conditional on assuming a sparse model, we want a prior that drives posterior inference as little as possible.

We illustrate the properties of the LASSOplus prior in Figure 1.Footnote ¹⁸ In each plot, the $y$ -axis contains the prior probability on a log scale and the $x$ -axis contains the magnitude of the effect, $|\unicode[STIX]{x1D6FD}|$ . The left figure plots the unconditional LASSOplus prior

(58) $$\begin{eqnarray}\displaystyle \text{Pr}(\unicode[STIX]{x1D6FD})=\int \Pr (\unicode[STIX]{x1D6FD}|\unicode[STIX]{x1D6FE})\Pr (\unicode[STIX]{x1D6FE})d\unicode[STIX]{x1D6FE} & & \displaystyle\end{eqnarray}$$

against the normal, LASSO (double exponential), Cauchy, and horseshoe priors.

The left-hand plot in Figure 1 illustrates how the different priors will handle large, intermediate, and small effects. The LASSOplus, LASSO, and horseshoe priors all concentrate at zero relative to the normal and Cauchy. The LASSOplus and horseshoe have the most pronounced spike at zero, and therefore will be the most aggressive in shrinking small effects to zero. This suggests that the LASSOplus and horseshoe should make the fewest false positive discoveries. LASSOplus is also relatively aggressive in shrinking intermediate effects to zero. For example, the normal prior places a high mass on effects less than 2 in magnitude, and will therefore shrink those the least. For larger values, say much larger than 4 in magnitude, the normal prior places a vanishingly low probability and will therefore shrink effect estimates quite a bit. The horseshoe and LASSO both place a higher mass on intermediate effects less than approximately 4 in magnitude. This suggests that the horseshoe and LASSO are more likely to discover intermediate effects, and that the LASSOplus will be relatively conservative in this range. We find this behavior in the simulations below. The horseshoe, LASSOplus, and Cauchy all have fatter tails than the LASSO and normal. Of these three, the LASSOplus has the fattest tails, and therefore will have the least impact on posterior inference for large effects.

The right-hand plot in Figure 1 shows how the decay parameter $\unicode[STIX]{x1D6FE}$ affects the conditional prior, $\text{Pr}(\unicode[STIX]{x1D6FD}|\unicode[STIX]{x1D6FE})$ . For a larger $\unicode[STIX]{x1D6FE}$ , the more aggressive the shrinkage is. For smaller $\unicode[STIX]{x1D6FE}$ , the more the model is adjusting to large effects. As $\unicode[STIX]{x1D6FE}$ increases, less mass is placed in the tails and more toward zero. The parameter $\unicode[STIX]{x1D6FE}$ is estimated endogenously within the model.

Another way to evaluate different prior distributions is to compare the amount of information they contribute to posterior inference. Practitioners are naturally wary of results driven by prior assumptions. A long strand of research has sought to identify reference priors such that the resultant posterior is dominated by the data (see particularly Bernardo Reference Bernardo1979; Jaynes Reference Jaynes1982; Berger and Bernardo Reference Berger and Bernardo1989; Bernardo Reference Bernardo, Dey and Rao2005; Berger et al. Reference Berger, Bernardo and Sun2009; Berger Reference Berger2006). A standard measure of the information in a prior $p$ with support ${\mathcal{B}}$ is its entropy, with higher values indicating less prior influence on the posterior:Footnote ¹⁹

(59) $$\begin{eqnarray}\displaystyle H(p)=-\int _{{\mathcal{B}}}p(\unicode[STIX]{x1D6FD})\log \left(p(\unicode[STIX]{x1D6FD})\right)d\unicode[STIX]{x1D6FD}. & & \displaystyle\end{eqnarray}$$

The reference prior is the prior from a class of densities that maximizes the entropy. For a single-parameter, asymptotically normal posterior, the reference prior is the Jeffreys’ prior.

We present the entropy of the priors in Table 1. The unconditional LASSOplus has the highest entropy among sparse priors (LASSO, horseshoe). The unconditional LASSOplus even has higher entropy than the Cauchy, which is proper but has no finite moments. Among the conditional LASSOplus priors, smaller values of $\unicode[STIX]{x1D6FE}$ correspond with a larger entropy. In the absence of information about expected effect size and the underlying level of sparsity, we prefer the conditional LASSOplus prior with the global decay parameter $\unicode[STIX]{x1D6FE}$ estimated from the data.

Table 1.

Comparison of Entropy Across Priors. From an objective Bayesian viewpoint, priors with a higher entropy are favored. Higher-entropy priors, intuitively, add less information to the model and have less impact on posterior inference. We present the entropy of the priors shown in Figure 1. The unconditional LASSO plus has the highest entropy among the unconditional priors. Among the conditional LASSOplus priors, smaller values of $\unicode[STIX]{x1D6FE}$ correspond with a larger entropy. In the absence of information about expected effect size, we prefer the conditional LASSOplus prior with the global decay parameter $\unicode[STIX]{x1D6FE}$ estimated from the data.

5.1 Subgroup analysis through repeatedly split sample analyses

In experimental analyses, subgroup analysis often consists of analysis of repeatedly split samples. For example, authors may conduct a long series of subgroup analyses by splitting the sample into different groups based on one dichotomous or dichotomized pretreatment covariate at a time. A regression model is fit to each group, and then the marginal effects of the different effects are evaluated within each subsample.

While showing these interactions between treatment conditions and pretreatment covariates is theoretically interesting, and often requested by audience members and reviewers, the methodology employed has a number of limitations. The decision to split subgroups one at a time introduces both conceptual and statistical concerns. Conceptually, the decision to conduct multiple separate subgroup analyses means that we implicitly acknowledge that we have estimated the wrong model. Take for example the simple moderating effects of two separate pretreatment variables on a treatment effect. The fact that the moderating effect of one of these variables is not included when estimating the moderating influence of another variable implies an awareness of potential misspecification. Statistically, the decision to conduct multiple separate subgroup analyses by discarding part of the data each time means the estimates are inefficient. This can be remedied by allowing for interactions between the treatment condition and the covariate in the same model. But the proliferation of parameters means that the standard regression framework is ill-equipped to deal with these situations. The basic reason is simple: the number of parameters quickly proliferates and hence some sort of stabilization, such as through sparsity, is necessary. We illustrate the advantages of sparsity-inducing priors below.

5.2 Interpreting interaction terms and the sparsity assumption

We use LASSOplus to estimate saturated models containing tens, hundreds, or even thousands of interaction terms. Common, and correct, practice proscribes including interaction terms in the absence of lower-order terms. We agree. Our implementation of LASSOplus, at its default, fits a model with all lower-order and interaction terms.

We also want to address the concern that the proposed methodology works only under the assumption that most of the effects are in truth zero. This assumption may arise from the frequentist implementation of the LASSO, whose proponents argue for using the “bet on sparsity principle” (Hastie et al. Reference Hastie, Tibshirani and Friedman2010, p. 611–613), arguing that sparse models are to be preferred over dense models. However, we do not assume that the true effects are zero in any sense. Our prior places mass 0 at the point 0. This differs from spike-and-slab priors, where the researcher implicitly places some prior mass on each parameter being zero (Mitchell and Beauchamp Reference Mitchell and Beauchamp1988; O’Hara and Silanapaa Reference O’Hara and Silanapaa2009; Gill Reference Gill2014, ch. 4.6). Instead, we generate a summary of the posterior that takes the value zero with some nonzero probability, allowing the data to tell us ex post that some effect is negligible. We are not assuming that the true value, even after the model is fit, is zero—instead, we seek the best sparse representation of the outcome in terms of main and subgroup effects.

Furthermore, researchers interested in characterizing every effect are able to do so. The pseudoposterior density of each effect is not sparse, as illustrated in our applied example. Researchers interested in evaluating nonselected effects are able to do so by evaluating this full posterior density. An important step in any analysis is to look to the data to determine which variables have a nonnegligible effect. LASSOplus provides an answer to this question by returning posterior median estimates of zero for negligible variables.

Lastly, the researcher may be nervous in interpreting higher-order interaction terms when the lower-order terms are not selected. The reason is that standard interaction terms cannot be interpreted without referencing their lower-order terms. The problem of interpreting interactions arises because of a correlation between the lower-order terms and the interaction term: the effect of one cannot be considered independently of the other (Esarey and Summer Reference Esarey and Summer2015). To solve this problem, we include interactions terms that are uncorrelated with their lower-order terms. To do so, first we construct the interaction term through elementwise multiplication, regress this term on its lower-order terms, and enter the residuals from this regression into LASSOplus. Under this construction, the effects of interaction terms can be interpreted as the effect above and beyond any lower-order effects. For a proof, see the online supplemental materials.

5.3 Scope of method for subgroup analysis

The development of the LASSOplus method to, in part, facilitate subgroup analysis had two goals. The first is to allow the experimentalist to uncover potentially relevant subgroup effects after implementing an experiment. Second, we designed the method to apply from the simplest of experiments to a conjoint analysis with repeated observations. The analysis of experimental data normally occurs in two steps. An experiment is designed to test a set of prespecified hypotheses. Upon completing an experiment, these hypotheses are then tested, and the point estimates and $p$ -values for each effect are reported.

Our method is designed for use in the subsequent analysis. After the inferential stage, researchers may be interested in higher-order effects, including treatment heterogeneity, when two treatments have an interactive effect; treatment effect heterogeneity, when the effect of a treatment varies with a covariate; or targeted treatment effect, when the effect of a treatment interaction varies with a covariate (Imai and Ratkovic Reference Imai and Ratkovic2013). Common practice involves repeatedly subsetting the data on pretreatment covariates and running separate regressions in each subset. For example, the researcher may find no treatment effect on average, but may find effects of opposite signs for males and females in her experimental data. The problem is the sheer number of effects quickly overwhelms, and researcher-driven repeated subsetting quickly devolves into interaction fishing.

The proposed method provides a means for considering all possible interaction terms and returning a sparse subset estimated as nonzero. The estimation is postinferential and descriptive, returning the effects that seem pronounced in this experimental data and might be considered for further study in the next experiment.

6.1 Setup

Our simulation is motivated by conjoint experiments, where the researcher wants to search through a large number of treatment/covariate interactions. The simulation mimics an experiment with 3 treatments, having 2, 3, and 4 levels, respectively. We also assume a set of $p$ pretreatment covariates, with $p\in \{10,25,50,75,100\}$ . The design matrix consists of the matrix of treatment indicators, $T$ , the pretreatment covariates, $X^{0}$ , and their residualized interactions. We designate the first level of each treatment as the baseline, dropping it and all of its interactions from $X$ . After creating all interactions, we are left with a design matrix $X$ with one of $\{76,181,356,531,706\}$ covariates. We run simulations with sample sizes $N\in \{50,100,250,500,1000,2500\}$ . For 16/30 of our simulation setups ( $53\%$ ), we have more observations than covariates.

We assume the following model:

(60) $$\begin{eqnarray}\displaystyle E(Y_{i}|X_{i},T_{i1},T_{i2},T_{i3}) & \propto & \displaystyle 3X_{i2}+3X_{i3}+3X_{i4}\end{eqnarray}$$

(61) $$\begin{eqnarray}\displaystyle & & \displaystyle +\,2\times \mathbf{1}(T_{i2}=b)+2\times \mathbf{1}(T_{i3}=b)+2\times \mathbf{1}(T_{i3}=c)\end{eqnarray}$$

(62) $$\begin{eqnarray}\displaystyle & & \displaystyle +\,X_{i2}\times \mathbf{1}(T_{i2}=b)+X_{i2}\times \mathbf{1}(T_{i3}=c)+X_{i2}\times \mathbf{1}(T_{i2}=a)\end{eqnarray}$$

(63) $$\begin{eqnarray}\displaystyle & & \displaystyle +\,X_{i4}\times \mathbf{1}(T_{i1}=b)+X_{i3}\times \mathbf{1}(T_{i3}=b)+X_{i4}\times \mathbf{1}(T_{i3}=c)\end{eqnarray}$$

with noise from a $t_{5}$ density and the systematic component is scaled to give a true $R^{2}$ of 0.5. The variables in $X_{i}$ are drawn from a $\mathit{Wishart}_{p}(N)$ , and each treatment condition is equiprobable. Each simulation setting was run 1000 times. The simulation design mimics a situation where there a few large main effects, some medium-sized average treatment effects, and the remaining interaction terms are small relative to the other effects.

6.2 Alternative estimators

We compare LASSOplus to three alternative sparse estimators: the frequentist LASSO, frequentist adaptive LASSO, and the horseshoe estimator. The LASSO and adaptive LASSO are fit using glmnet from R package glmnet. First-stage estimates come from ridge regression with the tuning parameter selected to minimize cross-validated error. The horseshoe prior is implemented in STAN.Footnote ²⁰ We coded up an implementation of LASSO $+$ OLS with all tuning parameters set at the defaults suggested by the original authors.

For details on the implementation of additional methods, see the online supplemental materials.