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.
35,796 characters · 13 sections · 11 citation commands
Unbiased Shrinkage Estimation
\begingroup \footnote{ Jann Spiess, Department of Economics, Harvard University, [email removed]}. I thank Gary Chamberlain, Maximilian Kasy, Carl Morris, and Jim Stock for insightful conversations, and seminar participants at Harvard for helpful comments. } \addtocounter{footnote}{-1} \endgroup
Many inference tasks have the following feature: the researcher wants to obtain a high-quality estimate of a set of target parameters (for example, a set of treatment effects in an RCT), but also estimates a number of nuisance parameters she does not care about separately (for example, coefficients on control variables). In these cases, can we reduce variance in the estimation of a target parameter without inducing bias by shrinking in the estimation of possibly high-dimensional nuisance parameters? In a linear regression model with homoscedastic, Normal noise, I show that a natural application of James--Stein shrinkage to the parameters associated with at least three control variables reduces loss in the possibly low-dimensional treatment effect parameter without producing bias provided that treatment is random.
The proposed estimator effectively averages between regression models with and without control variables, similar to the Hansen:2016ix model-averaging estimator and coinciding up to a degrees-of-freedom correction with the corresponding Mallows estimator from Hansen:2007cqa. For the specific choice of shrinkage, I contribute three finite-sample properties: First, I note that by averaging over the distribution of controls we obtain dominance of the shrinkage estimator even for low-dimensional target parameters, unlike other available results that require a loss function that is at least three-dimensional. Second, I establish that the resulting estimator remains unbiased under exogeneity of treatment. Third, I conceptualize it as a two-step estimator with a first-stage prediction component. Fourth, I show that it can be seen as a natural, invariant generalized Bayes estimator with respect to a partially improper prior corresponding to uninformativeness in the target parameter.
The linear regression model is set up in Section (ref). Section (ref) proposes the estimator and establishes loss improvement relative to a benchmark OLS estimator provided treatment is exogenous. Section (ref) motivates the estimator as an invariant generalized Bayes estimator (with respect to an improper prior) in a suitably transformed many-means problem.
I consider estimation of the structural parameter $\beta \in \mathbb{R}^k$ in the canonical linear regression model
from $n$ iid observations $(Y_i,X_i,W_i)$, where $X_i \in \mathbb{R}^m$ are the regressors of interest, $W_i \in \mathbb{R}^k$ control variables, and $U_i \in \mathbb{R}$ is homoscedastic, Normal noise. $\alpha$ is an intercept, \footnote{We could alternatively include a constant regressor in $X_i$ and subsume $\alpha$ in $\beta$. I choose to treat $\alpha$ separately since I will focus on the loss in estimating $\beta$, ignoring the performance in recovering the intercept $\alpha$.} and $\gamma$ is a nuisance parameter. To obtain identification of $\beta$ in (ref), I assume that $U_i$ is orthogonal to $X_i$ and $W_i$ (no omitted variables).
Throughout this document, I write upper-case letters for random variables (such as $Y_i$) and lower-case letters for fixed values (such as when I condition on $X_i = x_i$). When I suppress indices, I refer to the associated vector or matrix of observations, e.g. $Y \in \mathbb{R}^n$ is the vector of outcome variables $Y_i$ and $X \in \mathbb{R}^{n \times m}$ is the matrix with rows $X'_i$.
By assumption there are control variables $W$ available with
where $\sigma^2$ need not be known. We care about the (possibly high-dimensional) nuisance parameter $\gamma$ only in so far as it helps us to estimate the (typically low-dimensional) target parameter $\beta$, which is our object of interest.
Given $x \in \mathbb{R}^{n \times m}$ and $w \in \mathbb{R}^{n \times k}$, where we assume that $(\mathbf{1},x,w)$ has full rank $1 + m + k \leq n$, let $q = (q_{\mathbf{1}},q_x,q_w,q_r) \in R^{n \times n}$ orthonormal where $q_{\mathbf{1}} \in \mathbb{R}^n, q_x \in \mathbb{R}^{n \times m}, q_w \in \mathbb{R}^{n \times k}$ such that $\mathbf{1}$ is in the linear subspace of $\mathbb{R}^n$ spanned by $q_{\mathbf{1}} \in \mathbb{R}^{n}$ (that is, $q_{\mathbf{1}} \in \{\mathbf{1}/\sqrt{n},-\mathbf{1}/\sqrt{n}\}$), the columns of $(\mathbf{1},x)$ are in the space spanned by the columns of $(q_{\mathbf{1}},q_x)$, and the columns of $(\mathbf{1},x,w)$ are in the space spanned by the columns of $(q_{\mathbf{1}},q_x,q_w)$. (Such a basis exists, for example, by an iterated singular value decomposition.) Then,
Writing $Y^*_x$, $Y^*_w, Y^*_r$ for the appropriate subvectors of $Y^*$, we find, in particular, that
where $\mu_x = q'_{x} x \beta \in \mathbb{R}^m$, $\mu_w = q'_{w} w \gamma \in \mathbb{R}^k$, and $a = q'_x w (q'_w w)^{-1} \in \mathbb{R}^{m \times k}$. \footnote{Alternatively, we could have denoted by $\mu_x$ the mean of $Y_x^*$. However, by separating out $\mu_x$ from $a \mu_w$ I feel that the role of $\mu_w$ as a relevant nuisance parameter becomes more transparent.} In transforming linear regression to this Normal-means problem, as well as in partitioning the coefficient vector into two groups, for only one of which I will propose shrinkage, I follow Sclove:1968ja.
Conditional on $X{=}x, W{=}w$ and given an estimator $\hat{\mu}_w = \hat{\mu}_w(Y^*_w,Y^*_r)$ of $\mu_w$, a natural estimator of $\mu_x$ is $ \hat{\mu}_x = \hat{\mu}_x(Y^*_x,Y^*_w,Y^*_r) = Y^*_x - a \hat{\mu}_w $. An estimator of $\beta$ is obtained by setting $\hat{\beta} = (q'_x x)^{-1} \hat{\mu}_x$. (The linear least-squares estimator for $\beta$ is obtained from $\hat{\mu}_w = Y^*_w$.) A natural loss function for $\hat{\beta}$ that represents prediction loss units is the weighted loss $(\hat{\beta} - \beta)' (x' q_x q'_x x) (\hat{\beta} - \beta) = \| \hat{\mu}_x - \mu_x \|^2$. We can therefore focus on the (conditional) expected squared-error loss in estimating $\mu_x$, for which we find
with the seminorm $\| v \|_{a'a} = \sqrt{v' a'a v}$ on $\mathbb{R}^k$.
For high-dimensional $\mu_w$ ($k \geq 3$), a natural estimator $\hat{\mu}_w$ with low expected squared-error loss is a shrinkage estimator of the form $\hat{\mu}_w = C Y^*_w$ with scalar $C$, such as the James:1992jm estimator for which $C = 1 - \frac{(k - 2) \| Y^*_r \|^2}{(n - m - k + 1) \| Y^*_w \|^2}$ (or its positive part). While improving with respect to expected squared-error loss ($a'a = \text{const.} \cdot \mathbb{I}_k$), this specific estimator may yield higher (conditional) expected loss in $\mu_x$ when the implied loss function for $\mu_w$ deviates from squared-error loss ($a'a \neq \text{const.} \cdot \mathbb{I}_k$, so the loss function is not invariant under rotations). We will show below that it is still appropriate in the case of independence of treatment and control.
For conditional inference it is known that the least-squares estimator is admissible for estimating $\beta$ provided $m \leq 2$ and inadmissible provided $m \geq 3$ no matter what the dimensionality $k$ of the nuisance parameter $\gamma$ is James:1992jm, as the rank of the loss function is decisive. The above construction does not provide a counter-example to this result: the rank of $a'a =(w' q_w)^{-1} w' q_x q'_x w (q'_w w)^{-1}$ is at most $m$, so for $m \leq 2$, $\hat{\mu}_w = Y^*_w$ remains admissible for the loss function on the right. While we could achieve improvements for $m \geq 3$ -- through shrinkage in $\hat{\mu}_w$ and/or directly in $\hat{\mu}_x$ -- our interest is in the case where $m$ is low and $k$ is high. Conditional on $X{=}x, W{=}w$ we can thus not hope to achieve improvements that hold for any $(\beta,\gamma)$, but we can still hope that shrinkage estimation of $\mu_w$ yields better estimates of $\beta$ on average over draws of the data.
To this end, assume that
(that is, $W_i | X{=}x \stackrel{\text{iid}}{\sim} \mathcal{N}(1 \alpha_W + x_i \beta_W,\Sigma_W)$). Here, $\Sigma_W \in \mathbb{R}^{k \times k}$ is symmetric positive-definite (but not necessarily known). $\alpha_W \in \mathbb{R}^{1 \times k}, \beta_W \in \mathbb{R}^{m \times k}$ describe the conditional expectation of control variables given the regressors $X{=}x$. The case where $x$ and $W$ are orthogonal ($\beta_W = \mathbb{O}_{m \times k}$) and controls $W$ thus not required for identification will play a special role below.
Given $X{=}x$, assume $(q_{\mathbf{1}},q_x)$ is deterministic, and fix $q_{\perp}$ such that $\tilde{q} = (q_{\mathbf{1}},q_x,q_{\perp}) \in \mathbb{R}^{n \times n}$ is orthonormal. Note that
In particular, $q_x' W \rotatebox[origin=c]{90}{$\models$} q_\perp' W$. It follows with
that indeed $q_x' (Y,W) \rotatebox[origin=c]{90}{$\models$} q_\perp' (Y,W)$.
Conditional on $W{=}w$ in the above derivation, $ a \hat{\mu}_w - a \mu_w = q'_x w \hat{\gamma} - q'_x w \gamma $ for $\hat{\gamma} = (q'_w w)^{-1} \hat{\mu}_w$ a function of $q'_\perp w$ and $(Y^*_w,Y^*_r) = (q'_w q_\perp,q'_r q_\perp) (q'_\perp Y)$, so $\hat{\gamma} = \hat{\gamma}(q'_\perp Y,q'_\perp w)$. Assuming measurability, $\hat{\gamma}(q'_\perp Y,q'_\perp W) \rotatebox[origin=c]{90}{$\models$} (q'_x Y, q'_x W)$. Now writing $\hat{\gamma} = \hat{\gamma}(q'_\perp y,q'_\perp w)$ this implies that
with $ \mathop{}\!\textnormal{E}[ W' q_x q'_x W|X{=}x] = \beta'_W x' q_x q'_x x \beta_W + m \Sigma_W $ of full rank $k$. For the expectation of the implied $\hat{\beta}$, we find
We obtain the following characterization of conditional bias and squared-error loss of the implied estimator $\hat{\beta}$:
Note that this lemma does not rely on $n \geq 1 + m + k$, and indeed generalizes to the case $n > 1 + m$ for any $k \geq 1$, including $k > n$.
We consider the special case where treatment is exogenous, and thus $\beta_W = \mathbb{O}_{m \times k}$. This assumption could be justified, for example, in a randomized trial. Note that in this case in addition to the linear least-squares estimator in the “long” regression that includes controls $W$ another natural unbiased (conditional on $X{=}x$) estimator is available, namely the coefficient $(q'_x x)^{-1} q'_x Y$ in the “short” regression without controls. The “long” and “short” regression represent special (edge) cases in the class of two-step estimators introduced above, which are all unbiased in that sense under the exogeneity assumption:
This corollary clarifies that the class of natural estimators derived above are unbiased conditional on $X{=}x$ (but not necessarily on $X{=}x,W{=}w$ jointly), with expected loss equal to the expected out-of-sample prediction loss in a prediction problem where the prediction function $\tilde{w}_0 \mapsto \tilde{w}'_0 \hat{\gamma}$ is trained on $n-1-m$ iid draws, and evaluated on an additional, independent draw $(\tilde{Y}_0,\tilde{W}_0)$ from the same distribution. The “long” and “short” regressions are included as the special cases $\hat{\gamma}(\tilde{w},\tilde{y}) = (\tilde{w}'\tilde{w})^{-1} \tilde{w}' \tilde{y}$ and $\hat{\gamma} \equiv \mathbf{0}_k$, respectively.
The covariates in training and test sample follow the same distribution, which suggests an estimator that is invariant to rotations in the corresponding $k$-means problem. Indeed, the dominating estimator I construct in the following results is of the form
where the standard James:1992jm estimator (for unnknown $\sigma^2$) is recovered at $p = \frac{k - 2}{n - m - k + 1}$.
Note that the result extends to the positive-part analog for which the shrinkage factor is set to zero whenever the expression is negative. For $m=1$, the following dominance is immediate:
The assumption of exogenous treatment is essential for this result, as dropping conditioning on $W$ and restricting interest to $\beta$ would not suffice to break optimality of linear least-squares.
Starting with the transformations in (ref), we consider the decision problem of estimating $\beta$ (equivalently, $\mu_x$). Guided by the treatment of a linear panel-data model in Chamberlain:2009gx, I develop the specific estimator proposed in (ref) as (the empirical Bayes version of) an invariant Bayes estimator with respect to a partially uninformative (improper) Jeffreys prior.
In this section, we condition on $X$ throughout and assume that covariates $W$ are Normally distributed given $X$. Writing $W^*_x = q'_x W, W^*_\perp = q'_\perp W,Y^*_x = q'_x Y, Y^*_\perp = q'_\perp Y$, the transformation developed in (ref) yields the joint distribution
where $\Sigma_W^{1/2}$ is the unique symmetric positive-definite square-root of the symmetric positive-definite matrix $\Sigma_W$, and $V_W$ and $V_Y$ are independent. Here, in addition to $\mu_x = q'_x x \beta$, also $\mu_W = q'_x x \beta_W$, and $s = n - m - 1$. I write $\mathcal{Z} = \mathbb{R}^{m + s} \times \mathbb{R}^{(m + s) \times k}$ for the sample space from which $(Y^*,W^*)$ is drawn according to this $\mathop{}\!\textnormal{P}_\theta$, where I parametrize $\theta = (\mu_x,\gamma) \in \Theta = \mathbb{R}^m \times \mathbb{R}^k$. (I take $\sigma^2,\Sigma_W,\mu_W$ to be constants.)
The action space is $\mathcal{A}= \mathbb{R}^m$, from which an estimate of $\mu_x$ is chosen. As the loss function $L: \Theta \times \mathcal{A} \rightarrow \mathbb{R}$ I take squared-error loss $L(\theta,a) = \| \mu_x - a \|^2$. An estimator $\hat{\beta}: \mathcal{Z} \rightarrow \mathcal{A}$ from the previous section is a feasible decision rule in this decision problem.
For an element $g = (g_\mu,g_x,g_W,g_\perp)$ in the (product) group $G = \mathbb{R}^m \times O(m) \times O(k) \times O(s)$, where $\mathbb{R}^m$ denotes the group of real numbers with addition (neutral element $0$) and $O(k)$ the group of ortho-normal matrices in $\mathbb{R}^{k \times k}$ with matrix multiplication (neutral element $\mathbb{I}_k$), consider the following set of transformations (which are actions of $G$ on $\mathcal{Z}, \Theta, \mathcal{A}$):
For exogenous treatment, these transformations are tied together by leaving model and loss invariant. Indeed, the following is immediate from (ref): \footnote{Alternatively, we could have treated $\mu_W$ as an element of the parameter space and extend the analysis to the case of endogenous treatment. Adding $(g,\mu_W) \mapsto g_x \mu_W \Sigma_W^{-1/2} g'_W \Sigma_W^{1/2}$ to the action on the parameter space would have retained invariance.}
By Proposition (ref), a natural (generalized) Bayes estimator of $\mu_x$ is derived from an improper prior on $\theta$ that is invariant under the action of $G$ on $\Theta$, as this will yield a decision rule $d: \mathcal{Z} \rightarrow \mathcal{A}$ that is invariant in the sense that $d(m_\mathcal{Z}(g,(y,w))) = m_\mathcal{A}(g,d((y,w)))$ for all $(g,(y,w)) \in G \times \mathcal{Z}$. This implies for $\mu_x$ as an improper prior the Haar measure with respect to the translation action (i.e. up to a multiplicative constant the $\sigma$-finite Lebesgue measure on $\mathbb{R}^m$), and for $\gamma$ a prior that is uniform on ellipsoids $\gamma' \Sigma_W \gamma = \omega$. Taking $\frac{\omega}{\tau^2} \sim \chi^2_m$ with some $\tau > 0$ yields the prior $\gamma \sim \mathcal{N}(\mathbf{0},\tau^2 \Sigma_W^{-1})$. With a product prior for $\theta$, the resulting generalized Bayes estimator for $\mu_x$ -- which minimizes posterior loss conditional on the data -- is
Replacing $\Sigma_W$ by the specific sample analog $W_\perp' W_\perp / s$, we obtain the estimator $Y^*_x - \frac{s \tau^2}{s \tau^2 + \sigma^2} W^*_x ((W^*_\perp)' W^*_\perp)^{-1} (W^*_\perp)' Y^*_\perp$. Similarly assuming that $\gamma \sim \mathcal{N}(\mathbf{0},s \tau^2 W_\perp' W_\perp)$, an unbiased estimator of $\frac{s \tau^2}{s \tau^2 + \sigma^2}$ (given $W$) is
This estimator corresponds to the estimator from (ref) at $p = \frac{k - 2}{s - k} = \frac{k - 2}{n - m - k - 1}$. By construction, it retains the invariance of the associated generalized Bayes estimator. This is not specific to this value of $p$:
A natural application of James--Stein shrinkage to control variables in a Normal linear model consistently reduces expected prediction error without introducing bias in the treatment parameter of interest provided treatment is random. In this case, the linear least-squares estimator is thus inadmissible even among unbiased estimators.
In a companion paper JSIV, I show how shrinkage in at least four instrumental variables in a canonical structural form provides consistent bias improvement over the two-stage least-squares estimator. Together, these results suggests different roles of overfitting in control and instrumental variable coefficients, respectively: while overfitting to control variables induces variance, overfitting to instrumental variables in the first stage of a two-stage least-squares procedure induces bias.