Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
74,457 characters · 16 sections · 57 citation commands
The Fragility of Sparsity
\pagenumbering{Alph} \thispagestyle{empty}
\pagenumbering{arabic}
The linear regression model is the most common workhorse for estimating causal effects. The key assumption allowing for a causal interpretation of the coefficient on the treatment or intervention of interest is that conditional on the set of controls included in the model, the remaining variation in the treatment is as good as random. Since the plausibility of this assumption increases with additional control variables, researchers commonly estimate regression models with many controls.
To fix ideas, consider estimation of the coefficient $\beta$ on treatment $D_{i}$ in the regression model
based on $n$ observations indexed by $i$. Traditional estimation methods like \ac{OLS} limit the dimension $p$ of the control vector $W_{i}$ to be smaller than $n$. This limitation has fueled the development of alternative estimation methods that allow $p$ to exceed $n$, provided one is willing to make the (approximate) sparsity assumption that the control function $W_{i}'\gamma$ can be well-approximated, in a mean-squared error sense, by a linear combination of a small number $s$ of the controls, so that, up to log terms, $s^{2}/n\to 0$. In particular, several papers JaMo14,vbrd14,ZhZh14 have developed “debiased lasso” estimators of $\beta$ that are able to purge bias present when (ref) is estimated directly using the lasso or related methods. bch14 propose a corresponding post-double-selection procedure that first estimates the propensity score regression
by the lasso, and also runs a lasso regression of the outcome on the controls $W_{i}$. In a second step, $\beta$ is estimated by regressing the outcome on $D_{i}$ and a union of the controls selected by the two lasso regressions. A number of authors have developed variants of these methods that allow for treatment effect heterogeneity aiw18,bcfh17,ccddhnr18,farrell15. These developments have had a large influence on empirical practice in applied economics and related fields, and researchers routinely use these \acp{SBE} even in the regime $p<n$ as a complement to or replacement of \ac{OLS}.
This paper argues that this practice may not lead to reliable inference for two reasons. First, \acp{SBE}, like the sparsity assumption underlying them, are not invariant to linear reparametrizations of the control matrix $W$. For instance, while normalization choices such as the choice of baseline category for a set of dummy variables, or how and whether to center the variables prior to taking powers or interactions, are immaterial when running \ac{OLS}, they in general affect \acp{SBE}. We document that in three empirical applications of \acp{SBE}, which include the empirical illustration in bch14, seemingly innocuous normalization choices of this type induce variation in \acp{SBE} of the same magnitude as their standard error. It is thus impossible to interpret the resulting inference without substantively engaging with the question whether the particular normalization choice used is appropriate. Yet, papers using \acp{SBE} do not even tell readers what this choice was, let alone argue that it is the right one to obtain a sparse representation. While it is well-known that \acp{SBE} are not invariant to linear reparametrizations in the theoretical literature, this practice implies that the extent of their fragility is not currently appreciated, even by researchers in the field.
Second, one might be skeptical of the plausibility of the sparsity assumption more generally. In social science applications, it is often difficult to point to theoretical or institutional arguments for why a small number of controls should be able to soak up most of the confounding. We thus develop two statistical tests of the null hypothesis that (approximate) sparsity holds. Both tests tend to reject in the baseline specification of the three empirical applications. What is more, we are often unable to find any normalization choice for which the tests do not reject, suggesting that this second concern is not just a variation of the first.
Our conclusion is that applied researchers should be wary of using \acp{SBE} without a substantive defense of the sparsity assumption for the specific choice of the control matrix used in the analysis. Unless $p$ is close to or exceeds $n$, a robust default is to simply run \ac{OLS}.\footnote{One only needs to be careful in standard error construction, as \ac{OLS} no longer consistently estimates the residuals once $p\asymp n$. This renders the usual Eicker-Huber-White standard errors invalid; instead one needs to use a standard error formula that is robust to high-dimensional controls cjn18,DoSu18,dadamo19,kss20,jochmans22.}
Our analysis relates to several strands of literature. Many authors mention the \acp{SBE}' lack of invariance to linear reparametrization as an undesirable feature. In the specific context of including categorical controls, alternatives to standard lasso regression have been developed that are invariant to the choice of baseline category BoRe09,GeTu10,StShTi21. Similarly, one could perhaps incorporate the choice of a centering constant when taking powers as an additional tuning parameter when fitting \acp{SBE}. However, as we discuss further in (ref) below, such solutions come with their own challenges. Using a Bayesian framework, glp21 present empirical evidence suggesting sparsity may not be a compelling assumption in popular data sets used in economics. \Textcite{WuZh23} give simulation evidence and theoretical arguments showing that even if sparsity arguably holds, \acp{SBE} may still display substantial bias in finite samples. \Textcite{AnFr22} explore robustness of \acp{SBE} to tuning parameter choices.
The remainder of the paper is organized as follows. In (ref), we introduce the three empirical examples used to quantify fragility of \acp{SBE}: the empirical illustration in bch14 on the effect of abortion on crime, a ferrara22 study of occupational upgrading by Black southerners, and a study of the effect of moral values on voting behavior by enke20. Our analysis of fragility to linear reparametrizations focuses on two seemingly innocuous normalizations: first, we consider different ways of dropping collinear columns, such as which category to drop when including categorical variables. Second, when $W$ includes powers and interactions of a set of baseline controls, we consider different ways of centering the baseline controls, such as demeaning, subtracting the median or no centering. We find that these normalizations can move the estimates by two standard errors to more. We also consider alternative normalizations that further expand the set of possible specifications of $W$: expressing a set of categorical controls as indicators for subsets of the categories, rather than just as indicators for each individual category; and continuously varying where a variable is centered prior to taking powers, rather than just considering discrete values such as the mean or the median. These normalizations lead to an even greater variation in the estimates.
In (ref), we interpret these empirical results through the lens of a thought experiment: suppose a researcher picks one of the possible normalizations at random. How likely would such a random choice lead to an (approximately) sparse representation? We compute that probability in three stylized examples of normalizations and find that it rapidly decreases as function of $n$. Leaving default choices for how $W$ is constructed to statistical software or happenstance is very unlikely to lead to sparsity.
(ref) conducts an analysis of the relative efficiency of \acp{SBE} and \ac{OLS}. We focus on the regime where $p$ is smaller, but proportional to $n$. This is because if $p/n\to 0$, \ac{OLS} achieves the semiparametric efficiency bound under homoskedastic errors, which limits the arguments in favor of alternative estimators. On the other hand, if $p$ exceeds $n$, the \ac{OLS} benchmark is no longer available. We show that under homoskedasticity, the proportional variance reduction of \acp{SBE} relative to \ac{OLS} is capped at $1-p/n$ (effectively, we can at best avoid a degrees of freedom adjustment). This benchmark forms the basis for our recommendation to simply run \ac{OLS} unless $p$ is close to $n$: it yields robust estimates at little efficiency loss.
In (ref), we develop two tests of the sparsity assumption, again in the regime $p\asymp n$.\footnote{CaVe21 develop a test in the high-dimensional regime $p/n\to\infty$. However, their test requires explicit specification of the sparsity level under the null, as does the test developed by he20 in the regime $p<n$. In the context of factor models, BeSt24 consider testing the null of a zero vector against a sparse alternative.} We circumvent the difficulty that the sparsity assumption only restricts rates, rather than the actual number of non-zero coefficients for a given $n$, by testing it indirectly. The first test compares \ac{OLS} and lasso residuals: under (approximate) sparsity, the residual sum of squares of the lasso should exceed the \ac{OLS} residual sum of squares only by a small amount. The second test is an application of the hausman78 specification test: any difference between \ac{SBE} and \ac{OLS} estimates must be explained by differences in the relative efficiency of the estimators afforded by the sparsity assumption that \ac{SBE} exploits, but \ac{OLS} does not. Thus, the common practice of reporting \ac{SBE} estimates alongside \ac{OLS} should not be interpreted as a robustness check for the \ac{OLS} specification; rather, divergence between the estimates indicates a failure of the sparsity assumption.
(ref) concludes. All proofs are relegated to appendices.
In this section we revisit three applications that leveraged \acp{SBE}: the application in bch14 that probes the investigation by DoLe01 of the impact of abortion on crime, the ferrara22 study of employment opportunities for Black southerners in the aftermath of \ac{WW2}, and the examination of the relationship between moral values and voting behavior by enke20. In each application, we apply the original \ac{SBE} to (ref) after changing the specification of the control matrix $W$ (where rows correspond to the control vectors $W_{i}$) in seemingly innocuous ways that do not impact \ac{OLS} estimates. We show that the impact of these normalizations on \acp{SBE} is, on the other hand, substantial. To isolate the effect of the choice of control matrix from other implementation details, our analysis otherwise sticks to software defaults.\footnote{For estimates using post-double selection, we use the hdm package and its defaults, including the choice of the tuning parameter and normalizing the columns of $W$ to have unit variance.}
In (ref), we consider effects of two normalizations. First, we consider different ways of resolving multicollinearity in the control matrix by changing which columns are dropped to make the matrix $W$ full rank. For example, if multicollinearity arises due to the inclusion of categorical variables, we change which category is dropped. Original implementations of \acp{SBE} in each of the three applications remove multicollinearity in $W$ as a data-processing step, similar to standard implementations of least squares regression. Our exercise is therefore equivalent to changing the order of the columns of the control matrix, but otherwise conducting the analysis exactly as in the original. For lasso-based estimators, such a step is typically needed to ensure that so-called “compatibility conditions” or “restricted eigenvalue” assumptions hold; these assumptions ensure fast convergence rates for the lasso brt09.\footnote{A necessary condition for these assumptions is that submatrices of $W$ with $2s$ columns, where $s$ is the sparsity index, are full rank. It implies that if, say, a categorical variable has fewer than $2s$ categories, we cannot include all categories as well as the intercept.} Second, when $W$ includes powers and interactions of a baseline set of controls, we consider different normalizations of the baseline controls (such as demeaning vs subtracting the median) before taking powers and interactions. Appropriately centering the baseline variables puts them on the same scale, and may make the sparsity assumption more plausible and more easily interpretable. Relative to using raw polynomials, it typically also helps with numerical stability.
(ref) considers alternative normalizations of $W$ in the two scenarios above: expressing a set of categorical controls as indicators for subsets of the categories, rather than only as indicators for each individual category; and continuously varying where a variable is centered prior to taking powers and interactions, rather than just considering discrete values such as the mean or the median.
We begin by considering different ways of resolving multicollinearity issues that arise in each application.
The data in the first application consists of an annual panel of US states over the period 1985--97 originally analyzed by DoLe01, who ran a two-way fixed effects regression of crime rates on effective abortion rates, state and year fixed effects, and 8 baseline covariates.\footnote{These are: lags of the number of prisoners and police per capita, the unemployment rate, per-capita income and beer consumption, poverty rate, AFDC generosity lagged 15 years, and a dummy for a shall-issue concealed carry law.} BCH argue that this set of 8 controls may be insufficient to purge time-varying confounders. They consider a first-differences version of this specification, with 12 time effects and first-differences of the 8 baseline controls, which they augment with 136 further variables obtained by squaring the baseline controls and interacting them with each other and with linear and quadratic trends.\footnote{In particular, BCH add squares of the first differences as well as lags and squared lags of the baseline controls, and they interact the first differences and their squares with a linear and quadratic trend. They also add interactions of the first-differences, and interact them with a linear and a quadratic trend. Since the time series of the shall-carry dummy in each state never changes from $1$ to $0$, its first difference is also binary. These transformations therefore only yield 136 unique columns. } In addition, they include 49 variables that are time-invariant within each state, corresponding to initial values and averages of various transformations of the baseline controls, as well as interactions of these 49 variables with a time trend and its square. Since there are only 48 states in the data, these variables span the same column space as state fixed effects. The resulting control matrix has 303 columns, but rank of only 294: because time effects are included, 2 of the time-invariant variables are redundant, as are 2 of the interactions with a time trend and 2 with its square. In addition, one of the baseline controls, a shall-issue dummy, is binary and non-zero only 21 times, but it is interacted with 24 variables, so that 3 of the interactions are redundant. BCH first drop the collinear columns, and then estimate the model by the post-double lasso estimator that they develop.
Column 2 of Panel A in (ref) replicates the BCH estimates for each of the three versions of the crime rate outcome variable considered by BCH.\footnote{The replication differs slightly from the original (Table 2 in BCH) because the algorithm for dropping collinear columns in BCH had a coding error, yielding a matrix with 296 columns, and a rank of 293.} Comparing the standard error to the \ac{OLS} standard error in column 1, we see that the post-double lasso estimator appears to be substantially more precise than the \ac{OLS} estimate. It is also economically significant: the effective abortion rate is defined as the average legalized abortion rate among arrestee cohorts (the number of abortions per live birth in a cohort weighted by the cohort's share of arrestees, which is outcome-specific). For the violent crime outcome, its standard deviation across states in 1997 is about 1, so that the estimate implies a 16% reduction in crime rate per standard deviation increase in effective abortion rate. However, this result is very sensitive to how we resolve the collinearity. To resolve it, we may drop any 3 of the 24 interactions with the shall-issue dummy, and any 2 of the 49 time-invariant controls; the same holds when we interact these with a time trend and its square. This gives $\binom{24}{3}\binom{49}{2}^{3}\approx 3\times 10^{12}$ possible ways of obtaining a full rank control matrix. Columns 3 and 4 report a range of estimates we obtain by randomly choosing among these possibilities, showing that, depending on the outcome, the estimates move by 1.2 to 1.9 standard errors.
Similar collinearity issues arise in the second application that replicates the analysis in ferrara22, who studies to what extent post \ac{WW2} occupational upgrading of Black workers from low-skilled to semi-skilled can be attributed to war casualties among semi-skilled white soldiers. Using a decennial 1920--1960 unbalanced panel of county-level observations in 16 predominantly Southern US states, the study runs a two-way fixed effects specification with county and time fixed effects, interactions between state and time fixed effects, the share of semi-skilled Black workers as the outcome, and treatment given by the white casualty rate interacted with a post-war indicator. To purge time-varying confounders, the study also includes 24 baseline controls (including the county draft rate, \ac{WW2} spending, demographic and socioeconomic controls), their squares and interactions, as well as interactions between the 24 baseline controls and state and time effects. In addition, two baseline controls (number of slaves in 1860 and the unemployment rate in 1937) are also included in triple interactions with other controls, time and state effects. After dropping a reference state and a reference year in each interaction, and dropping zero and repeated columns, the resulting matrix has 2270 columns, but rank of only 2252, because the state of Delaware only contains 15 observations, but there are 33 Delaware-specific controls.\footnote{In particular, since the two baseline controls entering triple interactions are time-invariant and the specification includes county effects, we need to exclude the main effects of these controls, their squares and state-interactions, which yields a control matrix with 2588 columns. 154 columns are collinear since we need to drop a reference state and year in each interaction, and 18 Delaware-specific controls are also collinear. In addition, 164 columns are repeated or zero, yielding a rank equal to 2252.} ferrara22 first forms a full-rank control matrix, and then uses double-$t$ selection to estimate the model.\footnote{Specifically, ferrara22 first runs a selection step that regresses the outcome and the treatment on the control matrix using \ac{OLS}, and selects controls with a $t$-statistic that is larger than 2.575 in absolute value in either regression. In a second step, the outcome is regressed on the treatment and selected controls using \ac{OLS}.}
We replicate this specification in column 2 of Panel B in (ref).\footnote{The estimate is slightly higher than in the original because, to keep missing data and normalization issues separate, we restrict the sample to the 4,903 observations with no missing values in both estimation steps described in (ref). In addition, we include year fixed effects in both the first and second step, rather than exclude them from the first step.} However, similar to panel A, there are many ways of resolving multicollinearity in the control matrix, because there are many ways of specifying a reference state and a reference year in each interaction, and, if we keep 3 Delaware county effects, $\binom{30}{18}$ ways of dropping Delaware-specific controls. As shown in columns 3 and 4, depending on how the multicollinearity is resolved, the estimates vary by over three standard errors depending on how we order the state and time effects.\footnote{We also considered sensitivity of the results to ordering of the columns in the control matrix when the double-$t$ estimator is replaced by the post double lasso estimator. This yielded a range of variation in the estimates equal to about two thirds of a standard error.}
Our final empirical example uses data from enke20, who examines how voters' moral values affect their voting patterns. Enke uses survey data to construct an index of the relative importance of universalist moral values (individual rights, justice, fairness) vs communal values (loyalty, respect). He then runs a regression of three different measures of voting behavior on this index, controlling for 10 continuous or binary controls, as well as 5 sets of categorical variables.\footnote{The continuous or binary controls comprise: political liberalism, log of household income, education, log of population density of respondent's ZIP code, religiosity, gender, employment indicator, altruism, measure of trust, and the absolute value of the moral values index. The categorical variables are: county, year of birth, religious denomination and occupation fixed effects.} Columns 1 and 2 in panel C of (ref) replicate the \ac{OLS} and post-double lasso estimates reported in the paper. Here the inclusion of multiple sets of categorical variables necessitates a specification of a reference category for each set. Columns 3 and 4 report the range of estimates obtained by changing these reference categories, which varies between a third and a half of standard error depending on the outcome. This is a narrower range than in panels A and B, albeit still large enough to affect the economic interpretation of the estimates.
Columns 5 and 6 of (ref) show the range of estimates we obtain from considering different normalizations when taking powers and interactions between variables in the BCH and ferrara22 applications. In addition to no normalization (the choice in the original analyses), we consider demeaning, centering at the median, and setting the range to $[-1,1]$ and $[0,1]$. The table shows that such normalizations can again substantially move the estimates, by up to 1.3 standard errors.
We now consider two alternatives to the construction of $W$ in the presence of categorical variables and power and interactions.
In (ref), a categorical variable with $k$ categories was always specified as a set of $k-1$ dummies for individual categories, with a reference dummy dropped to prevent collinearity. We now consider expressing such variables as $k-1$ indicators for different subsets of the $k$ categories. As discussed in (ref), if the subsets are carefully chosen, such a specification may be more likely to lead to sparsity in certain contexts. Here we pick the subsets at random, subject to the constraint that they span the same column space. For the enke20 application, columns 3 and 4 of (ref) correspond to a special case of this exercise, when the subsets are constrained to be of cardinality one. Correspondingly, the range of variation in the estimates over all possible subsets, reported in columns 1 and 2 of (ref), is considerably greater, exceeding 0.9 standard errors in all regressions. For the ferrara22 application, the range exceeds two standard errors.\footnote{The range does not exceed that in columns 3 and 4 of (ref) because we do not vary which Delaware-specific controls are dropped.}
Similarly, the centering of baseline variables at the mean vs the median we considered in the previous subsection can be thought of as a special case of a more general problem: we could in principle choose the centering constant $\lambda$ to be any number. At a theoretical level, we are not aware of arguments for why $\lambda=0$ (no normalization), or setting it to the sample mean should more plausibly lead to sparsity than other choices. To further explore the sensitivity to the choice of $\lambda$, we standardize the baseline variables, and then vary the specification of $W$ by picking $\lambda$ from the range $[-1, 1]$ for each baseline variable. To align the exercise with the analytical calculations in (ref), we use Hermite rather than raw polynomials. Columns 3 and 4 of (ref) show that the resulting variation in the estimates is again substantial, exceeding 1.7 standard errors in all specifications. The sign of the estimated effect is also sensitive to the choice of $\lambda$ in two of the BCH specifications.
In this section we consider the thought experiment of picking a particular linear representation of the controls $W$ at random. In particular, we consider (i) a full rotation of the column space; (ii) different ways of controlling for categorical variables; and (iii) an offset for the scalar variable in a polynomial series regression. We find that sparse representations are exceedingly rare in each example. This suggests that in a given application, undiscerning default choices for expressing the column space will typically not satisfy the sparsity assumptions. This helps to explain the findings of the previous section: most normalization choices cannot yield sparse representations. As we discuss in detail in (ref) below, it also implies that the representation of the control matrix is a substantive choice and not just an implementation detail, like, say, the choice of the penalty parameter: it impacts not just the finite-sample estimate, but also the large-sample validity of \ac{SBE}-based inference.
For simplicity, the calculations in this section focus on sparsity in the outcome regression (ref). Let $A$ be a full-rank $p\times p$ matrix. Then an alternative expression for the column space of $W_{i}$ is $\tilde{W}_{i}=AW_{i}$ with coefficient $\tilde{\gamma}=A'^{-1}\gamma $, so that $W_{i}^{\prime}\gamma =\tilde{W} _{i}^{\prime}\tilde{\gamma}$. Note that any such transformation leaves the \ac{OLS} coefficient on $D_{i}$ numerically unaltered.
In our examples, $\gamma $ is exactly sparse, with just one non-zero element, $\norm{\gamma}_{0}=1$. The researcher uses the transformed regressors $\tilde{W}_{i}$ with coefficients $\tilde{\gamma}$, which typically is not exactly sparse. However, for sparsity-based inference to be asymptotically valid, it suffices for $\tilde{\gamma}$ to be approximately sparse. In particular, under our maintained assumption that $p\asymp n$, it suffices that for a sparsity index
the mean square approximation error of $W_{i}^{\prime}\gamma $ satisfies
bch14. Here, and throughout the paper, all limits are taken as $p\rightarrow \infty $ (or, equivalently, as $n\rightarrow \infty $).
By way of establishing a benchmark, we start with an extreme case, and consider the set of transformations
that is, $\tilde{W}_{i}$ is obtained by rotating $W_{i}$ by $R$, and $\tilde{\gamma}=R\gamma$. For the purposes of studying sparsity, this is a very large set of transformations, as for any $\gamma $ with $\norm{\gamma}_{2} >0$, there exists $R$ such that $\tilde{\gamma}$ has only one non-zero element. At the same time, there also exists an $R$ such that all elements of $\tilde{\gamma}$ are equal (and equal to $\norm{\gamma}_{2}/\sqrt{p}$).
Now suppose we take a random draw $\mathcal{R}$ from $\mathbf{R}$, with distribution equal to the Haar measure (so that for any $R\in \mathbf{R}$, $R\mathcal{R}\sim \mathcal{R}$). Recall that the multivariate standard normal distribution is spherically symmetrical, that is, for $\mathcal{Z}\sim \mathcal{N}(0,I_{p})$ and any $R\in \mathbf{R}$, $R\mathcal{Z}\sim \mathcal{Z}$. This implies that the induced distribution of $\tilde{\gamma}=\mathcal{R}\gamma $ satisfies
so the $p$ elements of $\mathcal{R}\gamma$ are, up to a common rescaling, i.i.d. standard normal. The normal distribution is thin-tailed, making it very unlikely that a few largest absolute values of $\mathcal{Z}$ dominate. This suggests that with very high probability, $\mathcal{R}\gamma$ is not approximately sparse, as formalized in the following result.
Given that a sparsity assumption allows for more informative inference, and that any coefficient vector can be rotated into an (extremely) sparse one, it is not surprising that a random rotation only rarely yields a sparse representation. The contribution of (ref) is to quantify just how exceedingly rare this event is: for $p\geq 50$, say, $p^{-p/4}<10^{-21}$. Of course, researchers rarely consider all rotations when specifying the regressor matrix. The relevance of this quantification lies in showing that even when the measure of plausible rotations is orders of magnitude smaller than the full set, the probability of arriving at a sparse representation is still minuscule.
An empirically common form of controls is dummies for categorical variables. If there are $p$ underlying categories, and a constant is included, then one must designate a reference category by dropping one dummy to avoid perfect collinearity. Suppose the model is exactly sparse with sparsity index $s$ when one chooses the right reference category. This means that $p-s$ categories have the same coefficient as the reference category. Thus, if one were to pick a reference category at random, one has a $1-s/p$ chance of inducing sparsity, which is a probability close to one under (ref).
However, an assumption that many categories have the same effect is not the only way to induce sparsity for a categorical variable. Suppose, for instance, that we want to nonparametrically control for a variable that measures age. Then a plausible form of sparsity arises under an assumption that the age effect is a step function, with some dividing line between young and old. If this threshold is not close to either very young or very old, then a reference category approach will not yield sparsity, but an appropriate re-expression of the column space will. Or maybe there are three distinct coefficients, one for the young, one for the middle-aged and one for the old, which again requires a different specification of the fixed effects to induce sparsity. Such specifications are less common in applied work. Arguably this is not because they are a priori less implausible, but rather, it is well understood that under \ac{OLS}, they are all equivalent to the reference category specification. However, when imposing a sparsity assumption on the coefficients, this equivalence breaks down.
A general specification of the column space with categorical data takes the form
where the $Z_{i}$ indicates the baseline categories (so each $Z_{i}$ is a column of $I_{p}$). By construction, elements of $W_{i}$ are all equal to zero or one in this specification, and since $A\in\mathbb{R}^{p\times p}$ is full rank, the constant vector is part of the column space spanned by $W_{i}$.
Suppose with $A_{0}\in \mathbf{A}$, the coefficient vector on $W_{i}$ is sparse. Now consider picking $\mathcal{A}\in\mathbf{A}$ at random by repeatedly drawing $p\times p$ matrices with i.i.d. Bernoulli$(q)$ entries, $0<q\leq 1/2$, until we find one that is full rank. By Theorem A of tikhomirov20, the probability of discarding a matrix in this process is vanishingly small, as it is smaller than $(1-q+\varepsilon)^{p}$ for all $\varepsilon >0$ and large enough $p$. This allows us to essentially ignore rank-deficient matrices, yielding the following result.
Our last example involves a large number of technical regressors. Specifically, suppose there is a scalar variable $z_{i}\sim_{\text{i.i.d.}}\mathcal{N}(0,1)$ which we want to control for nonparametrically by a series regression. We consider scaled Hermite polynomials
which are particularly convenient for an underlying standard Gaussian variable, since we have $E[\tilde{W}_{i}\tilde{W}_{i}^{\prime}] = I_{p}$ for $\tilde{W}_{i, j}=H_{j-1}(z_{i})$.
Suppose further that if $z_{i}$ is shifted by a constant $\lambda \in\mathbb{R}$, the expansion is extremely sparse, with only the highest order term being non-zero, that is, $Y_{i}=D_{i}\beta +H_{p-1}(z_{i}+\lambda)+U_{i}$. If we were to draw $\lambda $ uniformly on $[0,1]$, what is the probability that a researcher ignoring the shift obtains an approximately sparse representation? The following result shows that the probability will be bounded above by $\varepsilon/\log (p)$, for all $\varepsilon > 0$ and large enough $p$.
While the rate of $\log p$ in (ref) is quite slow, note that when sparsity fails, it does so dramatically in the sense that many coefficients diverge. This is not surprising. It is well-known that assuming that the first $s$ (say) terms of the polynomial expansion provide a good approximation to the regression function amounts to imposing a smoothness assumption on it. In this case, the exact basis approximation would matter little, and sparsity-based methods would perform well regardless of the choice of $\lambda$. However, under smoothness \acp{SBE} are not needed, since a low-order series expansion (i.e.\ only controlling for the first $s$ polynomial terms, with $s\ll p$) already performs well.
In the current context with technical regressors, the appeal of \acp{SBE} is that the assumption of sparsity they leverage is weaker than assuming smoothness, since we allow for the best $s$-term approximation to include high order polynomials. But if a good approximation to the regression function indeed requires high-order terms, it must be very non-smooth---and for non-smooth functions even small changes in the approximating basis imply wild swings in the coefficients. In other words, while sparsity is indeed weaker than smoothness, the particular class of non-smooth functions sparsity allows for is tightly linked to the basis used. If one suspects non-smoothness, and one wants to avoid running \ac{OLS} with many polynomial terms, one therefore needs a substantive argument for why the particular basis chosen is likely to capture it.
The fact that the sparsity assumption depends on normalization choices in the specification of $W$ is well-known in the theoretical literature. The three stylized thought-experiments above quantify the extent of this dependence: the vast majority of normalization choices do not induce an approximately sparse representation. The results in (ref) show that this lack of invariance induces variation in \ac{SBE} estimates that is of the same order as sampling uncertainty. This fragility appears underappreciated in the literature. We reviewed 50 most cited papers applying \acp{SBE}, and found that, with one exception, none made explicit what the normalization choices were, let alone discussed why they were appropriate to induce sparsity.\footnote{We ordered all papers citing one of the “debiased lasso” papers mentioned in the introduction by number of citations on Google Scholar. One paper mentioned that it included all dummy variables and did not pick a reference category; the paper did not discuss the implications of this choice for the restricted eigenvalue condition, or whether other normalization choices were made.}
A natural reaction to these results is to seek to modify \acp{SBE} in a way that they adapt to different forms of sparse representations. For instance, it may be possible to treat the offset $\lambda $ as another parameter to be estimated. Similarly, a number of proposals have been developed that modify the lasso to make it invariant to the choice of reference category, including the group lasso YuLi06, or variants of fused lasso tsrzk05, such as all-pairs penalties BoRe09,GeTu10 or the SCOPE estimator StShTi21. When combined with the debiasing techniques developed in the literature on debiased lasso JaMo14,vbrd14,ZhZh14,bch14, these approaches may yield algorithms that are less fragile.
However, it appears non-trivial to extend these modifications to the more complicated designs encountered in practice, such as when interactions are present, or to handle collinearity issues arising due to limited variability in the controls. Both the BCH and the ferrara22 applications of (ref) exhibit multicollinearity that goes beyond the issue of choosing baseline categories.
An alternative approach could be to exhaustively consider all plausible normalization choices, analogous to our analysis in (ref), and report the union of the confidence intervals. However, this is practically difficult, as the number of possible normalizations is often very large. Moreover, such an approach would only yield valid inference if at least one of the considered choices induces a sparse representation.
To ensure that the latter condition holds, one needs a substantive argument. In some instances, such arguments may necessitate further expanding the set of normalization choices. For instance, the grouped lasso approach in the references cited above require that the baseline categories can be partitioned into a small number of types, with all categories of the same type having (nearly) identical effect on the outcome variable. In the context of controlling for profession dummies, say, this amounts to an assumption that there are only a few profession types with heterogeneous effects, which is a strong assumption. It may be more reasonable to instead posit that professions are characterized by a combination of a small number of latent skills, and these skills affect the outcome in an additive manner. This then yields a different set of plausible normalizations that are not captured by grouped lasso-type approaches.
In absence of substantive arguments, the sparsity assumption may well fail to hold, even after considering a large set of potential normalization choices. We present evidence for this in our applications in (ref) below. Thus, fully data-driven methods to modify \acp{SBE} such that they adapt to an appropriate sparse representation regardless of the context does not seem like a feasible solution.
Applications of \acp{SBE} require researchers to take a stand on why a particular representation (or a set of representations) admits a sparse representation based on domain-specific knowledge. This is analogous to using problem-specific arguments to defend other substantive assumptions, such as assuming selection on observables. But it is at odds with the purported appeal of \acp{SBE} to be fully automatic, with the aim of reducing the researchers' degrees of freedom in model choice.
While the sensitivity of \acp{SBE} to the control matrix specification may seem similar to the usual issue that with flexible methods, implementation issues, such as tuning parameter choice, can matter a lot in finite samples, it is fundamentally different. It also differs from the well-known fact that the choice of basis matters in implementing series regression. Under smoothness, many basis choices for series regression are theoretically justified, as well as first-order equivalent in particular settings such as the partially linear model. Large sample validity of flexible methods can be justified under a wide range of tuning parameter choices, including the choice of penalty in implementing \acp{SBE}. In contrast, because the specification of the control matrix affects the validity of the sparsity assumption, and thus large-sample validity of inference, it is a substantive modeling choice, not merely an implementation detail. It needs to be defended as such.
One argument for using \acp{SBE} when $p$ is large, but still smaller than the sample size, is that \ac{OLS} estimates are too noisy. In this section, we quantify the potential efficiency gains of alternative estimators relative to \ac{OLS}, and we argue that the potential gains are limited unless $p$ is comparable to $n$.
We consider the model given in (ref) in the introduction. We maintain the assumption that the residuals $U_{i}$ and $\tilde{D}_{i}$ are conditionally mean zero to avoid technical complications; the assumption could be replaced by an assumption that restricts these conditional means to be small, so that the linear specifications in (ref) could be thought of as approximating non-linear conditional means cjn18. Using the Frisch-Waugh-Lowell theorem, we can write the \ac{OLS} estimator of $\beta$ as
where $\ddot{D}_{i}$ denotes the sample propensity score residual from estimating (ref) by \ac{OLS}. Specifically, employing the usual matrix notation, $\ddot{D}_{i}$ corresponds to the $i$th element of $\ddot{D}:=(I-P)D$, where $P=WW^{+}$ is the projection matrix onto the space spanned by the controls. Here $W^{+}$ denotes the pseudo-inverse (so that if $W$ is full rank, $P=W(W'W)^{-1}W'$).
To analyze $\hat{\beta}_{OLS}$, we need to impose more structure on the model in (ref).
Parts (i)--(iii) of (ref) are standard assumptions on sampling and the regression errors, ensuring that the errors are neither too thick-tailed nor degenerate. Part (iv) ensures that, when \ac{OLS} estimates of (ref) implicitly estimate the propensity score regression, we do not overfit so much that we eliminate any variation in the sample residuals $\ddot{D}_{i}$. This overfitting precludes using \ac{OLS} in settings with $p>n$. However, as long as the limit of $p/n$ is strictly smaller than $1$, the \ac{OLS} estimator remains asymptotically normal with the usual sandwich expression for its standard error:
The fact that \ac{OLS} remains asymptotically normal in the regime $p\asymp n$ is not new---results similar to (ref) appeared previously in, for instance, cjn18 who build on earlier work by mammen93. We state it here under a simpler set of assumptions to allow us to consider its implications for potential efficiency gains. To this end, observe that the asymptotic variance is in general larger than the semiparametric efficiency bound unless $p/n\to 0$. In particular, when the errors $U_{i}$ are homoskedastic, the efficient standard error is given by the square root of
in the sense that the semiparametric efficiency bound is given by the probability limit of $s_{*}^{2}/n$ robinson88. Semiparametrically efficient estimators $\hat{\beta}^{*}$ of $\beta$ thus need to satisfy
This lack of semiparametric efficiency is not a deficiency of \ac{OLS}: under Gaussian errors, one-sided $t$-tests based on \ac{OLS} are uniformly most powerful. Rather, it reflects the fact that when $p\asymp n$, achieving the semiparametric efficiency bound requires some type of restriction on $\gamma$. It turns out that restricting $\gamma$ to be sparse (see (ref)) is sufficient; indeed a number of \acp{SBE} satisfy (ref) bch14,JaMo14,vbrd14,ZhZh14.
The standard error ratio $s_{*}/s_{OLS}$ thus represents the potential efficiency gain from imposing sparsity. When $U_{i}$ is homoskedastic and (ref) holds, the ratio satisfies
where $\sigma^{2}_{\tilde{D}}$ is the variance of the error in the propensity score regression, and $\hat{\sigma}^{2}_{\tilde{D}}=\ddot{D}'\ddot{D}/(n-p)$, the mean squared error in the propensity score regression, can be thought of as an estimator of $\sigma^{2}_{\tilde{D}}$. The ratio $\kappa$ measures the relative bias of this estimator, with $\kappa=1$ when the estimator is unbiased. This is the case when $\tilde{D}_{i}$ is homoskedastic, and the standard error ratio $s_{*}/s_{OLS}$ then simplifies to a degrees-of-freedom correction, $\sqrt{1-p/n}$, as noted previously in cjn18et. When $\tilde{D}_{i}$ is heteroskedastic, $\hat{\sigma}^{2}_{\tilde{D}}$ may display downward bias if the leverages $P_{ii}$ are positively correlated with the conditional variances $E[\tilde{D}_{i}^{2}\mid W]$; this can be seen from writing $\kappa=E[\sum_{i}(1-P_{ii})E[\tilde{D}_{i}^{2}\mid W]]/ \sum_{i}(1-P_{ii})E[\tilde{D}_{i}^{2}]$. However, the correlation would need to be substantial to induce downward bias that is sizable enough to allow for large efficiency gains. When $p/n=0.2$ and $\kappa=0.9$, for instance, corresponding to 10% downward bias, the standard error ratio implies a $1-\sqrt{0.8\cdot 0.9}=15.1\%$ reduction in standard error. Consequently, the ratio $p/n$ needs to be much larger than $0.2$ to allow for sizable efficiency gains under sparsity: at $\kappa=0.9$, we need $p/n\geq 6/16\approx 0.38$ to reduce the standard errors by more than 25%.
While these efficiency calculations rely on homoskedastic errors $U_{i}$, similar conclusions are likely to hold under heteroskedasticity. The argument for this is that, in general,
where $s_{*, \hom}=\sigma_{U}/\sqrt{\tilde{D}'\tilde{D}}$ and $s_{OLS, \hom}=\sigma_{U}/\sqrt{\ddot{D}'\ddot{D}}$ are homoskedasticity-only standard error formulas. The ratios $s_{*}/s_{*, \hom}$ and $s_{OLS}/s_{OLS, \hom}$ measure the magnitude of the heteroskedasticity correction on the standard errors. Heteroskedasticity would thus need to have a large differential impact on the standard errors for \ac{OLS} versus $\hat{\beta}_{*}$ for the standard error ratio to deviate much from $\sqrt{(1-p/n)\kappa}$. But magnitudes of heteroskedasticity corrections tend to be modest in practice, so that substantive efficiency gains are unlikely when the controls number 20% or less of the sample size. In such cases, simply running \ac{OLS} offers a robust method of inference at little efficiency loss.
In this section, we develop two tests of the sparsity assumption under our maintained assumption that $p$ is smaller than, but proportional to the sample size. We then apply these tests to the empirical examples studied in (ref).
The first test we develop is a simple application of the idea popularized by hausman78 that if we have two estimators, one of which is more efficient, but requires stronger assumptions for its validity, we can indirectly test these stronger assumptions by checking whether the estimates are statistically significantly different. The next lemma formalizes the result.
The second part of (ref) allows us to construct a simple plug-in estimator of the standard error $s_{H}$ based on lasso or post-lasso residuals. In particular, when the regression functions admit a sparse approximation in the sense of (ref), the additional condition (i) in the second part of (ref) will hold with $f$ and $g$ given by the best sparse approximations to $W_{i}'\gamma$ and $W_{i}'\delta$, and condition (iv) will hold for the lasso or the post-lasso residuals. Conditions (ii) and (iii) are high-level conditions ensuring that if we include the approximation errors $r_{\delta}$ and $r_{\gamma}$ in the definition of the residuals, replacing $U_{i}$ with $U_{i}+r_{\gamma}(W_{i})$, and $\tilde{D}_{i}$ with $\tilde{D}_{i}+r_{\delta}(W_{i})$, this has negligible impact on the standard error $s_{H}$ in large samples; it is similar to the condition ASTE (P) (v) in bch14. We note that (ref) only imposes very weak conditions on the control matrix $W$, allowing it to be reduced-rank, so long as the rank is proportional to $p$, and allowing the rows to be dependent, and not identically distributed.
Sometimes, \acp{SBE} are used as a “robustness check” alongside a main specification based on \ac{OLS}. (ref) shows that such practice is in fact the opposite of a robustness check: if the two estimates are not close to one another this indicates failure of the sparsity assumption rather than lack of robustness in the \ac{OLS} estimates. When $U_{i}$ is homoskedastic, the Hausman standard error may be written as
so that when the efficiency gain $\sqrt{(1-p/n)\kappa}$ is small, the two estimates need to be tightly coupled, within a fraction of the \ac{SBE} standard error.
Our second approach to testing sparsity is based on the idea that if the identities $\mathcal{S}^{*}\subseteq \{1,\dotsc, p\}$ of the controls that give the best sparse approximation to the regression function are known, testing sparsity is equivalent to testing that the coefficients on the remaining controls are small, which can be gauged using a conventional $F$-statistic. Although the identities $\mathcal{S}^{*}$ are unknown in practice, it turns out that under the null hypothesis of sparsity, the residuals from the infeasible short regression that only includes controls in $\mathcal{S}^{*}$ are sufficiently well approximated by residuals from a lasso regression, so that the hypothesis can be tested by comparing the lasso and \ac{OLS} residuals.
To make the result precise, consider a linear regression
with $\epsilon_{i}$ independent across $i$, conditional on the regressors, and $p:=\dim(X_{i})<n$. In our setup, (ref) may correspond to one of three regressions. For testing sparsity of $\gamma$ in (ref), which all \acp{SBE} require, we can set $\mathsf{Y}_{i}=Y_{i}$ and $X_{i}=(D_{i}, W_{i}')'$. For testing sparsity in the propensity score regression, which is needed, for instance, for the validity of the post-double selection estimator of bch14, we can set $\mathsf{Y}_{i}=D_{i}$ and $X_{i}=W_{i}$. Sparsity in both regressions can be tested jointly in the reduced form regression of $\mathsf{Y}_{i}=Y_{i}$ onto $X_{i}=W_{i}$.
We wish to test the assumption that the regression function $X_{i}'\alpha$ admits a sparse approximation. To work out the implications of this hypothesis, we introduce some additional notation. For a subset $\mathcal{S}\subseteq\{1,\dotsc, p\}$ of the regressors, let $X_{S}$ denote the submatrix of $X$ that drops the columns corresponding to the complement of $\mathcal{S}$, let $P_{\mathcal{S}}=X_{\mathcal{S}}X_{\mathcal{S}}^{+}$ denote the projection matrix associated with $X_{\mathcal{S}}$, and let $P=XX^{+}$ denote the full projection matrix. In a slight departure from (ref), we gauge the quality of the approximation conditional on $X$, so that the approximation error from only using the regressors $X_{\mathcal{S}}$ is given by $(I-P_{\mathcal{S}})X\alpha$, the residual from projecting $X\alpha$ onto $X_{\mathcal{S}}$. The assumption that $X\alpha$ is sparse can then be stated as:
If we knew the identity of the subset $\mathcal{S}^{*}$, and we also assumed that the sparsity was exact, so that the $O_{p}(s)$ term in (ref) was zero, then a natural way of testing (ref) would be to compare the restricted and unrestricted sum of squared residuals,
The statistic $\mathcal{F}$ is the numerator of the homoskedastic $F$-statistic. While the $F$-statistic critical values are only valid under homoskedasticity, we can leverage the fact that when $p\to\infty$, as is the case under our asymptotics, $\mathcal{F}$ is asymptotically normal after centering and scaling to derive a critical value that is robust to heteroskedasticity. Furthermore, if we have an estimator $\hat{\alpha}$ such that $\mathsf{Y}-X\hat{\alpha}$ approximates the residuals $(I-P_{\mathcal{S}^{*}})\mathsf{Y}$ from the infeasible short regression sufficiently well, we can replace the infeasible sum of squared residuals $\mathsf{Y}'(I-P_{\mathcal{S}^{*}})\mathsf{Y}$ in $\mathcal{F}$ by $\norm{\mathsf{Y}-X\hat{\alpha}}_{2}^{2}$ to derive a feasible version of this test. Finally, it turns out that weakening exact sparsity to approximate sparsity in the sense of (ref) does not impact the null rejection probability of the test in large samples. The next lemma formalizes these arguments.
Conditions (a) and (d) are standard assumptions on the regression errors, while conditions (b) and (c) impose weak restrictions on the design matrix; condition (c) is analogous to the condition imposed in (ref), and (b) bounds leverage away from one. Conditions (i) and (ii) in the second part of the \namecref{theorem:non-normal-errors_epe} are standard mean-squared error and $\ell_{1}$ rate conditions that hold for the lasso and post-lasso estimators brt09,BeCh13. Conditions (iii) and (iv) are tail restrictions on the covariates and residuals analogous to those in bch14.
(ref) implies that we can test the sparsity assumption by calculating the lasso or post-lasso residuals $\hat{\epsilon}_{i}$, and then checking whether the residual sum of squares of the lasso is comparable to that of the \ac{OLS} residual sum of squares. If the difference satifies
where $z_{1-\alpha}$ is the $1-\alpha$ quantile of the standard normal distribution (1.645 for $\alpha=5\%$ level test), we reject (ref).
In general, testing hypotheses about $\alpha$ when $\lim_{n\to\infty}p/n>0$ can be quite involved, since plugging in regression residuals into asymptotic variance expressions leads to bias kss20,AnSo23,cjn18. As shown in the proof of (ref) (see (ref)), we avoid such difficulties here by virtue of the fact that under the assumptions of the lemma, the lasso residuals $\hat{\epsilon}_{i}$ are sufficiently accurate so that the variance part of the test statistic, the denominator in (ref), can be consistently estimated under the null hypothesis.
We now apply the tests developed in (ref) to the empirical illustrations considered in (ref). For each specification, we consider the Hausman test, which compares \ac{OLS} with \ac{SBE}, as well two versions of the residual test. The first version estimates the outcome regression in (ref) using the post-lasso, and compares the residuals to those based on \ac{OLS}. The second version compares post-lasso and \ac{OLS} residuals from estimating the propensity score regression in (ref).
(ref) reports the results. Column 1 applies the test to the original specification in each paper. We see that for 6 of the 7 outcomes, at least one of the tests rejects the assumption of sparsity.
To check that these results are not driven by finite-sample size distortions of these tests, (ref) conducts Monte Carlo simulations based on these applications; in these simulations, the size stays close to or below the nominal level.
One response to these findings is to seek alternative normalizations of the control matrix that are consistent with the sparsity assumptions. The remaining columns in (ref) report the range of $p$-values under the four different normalizations we considered in (ref). The table shows that in these applications, we were unable to find a normalization for five of the seven outcomes where at least one test did not reject. This is in spite of the large number of alternative specifications that these normalizations generate.
We have argued, using empirical evidence and theoretical arguments, that \acp{SBE} display a lack of robustness to the specification of the control matrix. In the three applications we have examined, the range of variation in the \ac{SBE} estimates under equally plausible alternative specifications is of the same order of magnitude as sampling uncertainty. Two reasons underlie this fragility. First, whether a small number of covariates can account for the bulk of the confounding depends on the particular specification of the control matrix: even if the sparsity assumption holds under a particular way of expressing the column space of the controls, most alternative plausible normalizations do not admit a sparse approximation. Second, it may be the case no control matrix in a large class of normalizations admits a sparse approximation.
What should a practitioner take away from these results? We have argued that unless $p$ is comparable to $n$, the potential efficiency gains of \acp{SBE} over \ac{OLS} are limited. Simply reporting \ac{OLS} with standard errors that are robust to the presence of many controls cjn18et,jochmans22,DoSu18,dadamo19 will deliver credible inference at little efficiency loss. When $p$ is comparable to $n$, \ac{OLS} becomes too noisy to be useful, and infeasible when the covariate dimension exceeds the sample size. Sparsity restrictions on the control vector, or other types of restrictions such as limiting the magnitude of the control coefficients LiMu21,akk20, are then necessary for informative inference. However, since these are substantive modeling restrictions, they need to be discussed and defended on substantive grounds, analogous to the discussion of other key assumptions, such as selection on observables. In particular, researchers who opt to leverage \acp{SBE} need to explain why sparsity should plausibly hold under the chosen specification of the control matrix, and not leave normalization choices to statistical software. The sparsity tests developed in this paper can serve as a complement to these arguments, provided they serve as a model specification check rather than a pretest.
We have focused on the sparsity assumption and \acp{SBE} because these estimators are used frequently, and their theory is well-developed. However, many other modern machine learning methods likewise lack invariance to linear reparametrization of the control matrix. When these methods are used for prediction, this lack of invariance is less important, for two reasons. First, one is typically interested in average performance over many predictions, and the overall prediction performance may be robust even if individual predictions are sensitive to normalizations. Second, one can gauge the performance of a given procedure directly using a test sample. When we incorporate these methods into econometric models, however, we are typically interested in inference on a single causal effect, and test sample benchmarking is unavailable. Understanding more generally when a lack of invariance leads to fragility is an interesting area for future research.