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.
72,834 characters · 11 sections · 50 citation commands
Many Proxy Controls
A recent, rapidly growing literature considers the problem of causal inference when a researcher observes only proxies for unobserved confounding factors. For example, one may wish to use test scores to proxy for ability, an unobserved confounder. The proxies are divided into two groups which are independent conditional on the unobserved confounders. One group is a set of negative control treatments: variables that have no direct causal effect on outcomes. The other group of proxies are negative control outcomes: variables that are not directly affected by the treatments.
Compared to standard factor-analytic methods, the proxy control approach has the advantage that the factor structure itself need not be identified. That is, neither the distribution of unobserved factors nor the causal effects of these factors must be identified. Thus proxy control methods may be applied even when the assumptions required for identification of the factor structure do not hold. However, this limits the application of proxy control methods to settings in which the factor structure itself is not of interest.
The proxy control approach is particularly apt for causal inference with high-dimensional data, that is, data that contain many covariates. In these settings one may hope to use the the rich covariates to adjust for confounding. Standard methods treat the covariates as controls, and are valid only if the confounders are non-random once we condition on the covariates. This condition is known as `unconfoundedness', `ignorability', or `selection on observables'. It holds, for example, if all of the confounders are included among the covariates.
By contrast, proxy control methods do not require that the covariates perfectly explain the confounders, only that they are informative about them. Proxy control methods can be valid when the covariates are only noisy proxies for the confounders. Thus proxy control methods allow us to exploit the richness of the covariates without having so assume unconfoundedness. However, to relax this assumption we must place other restrictions on the relationships between the relevant variables.
Nonparametric identification with proxy controls is achieved using the conditional independence assumptions listed in Sub-Figure 1.(a). $Y_i$ is an outcome of interest for an individual $i$ and $X_i$ is a vector of treatments assigned to $i$. $Y_i(\cdot)$ is individual $i$'s treatment response function, or equivalently, $Y_i(x)$ is individual $i$'s potential outcome from the counterfactual treatment level $x$. $W_i$ is a vector of unobserved confounding factors. $D_i$ is a vector of observed characteristics like age and gender. $Z_i$ and $V_i$ are two groups of proxies for $W_i$. $Z_i$ is a vector of negative control treatments and $V_i$ is a vector of negative control outcomes.
Sub-Figure 1.(b) contains two causal diagrams. The diagram on the left, taken from Tchetgen differs from that on the right in that it assumes $Z_i$ is a vector of pre-treatment variables. In the diagram on the right, $Z_i$ is a vector of post-treatment variables, as in Deaner2021. Each diagram is associated with a nonparametric structural equations model (Pearl2009) that implies the conditions in Sub-Figure 1.(a). Note that this is not the only causal structure that implies the conditions in Sub-Figure 1.(a), for example we could allow simultaneous causation between $V_i$ and $W_i$ and between $X_i$ and $Z_i$ (see Deaner2021).
1.(b) imposes that $Z_i$ and $V_i$ have no direct causal effect on each other. However, $Z_i$ and $V_i$ may both depend on $W_i$ and $D_i$. The dependence on $W_i$ means that we can understand these variables as noisy proxies for $W_i$. Correspondingly, condition 1.(a).3 implies that, after accounting for $W_i$ and $D_i$, the variables $V_i$ and $Z_i$ are independent.
the conditions in Figure 1 restrict the manner in which the treatments and outcomes are related to the proxies $Z_i$ and $V_i$. These restrictions are asymmetric. 1.(b) imposes that the negative control treatments $Z_i$ can both cause and be caused by the treatment $X_i$. However, the negative control outcomes $V_i$ cannot cause or be caused by the treatments $X_i$. On the other hand, the diagrams in 1.(b) allow $V_i$ to directly cause the outcomes $Y_i$ whereas $Z_i$ cannot directly cause $Y_i$.
Miao2018a show that if the conditional independence restrictions in Sub-Figure 1.(a) hold along with two statistical completeness assumptions and some regularity conditions, then the average structural function is nonparametrically identified. The average structural function is equal to $E[Y_i(x)]$ in potential outcomes notation or $ E[Y_i|do(x)]$ in the notation of Pearl2009. Deaner2021 shows that under the same conditional independence restrictions, related regularity conditions, and weaker completeness conditions, one can identify the conditional average structural function $E[Y_i(x_1)|X_i=x_2]$.\footnote{Identification of the conditional average structural function implies identification of the average structural function. However, the converse does not hold.} Note that Deaner2021 allows for additional covariates $D_i$ by assuming they are included in both $W_i$, $Z_i$, and $V_i$. Alternatively, one can apply the results in Miao2018a and Deaner2021 within one particular stratum of $D_i$.
Deaner2021 also provides an alternative characterization of the conditional average structural function and provides sufficient conditions for the conditional independence restrictions in the context of panel data.
Nonparametric estimation in this setting was considered by Deaner2019, and later by Tchetgen, Singh2020, Cui2020, and Kallus2021.\footnote{ Deaner2018 earlier considered nonparametric estimation in an equivalent setting that arises in the context of panel data. Deaner2019 was expanded to include a method for nonparametric inference in Deaner2020. } Miao2018 consider estimation in parametric models when a `confounding bridge' function is identified.
In this work we consider identification, estimation, and inference in linear models when the set of proxy controls and the set of confounding factors may be high-dimensional. In particular, we allow for the possibility that the number of proxies and confounders grows with the sample size. Existing work applies to low-dimensional settings in which the number of proxies and confounding factors is treated as fixed.
A key insight in this work is that if there are strictly fewer confounders than there are proxies in $V_i$ and $Z_i$, then a matrix of nuisance parameters has a low rank structure. In addition, a vector of nuisance parameters has a sparse structure. We exploit this low-rank structure and sparsity to reduce the number of free parameters to be estimated. This allows for more efficient estimation, particularly when the number of proxies is large. The number of confounders is generally unknown, so we propose a model selection methods that allows us to adapt to this quantity. The model selection methods are based on techniques from the literature on reduced-rank regression and $\ell_1$-penalized regression.
We present three different estimation methods. The first employs a known bound on the number of confounders, the second selects the number of confounders using the data. These first two methods are based on techniques from the literature on reduced-rank regression and have closed-forms. The third estimator is doubly-robust and is an example of a DML2 estimator of the kind analyzed in section 3.2 of Chernozhukov2018. Chernozhukov2018 shows that DML2 estimators are root-$n$ consistent, asymptotically unbiased, and asymptotically Gaussian, under relatively weak conditions on the nuisance parameter estimates. A disadvantage of the doubly-robust estimator (compared to our alternative methods) is that it does not have a closed-form, however the numerical part of the algorithm can be carried out using any standard Lasso procedure.
Linearity also allows us to weaken the identifying assumptions used in the nonparametric case. In particular, we need only assume variables are uncorrelated rather than independent, we can replace statistical completeness with more intuitive full rank conditions, and we avoid the need for regularity conditions like Assumption A3 in Miao2018.
The linear proxy control model has a long history, dating back to the work of Griliches1977 who considers two scalar proxies for a single confounding factor. However, to the best of our knowledge no existing work exploits the dimension reduction when there are fewer confounders than proxies in $V_i$ or $Z_i$, nor does any existing work allow for a growing number of proxies or confounders.
In sum, our contributions are as follows. We provide a set of identifying assumptions in the linear proxy control model. We present novel estimation methods that allow us to exploit the low rank structure and sparsity in the nuisance parameters when the number of unobserved confounders is less than the number of proxies in $Z_i$ and $V_i$. We develop asymptotic theory for the estimator and an associated inference method, and we provide simulation evidence of the efficacy of our methods.
Suppose that for each individual $i$, the potential outcome $Y_i(x)$ from treatment $x$ is linear in $x$ and differs between individuals only by an additive constant: \[ Y_i(x)=x'\beta_0+U_i \]
$\beta_0$ in the above is a vector of parameters that is the same for all individuals. $U_i$ is a random scalar that represents heterogeneity in potential outcomes. Note that the model implies treatment effects are the same for all individuals. By the definition of potential outcomes, the realized outcome $Y_i$ satisfies $Y_i=Y_i(X_i)$ and so:
Where $X_i$ is the realized treatment and is a column vector of length $d_X$. Let $W_i$ be a length-$d_W$ column vector of latent variables and let $V_i$ be a length-$d_V$ vector of proxies. Let $D_i$ be a length-$d_D$ vector of observed control variables whose first entry is $1$. We assume the following two relationships are linear:
Where $E[\varepsilon_i (W_i',D_i')]=E[\upsilon_i (W_i',D_i')]=0$. We can always find matrices $A_0$, $B_0$, $L_0$, and $R_0$ that satisfy these relationships, however we will place assumptions directly on $\varepsilon_i$ and $\upsilon_i$ and so the linear specifications do incur a loss of generality.
$V_i$ does not directly enter the expression for $U_i$, however this does not rule out the possibility that $V_i$ has a direct causal effect on $Y_i$. We could make this more explicit by including $V_i$ in the expression for $Y_i$ as follows:
Where $E[e_i (W_i',D_i')]=0$. If ((ref)) holds then the model for $Y_i$ above is equivalent to the combination of ((ref)) and ((ref)) with $A_0=F_0+\chi_0 B_0$, $L_0=K_0+\chi_0 R_0$, and $\varepsilon_i = e_i + \chi_0 \upsilon_i$.
We wish to identify and estimate $\beta_0$ in the model for $Y_i(\cdot)$. Identification of $\beta_0$ immediately implies identification of average treatment effects (ATEs): \[ E[Y_i(x_{1})-Y_i(x_{2})]=(x_{1}-x_{2})'\beta_0 \]
And also the conditional average structural function (CASF): \[ E[Y_i(x_{1})|X_i=x_{2}]=(x_{1}-x_{2})'\beta_0+E[Y_i|X_i=x_{2}] \]
Before we state our identifying assumptions let us introduce some additional notation. In particular we define the following objects for each of the variables $H_i=W_i,\, V_i,\, Z_i,\, Y_i,\, X_i$:
$M^{+}$ denotes the Moore-Penrose pseudo-inverse of a matrix $M$. For notational convenience we sometimes write $\tilde{H}_i=\tilde{H}_i(\gamma_{H,0})$ and $\bar{H}_i=\bar{H}_i(\omega_{H,0})$. So for example, $\tilde{X}_i$ is the residual from population linear regression of $X_i$ on $D_i$, or in other words, $X_i$ with $D_i$ partialled out. Similarly, $\bar{Z}_i$ is $Z_i$ with both $D_i$ and $X_i$ partialled out.
\theoremstyle{definition} \newtheorem*{A1.1}{Assumption 1.1 (Model and Exclusion restrictions)} \begin{A1.1} i. ((ref)), ((ref)), and ((ref)) hold. ii. $E[\varepsilon_i X_i]=0$. iii. $E[\varepsilon_i Z_i]=0$, $E[\upsilon_i Z_i]=0$, and $E[\upsilon_i X_i]=0$. \end{A1.1} \theoremstyle{definition}
\newtheorem*{A1.2}{Assumption 1.2 ($V_i$ is sufficiently informative about $W_i$)} \begin{A1.2} $E[{W}_i \bar{V}_i']$ has full row rank. \end{A1.2} \theoremstyle{definition} \newtheorem*{A1.3}{Assumption 1.3 ($Z_i$ is sufficiently informative about $W_i$)} \begin{A1.3} $E[{W}_i \bar{Z}_i']$ has full row rank. \end{A1.3} \newtheorem*{A1.4}{Assumption 1.4 (Full support)} \begin{A1.4} i. $E[\tilde{X}_i \tilde{X}_i']$ is non-singular. ii. $E[\bar{Z}_i \bar{Z}_i']$ and $E[D_i D_i']$ are non-singular. \end{A1.4}
\theoremstyle{definition} Assumptions 1.1.i, 1.1.ii, and 1.1.iii are analogous to conditions 1., 2., and 3., in Sub-Figure 1.(a). Because we impose linearity, we only require zero partial correlations rather than conditional independence.
We can formulate an equivalent set of conditions to Assumptions 1.1.ii and 1.1.iii in terms of the model ((ref)). Recall that $\varepsilon_i = e_i + \chi_0 \upsilon_i$. If $E[\upsilon X_i]=0$ and $E[\upsilon Z_i]=0$ then the remaining conditions in Assumptions 1.1.ii and 1.1.iii hold if and only if $E[e_i X_i]=0$ and $E[e_i Z_i]=0$.
Assumption 1.2 requires that once we have accounted for $X_i$ and $D_i$, $V_i$ is a sufficiently informative proxy for the confounders $W_i$. The assumption replaces the statistical completeness condition on $V_i$ required in the nonparametric setting. Note that the assumption is equivalent to the rank condition for identification in linear instrumental variables (IV) estimation: $W_i$ takes the role of the endogenous regressors, $D_i$ and $X_i$ take the role of the exogenous regressors, and $V_i$ acts as a vector of instruments.
Similarly, Assumption 1.3 requires that after accounting for the treatments and observed controls, $Z_i$ is sufficiently informative about the confounders. Again, this is the same condition required for identification in a linear IV model in which $Z_i$ is a vector of instruments for $W_i$, and the variables $X_i$ and $D_i$ are exogenous regressors.
Note that Assumptions 1.2 and 1.3 can only hold if the vectors $Z_i$ and $V_i$ each have weakly larger dimension than $W_i$.
Assumption 1.4.i is a very mild condition that after partialling out the observed controls $D_i$, the treatments are not perfectly colinear. Assumption 1.4.ii is without loss of generality for the identification results because we can ignore any linearly dependent components of $\bar{Z}_i$ and $D_i$.
\theoremstyle{plain} \newtheorem*{Th1}{Theorem 1} \begin{Th1}
Under Assumptions 1.1-1.4 $\beta_0$ and $d_W$ are identified. In particular, let $d_W\leq r\leq d_V$, then there exists a vector $\beta$ and matrices $A$, $B$, $C$, and $G$ of respective dimensions $1\times r$, $d_{V}\times r$, $r\times d_{Z}$, and $r\times d_{X}$, with $rank(B)\leq r$ so that:
For any such a solution, $\beta=\beta_0$, and $d_W$ is the smallest value of $r$ such that a solution exists. \end{Th1}
Theorem 1 provides a characterization of $\beta_0$ and $d_W$ in terms of a set of moment conditions. The moment conditions involve nuisance parameters $A$, $B$, $C$, and $G$, which are not uniquely determined. One solution to the moment conditions is $\beta=\beta_0$, $A=A_0$, $B=B_0$, $C=C_0$, and $G=G_0$, where $\beta_0$, $A_0$, and $B_0$ are as defined ((ref)), ((ref)), and ((ref)), and $C_0$ and $G_0$ are respectively the population coefficients from linear regression of $\tilde{W}_i$ on $\bar{Z}_i$ and $\tilde{X}_i$.
$\beta_0$ could be estimated directly from the moment conditions in Theorem 1 using the Generalized Method of Moments (GMM) (Hansen1982). However, such an estimation method is computationally challenging. It involves minimizing a GMM objective jointly over $\beta$ and matrices $A$, $B$, $C$, and $G$. Further complicating matters, the moment conditions are non-linear in parameters, the solutions $A$, $B$, $C$, and $G$ are only unique up to non-singular transformations, and (for $r\leq d_V$) the rank constraint on $B$ is non-linear.
As we discuss in Section 3, a more computationally expedient estimator can be attained using an alternative characterization of $\beta_0$ and $d_W$ that is equivalent to that in Theorem 1. The equivalent characterization is captured in Corollary 1 below.
\theoremstyle{plain} \newtheorem*{C1}{Corollary 1} \begin{C1} Under Assumptions 1.1-1.4 $\beta_0$ and $d_W$ are identified. In particular, if and only if $\beta=\beta_0$ then there exists $M\in\mathbb{R}^{d_V\times (d_Z +d_X)}$ and $\xi\in\mathbb{R}^{d_V}$ so that the moment conditions below hold:
The unique solution $M_0$ to ((ref)) has $rank(M_0)=d_W$. Moreover, there exists a $\xi$ with $||\xi||_0\leq d_W$ so that $\beta_0$, $M_0$, and $\xi$ satisfy ((ref)).\footnote{$||v||_0$ is the number of non-zero entries in the vector $v$.} \end{C1}
$\beta_0$ is the unique value of $\beta$ that satisfies the moment conditions in Corollary 1 for some $M$ and $\xi$. There is a unique matrix $M$ that satisfies the moment conditions, and we refer to this solution as $M_0$, and note that $M_0=B_0(C_0,G_0)$. If $d_W\leq d_V$ then there is not a unique choice of $\xi$ with $||\xi_0||_0\leq d_W$ that satisfies the moment conditions. However, it is useful to refer to one particular solution $\xi_0$ which satisfies the moment conditions, we let $\xi_0$ be the (generically unique) solution with minimal $\ell_1$ norm. Note that if the minimal $\ell_1$ solution is unique then $||\xi_0||_0\leq d_W$. $M_0$ and $\xi_0$ are nuisance parameters and are not of direct interest.
Corollary 1 suggest two different means of adapting to the number of confounding factors $d_W$. Firstly, $d_W$ is the rank of $M_0$. Secondly, there is a solution $\beta,M,\xi$ to the moment conditions where $\xi$ has at most $d_W$ non-zero entries. We exploit these results in Subsection 3.2.
To the best of our knowledge, the characterization of $\beta_0$ and the number of confounders $d_W$ using the moment condition in ((ref)) is original, and similarly for ((ref)) and ((ref)).
The characterization of $\beta_0$ is distinct from the characterization using instrumental variables-type moment conditions in Miao2018 and implicit in Griliches1977. As we discuss below, when our assumptions hold and $r<d_V$, our characterization provides additional restrictions that help identify $\beta_0$.
First let us compare with Miao2018. For simplicity let us assume there are no additional controls $D_i$. Miao2018 assume the existence of a function called a `confounding bridge' which then plays a key role in their analysis. A confounding bridge is a function $b$ with the property that for each $x$ in the support of $X_i$, with probability $1$: \[E[Y_i|W_i,X_i=x]=E[b(V_i,x)|W_i,X_i=x]\] Suppose our Assumptions 1.1-1.4 hold and $\epsilon_i$ and $\upsilon_i$ are mean independent of $W_i$ (rather than just uncorrelated with $W_i$), then our model admits a confounding bridge of the form $b(v,x)= \beta_0'x + A_0(B_0' Q B_0)^{-1} B_0' Q v$, where $Q$ is any non-singular matrix and $A_0$ and $B_0$ are as defined in ((ref)) and ((ref)).\footnote{Under Assumptions 1.1-1.4 $B_0$ has full column rank and so $B_0' Q B_0$ is non-singular. See the proof of Theorem 1.}
Miao2018 impose assumptions that imply the confounding bridge is unique and point identified. In our model it may be neither unique nor point identified. In fact, under Assumptions 1.1-1.4 the confounding bridge is generally not unique unless $ d_V = d_W $, otherwise it may depend on the matrix $Q$.\footnote{Under Assumptions 1.1-1.4 $B_0'B_0$ is non-singular. If $B_0'B_0$ is non-singular then $(B_0' Q B_0)^{-1}B_0' Q=(B_0'B_0)^{-1} B_0'$ for all non-singular $Q$ if and only if $B_0$ is a square matrix, i.e., $d_W=d_V$. Strictly speaking, even if $d_W\neq d_V$ the confounding bridge may be unique for certain values of $A_0$ (for example, if $A_0$ is a matrix of zeros).} Even if the confounding bridge is unique, in order to identify the bridge, $Z_i$ must be relevant instruments for $V_i$ after controlling for $X_i$ (see Assumption 5 in Miao2018). Again, under our assumptions this is only possible when $d_V=d_W$.
Applying Miao2018 in our model amounts to using GMM to estimate solutions $\beta$ and $\gamma$ to the following moment condition:\footnote{Miao2018 allow the instruments $(Z_i',X_i')$ to be replaced with any vector of transformations $q(Z_i,X_i)$ with finite variance. However, if $q$ is nonlinear then the resulting moment conditions are valid only when $\epsilon_i$ and $\upsilon_i$ are mean independent of $W_i$ rather than just uncorrelated with $W_i$.}
Under Assumptions 1.1-1.4, the condition above is satisfied when $\beta=\beta_0$ and $\xi=A_0(B_0' Q B_0)^{-1} B_0' Q$ for any non-singular $Q$. When $ d_W < d_V $, the solution $\xi$ is generally not unique, however $\beta_0$ may still be point identified by the moment condition.
((ref)) is a standard instrumental variables (IV) moment condition. Griliches1977 suggests a method for estimation with scalar proxy controls that amounts to performing a standard IV procedure to empirically solve ((ref)).
Corollary 1 imposes additional structure compared to ((ref)). Firstly, Corollary 1 shows that there is a $\xi$ with $d_W$ non-zero entries that satisfies the moment condition. Secondly, if ((ref)) in Corollary 1 holds, then ((ref)) is equivalent to ((ref)). Thus Corollary 1 supplements ((ref)) with an additional set of moment conditions ((ref)). ((ref)) provides additional identifying power when $d_W<\min\{ d_V,d_Z+d_X \} $ (there are fewer confounding factors than proxies). If $d_W$ is unknown, then ((ref)) helps to identify $d_W$, because $d_W$ is rank of the lowest-rank solution $M$ to ((ref)).
Conversely, if $d_W = d_V $ then the additional moment conditions ((ref)) do not help identify $\beta_0$ because there is no reduced-rank restriction on the solution $M$ to ((ref)). Further, if $d_W = d_V $, then the restriction that $\xi$ have only $d_W$ non-zero entries is trivial.
In sum, the results in Theorem 1 and Corollary 1 provide additional identifying power over existing results whenever there are fewer confounding factors than proxies, and also help identify the number of confounding factors. Our results show precisely how the number of latent factors implies sparsity and low-rank structure within the nuisance parameters.
For a given $r$, we could estimate $\beta_0$ by GMM from the moment condition ((ref)) in Theorem 1. However, the moment condition is non-linear in parameters, the objective is generally non-convex, the rank restriction on $B$ is non-linear when $r\leq d_V$, and the solutions $A$, $B$, $C$, and $D$ are (at most) unique up to a normalization. Thus direct minimization of the GMM objective by standard numerical methods may be computationally infeasible, particularly in high dimensions. This problem also applies when $r$ is chosen using model selection (in which case we select an $r$ that estimates $d_W$).
In order to avoid the computational difficulty of joint GMM estimation using ((ref)), we suggest a sequential method of moments estimator (see Newey1994) based on the moment conditions in Corollary 1. We use the first set of moment conditions ((ref)) to estimate the solution $M_0$. We plug the estimate into the moment conditions ((ref)) and use these to estimate the remaining coefficients, including $\beta_0$.
If $d_W\leq \min\{d_V,d_Z+d_X\}$ then the solution $M_0$ to ((ref)) has a reduced-rank structure. Imposing this reduced-rank structure in estimation reduces the number of free nuisance parameters that must be estimated. For estimation of $M_0$, we apply methods from reduced-rank regression. These methods provide closed-form solutions, even when the rank is unknown and must be chosen by model selection. We can also impose a sparse structure on an estimate of $\xi_0$ (recall $\xi=\xi_0$ satisfies the moment conditions in Corollary 1 and $||\xi_0||_0\leq d_W$). To induce sparsity we use $\ell_1$-penalization. In the sub-section on adaptive estimation we propose an estimator with a closed-form that does not induce sparsity in the estimate of $\xi_0$, and also a doubly-robust (but more computationally burdensome) estimator that does impose that $\xi_0$ is sparse.
Let us introduce some notation. Let $n$ be the number of available observations, let $X$ be the matrix with $n$ rows whose $i^{th}$ row is equal to $X_i'$ and similarly for $Y$, $V$, $Z$, and $D$. Recall that the partialled out variables $\tilde{X}_i$, $\tilde{Z}_i$, $\bar{V}_i$, etc., depend on population linear regression parameters. In estimation we must replace these population coefficients with sample analogues. For each $H=V,\, Z,\, Y,\, X$ we define the objects below:
We also let $\hat{H}_i$ be the transpose of the $i^{th}$ column vector of $\hat{H}$ and $\check{H}_i$ the transpose of the $i^{th}$ column vector of $\check{H}$. Thus $\hat{\gamma}_X$ is a sample estimate of $\gamma_{X,0}$, $\hat{X}_i$ is an estimate of $\tilde{X}_i$, and similarly for the other variables.
Corollary 1 states that $M_0$ has rank equal to $d_W$, the number of latent confounding factors. If $r$ upper bounds $d_W$ then $rank(M_0)\leq r$. In this sub-section we consider estimation under this restriction on the rank of $M$ for a fixed choice of $r$. Of particular interest is the case of $d_W$ known, then we can set $r=d_W$..
To estimate $\beta_0$ we first find a rank $\leq r$ matrix $\hat{M}_r$ that approximately solves an empirical analogue of ((ref)) in Corollary 1. In particular, we let $\hat{M}_r$ solve the minimization problem below:
Where $||\cdot||_F^2$ is the squared Frobenius norm (the sum of the squared entries of the matrix).
$\hat{M}_r$ has a closed-form solution described in Reinsel1998 and originally due to Izenman1975. To describe the solution, let $\hat{\Omega}=(\hat{Z},\hat{X} )'(\hat{Z},\hat{X} )$ and define a matrix $\hat{Q}$ given by:
Let $\hat{E}=eigen(\hat{Q})$ be the matrix whose columns are the right eigenvectors of $\hat{Q}$ normalized so that $\hat{E}\hat{E}'$ is the identity and ordered so that the $k^{th}$ column of $\hat{E}$ corresponds to the $k^{th}$ largest eigenvalue. Let $\hat{E}_{[:,1:r]}$ denote the sub-matrix of the first $r$ columns of $\hat{E}$. Then we have:
Having evaluated $\hat{M}_r$ we solve an empirical analogue of ((ref)) with $M$ in the moment condition replaced by $\hat{M}_r$. In particular, our estimate of $\beta_0$ is the vector $\hat{\beta}_r$ that minimizes the least-squares objective below:
Where $||\cdot||$ is the Euclidean norm.
Corollary 1 states that there exists $\xi$ that satisfies the moment conditions with $||\xi||_0\leq r$. We could impose this restriction by adding $||\xi||_0\leq r$ as a constraint in the minimization problem ((ref)). That is, instead of taking the minimum over $\xi\in \mathbb{R}^{d_V}$ we could take the minimum over $\xi\in \mathbb{R}^{d_V}:\,||\xi||_0\leq r$. However, if $\hat{M}_r$ is of rank $r$, then adding this restriction has no effect on the resulting estimator of $\beta_0$. This is because for any vector $\xi$ and rank-$r$ matrix $M$, $M'\xi$ always equals $M' \xi^*$ for some vector $\xi^*$ with only $r$ non-zero components and vice-versa.
The minimizer $\hat{\beta}_r$ in ((ref)) has a closed-form given below. \[ \hat{\beta}_r=(I_{d_{X}},0_{d_X\times d_V})\big(\hat{J}_r'(\hat{Z},\hat{X})'(\hat{Z},\hat{X})\hat{J}_r\big)^{+}\hat{J}_r'(\hat{Z},\hat{X})'\hat{Y} \] $I_{d_{X}}$ is the $d_X\times d_X$ identity matrix, $0_{d_X\times d_V}$ is a $d_X\times d_V$ matrix of zeros, and $\hat{J}_r=\bigg(
,\hat{M}_r\bigg)$. The complete procedure is detailed step-by-step in Algorithm 1 below.
Algorithm (ref) applies for a specific choice of rank restriction $r$ on the estimate of $M_0$. A smaller value of $r$ results in a greater dimension reduction in the nuisance parameters one must estimate. However, if $r< d_W$ (where $d_W$ is the number of unobserved confounders) then the moment conditions are misspecified, i.e., they have no solution. Ideally, $r$ would be chosen to equal $d_W$ which ensures the number of free parameters is minimized and the moment condition correctly specified. However, $d_W$ is generally unknown. We can choose the rank restriction $r$ by model selection and thus adapt to the unknown quantity $d_W$. These methods also result in an estimate of $d_W$, which may be of interest in itself.
We perform model selection by adding penalty terms in estimation so as to induce a low-rank structure on our estimate of $M_0$ and/or a sparse structure on our estimate of $\xi_0$ (recall $\xi_0$ is a choice of $\xi$ that satisfies the moment conditions and has $d_W$ non-zero entries). If we only apply a penalty in the estimation of $M_0$ and use a particular penalty described below, the resulting estimator has a closed-form. If we penalize in both the estimation of $M_0$ and $\xi_0$ this generally results in an estimator without a closed-form solution, but the required numerical minimization can be performed using a standard Lasso estimator (Tibshirani1996). The combination of both penalties results in a doubly-robust estimation method at the cost of additional computation.
In this sub-section we present both an estimator with a closed form that only applies penalization in the estimation of $M_0$, and a doubly-robust estimator without a closed-form that applies penalization in estimation of both $M_0$ and $\xi_0$.
Doubly-robust estimators are insensitive to estimation error in any one of a number of nuisance parameters. They are Neyman Orthogonal/Locally Robust, and generally allow for valid inference using a standard root-$n$ Gaussian approximation under weaker conditions than estimators that are not robust (see for example Chernozhukov2016 and Chernozhukov2018).
Let us first describe the estimator of $M_0$ with an adaptive rank restriction. We replace the first-stage objective ((ref)) with a penalized objective which we minimize over all $d_V\times (d_Z+d_X)$ matrices $M$ (rather than those with rank weakly less than $r$). The penalized objective is:
Where $pen(\cdot)$ is a penalty function and $\lambda_n$ is a scalar penalty parameter that may depend on the data. The penalty function should induce a low-rank structure in the minimizer, for example one could let $pen(M)$ return the nuclear norm of $M$, which is the sum of the singular values of $M$. We focus on the case in which $pen(M)$ simply returns the rank of the matrix $M$. The corresponding objective is given below:
Minimization of the objective above is a penalized reduced-rank regression problem as considered in Bunea2011, and admits a closed-form solution. Bunea2011 show that the solution $\hat{M}$ is a matrix of rank $\hat{r}$, where $\hat{r}$ is the number of eigenvalues of $\hat{Q}$ that exceed $\lambda_n$. Recall that $\hat{Q}$ is defined in ((ref)). $\hat{M}$ is equal to $\hat{M}_r$ defined in ((ref)) but with $r$ set to the data-dependent $\hat{r}$. Note that $\hat{r}$ can be understood as an estimator of the number of unobserved confounders $d_W$.
Having obtained the matrix $\hat{M}$ that minimizes the penalized objective, we may proceed as in the previous subsection to attain an estimate $\hat{\beta}$ of $\beta_0$. The full algorithm is given in detail Algorithm (ref) below. The first five steps provide the estimates $\hat{r}$ and $\hat{M}$. The procedure does not require any numerical optimization nor simulation.
In our simulations we select the penalty parameter $\lambda_n$ by cross-validation. More precisely, let $\{\mathcal{I}_j\}_{j=1}^J$ be a partition of the indices $\{1,...,n\}$. Let $\mathcal{I}_{-j}$ denote all the elements of $\{1,...,n\}$ that are not in $\mathcal{I}_j$. For a given value of $\lambda$ of the penalty $\lambda_n$, let $\hat{M}_{j,\lambda}$ be the result from carrying out steps 2 to 9 in Algorithm 2 using only the data with indices in $\mathcal{I}_{-j}$. We choose $\lambda_n$ to be the value of $\lambda$ that minimizes the cross validation objective below: \[ \sum_{j=1}^J\sum_{i\in \mathcal{I}_j}||\hat{V}_i-\hat{M}_{j,\lambda_n}'(\hat{Z}_i',\hat{X}_i')'||^2 \]
Corollary 1 states that there exists a vector $\xi_0$ with at most $d_W$ non-zero entries so that the moment conditions are satisfied when $\xi=\xi_0$, $\beta=\beta_0$, and $M=M_0$.
The estimator in Algorithm (ref) does not impose sparsity on an estimate of $\xi_0$. We propose a doubly-robust estimator that exploits both the sparsity in $\xi_0$ and the low-rank structure of $M_0$.
Before we describe the procedure as a whole, let us motivate an estimator of $\xi_0$ that induces sparsity. Using ((ref)) to substitute out $M$ from ((ref)) and then partialling out $\tilde{X}_i$ from the resulting equation, we get the following:
Our estimate of $\xi_0$ is indirectly based on the moment condition above. However, it is useful to adapt the moment condition so that $\xi_0$ may take the form of a penalized least-squares estimator. Let $M_{0,[:,1:d_Z]}$ be the sub-matrix of $M_0$ that contains its first $d_Z$ columns. If we substitute the definition of $M_0$ into the condition above and multiply both sides by $M_{0,[:,1:d_Z]}$ then we get the following moment condition:
$\xi_0$ satisfies the condition above if and only if it minimizes the following least squares criterion: \[ E\big[(\bar{Y}_i-\xi_0'M_{0,[:,1:d_Z]}\bar{Z}_i )^2\big] \] To estimate $\xi_0$, we minimize a penalized empirical analogue of the criterion above. In particular, our estimate of $\xi_0$ is the vector $\xi$ that minimizes the empirical objective below:
Where $||\cdot||_1$ is the $\ell_1$ norm and $\delta_n$ is a penalty parameter. $\hat{M}_{d_V,[:,1:d_Z]}$ is the sub-matrix consisting of the first $d_Z$ columns of $\hat{M}_{d_V}$, which is the estimator of $M_0$ defined in Algorithm 1 for $r=d_V$. Note that the estimate $\hat{M}_{d_V}$ does not have a restricted rank.
Minimization of ((ref)) is an $\ell_1$-penalized least squares problem and can be solved using any standard Lasso algorithm. A number of methods are available for selecting the penalty parameter in Lasso regression. For example, $\delta_n$ could be chosen using cross-validation. In our simulations we normalize each of the regressors and the outcomes so that they have mean zero and unit variance and then perform Lasso with the penalty parameter equal to $d_V /n$.
We now turn to the definition of the doubly-robust moment condition. $M_0$ only enters the doubly-robust moment condition indirectly through a parameter $\mu_0$ defined below: \[ \mu_0=E[\tilde{X}_i(\tilde{Z}_i',\tilde{X}_i')]M_0'\big(M_{0,[:,1:d_Z]}E[\bar{Z}_i\bar{Z}_i']M_{0,[:,1:d_Z]}'\big)^{+}M_{0,[:,1:d_Z]} \]
Our estimate $\mu_0$ is a sample analogue of the above. We replace the expectations with sample averages, we replace the partialled out variables with their sample counterparts, and we replace $M_0$ with the adaptive reduced-rank estimate $\hat{M}$ from Algorithm (ref). The estimate $\hat{\mu}$ is given below: \[ \hat{\mu}=\hat{X}'(\hat{Z},\hat{X})\hat{M}'\big(\hat{M}_{[:,1:d_Z]}\check{Z}'\check{Z}\hat{M}_{[:,1:d_Z]}'\big)^{+}\hat{M}_{0,[:,1:d_Z]} \]
In the above $\hat{M}_{[:,1:d_Z]}$ is the sub-matrix that contains the first $d_Z$ columns of $\hat{M}$.
We can now define the doubly-robust moment condition. The condition involves a number of nuisance parameters. These nuisance parameters are $\mu_0$, $\xi_0$, and all the nuisance parameters involved in partialling out $D_i$ and $X_i$. The nuisance parameters involved in partialling out are $\gamma_{Y,0}$, $\gamma_{V,0}$, $\gamma_{X,0}$, $\omega_{Z,0}$, $\omega_{Y,0}$, and $\omega_{V,0}$.
The doubly-robust moment condition is as follows:\footnote{Note that $\gamma_{X,0}$ appears twice as an argument in the moment condition, this is not an error. $\gamma_{X,0}$ enters the score function in two different places, and for analytical purposes it is useful to think of these two occurrences as separate parameters.} \[ E[g_i(\beta_0;\xi_0,\mu_0,\gamma_{Y,0},\gamma_{V,0},\gamma_{X,0},\gamma_{X,0},\omega_{Z,0},\omega_{Y,0},\omega_{V,0})]=0 \] The score function $g_i(\cdot)$ in the moment condition is defined below:
The moment condition is doubly-robust because it holds when we replace a single one of its nuisance parameter arguments with any alternative value. For example, for any vector $\xi$ (not just the correct value $\xi_0$), the moment condition holds: \[ E[g_i(\beta_0;\xi,\mu_0,\gamma_{Y,0},\gamma_{V,0},\gamma_{X,0},\gamma_{X,0},\omega_{Z,0},\omega_{Y,0},\omega_{V,0})]=0 \]
We prove the validity of the moment condition (under Assumptions 1.1-1.4) and establish that it is doubly-robust in Appendix B.
The doubly robust estimate $\hat{\beta}_{DR}$ of $\beta$, sets an empirical expectation of the score function to zero when the nuisance parameters are replaced with first-stage estimates. As we discuss below, it may be advantageous to employ sample-splitting in doubly-robust estimation, but for ease of exposition we begin by defining a doubly robust estimator without sample-splitting.
Recall the definitions of the estimates $\hat{\xi}$, $\hat{\mu}$, $\hat{\gamma}_Y$, $\hat{\gamma}_V$, $\hat{\gamma}_X$, $\hat{\omega}_Z$, $\hat{\omega}_Y$, and $\hat{\omega}_V$. $\hat{\beta}_{DR}$ solves: \[ \frac{1}{n}\sum_{i=1}^n g_i(\hat{\beta}_{DR};\hat{\xi},\hat{\mu},\hat{\gamma}_Y,\hat{\gamma}_V, \hat{\gamma}_X, \hat{\gamma}_X, \hat{\omega}_Z,\hat{\omega}_Y,\hat{\omega}_V)=0 \] We can write the solution succinctly as follows: \[ \hat{\beta}_{DR}=\hat{\Sigma}^{+}\big(\hat{X}' \big(\hat{Y} - \hat{V} \hat{\xi} \big) -\hat{\mu} \check{Z}' (\check{Y} -\check{V} \hat{\xi})\big)/n \] In the above $\hat{\Sigma}=\hat{X}'\hat{X}/n$. Recent work including Chernozhukov2018 and Chernozhukov2016 shows that there may be advantages to sample-splitting in doubly robust and locally robust estimators. Below we describe a version of the estimator above that employs sample splitting.
We partition the data into $J$ sub-samples. In particular, let $\{\mathcal{I}_j\}_{j=1}^J$ be a partition of $\{1,...,n\}$ and let $n_j$ be the number of entries in $\mathcal{I}_j$. Thus each index $i=1,...,n$ is a member of precisely one subset $\mathcal{I}_j$ and $\sum_{j=1}^J n_j =n$. We will use the shorthand $\mathcal{I}_{-j}$ to denote all the elements of $\{1,...,n\}$ that are not in $\mathcal{I}_j$ (i.e., the complement of $\mathcal{I}_j$).
For each $j=1,...,J$ the researcher evaluates each of the nuisance parameter estimates using only the observations with indices in $\mathcal{I}_{-j}$, that is, the data outside of the $j^{th}$ subsample. Thus, for each $j$, one evaluates estimates $\hat{\xi}_j$, $\hat{\mu}_j$, $\hat{\gamma}_{Y,j}$, $\hat{\gamma}_{V,j}$, $\hat{\gamma}_{X,j}$, $\hat{\omega}_{Z,j}$, $\hat{\omega}_{Y,j}$, and $\hat{\omega}_{V,j}$. These estimates are calculated using only data outside the $j^{th}$ subsample, but are otherwise identical to $\hat{\xi}$, $\hat{\mu}$, $\hat{\gamma}_{Y}$, $\hat{\gamma}_{V}$, $\hat{\gamma}_{X}$, $\hat{\omega}_{Z}$, $\hat{\omega}_{Y}$, and $\hat{\omega}_{V}$ respectively.
The estimate $\hat{\beta}_{DR}$ with sample-splitting satisfies the formula:
The solution to the above is as follows: \[ \hat{\beta}_{DR}=\hat{\Sigma}_X^{+}\frac{1}{n}\sum_{j=1}^J\sum_{i\in \mathcal{I}_j}\big(\hat{X}_{i,j} \big(\hat{Y}_{i,j} - \hat{V}_{i,j}' \hat{\xi}_j \big) +\hat{\mu}_j \check{Z}_{i,j} (\check{Y}_{i,j} -\check{V}_{i,j}' \hat{\xi}_j)\big) \] In the above $\hat{\Sigma}_X = \sum_{j=1}^J\sum_{i\in\mathcal{I}_j}\hat{X}_{j,i}\hat{X}_{j,i}'/n$. The full doubly-robust procedure, with sample-splitting, is detailed in Algorithm (ref) in Appendix A.
In some cases the vector of additional covariates $D_i$ may be high-dimensional, that is, there may be many available additional covariates. In addition, we may believe that only a subset of these covariates are linearly predictive of $V_i$, $Z_i$, $X_i$, and $Y_i$. In this case $\gamma_{H,0}$ and/or $\omega_{H,0}$ may be sparse or approximately sparse for some $H\in \{V,Z,X,Y\}$.
To exploit this sparsity or approximate sparsity, we can replace the linear regression estimators $\hat{\gamma}_{H}$ and $\hat{\omega}_{H}$ with Lasso estimators. Let $\eta_{H,n,\gamma}$ and $\eta_{H,n,\omega}$ be a scalar penalty parameter. For $H=V, \, Z,\, X,\, Y$ we can define the following alternative estimators $\hat{\gamma}_{H}$ and $\hat{\omega}_{H}$:
We can then use these estimates in place of the linear regression estimates in the algorithms detailed in this section.
A number of methods exist for choosing the penalty parameters in Lasso regression, for example cross-validation.
The methods in the previous section estimate a low-dimensional parameter of interest $\beta_0$ (and perhaps $d_W$) in the presence of possibly high-dimensional nuisance parameters. We take the standard approach to asymptotic analysis in such settings which is to find conditions under which the estimates are root-$n$ consistent and admit an asymptotic Gaussian approximation. We focus on the doubly-robust estimator with sample-splitting $\hat{\beta}_{DR}$ detailed in the previous section.
The doubly-robust estimator is a Double-Machine Learning 2 (DML2) estimator of the kind analyzed in section 3.2 in Chernozhukov2018. DML2 estimators (along with the DML1 estimators in Chernozhukov2018) have the advantage that they are root-$n$ consistent and asymptotically normal centered at the true parameter under relatively weak conditions on the rates at which the nuisance parameters converge.
Chernozhukov2018 suggest a variance estimator for DML2 estimators. In the case of $\hat{\beta}_{DR}$ this simplifies to: \[ \hat{\sigma}^2=\frac{1}{n}\sum_{j=1}^J\sum_{i\in \mathcal{I}_j}\hat{\Sigma}_X^{+}\hat{g}_i\hat{g}_i '\hat{\Sigma}_X^{+} \] For each $j\in1,...,J$ and $i\in\mathcal{I}_j$, $\hat{g}_i$ is defined as follows: \[ \hat{g}_i= g_i(\hat{\beta}_{DR};\hat{\xi}_j,\hat{\mu}_j,\hat{\gamma}_{Y,j},\hat{\gamma}_{V,j},\hat{\gamma}_{X,j},\hat{\gamma}_{X,j},\hat{\omega}_{Z,j},\hat{\omega}_{Y,j},\hat{\omega}_{V,j}) \] $\hat{g}_i$ is a sample analogue of $g_i$, which is defined by: \[ {g}_i= g_i(\beta_0;\xi_0,\mu_0,\gamma_{Y,0},\gamma_{V,0},\gamma_{X,0},\gamma_{X,0},\omega_{Z,0},\omega_{Y,0},\omega_{V,0}) \] If the variance estimator is consistent and $\hat{\beta}_{DR}$ is asymptotically Gaussian centered at $\beta_0$, then a confidence interval for $l'\beta_0$ (where $l$ is some vector) can be obtained as follows: \[ CI=\big[l'\hat{\beta}\pm \Phi^{-1}(1-\alpha /2)\sqrt{l'\hat{\sigma}^2 l/n}\big] \] The formula above is suggested in Chernozhukov2018. $\Phi$ is the cumulative distribution function of a standard Gaussian random variable.
We now present high-level assumptions for root-$n$ consistency and asymptotic normality of the doubly-robust estimator with sample splitting as defined in ((ref)). Note that our results apply for any choice of estimators for the nuisance parameters not just those specified in Section 3.
Our asymptotic analysis is based on Theorems 3.1 and 3.2 in Chernozhukov2018. The Assumptions 1.1-1.4 and 4.1-4.3 (stated below) act as primitive conditions for the assumptions in that paper.
In order to derive results that are uniform over some parameter space, we suppose that for each sample size $n$, the data generating process, denoted by $P$, belongs to some set $\mathcal{P}_n$. The Assumptions below then restrict $\mathcal{P}_n$.
For notational convenience, for any random column vector $H_i$ we let $\Sigma_H=E[H_i H_i']$, however in the case of $H_i=(X_i',D_i')'$ we write $\Sigma_{XD}$. As in previous sections, if $b$ is a vector then $||b||$ is the Euclidean norm of $b$. If $A$ is a matrix then $||A||$ is the Euclidean matrix norm of $A$, so if $A$ has $d$ rows then $||A||=\sup_{b\in\mathbb{R}^d: ||b||=1}||Ab||$. For a positive semi-definite matrix $A$, $A^{1/2}$ is the unique positive semi-definite matrix $B$ so that $BB=A$ and if $A$ is strictly positive definite then $A^{-1/2}$ is the inverse of $A^{1/2}$.
\theoremstyle{definition} \newtheorem*{A4.1}{Assumption 4.1 (Restrictions on the DGP)} \begin{A4.1} If $P\in\mathcal{P}_n$ the following hold with all scalars on the right-hand sides finite. i. $\Sigma_X$, $\Sigma_{\tilde{X}}$ ,$\Sigma_{\tilde{V}}$, $\Sigma_{\bar{V}}$, $\Sigma_{\bar{Z}}$, $\Sigma_{D}$, and $\Sigma_{XD}$, are finite, and strictly positive definite. $E[\tilde{Y}_{i}^{2}]>0$ and $E[\bar{Y}_{i}^{2}]>0$. ii. $E[{Y}_{i}^{2}]\leq\sigma_{Y}^2$, $||\Sigma_X^{1/2} \beta_{0}||\leq\bar{\beta}$, $||\Sigma_{X}||\leq \sigma_{X}^2$, $||\Sigma_{\tilde{V}}^{1/2} \xi_{0}||\leq\bar{\xi}$, $||\mu_{0}\Sigma_{\bar{Z}}^{1/2}||\leq\bar{\mu}$. iii. With probability 1 for $H=X,V,Y$ $||\Sigma_{\tilde{H}}^{-1/2}E[\tilde{H}_{i}\tilde{H}_{i}'|{D}_{i}]^{1/2}||\leq\bar{\sigma}_{\tilde{H}|D}$, for $H=Z,V,Y$ $ ||\Sigma_{\bar{H}}^{-1/2}E[\bar{H}_i \bar{H}_i '|X_{i},D_i]^{1/2}||\leq\bar{\sigma}_{\bar{H}|XD} $, $||\Sigma_{\tilde{V}}^{-1/2}E[\tilde{V}_{i}\tilde{V}_{i}'|\tilde{X}_{i}]^{1/2}||\leq \bar{\sigma}_{\tilde{V}|\tilde{X}}$, $||\Sigma_{\bar{V}}^{-1/2}E[\bar{V}_{i}\bar{V}_{i}'|\bar{Z}_{i}]^{1/2}||\leq \bar{\sigma}_{\bar{V}|\bar{Z}}$, and $E[\bar{Y}_i^2|\bar{Z}_i]/E[\bar{Y}_i^2]\leq \bar{\sigma}_{\bar{Y}|\bar{Z}}^2$. iv. $E[||\Sigma_D^{-1/2} D_i||^4]\leq d_D s_D^2$ and $E[||\Sigma_{XD}^{-1/2} (X_i',D_i')'||^4]\leq (d_D+d_X) s_{XD}^2$.
\end{A4.1}
\theoremstyle{definition} \newtheorem*{A4.2}{Assumption 4.2 (Convergence rates of the nuisance parameter estimates)} \begin{A4.2} There is a sequence $\alpha_n$ with $\alpha_n\to 0$ so that if $P\in\mathcal{P}_n$ then with probability at least $1-\alpha_n$ the following hold for $j=1,...,J$. i. $||(\hat{\mu}_j-\mu_{0})\Sigma_{\bar{Z}}^{1/2}||\leq\delta_{\mu,n}$. ii. $||\Sigma_V^{1/2}(\hat{\xi}_j-\xi_{0})||\leq\delta_{\xi,n}$. iii. For $H=X,V,Y$, $||\Sigma_{\tilde{H}}^{-1/2}(\hat{\gamma}_{H,j}-\gamma_{H,0})\Sigma_D^{1/2}||\leq\delta_{\gamma,H,n}$. iv. For $H=Z,V,Y$, $ ||\Sigma_{\bar{H}}^{-1/2}(\hat{\omega}_{H,j}-\omega_{H,0})\Sigma_{XD}^{1/2}||\leq\delta_{\omega,H,n}$. \end{A4.2}
\theoremstyle{definition} \newtheorem*{A4.3}{Assumption 4.3 (Higher order moments)} \begin{A4.3} For some $q>2$, if $P\in\mathcal{P}_n$ then the following inequalities hold with the scalars on the right-hand sides finite. i. $E\big[||\tilde{X}_{i}\tilde{V}_{i}'\xi_{0}||^{q}\big]^{1/q}\leq\bar{S}_{q}$, $E\big[||\mu_{0}\bar{Z}_{i}\bar{Y}_{i}||^{q}\big]^{1/q}\leq\bar{S}_{q}$, and $E\big[||\mu_{0}\bar{Z}_{i}\bar{V}_{i}'\xi_{0}||^{q}\big]^{1/q}\leq\bar{S}_{q}$ ii. $E\big[||\Sigma_{D}^{-1/2}D_{i}D_{i}\Sigma_{D}^{-1/2}||^{q}\big]\leq d_{D}^{2}S_{q,D,D}$, for $H=X,Y,V$, $E\big[||\Sigma_{\tilde{X}}^{-1/2}\tilde{X}_{i}\tilde{H}_{i}\Sigma_{\tilde{H}}^{-1/2}||^{q}\big]\leq d_{X}d_{H}S_{q,\tilde{X},\tilde{H}}$ and for $H=X,Y,V$ $E\big[||\Sigma_{D}^{-1/2}D_{i}\tilde{H}_{i}\Sigma_{\tilde{H}}^{-1/2}||^{q}\big]\leq d_{D}d_{H}S_{q,D,\tilde{H}}$. iii. For $H=Z,Y,V$, $E\big[||\Sigma_{\bar{Z}}^{-1/2}\bar{Z}_{i}\bar{H}_{i}\Sigma_{\bar{H}}^{-1/2}||^{q}\big]\leq d_{Z}d_{H}S_{q,\bar{Z},\bar{Z}}$ and $E\big[||\Sigma_{XD}^{-1/2}(X_{i}',D_{i}')'\bar{H}_{i}\Sigma_{\bar{H}}^{-1/2}||^{q}\big]\leq(d_{X}+d_{D})d_{H}S_{q,XD,\bar{Z}}$.
\end{A4.3}
Assumption 4.1 imposes bounds on the magnitudes of some population objects. The bounds must apply for all sample sizes $n$. The bounds given in the assumption do not change with the sample size apart from in the case of 4.1.iv. The upper bound $d_D s_D^2$ in 4.1.iv may grow with the sample size at because $d_D$, the number of additional controls, may grow with $n$, likewise for $(d_D + d_X ) s_{XD}^2$. This flexibility is important because $E[||\Sigma_D^{1/2} D_i||^4]$ is equal to the sum of fourth moments of each of the $d_D$ components of $\Sigma_D^{1/2} D_i$. Suppose the bound on $E[||\Sigma_D^{1/2} D_i||^4]$ were fixed, this would imply that as $d_D$ grows, the average fourth moment of each component of $\Sigma_D^{1/2} D_i$ would have to go to zero. Note that the bounds in 4.1.iii are trivially satisfied under conditional homoskedasticity. In that case each of the upper bounds in 4.1.iii simply equals $1$.
Assumption 4.2 simply imposes convergence rates for each of the nuisance parameter estimates. Note that the convergence rates are required to hold uniformly over sequences of DGPs in $\{\mathcal{P}_n\}_{n=1}^\infty$.
Assumption 4.3 imposes some bounds on higher-order moments. This condition is required to hold for consistency of the variance estimate $\hat{\sigma}$ and thus for asymptotic validity of the confidence intervals.
\theoremstyle{plain} \newtheorem*{Th2}{Theorem 2} \begin{Th2}
Suppose that for each $n$, $P\in\mathcal{P}_{n}$ so that Assumptions 1.1-1.4, Assumptions 4.1, 4.2, and 4.3 all hold, $d_{X}$ is fixed, and all of the eigenvalues of $E[g_{i}g_{i}']$ are bounded below away from zero. Suppose that the nuisance parameter estimates are all consistent. That is, $\delta_{\mu,n}\prec1$, $\delta_{\xi,n}\prec1$, for $H=X,Y,V$ we have $\delta_{\gamma,H,n}\prec1$, and for $H=Z,Y,V$ $\delta_{\omega,H,n}\prec1$. Moreover, suppose that:
Finally suppose that there exits a constant $c$ so that:
Where $q>2$ is the constant in Assumption 4.3. Then uniformly over all $P\in\mathcal{P}_{n}$, $\hat{\beta}$ is root-$n$ consistent and asymptotically normal: \[ \sqrt{n}\sigma^{-1}(\beta_{0}-\hat{\beta})\rightsquigarrow N(0,I) \]
Where the asymptotic variance $\sigma$ is given by: $ \sigma=\Sigma_{X}^{-1}E[g_{i}g_{i}']\Sigma_{X}^{-1} $. Moreover, the variance estimator $\hat{\sigma}$ is consistent for $\sigma$ and the confidence described earlier in this section have asymptotically correct coverage.
\end{Th2}
Theorem 2 establishes uniform root-$n$ consistency of the estimator and asymptotic validity of the confidence intervals. The theorem requires conditions ((ref)), ((ref)), and ((ref)). One set of primitive conditions for ((ref)) and ((ref)) is that $d_D$ grows more slowly than $n$ and all of the nuisance parameter estimates go to zero at a strictly faster rate than $n^{-1/4}$.
The condition ((ref)) directly restricts the rate at which $d_D$, $d_V$, and $d_Z$ may grow with the sample size. The condition is weaker when $q$ is large. Indeed, if the nuisance parameter estimates are consistent then ((ref)) must hold for a sufficiently large value of $q$. However when $q$ is larger then Assumption 4.3 requires the existence of even higher order moments.
In order to assess the efficacy of the methods we present in Section 3 we carry out a Monte Carlo simulation. We implement our methods on a number of simulated datasets. For each simulation, we draw observations independently and identically from the following model:
The residuals $\upsilon_i$, $\epsilon_i$, $\eta_i$, and $e_i$ are drawn independently of each other from zero mean Gaussian distributions: $W_i\sim N(0,I)$, $\upsilon_i\sim N(0,\Sigma_V)$, $\epsilon_i\sim N(0,\Sigma_X)$, $\eta_i\sim N(0,\Sigma_Z)$, and $e_i\sim N(0,\Sigma_Y)$. Note that we do not include additional controls $D_i$ in our simulations.
In each simulation we must choose parameters $\beta_0$, $B_0$, $C_0$, $G_0$, $T_0$, $F_0$, $\chi_0$, $\Sigma_Y$, $\Sigma_V$, $\Sigma_X$, and $\Sigma_Z$. Rather than use a fixed value of each parameter in all of our simulations, we draw the parameters at random in each simulation. Thus our simulation results show the weighted average performance of our estimators over a parameter space.
We draw the parameters as follows. The elements of the coefficient matrices $\beta_0$, $B_0$, $C_0$, $G_0$, $T_0$, $F_0$, and $\chi_0$ are all independently mean-zero normal with variance equal to the square root of the number of columns of the matrix. For example, the elements of $F_0$ are all independent with distribution $N(0,1/\sqrt{d_W})$. This choice of the variances of the normal distributions ensures that the ratio of the variance in each variable to the residual variance remains roughly constant as the dimension changes.
The covariance matrices have a re-scaled inverse Wishart distribution, for example $d_V p\Sigma_V^{-1} \sim W_{d_V}(I,d_V p)$. The natural number $p$ is a hyper-parameter that determines the degrees of freedom of the Wishart distribution.
We are left with hyperparameters $p$, $d_W$, $d_X$, $d_V$, $d_Z$, and the sample size $n$. In all of our simulations we let $d_X=1$ so that there is a single treatment of interest. We set $p=2$ which means the covariance matrices are concentrated around the identity. In all of our simulations $d_Z=d_V$ so there are the same number of proxies in $Z_i$ as in $V_i$. We carry out simulations for a range of choices for the remaining hyperparameters $d_W$, $d_V$, and $n$.
Figure 5.1 shows the mean-squared errors of alternative estimators for a variety of different hyperparameters. The estimators that are compared are: a naive least-squares estimator that simply treats $V_i$ as a set of controls, and three different proxy control estimators. These are the proxy control estimator in Algorithm 1 with $r=d_V$ (i.e, no rank restriction), the computationally expedient adaptive estimator in Algorithm 2, and the doubly-robust estimator with sample-splitting in Algorithm 3.
The proxy control estimator with no rank restriction is equivalent to the two-stage least squares strategy of Griliches1977 in which ${V}_i$ is a vector of endogenous regressors, ${X}_i$ is a vector of exogenous regressors, and ${Z}_i$ is a vector of instruments. In all cases shows in Figure 5.1 this estimator performs worse (in terms of median squared error) than the other two proxy control estimators. The difference is particularly stark when the ratio of the number of proxies to the number of confounders is large.
In the left-most sub-figures in Figure 5.1, the number of proxies in each group is the same as the number of unobserved confounders. Therefore there is no reduced-rank structure or sparsity in the true parameters. Nonetheless, the model selection estimators out-perform the unrestricted estimator. This may be due to the regularizing effect of the reduced-rank regression and $\ell_1$ penalization.
The rank-selection estimator in Algorithm 2 and doubly-robust estimator perform similarly in all simulations, particularly in large samples. The double-robustness of the estimator in Algorithm 3 is intended to reduce the bias, but may come at the cost of increasing the variance.
The naive estimator that treats $V_i$ as a set of controls is necessarily inconsistent, and indeed the mean squared error of this estimator is stable for all but the smallest sample sizes shown in the figures. Nonetheless, the estimator regularly out-performs the unrestricted proxy control estimator in small samples, and it outperforms this estimator even in the largest samples in the most high dimensional case where $d_V=d_Z=100$. The naive method involves estimation of a $d_V$-dimensional nuisance parameter, whereas in the unrestricted proxy control method the nuisance parameters have dimension $d_V(1+d_Z + d_X)$. The larger dimension of the nuisance parameters may lead to high variance of the unrestricted proxy control method, which may dominate the bias in the naive estimator in finite samples.
The naive estimator out-performs all the other methods when $d_W=d_V=d_Z=20$, even in large samples. When $d_W=d_V=d_Z=10$ the naive estimator out-performs the other methods for sample-sizes in the bottom half of the range we examine, but performs worse that all three proxy methods in the largest samples. When the number of proxies in each group exceeds the number of confounders, the adaptive proxy control methods substantially outperformed the naive estimator in all but the very smallest sample sizes we examine.
Figure 5.2 shows the percentage of simulations in which $99\%$, $95\%$, and $90\%$ confidence intervals cover the true parameter $\beta_0$ (recall $\beta_0$ is drawn at random in each simulation). The confidence intervals are those based on a Gaussian approximation for the doubly-robust estimator as described in Section 4.
We see from the left-most sub-figures that the confidence intervals greatly under-cover when the number of proxies in $V_i$ and $Z_i$ equal the number of confounding factors. When there are more proxies than confounders the coverage is close to the desired levels, particularly in large samples.
In Table 1 we give the proportion of simulations in which the rank $r$ selected according to Algorithm 2, is equal to the number of confounders $d_W$ (which is the rank of the matrix $M_0$). In each row we see that the probability of correctly selecting the rank increases with the sample size. For all choices of $d_W$, $d_V$, and $d_Z$ other than those in the first two rows of the table, the algorithm selects the correct rank in at least 87% of simulations when the sample size is $10,000$ or higher, and at least 69% of the time when the sample size is $5000$ or higher.
However, when $d_V=d_Z=d_W=10$ the algorithm selects the incorrect rank in over 80% of simulations, even with $100,000$ observations. When $d_V=d_Z=d_W=20$ the incorrect rank is chosen over 90% of the time. This may explain the under-coverage when $d_V=d_Z=d_W$ as seen in Figure 5.2. In these cases the algorithm regularly selects a rank that is smaller than the true rank (the algorithm cannot select a rank smaller than the minimum of $d_V$ and $d_Z$), and this leads to bias in the resulting coefficient estimates. This bias then leads to under-coverage of the confidence intervals, and because the performance of the rank selection only gradually increases with the sample size, this bias decreases more slowly than the variance, leading to increasingly poor coverage as the sample size grows.
We present novel identification results for the linear model with proxy controls. Our identification results suggest method of moments estimators that can take advantage of the dimension reduction when the number of unobserved confounding factors is smaller than the number of proxies. We present model selection methods that adapt to the unknown number of confounding factors. We provide conditions for uniform root-$n$ consistency of our estimates and asymptotic validity of an inference procedure. Our simulation results suggest that our estimators are more effective than proxy control methods that do not exploit the dimension reduction, particularly when the the number of proxies substantially exceeds the number of unobserved confounders. In the latter case inference based on our doubly-robust adaptive proxy control method performs well.