Bias-Aware Inference in Regularized Regression Models

Bias-Aware Inference in Regularized Regression Models

Introduction

We are interested in estimation and inference on a scalar coefficient $\beta$ in a linear regression model

when the $k$-vector $z_i$ of controls is large. In such settings, the classic \ac{OLS} estimator is often uninformative, exhibiting variance that is too large; the estimator is not even defined when $k>n$. This motivates modifying the \ac{OLS} objective function to penalize large values of $\gamma$, thereby lowering variance at the cost of introducing bias.

The most popular of these approaches is the lasso tibshirani_regression_1996 or other variants of $\ell_{1}$ penalization \citep*[e.g.][]{CaTa07,BeChWa11}. There is a large literature buhlmann_statistics_2011 showing favorable \ac{MSE} properties of these estimators when $\gamma$ is sparse. For inference, several papers have proposed \acp{CI} based on “double lasso” estimators belloni_inference_2014,javanmard_confidence_2014,van_de_geer_asymptotically_2014,zhang_confidence_2014, with asymptotic justification relying on rate conditions for the sparsity of $\gamma$. However, in many applications in economics, the sparsity assumption is not compelling and may be hard to motivate. Furthermore, it is unclear what sparsity level this approach implicitly imposes in a given finite sample.

We propose bounding the magnitude of the control coefficients, rather than their sparsity level, by assuming that $\operatorname{Pen}(\gamma)\leq C$. The penalty function $\operatorname{Pen}(\cdot)$ formalizes the notion of magnitude, and it can incorporate any restrictions on $\gamma$ that place it in a convex symmetric set. Such restrictions arise naturally in a plethora of applications. For instance, the dimensionality of the control vector is often large due to the inclusion of additional controls that are collectively believed to only be weakly associated with the outcome, but are nonetheless included to purge any possible confounding. One can then take the penalty to be an $\ell_{p}$ norm for the additional controls. If $z_i'\gamma$ is a basis approximation to some smooth function, we can define $\operatorname{Pen}(\gamma)$ to incorporate bounds on the derivatives of this function. The regularity parameter $C$ plays a role analogous to a sparsity bound.

We obtain sharp finite-sample results deriving near-optimal estimators and \acp{CI} under this penalty constraint and the idealized assumption that the regression errors $\varepsilon_{i}$ are Gaussian with a known homoskedastic variance. We show that the class of estimators that exactly resolves the trade-off between worst-case bias and variance can be obtained by (1) running a penalized propensity score regression of $w_i$ on $z_i$ using $\operatorname{Pen}(\cdot)$ as the penalty function; and then (2) using the residuals from this regression as an instrument in the univariate regression of $Y_{i}$ on $w_{i}$. \acp{CI} based on these estimators can be constructed by using a critical value that incorporates the worst-case bias of the estimator, which we show obtains automatically as a byproduct of the regularized regression in step (1). Because these \acp{CI} are bias-aware---they account for the potential finite-sample bias of the estimator---they are valid in finite samples in this idealized Gaussian setup. We show how to choose the weight on the penalty function in step (1) to optimize the \ac{MSE} of the resulting estimator, or the length of the resulting \ac{CI}.

In the more realistic setting with heteroskedasticity and an unknown error distribution, bias-aware \acp{CI} can be formed using heteroskedasticity-robust variance estimators, and we give conditions for their asymptotic validity, allowing for high-dimensional asymptotics with $k\gg n$. Our setup can allow for the effect of $w_{i}$ on the outcome to be heterogeneous, either by including interactions of the treatment and demeaned covariates among the controls, or by reinterpreting $\beta$ in (ref) as a weighted average treatment effect. We show that the treatment weights solve a bias-variance tradeoff in a problem where we can pick the estimand to make the estimation problem as easy as possible.

We also employ the high-dimensional asymptotics to study rates of convergence of the bias-aware \acp{CI} when $\operatorname{Pen}(\gamma)$ is an $\ell_p$ norm. We show that if $k\gg n$ and $C$ does not shrink with $n$, the optimal \ac{CI} shrinks more slowly than $n^{-1/2}$, so that the bias term asymptotically dominates; accounting for bias in \ac{CI} construction thus cannot be avoided even in large samples. Furthermore, we show that, in the $\ell_1$ case, this rate cannot be improved even if one additionally imposes the same $\ell_1$ bound in the propensity score regression of $w_i$ on $z_i$, as well as a certain degree of sparsity in both regressions.

Explicit specification of the regularity parameter $C$ that bounds the magnitude of $\gamma$ is a key input for our approach. Our efficiency bounds show that it is impossible to automate the choice of $C$ when forming \acp{CI}. We discuss how relating the magnitude of $\gamma$ to other quantities, such as the magnitude of the control coefficients in a short regression that only includes baseline controls, can help guide its choice. We develop a rule of thumb specification for $C$ based on this idea that we use in our simulations and empirical application. Robustness of the results can be assessed by computing a breakdown value of $C$, its largest value such that the empirical finding of interest, such as rejecting a particular null hypothesis, holds. Selection of $C$ cannot be automated due to the impossibility of getting a sufficiently informative data-driven upper bound for it. We show, however, that it is possible to obtain a lower \ac{CI} for $C$, which can be used as a specification check to ensure that the chosen value is not too low.

The requirement to explicitly choose $C$ may seem like a limitation of our approach relative to sparsity-based approaches, where the analogous tuning parameter, the degree of sparsity, does not need to be explicitly specified. However, good finite-sample performance of such methods relies on bounding these tuning parameters implicitly, and such implicit bounds are hard to calculate or evaluate in a given problem. We demonstrate this issue in a Monte Carlo analysis, where we show that double-lasso \acp{CI} suffer from moderate to severe undercoverage even in designs that are apparently sparse. Indeed, we view the explicit specification of $C$ as an advantage of our approach, because our coverage guarantees and efficiency bounds are based on transparent assumptions rather than “asymptotic promises” about tuning parameters that are hard to evaluate in a particular sample.

Our results relate to several strands of literature. Our procedures and efficiency bounds apply the general theory of estimation and inference on linear functionals in convex Gaussian models developed in IbKh85, donoho94, low95 and ArKo18optimal, and add to a growing literature applying this approach to various settings, including ArKo20ate,ArKo20sensitivity, kolesar_inference_2018, ImWa19, RaRo19, noack_bias-aware_2019, and kwon_inference_2020. muralidharan_factorial_2020 apply the approach in the present paper to experiments with factorial designs and bounds on interaction effects.

The idea of using propensity score residuals to estimate $\beta$ goes back at least to the work of robinson_root-n-consistent_1988 on the partly linear model. We provide a novel finite-sample justification for this idea, as well as an exact result giving the optimal penalization of this regression. Our setup allows for a general form of $\operatorname{Pen}(\cdot)$, and yields existing estimators in a few special cases; the bias-aware \acp{CI} to accompany such estimators are novel. First, we recover the optimal linear estimators in heckman_minimax_1988, who considered the partly linear model with a penalty function bounding the first or second derivative of a univariate nonparametric regression function. Next, we reproduce the result in li82 that the optimal estimator uses ridge regression when the penalty corresponds to an $\ell_{2}$ norm. Finally, li_linear_2020 consider the weighted $\ell_2$ norm $\operatorname{Pen}(\gamma)=\left(\sum_{i=1}^n (z_i'\gamma)^{2} \right)^{1/2}$. They develop bias-aware \acp{CI} under this penalty based on a likelihood ratio statistic, which are numerically shown to be close to optimal under homoskedasticity and a particular weighted average length criterion. However, unlike our \ac{CI}, the li_linear_2020 \ac{CI} may end up being longer than the long regression \ac{CI}, as we illustrate in our empirical application in (ref).

The next section presents our finite-sample results in the idealized model with Gaussian errors. (ref) discusses implementation in the more realistic setting with unknown error distribution. (ref) derives rates of convergence under high-dimensional asymptotics and bounds on an $\ell_p$ norm. (ref) compares our approach to \acp{CI} motivated by sparsity constraints. The performance of our methods is evaluated in a Monte Carlo study in (ref), while (ref) illustrates them in an empirical application. Proofs and auxiliary results appear in appendices.

Finite-sample results

This section sets up an idealized version of our model with Gaussian homoskedastic errors. We then show how to construct estimators and \acp{CI} in this model that are near-optimal in finite samples.

Setup

We write the model in (ref) in vector form as

where $w=(w_1,\dotsc, w_n)'\in\mathbb{R}^n$ is the variable of interest with coefficient $\beta\in\mathbb{R}$ and $Z=(z_{1}, \dotsc, z_{n})'\in \mathbb{R}^{n\times k}$ is a matrix of control variables. The design matrix $X=(w, Z)$ is treated as fixed. To obtain finite-sample results, we further assume that the errors are normal and homoskedastic $\varepsilon\sim\mathcal{N}(0, \sigma^{2}I_{n})$, with $\sigma^{2}$ known. To ensure informative inference on $\beta$ when $k$ is large relative to $n$ (including the case $k>n$), the researcher needs to make a priori restrictions on the control coefficients $\gamma$. We assume that these restrictions can be formalized by restricting the parameter space for $(\beta, \gamma')'$ to be $\mathbb{R}\times \Gamma$ where, for some linear subspace $\mathcal{G}$ of $\mathbb{R}^{k}$ and some seminorm $\operatorname{Pen}(\cdot)$ on $\mathcal{G}$,

The requirement that $\operatorname{Pen}(\cdot)$ be a seminorm means that it satisfies the triangle inequality ($\operatorname{Pen}(\gamma+\tilde\gamma)\le \operatorname{Pen}(\gamma)+\operatorname{Pen}(\tilde\gamma)$), and homogeneity ($\operatorname{Pen}(c\gamma)=\abs{c}\operatorname{Pen}(\gamma)$ for any scalar $c$), but, unlike a norm, it is not necessarily positive definite ($\operatorname{Pen}(\gamma)=0$ does not imply $\gamma=0$). This allows us to cover settings where only a subset of the control coefficients is restricted.

A common class of restrictions arises when $\operatorname{Pen}(\gamma)$ is a weighted $\ell_{p}$ norm on a subset of the coefficients. To describe two examples in this class of restrictions, partition the controls into a set of $k_{1}\geq 0$ unrestricted baseline controls and a set of $k_{2}=k-k_{1}$ additional controls, $Z=(Z_{1}, Z_{2})$. Partition $\gamma=(\gamma_{1}', \gamma_{2}')'$ accordingly. Let $H_{A}$ denote the projection matrix onto the column space of a matrix $A$. Let $\|\cdot\|_p$ denote the $\ell_p$ norm.

In addition to selecting the penalty, the specification of $\Gamma$ also requires the researcher to pick the regularity parameter $C$; here we take it as given, and defer a discussion of its choice to (ref).

Formulating the parameter space $\Gamma$ in terms of a seminorm is not restrictive in the sense that essentially any convex set $\Gamma$ that is symmetric ($\gamma\in\Gamma$ implies $-\gamma\in\Gamma$) can be defined in this way yosida_functional_1995. Although we rule out non-convex constraints on $\Gamma$, such as sparsity, our results nonetheless have implications for such settings, as we discuss in (ref).

Our goal is to construct estimators and \acp{CI} for $\beta$. To evaluate estimators $\hat{\beta}$ of $\beta$, we consider their worst-case performance over the parameter space $\mathbb{R}\times \Gamma$ under the \ac{MSE} criterion,

where $E_{\beta, \gamma}$ denotes expectation under $(\beta, \gamma')'$. An interval $\{\hat\beta\pm \hat\chi\}$ with half-length $\hat\chi=\hat\chi(Y, X)$ is a \ac{CI} with level $1-\alpha$ if it satisfies the coverage requirement

where $P_{\beta, \gamma}$ denotes probability under $(\beta, \gamma')'$. To compare two CIs under a particular parameter vector $(\beta, \gamma')'$, we prefer the one with shorter expected length $E_{\beta, \gamma}[2\hat\chi]$. Note that optimizing expected length will not necessarily lead to \acp{CI} centered at an estimator $\hat\beta$ that is optimal under the \ac{MSE} criterion.

Linear estimators and CIs

We start by considering estimators that are linear in the outcomes $Y$, $\hat{\beta}=a'Y$, and derive \acp{CI} based on such estimators. The $n$-vector of weights $a$ may depend on the design matrix $X$ or the known variance $\sigma^{2}$. In (ref) below, we show how to choose the weights $a$ optimally, and in (ref) we show that when $a$ is optimally chosen, the resulting estimators and \acp{CI} are optimal or near-optimal among all procedures, not just linear ones.

Under a given parameter vector $(\beta, \gamma')'$, the bias of $\hat{\beta}=a'Y$ is given by $a'(w\beta+Z\gamma) - \beta$. As $(\beta, \gamma')'$ ranges over the parameter space $\mathbb{R}\times \Gamma$, the bias ranges over the interval $[-\operatorname{\overline{bias}}_{\Gamma}(\hat\beta), \operatorname{\overline{bias}}_{\Gamma}(\hat\beta)]$, where

denotes the worst-case bias. The variance of $\hat\beta$ does not depend on $(\beta, \gamma')'$, and is given by $\operatorname{var}(\hat\beta) = \sigma^2 a'a$.

To form a CI centered at $\hat\beta$, note that the $z$-statistic $(\hat\beta-\beta)/\operatorname{var}(\hat\beta)^{1/2}$ follows a $\mathcal{N}(b,1)$ distribution with mean bounded by $\abs{b}\le \operatorname{\overline{bias}}_{\Gamma}(\hat\beta)/\operatorname{var}(\hat\beta)^{1/2}$. Thus, a two-sided \ac{CI} can be formed as

and $\operatorname{cv}_{\alpha}(B)$ denotes the $1-\alpha$ quantile of the folded normal distribution, $\abs{\mathcal{N}(B,1)}$.\footnote{The critical value $\operatorname{cv}_{1-\alpha}(B)$ can be computed as the square root of the $1-\alpha$ quantile of a non-central $\chi^{2}$ distribution with $1$ degree of freedom and non-centrality parameter $B^{2}$.}

This \ac{CI} is bias-aware in that the critical value $\operatorname{cv}_{\alpha}(\cdot)$ reflects the potential finite-sample bias of $\hat{\beta}$. Following the terminology in donoho94, we refer to the \ac{CI} as a \acf{FLCI}, since its length $2\chi$ is fixed: it depends only on the non-random design matrix $X$, and known variance $\sigma^{2}$, but not on $Y$ or the parameter vector $(\beta, \gamma')'$.

Optimal weights

Both the \ac{MSE} $R(\hat\beta;\Gamma)= \operatorname{\overline{bias}}_{\Gamma}(\hat\beta)^2 + \operatorname{var}(\hat\beta)$ and the \ac{FLCI} length $2\chi$ given in (ref) are increasing in the variance of $\hat{\beta}$ and in its worst-case bias $\operatorname{\overline{bias}}_{\Gamma}(\hat{\beta})$. Therefore, to find the optimal weights, we first minimize variance subject to a bound $B$ on worst-case bias,

We then vary the bound $B$ to find the optimal bias-variance tradeoff for a given criterion (\ac{MSE} or \ac{FLCI} length). Since this optimization does not depend on the outcome data $Y$, optimizing the weights does not affect the coverage properties of the resulting \ac{CI}.

Our main computational result, in (ref) below, shows that the estimator solving the optimization problem in (ref) is given by a simple two-step procedure. In the first step, we estimate a penalized regression of $w$ on $Z$ with penalty $\operatorname{Pen}(\pi)$, so that the coefficient estimate on $Z$, $\pi_\lambda$, solves

where $t_\lambda$ is a bound on the penalty term. We refer to (ref) as a (regularized) propensity score regression, even though we don't require $w_{i}$ to be binary. In the second step, we use the residuals $\tilde{w}_{\lambda}:=w-Z\pi_{\lambda}$ from the propensity score regression as instruments in a univariate regression of $Y$ on $w$. The tuning parameter $\lambda$ indexes the weight placed on the constraint in (ref), and its selection depends on the criterion we are optimizing. It may correspond to the Lagrange multiplier in a Lagrangian formulation of (ref), or, if we can solve (ref) directly, we may take $t_\lambda=\lambda$.

The result follows by applying the general theory of IbKh85, donoho94, low95, and ArKo18optimal to our setting, which allows us to rewrite (ref) as a convex optimization problem. Solving it then yields the result.

With the solution to (ref) in hand, for estimation and \ac{CI} construction, we select penalties $\lambda^{*}_{\textnormal{MSE}}$ and $\lambda^{*}_{\textnormal{FLCI}}$ that optimize the \ac{MSE} and \ac{CI} length, respectively. Specifically, the penalties solve the univariate optimization problems

with $V_\lambda$ and $\overline B_\lambda$ given in (ref). The optimal linear estimator is then given by $\hat{\beta}_{\lambda^{*}_{\textnormal{MSE}}}$, and the optimal \ac{FLCI} takes the form $\hat{\beta}_{\lambda^{*}_{\textnormal{FLCI}}} \pm \sigma\norm{{a_{\lambda^{*}_{\textnormal{FLCI}}}}}_{2} \cdot \operatorname{cv}_{\alpha}\left( \frac{C}{\sigma}\frac{{a_{\lambda^{*}_{\textnormal{FLCI}}}} 'Z\pi_{\lambda^{*}_{\textnormal{FLCI}}} }{ \operatorname{Pen}(\pi_{\lambda^{*}_{\textnormal{FLCI}}})\norm{{a_{\lambda^{*}_{\textnormal{FLCI}}}}}_{2}} \right)$.

As $t_{\lambda}\to 0$, provided that $\operatorname{Pen}(\cdot)$ is a norm on $Z_{2}$, $\hat{\beta}_{\lambda}$ converges to the short regression estimate $\hat{\beta}_{\textnormal{short}}=\frac{w'(I-H_{Z_{1}})Y}{w'(I-H_{Z_{1}})w}$ that only includes the unrestricted controls $Z_{1}$. This estimator minimizes variance among all linear estimators with finite worst-case bias. In the other direction, as $t_{\lambda}\to\infty$, $\hat{\beta}_{\lambda}$ converges to the long regression estimate $\hat{\beta}_{\textnormal{long}}=\frac{w'(I-H_{Z})Y}{w'(I-H_{Z})w}$, provided that $w$ is not in the column space of $Z$ (which ensures that the condition $\norm{\tilde{w}_{\lambda}}_{2}>0$ in (ref) holds for all $\lambda$). This estimator minimizes variance among all linear estimators that are unbiased, so (ref) reduces to the Gauss-Markov theorem in this case. In other words, the short and long regressions are corner solutions of the bias-variance tradeoff, in which weight is entirely placed on variance, or on bias.

Efficiency among non-linear procedures

So far, we have restricted attention to procedures that are linear in the outcomes $Y$. We now show that the estimator $\hat{\beta}_{\lambda^{*}_{\textnormal{MSE}}}$, and the \ac{CI} based on the estimator $\hat{\beta}_{\lambda^{*}_{\textnormal{FLCI}}}$ are in fact highly efficient among all procedures, not just linear ones. This is due to the convexity and symmetry of the parameter space $\Gamma$, and follows from the general results in donoho94, low95 and ArKo18optimal for estimation of linear functionals in Gaussian models with convex parameter spaces.

By construction, the estimator $\hat{\beta}_{\lambda}$ minimizes variance among all linear estimators with a bound $C\overline{B}_{\lambda}$ on the bias (or equivalently, it minimizes bias among all linear estimators with a bound $V_{\lambda}$ on the variance). (ref)(ref) shows that this optimality property is retained if we enlarge the class of estimators to all estimators, including non-linear ones. As a result, the minimax linear estimator $\hat{\beta}_{\lambda^{*}_{\textnormal{MSE}}}$ (i.e.\ the estimator attaining the lowest worst-case \ac{MSE} in the class of linear estimators) continues to perform well among all estimators, including non-linear ones: by (ref)(ref), its worst-case \ac{MSE} efficiency is at least 80%. The exact efficiency bound $\kappa^*_{\textnormal{MSE}}(X, \sigma, \Gamma)$ depends on the design matrix, noise level, and particular choice of the parameter space, and can be computed explicitly in particular applications. We have found that typically the efficiency is considerably higher than the 80% lower bound.

Finally, (ref)(ref) shows that it is not possible to substantively improve upon the \ac{FLCI} based on $\hat{\beta}_{\lambda^{*}_{\textnormal{FLCI}}}$ in terms of expected length when $\gamma=0$, even if we consider variable length \acp{CI} that “direct power” at $\gamma=0$ (potentially at the expense of longer expected length when $\gamma\neq 0$). The construction of the \ac{FLCI} may appear conservative: its length depends on the worst-case bias over the parameter space for $(\beta, \gamma')'$, which, as the proof of (ref) shows, attains at $\gamma=Ct_{\lambda^{*}_{\textnormal{FLCI}}}^{-1}\pi_{\lambda^{*}_{\textnormal{FLCI}}}$, with $\operatorname{Pen}(\gamma)=C$. Therefore, one may be concerned that when the magnitude of $\gamma$ is much smaller than $C$, the \ac{FLCI} is too long. (ref)(ref) shows that this is not the case, and the efficiency of the \ac{FLCI} is at least 71.7% relative to variable-length \acp{CI} that optimize their expected length when $\gamma=0$. The exact efficiency bound $\kappa^*_{\textnormal{MSE}}(X, \sigma, \Gamma)$ can be computed explicitly in particular applications, and we have found that it is typically considerably higher than $71.7\%$.

A consequence of (ref)(ref) is that it is impossible to form a CI that is adaptive with respect to the regularity parameter $C$ that bounds $\operatorname{Pen}(\gamma)$. In the present setting, an adaptive CI would have length that automatically reflects the true regularity $\operatorname{Pen}(\gamma)$ while maintaining coverage under a conservative a priori bound on $\operatorname{Pen}(\gamma)$. However, according to (ref)(ref), any CI must have expected width that reflects the conservative a priori bound $C$ rather than the true regularity $\operatorname{Pen}(\gamma)$, even when $\operatorname{Pen}(\gamma)$ is much smaller than the conservative a priori bound $C$. In particular, it is impossible to automate the choice of the regularity parameter $C$ when forming a CI\@. We therefore recommend varying $C$ as a form of sensitivity analysis, or using auxiliary information to choose $C$; see (ref).

Heterogeneous treatment effects

If $w_{i}$ is as good as randomly assigned conditional on $z_{i}$, the coefficient $\beta$ in (ref) can be interpreted as the \ac{ATE} of a one-unit increase in the variable $w$. This interpretation requires that the individual \acp{TE} are mean independent of $z_{i}$. To relax this assumption, we replace $\beta$ in (ref) with a covariate-specific coefficient $\beta(z)$ that represents the conditional \ac{ATE} for units $i$ with $z_{i}=z$, obtaining the model

Suppose the parameter of interest takes the form $\int \beta(z)\, d\mu(w, z)$ where $\mu$ is a signed measure defined on some set that includes the empirical support $\{(w_{i}, z_{i})\}_{i=1}^{n}$. Allowing $\mu$ to be signed allows for inference on non-convex averages of $\beta(z)$. We allow $\mu$ to place nonzero mass outside the empirical support $\{(w_{i}, z_{i})\}_{i=1}^{n}$, thereby allowing for extrapolation. The general theory of estimation and inference on linear functionals developed in donoho94 and ArKo18optimal underlying (ref) can be applied to inference on this parameter under any convex and symmetric restriction on the function $\beta(\cdot)$ and the parameter $\gamma$---only the worst-case bias calculation in (ref) changes. We now discuss two particular specifications for $\mu$ and $\beta(\cdot)$.

For the first approach, let $\mu$ correspond to a weighted empirical measure with weights $c_{i}$ that sum to one, so the parameter of interest is given by $\tilde{\beta}=\sum_{i}c_{i}\beta(z_{i})$. For example, the unweighted case $c_{i}=1/n$ gives the (conditional on the sample) \ac{ATE}, while setting $c_{i}=w_{i}/\sum_{j}w_{j}$ gives the \ac{ATE} for the treated. Assume also that the function $\beta(z)$ is linear, $\beta(z)=z'\delta$, and that the first element of $z_{i}$ is a constant. Consider a parameter space for $(\delta',\gamma')'$ given by $\{(\delta',\gamma')'\in\mathcal{B}\colon \operatorname{Pen}(\delta,\gamma)\leq C\}$, where $\mathcal{B}$ is a subspace of a Euclidean space with a seminorm $\operatorname{Pen}$. Allowing $\mathcal{B}$ to be a subspace allows us to restrict $\beta(z)$ to only depend on a subset of the controls. As noted by imbens_recent_2009 in the unpenalized case, we can map this problem back into the model in (ref) by rewriting (ref) under these assumptions as

where $z_{i,-1}$ is the vector of controls excluding a constant and $\delta_{-1}$ is the corresponding subvector of $\delta$. This is exactly our problem in (ref), with the control vector consisting of the original controls $z_{i}$ as well as the interaction of the treatment with the demeaned controls, $w_{i}(z_{i,-1}-\sum_{i}c_{i}z_{i,-1})$. We use this approach in our empirical application in (ref).

For the second approach, we compute the same linear estimator and bias-aware \acp{CI} as in the homogeneous \ac{TE} model in (ref); we only change the interpretation of the estimand as targeting a particular weighted average of \acp{TE} given in the next theorem. To set the stage for the theorem, let us denote the worst-case bias of a linear estimator $\hat{\beta}$ relative to an estimand $\int \beta(z)d\mu(w, z)$ when the heterogeneity is completely unrestricted by $\widetilde{\operatorname{\overline{bias}}}_{\Gamma}(\hat{\beta};\mu)=\sup_{\beta(\cdot), \gamma} \left[ \sum_{i=1}^n a_{i}w_i \beta(z_i)+a'Z\gamma - \int \beta(z)\, d\mu(w, z) \right]$, where the supremum is over $\gamma\in\Gamma$ and all functions $\beta(\cdot)$. For an $n$-vector $a$, let $\mu_{a, w}^{*}(w, z)$ denote a weighted empirical measure with (possibly negative) weights $a_{i}w_{i}$, so that the parameter of interest becomes $\tilde{\beta}_{a, w}=\sum_{i}a_{i}w_{i}\beta(z_{i})$.

The theorem gives three results for inference when $\beta(\cdot)$ is unrestricted. First, given a linear estimator $a'Y$, the only estimand for which the bias is finite is $\tilde{\beta}_{a, w}$. Second, for this estimand, the bias is the same as in the homogeneous \ac{TE} model with $\beta(z)=\beta$. Thus, assuming homogeneous \acp{TE} leads to valid inference for $\tilde{\beta}_{a, w}$ when \acp{TE} are in fact heterogeneous. In particular, the bias-aware \ac{CI} based on the estimator $\hat{\beta}_{\lambda}$ provides inference on the weighted average $\tilde{\beta}_{a_{\lambda}, w}= \sum_{i}\tilde{w}_{\lambda, i}w_{i}\beta(z_{i})/\sum_{i}\tilde{w}_{\lambda, i}w_{i}$. Observe that, when the treatment $w_{i}$ is binary, the weights are non-negative if and only if the residual in the propensity score regression $\tilde{w}_{\lambda, i}$ is positive whenever $w_{i}=1$---this can be easily verified in a given application. Equivalently, the weights are positive if the fitted values $z_{i}'\pi_{\lambda}$ are smaller than one: the fitted values in the propensity score regression must respect the population constraint that treatment probabilities must be smaller than one. This finding is a finite-sample analog of the identification result in ghk21, who show that the estimand in the partly linear model has an analogous weighted \ac{ATE} interpretation under heterogeneous \acp{TE}. An analogous identification result in a random design setting dates to at least angrist_estimating_1998, who gave a weighted \ac{ATE} interpretation to the \ac{OLS} estimand.

Third, the estimator $\hat{\beta}_{\lambda}$ remains optimal in the heterogeneous \ac{TE} model in that it solves a bias-variance tradeoff in a problem where we can pick the estimand to make the estimation problem as easy as possible. The problem (ref) is a finite-sample version of the “moving the goalposts” problem considered by crump_moving_2006.\footnote{The optimization problem (ref) was also used by ImWa19 in the context of a regression discontinuity design with multiple cutoffs, although they did not explicitly note the optimality properties of the resulting estimator.} crump_moving_2006 derived the measure $\mu$ that minimizes the asymptotic variance of a particular class of inverse propensity score weighted estimators of weighted \acp{ATE} in a random design setup. ghk21 show this measure in fact minimizes the variance among all regular estimators; the measure coincides with the weighting of the treatment effects given in angrist_estimating_1998 under homoskedasticity, giving an optimality property to the \ac{OLS} estimator.

Implementation with non-Gaussian and heteroske\-dastic errors

We now discuss practical implementation issues, allowing $\varepsilon$ to be non-Gaussian and heteroskedastic. As a baseline, we propose the following implementation:

The following remarks discuss the implementation choices, and the optimality and validity of the baseline procedure.

Rates of convergence

We now derive the rates of convergence for the optimal linear \acp{FLCI} as $n\to\infty$. For ease of notation, we assume all coefficients are constrained, and focus on the case $\operatorname{Pen}(\gamma)=\norm{\gamma}_{p}$ for some $p\ge 1$, and the case $\operatorname{Pen}(\gamma)=\norm{Z\gamma/\sqrt{n}}_{2}$ (see (ref)). We allow the regularity parameter $C=C_{n}$ go to $0$ or $\infty$ with the sample size, and consider high dimensional asymptotics where $k=k_n\gg n$. We consider a standard “high dimensional” setting, placing conditions on the design matrix $X$ that hold with high probability when $w_i, z_i$ are drawn i.i.d.\ over $i$, with the eigenvalues of $\operatorname{var}((w_i, z_i')')$ bounded away from zero and infinity.

Let $q\in[0,\infty]$ denote the Hölder conjugate of $p$, satisfying $1/p+1/q=1$. We will show that when $\operatorname{Pen}(\gamma)=\norm{\gamma}_{p}$, the optimal linear \ac{FLCI} shrinks at the rate

Furthermore, for $p=1$ and $p=2$, we will show that no other \ac{CI} can shrink at a faster rate. For $p=1$, we will in fact prove a stronger result showing that imposing sparsity bounds on the outcome and propensity score regressions, in addition to the bound on $\operatorname{Pen}(\gamma)$, does not help achieve a faster rate, unless one assumes sparsity of order greater than $C_n\sqrt{n/\log(k)}$ (termed the “ultra sparse” case in cai_lasso_2017). For the case $\operatorname{Pen}(\gamma)=\norm{Z\gamma/\sqrt{n}}_{2}$, we will show that the optimal rate is given by $n^{-1/2}+C$ when $k\gg n$.

If $k\gg n$ and $C=C_{n}$ does not decrease to zero with $n$, these rates require $p<2$ (so that $q>2$) for consistency. When $p=1$, we can then allow $k$ to grow exponentially with $n$, whereas setting $1<p<2$ allows for $k$ to grow at a polynomial rate in $n$ that depends on $p$. Since taking $C_n\to 0$ rules out even a single coefficient being bounded away from zero, these bounds imply that taking $p<2$ in “high dimensional” settings is necessary for consistency, with $p=1$ offering the best rate conditions. It also follows from these rate results that if $C_{n}=C$ does not decrease to zero with $n$, the bias term can dominate asymptotically, making it necessary to explicitly account for bias in \ac{CI} construction even in large samples.

Upper bounds

To state the result, given $\eta>0$, let $\mathcal{E}_n(\eta)$ denote the set of design matrices $X$ for which there exists $\delta\in\mathbb{R}^k$ such that

Let $R^{*}_{\textnormal{FLCI}}(X, C)=2\operatorname{cv}_\alpha(C \overline{B}_{\lambda^{*}_{\textnormal{FLCI}}}/V_{\lambda^{*}_\textnormal{FLCI}}^{1/2})\cdot V_{\lambda^{*}_\textnormal{FLCI}}^{1/2}$ denote the length of the optimal linear \ac{FLCI}.

The second part of the \namecref{thm:upper_rate_bound} follows since the short regression without any controls achieves a bias that is of the order $C$. The first part shows that the upper bounds on the rate of convergence match those in (ref) if the high-level condition $X\in \mathcal{E}_n(\eta)$ holds. The next lemma shows that this high-level condition holds with high probability when $w_{i}, Z_{i}$ are drawn i.i.d.\ from a distribution satisfying mild conditions on moments and covariances.

Lower bounds

We now show that the rates in (ref) are sharp when $p=2$, or $p=1$.

\texorpdfstring{$p=2$}{p=2}

As with the upper bound in (ref), we derive a bound that holds when the design matrix $X$ is in some set, and then show that this set has high probability when $w_i, z_i$ are drawn i.i.d.\ from a sequence of distributions satisfying certain conditions. We focus on the case $k\geq n$. Let $\widetilde{\mathcal{E}}_n(\eta)$ denote the set of design matrices $X$ such that

where $\operatorname{eig}(A)$ denotes the set of eigenvalues of a square matrix $A$.

If $z_{i}$ is i.i.d.\ over $i$, then $EZZ'/k $ is equal to the $n\times n$ identity matrix times the scalar $\frac{1}{k}\sum_{j=1}^k E[z_{ij}^2]$. Thus, the condition on the minimum eigenvalue of $ZZ'/k$ will hold under concentration conditions on the matrix $Z'Z$ so long as the second moments of the covariates are bounded from below. Here, we state a result for a special case where the $z_{ij}$'s are i.i.d.\ normal, which is immediate from donoho_for_2006.

\texorpdfstring{$p=1$}{p=1}

We now consider the case where $p=1$, as in (ref). Rather than imposing conditions on $X$ in a fixed design setting that hold with high probability (as in (ref) and (ref)), we directly consider a random design setting, and we do not condition on $X$ when requiring coverage of \acp{CI}. This allows us to strengthen the conclusion of our theorem by showing that the rate in (ref) is sharp even if one imposes a linear model for $w_i$ given $z_i$ along with sparsity and $\ell_1$ bounds on the coefficients in this model.

We introduce some additional notation to cover the random design setting, which we use only in this section. We consider a random design model

We use $P_{\vartheta}$ and $E_{\vartheta}$ for probability and expectation when $Y, X$ follow this model with parameters $\vartheta=(\beta, \gamma', \delta', \sigma^2,\sigma^2_{v})'$. Let $\sigma_0^2>0$ and $\sigma^2_{v, 0}>0$ be given and let $\Theta(C, s, \eta)$ denote the set of parameters $\vartheta=(\beta, \gamma', \delta', \sigma^2, \sigma^2_{v})$ where $\abs{\sigma^2-\sigma^2_0}\le \eta$, $\abs{\sigma_{v}^2-\sigma^2_{v, 0}}\le \eta$, $\norm{\gamma}_{1}\le C$, $\norm{\delta}_{1}\le C$, $\norm{\gamma}_0\le s$ and $\norm{\delta}_0\le s$.

(ref) follows from similar arguments to cai_lasso_2017 and javanmard_debiasing_2018, who provide similar bounds for the case where only a sparsity bound is imposed. According to (ref), imposing sparsity does not allow one to improve upon the \acp{CI} that uses only the $\ell_1$ bound $\norm{\gamma}_1\le C_n$ (thereby attaining the rate in (ref)), unless one imposes sparsity of order greater than $C_n\sqrt{n/\log k}$. We provide further comparison with \acp{CI} that impose sparsity in the next section.

Comparison with sparsity constraints

Several authors have considered \acp{CI} for $\beta$ using “double lasso” estimators belloni_inference_2014,javanmard_confidence_2014,van_de_geer_asymptotically_2014,zhang_confidence_2014. These \acp{CI} are valid under the parameter space

where $\norm{\gamma}_0=\#\{j\colon \gamma_j\ne 0\}$ is the $\ell_0$ “norm,” which indexes the sparsity of $\gamma$, and with $s$ increasing slowly enough relative to $n$ and $k$. Since $\norm{\gamma}_0$ is not a true norm or seminorm (it is non-convex), this parameter space is not covered by our setup. Nonetheless, as we show in (ref), if the sparsity assumption is used to bound the $\ell_{1}$ loss of a preliminary lasso estimator, arguments from (ref) lead to estimators and \acp{CI} that are analogous to those proposed in the double lasso literature. In (ref), we provide a comparison of our approach to these double lasso \acp{CI}.

Connection between double lasso and optimal estimator under \texorpdfstring{$\ell_1$}{ell1} constraints

When $\operatorname{Pen}(\gamma)=\norm{\gamma}_1$ ((ref)), the solution $\pi_\lambda$ to (ref) is the lasso estimate in the propensity score regression of $w$ on $Z$, and our estimator (ref) uses residuals from this lasso regression. This is related to “double lasso” estimators used to form \acp{CI} for $\beta$ under sparsity constraints on $\gamma$ belloni_inference_2014,javanmard_confidence_2014,van_de_geer_asymptotically_2014,zhang_confidence_2014. For concreteness, we focus on the estimator in zhang_confidence_2014, which is given by

where $\hat\beta_{\textnormal{lasso}}, \hat\gamma_{\textnormal{lasso}}$ are the lasso estimates from regressing $Y$ on $X$:

for some penalty parameter $\tilde \lambda>0$.

To further understand the connection between these estimators, we note that zhang_confidence_2014 motivate their approach by bounds of the form

which hold with high probability with the constant depending on certain “compatibility constants” that describe the regularity of the design matrix $X$ buhlmann_statistics_2011. This suggests correcting the initial estimate $\hat\beta_{\textnormal{lasso}}$ by estimating $\tilde\beta=\beta-\hat\beta_{\textnormal{lasso}}$ in the regression

where $\tilde Y=Y-\hat\beta_{\textnormal{lasso}}-Z\hat\gamma_{\textnormal{lasso}}$. Heuristically, we can treat the bound in (ref) as a constraint $\norm{\tilde\gamma}_{1}\le \tilde C$ on the unknown parameter $\tilde\gamma=\gamma-\hat\gamma_{\textnormal{lasso}}$ and search for an optimal estimator of $\tilde\beta=\beta-\hat\beta_{\textnormal{lasso}}$ under this constraint. Applying the optimal estimator derived in (ref) then suggests estimating $\beta-\hat\beta_{\textnormal{lasso}}$ with $\frac{\tilde{w}_{\lambda}'\tilde Y}{\tilde{w}_{\lambda}'w}$. Adding this estimate to $\hat\beta_{\textnormal{lasso}}$ gives the estimate $\hat\beta_{\textnormal{ZZ}}$ proposed by zhang_confidence_2014. Whereas zhang_confidence_2014 motivate their approach as one possible way of correcting the initial estimate $\hat\beta_{\textnormal{lasso}}$ using the bound in (ref), the above analysis shows that their correction is in fact identical to an approach in which one optimizes this correction numerically.\footnote{The estimator proposed by javanmard_confidence_2014 performs a numerical optimization of this form, but with the constraint (ref) replaced by a constraint on $\abs{\hat\beta_{\textnormal{lasso}}-\beta}+\norm{\hat\gamma_{\textnormal{lasso}}-\gamma}_{1}$. Thus, (ref) shows that a modification of the constraint used in javanmard_confidence_2014 yields the same estimator as zhang_confidence_2014.}

Under the bound in (ref) it follows that $\hat\beta_{\textnormal{ZZ}}-\beta =\tilde b + {a_\lambda}'\varepsilon$ where $a_\lambda=\frac{\tilde{w}_\lambda}{\tilde{w}_{\lambda}'w}$ are the optimal weights under the $\ell_1$ constraint $\norm{\tilde\gamma}_{1}\le \tilde C$, given in (ref). Furthermore, $\abs{\tilde b}\le \tilde C\overline B_\lambda$, with $\overline B_\lambda$ given in (ref) and $\tilde C$ given in (ref), and the variance of the random term ${a_\lambda}'\varepsilon$ is given by $V_\lambda$ in (ref). Using arguments similar to those used to prove (ref), it follows that $\tilde C\overline B_\lambda/\sqrt{V_\lambda}$ is bounded by a constant times $s(\log k)/\sqrt{n}$, so that one can ignore bias in large samples as long as this term converges to zero. This leads to the \ac{CI} proposed by zhang_confidence_2014, which takes the form

where $\hat V_\lambda$ is an estimate of the variance $V_\lambda$. We use the term “double lasso \ac{CI}” to refer to this \ac{CI}, and to related \acp{CI} such as those proposed in belloni_inference_2014,javanmard_confidence_2014,van_de_geer_asymptotically_2014.

Comparison of our approach with CIs based on double lasso estimators

When should one use a double lasso \ac{CI}, and when should one use the approach in the present paper? In principle, this depends on the a priori assumptions one is willing to make, and whether they are best captured by a sparsity bound or a bound on convex penalty function, such as the $\ell_{1}$ or $\ell_{2}$ norm. In many settings, it may be difficult to motivate the assumption that a regression function has a sparse approximation, whereas upper bounds on the magnitude of the coefficients may be more plausible.

A key advantage of the \acp{CI} and estimators we propose is that they have sharp finite-sample optimality properties and coverage guarantees in the fixed design Gaussian model with known error variance. While this is an idealized setting, the worst-case bias calculations do not depend on the error distribution, and remain the same under non-Gaussian, heteroskedastic errors. Our approach directly accounts for the potential finite-sample bias of the estimator, rather than relying on “asymptotic promises” about rates at which certain constants involved in bias terms converge to zero.

On the flip side, our \acp{CI} require an explicit choice of the regularity parameter $C$ in order to form a “bias-aware” \ac{CI}. In contrast, \acp{CI} based on double lasso estimators do not require explicitly choosing the regularity (in this case, the sparsity $s$), since they ignore bias. This is justified under asymptotics in which $s$ increases more slowly than $\sqrt{n}/\log k$, which lead to the bias of $\hat\beta_{\textnormal{ZZ}}$ decreasing more quickly than its standard deviation. Thus, the \ac{CI} in (ref) is “asymptotically valid” without the need to explicitly specify the sparsity index $s$: one need only make an “asymptotic promise” that $s$ increases slowly enough. However, such asymptotic promises are difficult to evaluate in a given finite-sample setting. Indeed, as shown by wuthrich_omitted_2021 and confirmed in our Monte Carlos in (ref) below, the double lasso CI leads to undercoverage in finite samples even in relatively sparse settings. To ensure good finite-sample coverage of the \ac{CI} in (ref), one needs to ensure that the actual finite-sample bias is negligible relative to the standard deviation of the estimator. But since any bias bound depends on the sparsity index $s$ (as in the bound in (ref)), this gets us back to having to explicitly specify $s$.

Thus, \acp{CI} that ignore bias such as conventional \acp{CI} based on double lasso estimators do not avoid the problem of specifying $s$ or $C$: they merely make such choices implicit in their asymptotic promises. These issues show up formally in the asymptotic analysis of such \acp{CI}. In particular, double lasso \acp{CI} require the “ultra sparse” asymptotic regime $s=o(\sqrt{n}/\log k)$, and they undercover asymptotically in the “moderately sparse” regime where $s$ increases more slowly than $n$ with $s\gg \sqrt{n}/\log k$. Indeed, (ref) above, as well as the results of cai_lasso_2017 and javanmard_debiasing_2018 show that it is impossible to avoid explicitly specifying $s$ if one allows for the moderately sparse regime.

On the other end of the spectrum, in the “low dimensional” regime where $k\ll n$, the double lasso \ac{CI} is asymptotically equivalent to the usual \ac{CI} based on the long regression. Thus, the double lasso \ac{CI} cannot be used when the goal is to use a priori information on $\gamma$ to improve upon the \ac{CI} based on the long regression muralidharan_factorial_2020. In contrast, our approach optimally incorporates the bound $C$ regardless of the asymptotic regime.

Simulation results

We now illustrate the performance of our methods when the penalty takes the form of an $\ell_{1}$ norm on a subset of $k_{2}$ controls, as in (ref). We consider a design taken from belloni_inference_2014, with data generated from a random regressor model that supplements (ref) with a propensity score regression

with $\tilde{w}_{i}$ and $\varepsilon_{i}$ independent standard Gaussian, and independent of $z_{i}$, which are distributed i.i.d.\ $\mathcal{N} (0, \Sigma)$ with $\Sigma_{ij}=2^{-\abs{i-j}}$. Similar to wuthrich_omitted_2021, we tweak the belloni_inference_2014 design by considering regression coefficients that are of similar magnitude rather than decaying. This allows us to separately vary the degree of sparsity and the signal-to-noise ratio. Specifically, we set

We consider three methods for constructing \acp{CI} for $\beta$ with nominal level 95%. The first two methods implement (ref), with the penalty given by the $\ell_{1}$ norm of $\gamma_{2}$, the last $k_{2}$ regression coefficients. The first method, which we refer to as “oracle,” sets the penalty parameter $C$ to the actual value of $\norm{\gamma_{2}}_{1}$, and uses knowledge of the variance of the error term $\varepsilon_{i}$. The second \ac{CI}, termed “AKK,” uses initial residual estimates based on the lasso estimator (that only penalizes $\gamma_{2}$), with the penalty chosen via 10-fold cross-validation. The \ac{CI} uses the rule of thumb calibration $C^{rot}=\norm{\widehat{\gamma}_{short}}_{1}$ from (ref), where $\widehat{\gamma}_{short}$ are \ac{OLS} estimates from a short regression that only includes the first $k_{1}$ controls (the “baseline” controls). The final method, termed “BCH”, implements the double lasso procedure by belloni_inference_2014, using the R package hdm, without penalizing the $k_{1}$ baseline controls and including an intercept.

The \ac{DGP} in our random regressor model depends on 8 parameters: $n$, $k_{1}$, $k_{2}$, $s$, $\beta$, $\sigma_{\tilde{w}}$, $c_{1}$ and $c_{2}$. We consider $n \in \{500, 1000\}$, $k_{1} \in \{5, 10\}$, $k_{2} \in \{100, 200, 500, 1000\}$, $s\in \{10, 20, 100\}$, $\beta \in \{0, 2\}$, and $\sigma_{\tilde{w}} \in \{0.5, 1\}$. We calibrate $c_{1}$ and $c_{2}$ by fixing the population $R^{2}$ from the regression of $Y$ on $Z$, and fixing the ratio $\nu_{rot} = \norm{\gamma_{2}}_{1} / \norm{\tilde{\gamma}_{short}}_{1}$. This allows us to directly control the signal-to-noise ratio, and the validity of the rule-of-thumb calibration: if $\nu_{rot}\leq 1$, then the population restriction underlying our rule of thumb is valid. We consider 4 values for the population $R^{2}$, $\{0.01, 0.1, 0.25, 0.5\}$, and 12 values for $\nu_{rot}$, $\{0.2, 0.4, \dotsc, 2.4 \}$. This gives a total of 9,216 \acp{DGP}.

(ref) reports the simulation results for $n=500$. The results for $n=1000$ are reported in (ref). In line with the theory, the coverage of the oracle \ac{CI} is close to nominal across all designs.\footnote{The slight undercoverage reported in the tables is due to Monte Carlo error: with 1000 simulation draws, the expected worst-case coverage over 160 DGPs is 93% if the true coverage for each DGP is 95%.} When $\nu_{rot}\leq 1$, coverage of the AKK \ac{CI} is likewise close to nominal. Under mild violations of the population constraint, the \acp{CI} display moderate undercoverage: when $\nu_{rot}\leq 1.5$, coverage remains over 86.6% across all designs, and over 90.7% when $k_{2} \leq 200$. Only when $\nu_{rot} > 1.5 $ and $k_{2} \geq 500$, the undercoverage becomes more severe. In contrast, the BCH method displays moderate undercoverage even in sparse designs with $s=10$, with coverage at about 85% when $k_{2}=1000$ and $n=500$. The undercoverage gets more severe, with coverage dipping below 60% once $s=20$, and the \acp{CI} almost entirely miss the true parameter in dense designs with $s=100$. These results illustrate the concern discussed in (ref) that asymptotic sparsity requirements may be difficult to evaluate in finite samples.

The favorable coverage of the AKK \acp{CI} relies heavily on using the bias-aware critical value. Unreported simulations show that the coverage of \acp{CI} constructed using the same estimators as the AKK \acp{CI} but with standard critical values (i.e., 1.96 for 95% coverage), rather than our bias-aware critical values, can be as low as 74.8% for \acp{DGP} with $\nu_{rot}\leq 1$.

The AKK \acp{CI} display a mild increase in average length relative to the oracle, with the length penalty ranging between 0 and 16%. The length penalty relative to the BCH method is also in this range for designs where both methods achieve good coverage. This is a bargain price to pay for the much more reliable and transparent coverage performance.

Empirical application

This section shows the performance of our methods using survey data on $n=496$ winners of major and minor prizes in the Massachusetts lottery in 1984--88 from imbens2001lottery to estimate the \ac{MPE} out of unearned income, a key structural parameter in labor and public economics. While unearned income is typically endogenous, imbens2001lottery argue that in this sample, observable individual characteristics proxy well enough for the frequency of lottery ticket purchases that the magnitude of winnings is as good as random. The lottery winnings are paid out over 20 years, so that in a regression of the average social security earnings in the 6 years after the lottery, $Y_{i}$, onto yearly lottery payments, $X_{i}$, and individual controls, the coefficient on $X_{i}$ may be interpreted as the \ac{MPE}.

We focus on a specification taken from li_linear_2020, who augment a baseline set of $k_{1}=7$ individual controls $Z_{1}$ consisting of the intercept, two continuous controls (years of education and age), and 4 binary controls (indicators for male, college, age over 55, and age over 65) with $k_{2}=25$ additional controls $Z_{2}$ that are constructed by taking demeaned cross-products of 4 the baseline binary controls and their interactions with $X_{i}$ and dropping collinear terms. Both $Z_{1}$ and $Z_{2}$ are standardized. Following the discussion in (ref), the coefficient on $X_{i}$ in this specification can be interpreted as the average \ac{MPE}, allowing for heterogeneity in the \ac{MPE} with respect to the binary controls. In contrast, the short regression estimand in a regression that only includes $Z_{1}$ is biased for the average \ac{MPE} in presence of such heterogeneity.

The \ac{MPE} estimate in the long regression equals $-0.049$, close to the short regression estimate $-0.052$ that only includes the baseline controls $Z_{1}$ and corresponds to the specification in Table 4, column II row 1 in imbens2001lottery. However, the long regression estimate is very noisy: the 95% confidence interval $(-0.115, 0.016)$ includes positive values for the average \ac{MPE} which economic theory rules out, and it is over 3 times longer than the short regression \ac{CI} $(-0.073, -0.032)$. To increase precision of inference, li_linear_2020 restrict the average squared mean effects $z_{2i}'\gamma_{2}$ using an $\ell_{2}$ penalty given in (ref). Calibrating $C$ to the rule of thumb value from (ref), $C^{rot}=7.2$, yields the \ac{CI} $(-0.116, 0.018)$ using the li_linear_2020 method, which is even longer than the long regression \ac{CI}.\footnote{li_linear_2020 show that their method is close to optimal in terms of weighted average length under a homoskedastic benchmark. This may no longer be the case under heteroskedasticity. Their \ac{CI} is variable length, and may be longer than the long regression \ac{CI} in some samples even in the homoskedastic case. In contrast, our construction guarantees length improvements over the long regression in all samples under homoskedasticity.} The \ac{CI} constructed using our method, $(-0.114, 0.015)$, improves slightly upon the long \ac{CI}, but it is still too wide to be informative.\footnote{To make the methods more comparable and not conflate the comparison with differences in standard error construction, variance estimates underlying \acp{CI} for all methods use residual estimates based on a lasso estimator that penalizes only $\gamma_{2}$, with penalty chosen by 10-fold cross-validation.} The li_linear_2020 penalty affords only marginal precision gains because the penalty limits the average influence of the additional regressors---but these regressors only marginally increase the regression $R^{2}$ in the long regression: the adjusted $R^{2}$ increases from 0.233 in the short regression to 0.236 in the long regression.