Doubly Robust Inference in Causal Latent Factor Models

\setcounter{page}{1} \pagestyle{plain} \vskip 80pt \centerline{\bf Doubly Robust Inference in Causal Latent Factor Models}

\advance\leftmargini -.8em

\@startsection{section}{1}{0mm}{-\baselineskip}{0.5\baselineskip}{\center\normalfont\bf}{Introduction}

This article presents a novel framework for the estimation of average treatment effects in modern data-rich environments in the presence of unobserved confounding. We define modern data-rich environments as those featuring many outcome measurements across a wide range of units. Our interest in data-rich environments stems from the emergence of digital platforms (e.g., internet retailers, social media companies, and ride-sharing companies), electronic medical records systems, IoT devices, and other real-time digitized data systems, which gather economic and social behavior data with unprecedented scope and granularity.

Take the example of an internet retailer. The platform collects not only information on purchases of many customers across many products or product categories, but also on glance views, impressions, conversions, engagement metrics, navigation paths, shipping choices, payment methods, returns, reviews, and more. While some variables, such as geo-location and type of device or browser, can be safely treated as pre-determined relative to the platform's treatments (advertisements, discounts, web-page design, etc.), most are outcomes affected by the treatments, latent customer preferences, and unobserved product features. We leverage the availability of many outcome measures in modern data-rich environments to estimate average treatment effects in the presence of unobserved confounding. The core identification concept is that if each element of a high-dimensional outcome vector is influenced by a common low-dimensional vector of unobserved confounders, it becomes possible to remove the influence of the confounders and identify treatment effects.

Two primary approaches to the estimation of treatment effects are outcome-based and assignment-based methods. Consider again the example of an internet-retail platform where customers interact with various product categories. For each consumer-category pair, the platform makes decisions to either offer a discount or not, and records whether the consumer purchased a product in the category. Outcome-based methods operate by imputing the missing potential outcomes for each consumer-product category pair. This process involves predicting whether a consumer, who received a discount, would have made the purchase without the discount (i.e., the potential outcome without discount), and conversely, if a consumer who did not receive the discount would have purchased the product had they received the discount (i.e., the potential outcome with discount). In contrast, assignment-based methods estimate the probabilities of consumers receiving discounts in each product category and adjust for missing potential outcomes by weighting observed outcomes inversely to the probability of missingness.

A substantial body of literature has explored outcome-based methods, particularly in settings where all confounding factors are measured cochran1968effectiveness,rosenbaum1983central,angrist1998estimating,abadie2006large. Imputing potential outcomes in the presence of unobserved confounders poses a more complex challenge. In this context, a commonly adopted framework is the synthetic control method and its variants abadie2003economic,abadie2010synthetic,cattaneo2021prediction,arkhangelsky2021synthetic. An alternative but related approach to outcome imputation under unobserved confounding is the latent factor framework bai2002determining,bai2009panel,xiong2023large, wherein each element of the large-dimensional outcome vector is influenced by the same low-dimensional vector of unobserved confounders. Matrix completion methods chatterjee2015matrix, athey2021matrix, bai2021matrix, dwivedi2022counterfactual, agarwal2023causal which have found widespread applications in recommendation systems and panel data models, are closely related to latent factor models. Similarly, existing assignment-based procedures to estimate average treatment effects rely on the assumption of no unmeasured confounding robins2000marginal, hirano2003efficient, wooldridge2007inverse, common trends restrictions abadie2005, or the availability of an instrumental variable abadie2003semiparametric,tymon2024kappa.

In this article, we propose a doubly-robust estimator robins1994estimation, bang2005doubly, chernozhukov2018double of average treatment effects in the presence of unobserved confounding. This estimator leverages information on both the outcome process and the treatment assignment mechanism under a latent factor framework. It combines outcome imputation and inverse probability weighting with a new cross-fitting approach for matrix completion. We show that the proposed doubly-robust estimator has better finite-sample guarantees than alternative outcome-based and assignment-based estimators. Furthermore, the doubly-robust estimator is approximately Gaussian, asymptotically unbiased, and converges at a parametric rate, under provably valid error rates for matrix completion, irrespective of other properties of the matrix completion algorithm used for estimation.

To our knowledge, this is the first article that leverages latent structures in both the assignment and the outcome processes to obtain a doubly-robust estimator of average treatment effects in the presence of unobserved confounding. arkhangelsky2022doubly study doubly-robust identification with longitudinal data under the assumption that conditioning of a function of the treatment assignments over time (e.g., the fraction of times an individual is exposed to treatment) is enough to remove confounding. athey2021matrix, bai2021matrix, dwivedi2022counterfactual, agarwal2023causal, and xiong2023large propose estimators that apply matrix completion techniques to impute potential outcomes. Although these studies utilize low-rank restrictions in the outcome process, they do not investigate the possibility of similar latent structures in the treatment assignment process. Our article addresses this question, and demonstrate substantial benefits from incorporating knowledge about the structure of the assignment mechanism.

{\bf Terminology and notation.} For any real number $b\in \mathbb{R}$, $\left\lfloor b \right \rfloor$ is the greatest integer less than or equal to $b$. For any positive integer $b$, $[b]$ denotes the set of integers from $1$ to $b$, i.e., $[b] \triangleq \{ 1,\cdots, b \}$. We use $c$ to denote any generic universal constant, whose value may change between instances. For any $c > 0$, $m(c)=\max\{ c, \sqrt{c} \}$ and $\ell_c=\log(2/c)$. For any two deterministic sequences $a_n$ and $b_n$ where $b_n$ is positive, $a_n = O(b_n)$ means that there exist a finite $c > 0$ and a finite $n_0 > 0$ such that $| a_n | \leq c\,b_n$ for all $n \geq n_0$. Similarly, $a_n = o(b_n)$ means that for every $c > 0$, there exists a finite $n_0 > 0$ such that $| a_n | < c\,b_n$ for all $n \geq n_0$. Further, $a_n = \Omega(b_n)$ means that there exist a finite $c > 0$ and a finite $n_0 > 0$ such that $| a_n | \geq c\,b_n$ for all $n \geq n_0$. For a sequence of random variables, $x_n = O_p(1)$ means that the sequence $| x_n |$ is stochastically bounded, i.e., for every $\varepsilon > 0$, there exists a finite $\delta >0$ and a finite $n_0 > 0$ such that $\mathbb{P}\big( | x_n | > \delta \big) < \varepsilon$ for all $n \geq n_0$. Similarly, $x_n = o_p(1)$ means that the sequence $| x_n |$ converges to zero in probability, i.e., for every $\varepsilon > 0$ and $\delta >0$, there exists a finite $n_0 > 0$ such that $\mathbb{P}\big( | x_n | > \delta \big) < \varepsilon$ for all $n \geq n_0$. For sequences of random variables $x_n$ and $b_n$, $x_n = O_p(b_n)$ means $x_n = \wbar{x}_n b_n$ where the sequence $\wbar{x}_n = O_p(1)$. Likewise, $x_n = o_p(b_n)$ means $x_n = \wbar{x}_n b_n$ where the sequence $\wbar{x}_n = o_p(1)$.

A mean-zero random variable $x$ is $\textrm{subGaussian}$ if there exists some $b >0$ such that $\mathbb{E}[\exp(sx)] \leq \exp(b^2s^2/2)$ for all $s\in\mathbb{R}$. Then, the $\textrm{subGaussian}$ norm of $x$ is given by $\snorm{x}_{\psi_2} = \inf\{ t > 0: \mathbb{E}[\exp(x^2/t^2)] \leq 2 \}$. A mean-zero random variable $x$ is $\textrm{subExponential}$ if there exist some $b_1, b_2 >0$ such that $\mathbb{E}[\exp(sx)] \leq \exp(b_1^2s^2/2)$ for all $-1/b_2 < s < 1/b_2$. Then, the $\textrm{subExponential}$ norm of $x$ is given by $\snorm{x}_{\psi_1} = \inf\{ t > 0: \mathbb{E}[\exp(| x |/t)] \leq 2 \}$. $\texttt{Uniform}(a,b)$ denotes the uniform distribution over the interval $[a,b]$ for $a,b \in \mathbb{R}$ such that $a<b$. $\mathcal{N}(\mu,\sigma^2)$ denotes the Gaussian distribution with mean $\mu$ and variance $\sigma^2$.

For a vector $u \in \mathbb{R}^n$, we denote its $t^{th}$ coordinate by $u_t$ and its $2$-norm $\snorm{u}_2$. For a matrix $U \in \mathbb{R}^{n_1 \times n_2}$, we denote the element in $i^{th}$ row and $j^{th}$ column by $u_{i,j}$, the $i^{th}$ row by $U_{i, \cdot}$, the $j^{th}$ column by $U_{\cdot, j}$, the largest eigenvalue by $\lambda_{\max}(U)$, and the smallest by $\lambda_{\min}(U)$. Given a set of indices $\mathcal{R} \subseteq [n_1]$ and $\mathcal{C} \subseteq [n_2]$, $U_{\mathcal{I}} \in \mathbb{R}^{| \mathcal{R} | \times | \mathcal{C} |}$ is a sub-matrix of $U$ corresponding to the entries in $\mathcal{I} \triangleq \mathcal{R} \times \mathcal{C}$, and $U_{-\mathcal{I}}= \{u_{i,j} : (i,j)\in\{[n_1]\times [n_2]\}\setminus \mathcal{I}\}$. Further, we denote the Frobenius norm by $|\!| U | \! |_{\mrm{F}} \triangleq \big( \sum_{i \in [n_1], j \in [n_2]} u_{i,j}^2 \big)^{1/2}$, the $(1,2)$ operator norm by $|\!| U | \! |_{1, 2} \triangleq \max_{j \in [n_2]} \big( \sum_{i \in [n_1]} u_{i,j}^2 \big)^{1/2}$, the $(2,\infty)$ operator norm by $\snorm{U}_{2, \infty}\triangleq \max_{i \in [n_1]} \big( \sum_{j \in [n_2]} u_{i,j}^2 \big)^{1/2}$, and the maximum norm by $|\!| U | \! |_{\max} \triangleq \max_{i \in [n_1], j \in [n_2]} | u_{i,j} |$. Given two matrices $U, V \in \mathbb{R}^{n_1 \times n_2}$, the operators $\odot$ and $\mathrel{ \vphantom{\oplus}\mathpalette\o@generic{{0.4}{\boldsymbol{/}}} }$ denote element-wise multiplication and division, respectively, i.e., $t_{i,j} = u_{i,j} \cdot v_{i,j}$ when $T = U \odot V$, and $t_{i,j} = u_{i,j} / v_{i,j}$ when $T = U \mathrel{ \vphantom{\oplus}\mathpalette\o@generic{{0.4}{\boldsymbol{/}}} } V$. When $V$ is a binary matrix, i.e., $V \in \{ 0,1 \}^{n_1 \times n_2}$, the operator $\otimes$ is defined such that $t_{i,j} = u_{i,j}$ if $v_{i,j} = 1$ and $t_{i,j} = \,?$ if $v_{i,j} = 0$ for $T = U \otimes V$. Given two matrices $U \in \mathbb{R}^{n_1 \times n_2}$ and $V \in \mathbb{R}^{n_1 \times n_3}$, the operator $*$ denotes the (transposed column-wise) Khatri-Rao product of $U$ and $V$, i.e., $T = U * V \in \mathbb{R}^{n_1 \times n_2 n_3}$ such that $t_{i,j} = u_{i,j-n_2 \bar{j}} \cdot v_{i,1+\bar{j}}$ where $\bar{j} = \left\lfloor (j-1)/n_2 \right \rfloor$. For random objects $U$ and $V$, $U \hspace{1mm}{\perp \!\!\! \perp}~ V$ means that $U$ is independent of $V$.

\@startsection{section}{1}{0mm}{-\baselineskip}{0.5\baselineskip}{\center\normalfont\bf}{Setup}

Consider a setting with $N$ units and $M$ measurements per unit. For each unit-measurement pair $i \in [N] \ \text{and}\ j \in [M]$, we observe a treatment assignment $a_{i,j} \in \{ 0,1 \}$ and the value of the outcome $y_{i,j} \in \mathbb{R}$. Although our results can be easily generalized to multi-ary treatments, for the ease of exposition, we focus on binary treatments.

We operate within the Neyman-Rubin potential outcomes framework and denote the potential outcome for unit $i \in [N]$ and measurement $j \in [M]$ under treatment $a \in \{ 0,1 \}$ by $y_{i,j}^{(a)} \in \mathbb{R}$. A no-spillover assumption is implicit in the notation, i.e., the potential outcome $y_{i,j}^{(a)}$ does not depend on the treatment assignment for any other unit-measurement pair. In the context of online retail data, the assumption of no spillovers across measurements is justified if the cross-elasticity of demand across product categories, $j$, is low. Our framework allows for the possibility that the same treatment affects multiple outcomes (e.g., $a_{i,j}=a_{i,j'}$ with probability one, for some $j$ and $j'$ in $[M]$). Realized outcomes, $y_{i,j}$, depend on potential outcomes and treatment assignments,

for all $i \in [N] \ \text{and}\ j \in [M]$. (ref) and (ref) extend the framework proposed in this article to a panel data setting with lagged treatment effects.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Sources of stochastic variation}

In the setup of this article, each unit $i \in [N]$ is characterized by a set of unknown parameters, $\{(\theta_{i,j}^{(0)},\theta_{i,j}^{(1)}, p_{i,j})\in \mathbb R^2\times [0, 1] \}_{j\in [M]}$, which we treat as fixed. Potential outcomes and treatment assignments are generated as follows: for all $i \in [N], j \in [M]$, and $a \in \{ 0,1 \}$,

where $\varepsilon_{i,j}^{(a)}$ and $\eta_{i,j}$ are mean-zero random variables, and

It follows that $\theta_{i,j}^{(a)}$ is the mean of the potential outcome $y_{i,j}^{(a)}$, and $p_{i,j}$ is the unknown assignment probability or latent propensity score. The matrices ${\Theta}^{(0)} \triangleq \{ \theta_{i,j}^{(0)} \}_{i \in [N], j \in [M]}$, ${\Theta}^{(1)} \triangleq \{ \theta_{i,j}^{(1)} \}_{i \in [N], j \in [M]}$, and ${P} \triangleq \{ p_{i,j} \}_{i \in [N], j \in [M]}$ collect mean potential outcomes and assignment probabilities. Then, the matrices ${E}^{(0)}\triangleq \{ \varepsilon_{i,j}^{(0)} \}_{i \in [N], j \in [M]}, {E}^{(1)}\triangleq \{ \varepsilon_{i,j}^{(1)} \}_{i \in [N], j \in [M]}$, and $ {W}\triangleq \{ \eta_{i,j} \}_{i \in [N], j \in [M]}$ capture all sources of randomness in potential outcomes and treatment assignments.

Our setup allows ${\Theta}^{(0)}, {\Theta}^{(1)}$ to be arbitrarily associated with ${P}$, inducing unobserved confounding. The assumptions in (ref) imply that ${\Theta}^{(0)}, {\Theta}^{(1)}$, and ${P}$ include all confounding factors, and require $(\varepsilon^{(0)}_{i, j}, \varepsilon^{(1)}_{i, j}) \hspace{1mm}{\perp \!\!\! \perp}~ \eta_{i, j}$ for every $i \in [N]$ and $j \in [M]$.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Target causal estimand}

For any given measurement $j \in [M]$, we aim to estimate the effect of the treatment averaged over all units,

$\mathrm{ATE}_{\cdot,j}$ akin to the conditional average treatment effect of abadie2006large, but based on the latent means, $\theta_{i,j}^{(a)}$, in (ref) rather than on conditional means that depend on observed covariates only. It is straightforward to adapt the methods in this article to the estimation of alternative parameters, like the average treatment effect across measurements for each unit $i$, or the estimation of treatment effects over a subset of the units, $S \subset [N]$.

\@startsection{section}{1}{0mm}{-\baselineskip}{0.5\baselineskip}{\center\normalfont\bf}{Estimation}

In this section, we propose a procedure that uses the treatment assignment matrix ${A}$ and the observed outcomes matrix ${Y}$ to estimate $\mathrm{ATE}_{\cdot,j}$, where

The estimator proposed in this section leverages matrix completion as a key subroutine. We start the section with a brief overview of matrix completion methods. \@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Matrix completion: A primer} Consider a matrix of parameters $T \in \mathbb{R}^{N \times M}$. While $T$ is unobserved, we observe the matrix $S \in \{\mathbb{R} \cup \{\,?\}\}^{N \times M}$ where $\,?$ denotes a missing value. The relationship between $S$ and $T$ is given by

Here, $H \in \mathbb{R}^{N \times M}$ is a noise matrix, and ${F} \in \{0, 1\}^{N \times M}$ is a masking matrix with ones for the recorded entries of $S$ and zeros for the missing entries.

A matrix completion algorithm, denoted by $\texttt{MC}$, takes the $S$ as its input, and returns an estimate of $T$, which we denote by $\what {T}$ or $\texttt{MC}(S)$. In other words, $\texttt{MC}$ produces an estimate of a matrix from noisy observations of a subset of all the elements of the matrix.

The matrix completion literature is rich with algorithms $\texttt{MC}$ that provide error guarantees, namely bounds on $\snorm{\texttt{MC}(S)-T}$ for a suitably chosen norm/metric $\snorm{\cdot}$, under a variety of assumptions on the triplet $(T, H, {F})$. Typical assumptions are $(i)$ $T$ is low-rank, $(ii)$ the entries of $H$ are independent, mean-zero and sub-Gaussian random variables, and $(iii)$ the entries of ${F}$ are independent Bernoulli random variables. Though matrix completion is commonly associated with the imputation of missing values, a typically underappreciated aspect is that it also denoises the observed matrix. Even when each entry of $S$ is observed, $\texttt{MC}(S)$ subtracts the effects of $H$ from $S$, i.e., it performs matrix denoising. nguyen2019low provide a survey of various matrix completion algorithms.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Key building blocks}

We now define and express matrices that are related to the quantities of interest ${\Theta}^{(0)}, {\Theta}^{(1)}$, and ${P}$ in a form similar to (ref). See Figure (ref) for a visual representation of these matrices.

Building on the earlier discussion, the application of matrix completion yields the following estimates:

where the algorithm $\texttt{MC}$ may vary for ${\what{\Theta}}^{(0)}$, ${\what{\Theta}}^{(1)}$, and ${\what{P}}$. Because all entries of ${A}$ are observed, $\texttt{MC}({A})$ denoises ${A}$ but does not need to impute missing entries. From (ref) and (ref), it follows that ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$ depend on ${A}$ and ${Y}$, whereas ${\what{P}}$ depends only on ${A}$.

In this section, we deliberately leave the matrix completion algorithm $\texttt{MC}$ as a “black-box”. In (ref), we establish finite-sample and asymptotic guarantees for our proposed estimator, contingent on specific properties for $\texttt{MC}$. In (ref), we propose a novel end-to-end matrix completion algorithm that satifies these properties.

Given matrix completion estimates of $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)}, {\what{P}})$, we formulate two preliminary estimators for $\mathrm{ATE}_{\cdot,j}$: $(i)$ an outcome imputation estimator, which uses ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$ only, and $(ii)$ an inverse probability weighting estimator, which uses ${\what{P}}$ only. Then, we combine these to obtain a doubly-robust estimator of $\mathrm{ATE}_{\cdot,j}$.

{\bf Outcome imputation (OI) estimator.} Let $\what{\theta}_{i,j}^{(a)}$ denote the $(i,j)$-th entry of ${\what{\Theta}}^{(a)}$ for $i \in [N], j \in [M]$, and $a \in \{0, 1\}$. The OI estimator for $\mathrm{ATE}_{\cdot,j}$ is defined as follows:

That is, the OI estimator is obtained by taking the difference of the average value of the $j$-th column of the estimates ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$. The quality of the OI estimator depends on how well ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$ approximate the mean potential outcome matrices ${\Theta}^{(0)}$ and ${\Theta}^{(1)}$, respectively.

{\bf Inverse probability weighting (IPW) estimator.} Let $\what{p}_{i,j}$ denote the $(i,j)$-th entry of ${\what{P}}$ for $i \in [N]$ and $j \in [M]$. The IPW estimate for $\mathrm{ATE}_{\cdot,j}$ is defined as follows:

That is, the IPW estimator is obtained by taking the difference of the average value of the $j$-th column of the matrices ${{Y}}^{(0),\mathrm{obs}}_{}$ and ${{Y}}^{(1),\mathrm{obs}}_{}$, replacing unobserved entries with zeros, and weighting each outcome by the inverse of the estimated assignment probability to account for confounding. The quality of the IPW estimate depends on how well ${\what{P}}$ approximates the probability matrix ${P}$.

The matrix completion-based OI and IPW estimators in (ref) and (ref) have the same form as the classical OI and IPW estimators, which are derived for settings where all confounders are observed imbens_rubin_2015. In contrast to the classical setting, our framework is one with unmeasured confounding.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Doubly-robust (DR) estimator}

The DR estimator of $\mathrm{ATE}_{\cdot,j}$ combines the estimates ${\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)}$, and ${\what{P}}$ from (ref). It is defined as follows:

In (ref), we prove that $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$ consistently estimates $\mathrm{ATE}_{\cdot,j}$ as long as either $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)})$ is consistent for $({\Theta}^{(0)}, {\Theta}^{(1)})$ or ${\what{P}}$ is consistent for ${P}$, i.e., it is doubly-robust. Furthermore, we show that the DR estimator provides superior finite sample guarantees than the OI and IPW estimators, and that it satisfies a central limit theorem at a parametric rate under weak conditions on the convergence rate of the matrix completion routine. Using simulated data, (ref) demonstrates the improved performance of DR, relative to OI and IPW. Despite substantial biases observed in both OI and IPW estimates, the error of the DR estimate closely follows a mean-zero Gaussian distribution. We provide a detailed description of the simulation setup in (ref).

\@startsection{section}{1}{0mm}{-\baselineskip}{0.5\baselineskip}{\center\normalfont\bf}{Main Results}

This section presents the formal results of the article. (ref) details assumptions, (ref) discusses finite-sample guarantees, and (ref) presents a central limit theorem for $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Assumptions}

{\bf Requirements on data generating process.} We make two assumptions on how the data is generated. First, we impose a positivity condition on the assignment probabilities.

(ref) requires that the propensity score for each unit-outcome pair is bounded away from $0$ and $1$, implying that any unit-item pair can be assigned either of the two treatments. An analogous assumption is pervasive in causal inference models with no-unmeasured confounding. For simplicity of exposition and to avoid notational clutter, (ref) requires (ref) for all outcomes, $j\in [M]$. In practical applications, however, $\mathrm{ATE}_{\cdot,j}$ may be estimated for a select group of those outcomes. In that case, the positivity assumption applies only for the selected subset of outcomes for which $\mathrm{ATE}_{\cdot,j}$ is estimated.

Next, we formalize the requirements on the noise variables.

(ref)(ref) defines $({\Theta}^{(0)}, {\Theta}^{(1)}, {P})$ as matrices collecting the means of the potential outcomes and treatment assignments in (ref). Further, for every measurement, it imposes independence across units in the noise variables. (ref)(ref) imposes independence between the noise in the potential outcomes and noise in treatment assignment, and implies that for each particular unit $i$ and measurement $j$, confounding emerges only from the interplay between $(\theta_{i,j}^{(0)}, \theta_{i,j}^{(1)})$ and $p_{i,j}$. Finally, (ref)(ref) is mild and useful to derive finite-sample guarantees. For the central limit theorem in (ref), subGaussianity could be disposed of by restricting the moments of $\varepsilon_{i,j}^{(a)}$. (ref) does not restrict the dependence between $\varepsilon_{i,j}^{(0)}$ and $\varepsilon_{i,j}^{(1)}$. Neither (ref) restricts the dependence of $\eta_{i,j}$ across outcomes. In particular, (ref) allows for the existence of pairs of outcomes $(j,j')$ such that $\mathbb{E}[\eta_{i,j}^2]=\mathbb{E}[\eta_{i,j'}^2]=\mathbb{E}[\eta_{i,j}\eta_{i,j'}]$, in which case $a_{i,j}=a_{i,j'}$ with probability one.

{\bf Requirements on matrix completion estimators.} First, we assume the estimate ${\what{P}}$ is consistent with (ref).

(ref) holds when the entries of ${\what{P}}$ are truncated to the range $[{\bar\lambda}, 1-{\bar\lambda}]$, provided ${\bar\lambda}$ is not greater than $\lambda$. Second, our theoretical analysis requires independence between certain elements of the estimates $({\what{P}}, {\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)})$ from (ref), and the noise matrices $({W}, {E}^{(0)},{E}^{(1)})$. We formally state this independence condition as an assumption below.

(ref) requires that within each of the two partitions of the units, estimated mean potential outcomes and estimated assignment probabilities are jointly independent of the error in assignment probabilities, for every measurement. Similarly, (ref) requires that within each of the two partitions of the units, estimated assignment probabilities are independent jointly of the noise in assignment probabilities and potential outcomes, for every measurement. Conditions like (ref) and (ref) are familiar in the doubly-robust estimation literature. chernozhukov2018double employ a cross-fitting device to enforce an assumption similar to (ref) in a context with no unmeasured confounders. (ref) provides a novel cross-fitting procedure for matrix estimation under which (ref) holds for any $\texttt{MC}$ algorithm (under additional assumptions on the noise variables).

{\bf Matrix completion error rates.} The formal guarantees in this section depend on the normalized $(1,2)$-norms of the errors in estimating the unknown parameters $({\Theta}^{(0)}, {\Theta}^{(1)}, {P})$. We use the following notation for these errors:

A variety of matrix completion algorithms deliver $\mc{E}\big( {\what{P}} \big)=O_p(\min\{ N,M \}^{-\alpha})$ and $\mc{E}\big( {\what{\Theta}}^{} \big)=O_p(\min\{ N,M \}^{-\beta})$, where $0<\alpha, \beta \leq 1/2$. The conditions in this section track dependence on $N$ only. We say that the normalized errors $\mc{E}\big( {\what{P}} \big)$ and $\mc{E}\big( {\what{\Theta}}^{} \big)$ achieve the parametric rate when they have the same rate as $O_p(N^{-1/2})$. (ref) explicitly characterizes how the rates of convergence $\mc{E}\big( {\what{P}} \big)$ and $\mc{E}\big( {\what{\Theta}}^{} \big)$ depend on $N$ and $M$ for a particular matrix completion algorithm based on bai2021matrix.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Non-asymptotic guarantees}

The first main result of this section provides a non-asymptotic error bound for $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}- \mathrm{ATE}_{\cdot,j}$ in terms of the errors $\mc{E}\big( {\what{P}} \big)$ and $\mc{E}\big( {\what{\Theta}}^{} \big)$ defined in (ref).

The proof of (ref) is given in (ref). (ref) bound the absolute error of the DR estimator by {the rate of $\mc{E}\big( {\what{\Theta}}^{} \big) (\mc{E}\big( {\what{P}} \big) + N^{-0.5}) + N^{-0.5}$.} When $\mc{E}\big( {\what{P}} \big)$ is lower bounded at the parametric rate of $N^{-0.5}$, $\mrm{Err}^{\mathrm{DR}}_{N, \delta}$ has the same rate as $\mc{E}\big( {\what{P}} \big) \mc{E}\big( {\what{\Theta}}^{} \big) + N^{-0.5}$.

{\bf Doubly-robust behavior of \texorpdfstring{$\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$}.} The error rate of $\mc{E}\big( {\what{P}} \big) \mc{E}\big( {\what{\Theta}}^{} \big) + N^{-0.5}$ immediately reveals that the DR estimate is doubly-robust with respect to the error in estimating the mean potential outcomes $({\Theta}^{(0)},{\Theta}^{(1)})$ and the assignment probabilities ${P}$. First, the error $\mrm{Err}^{\mathrm{DR}}_{N, \delta}$ decays at a parametric rate of $O_p(N^{-0.5})$ as long as the product of error rates, $\mc{E}\big( {\what{P}} \big) \mc{E}\big( {\what{\Theta}}^{} \big)$, decays as $O_p(N^{-0.5})$. As a result, $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$ can exhibit a parametric error rate even when neither the mean potential outcomes nor the assignment probabilities are estimated at a parametric rate. Second, $\mrm{Err}^{\mathrm{DR}}_{N, \delta}$ decays to zero as long as either of $\mc{E}\big( {\what{P}} \big)$ or $\mc{E}\big( {\what{\Theta}}^{} \big)$ decays to zero, provided both errors are $O_p(1)$.

We next compare the performance of DR estimator with the OI and IPW estimators from (ref), respectively. Towards this goal, we characterize the $\mathrm{ATE}_{\cdot,j}$ estimation error of $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{OI}}$ in terms of $\mc{E}\big( {\what{\Theta}}^{} \big)$ and of $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{IPW}}$ in terms of $\mc{E}\big( {\what{P}} \big)$.

The proofs of (ref) and (ref) are given in (ref), respectively. (ref) implies that \ifxarxivarxiv in an asymptotic sequence \fi with bounded $\theta_{\max}$, OI and IPW attain the parametric rate $O_p(N^{-0.5})$ provided $\mc{E}\big( {\what{\Theta}}^{} \big)$ and $\mc{E}\big( {\what{P}} \big)$ are $O_p(N^{-0.5})$, respectively. The next corollary, proven in (ref), compares these error rates with those obtained for the DR estimator in (ref).

(ref) shows that the DR estimate's error decay rate is consistently superior to that of the OI and IPW estimates across a variety of regimes for $\alpha, \beta$. Specifically, the error $\mrm{Err}^{\mathrm{DR}}_{N, \delta}$ scales strictly faster than both $\mrm{Err}^{\mathrm{OI}}_N$ and $\mrm{Err}^{\mathrm{IPW}}_{N, \delta}$ if the estimation errors of ${\what{\Theta}}^{(0)}$, ${\what{\Theta}}^{(1)}$, and ${\what{P}}$ converge slower than at the parametric rate $O_p(N^{-1/2})$. When the estimation errors of ${\what{\Theta}}^{(0)}$, ${\what{\Theta}}^{(1)}$, and ${\what{P}}$ all decay at a parametric rate, OI, IPW, and DR estimation errors decay also at a parametric rate.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Asymptotic guarantees}

The next result, proven in (ref) as a corollary of (ref), provides conditions on $\mc{E}\big( {\what{P}} \big)$ and $\mc{E}\big( {\what{\Theta}}^{} \big)$ for consistency of $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$.

(ref) states that $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$ is a consistent estimator for $\mathrm{ATE}_{\cdot,j}$ as long as either the mean potential outcomes or the assignment probabilities are estimated consistently.

The next theorem, proven in (ref), establishes a Gaussian approximation for $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$ under mild conditions on error rates $\mc{E}\big( {\what{P}} \big)$ and $\mc{E}\big( {\what{\Theta}}^{} \big)$.

(ref) describes two simple requirements on the estimated matrices ${\what{P}}$ and $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)})$, under which $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$ exhibits an asymptotic Gaussian distribution centered at $\mathrm{ATE}_{\cdot,j}$. Condition (ref) requires that the estimation errors of ${\what{P}}$ and $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)})$ converge to zero in probability. Condition (ref) requires that the product of the errors decays sufficiently fast, at a rate $o_p\big( N^{-1/2} \big)$, ensuring that the bias of the normalized estimator in (ref) converges to zero. Condition (ref) is similar to conditions in the literature on doubly-robust estimation of average treatment effects under observed confounding chernozhukov2018double. Specifically, in that context, chernozhukov2018double assume that the product of propensity estimation error and outcome regression error decays faster than $N^{-1/2}$.

{\bf Black-box asymptotic normality.} We emphasize that (ref) applies to any matrix completion algorithm $\texttt{MC}$, provided conditions (ref) and (ref) hold. This level of generality is useful because the product of $\mc{E}\big( {\what{P}} \big)$ and $\mc{E}\big( {\what{\Theta}}^{} \big)$ is $o_p\big( N^{-1/2} \big)$ for a wide range of $\texttt{MC}$ algorithms, under mild assumptions on $({\Theta}^{(0)},{\Theta}^{(1)}, {P})$. In contrast, achieving such black-box asymptotic normality for OI or IPW estimates is challenging. Their biases are tied to the individual error rates, $\mc{E}\big( {\what{\Theta}}^{} \big)$ and $\mc{E}\big( {\what{P}} \big)$, which are typically lower-bouded at the parametric rate of $N^{-0.5}$.

The next result, proven in (ref), provides a consistent estimator for the asymptotic variance $\wbar{\sigma}_j^2$ from (ref).

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Application to panel data with lagged treatment effects} (ref) considered a model where the outcome $y_{i,j}$ for unit $i \in [N]$ and measurement $j \in [M]$ depends on treatment assignment only for unit $i$ and measurement $j$, i.e., $a_{i,j}$. (ref) discusses how to extend the results of this section to a setting of panel data with lagged treatment effects. In a panel data setting, the $M$ measurements correspond to $T$ time periods, and $t$ denotes the time index. Then, (ref) considers an auto-regressive setting, where the potential outcomes at time $t$ depends on the treatment assignment at time $t$ and the realized outcome at time $t-1$, i.e., for all $i \in [N], t \in [T]$, and $a \in \{ 0,1 \}$,

and observed outcomes satisfy

The presence of lagged treatment effects in this model makes it crucial to define causal estimands for entire sequences of treatments. (ref) describes how the proposed doubly-robust estimation can be extended to treatment sequences and derives a generalization of (ref).

\@startsection{section}{1}{0mm}{-\baselineskip}{0.5\baselineskip}{\center\normalfont\bf}{Matrix Completion with Cross-Fitting}

In this section, we introduce a novel algorithm designed to construct estimates $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)}, {\what{P}})$ that adhere to (ref) and satisfy conditions (ref) in (ref). We first explain why traditional matrix completion algorithms fail to deliver the properties required by (ref). We then present $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$, a meta-algorithm that takes any matrix completion algorithm and uses it to construct $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)}, {\what{P}})$ that satisfy (ref), and the stronger independence condition in (ref). Finally, we describe $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$, an end-to-end algorithm obtained by combining $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ with the singular value decomposition (SVD)-based algorithm of bai2021matrix, and establish that it also satisfies conditions (ref) in (ref).

{\bf Traditional matrix completion.} Estimates $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)}, {\what{P}})$ obtained from existing matrix completion algorithms need not satisfy (ref). In particular, using the entire assignment matrix ${A}$ to estimate each element of ${P}$ typically results in a violation of $\big\{ \what{p}_{i, j} \big \}_{i\in\mathcal{R}_s} \hspace{1mm}{\perp \!\!\! \perp}~ \big\{ \eta_{i, j} \big \}_{i\in\mathcal{R}_s}$ in (ref), as each entry of $\what{P}$ is allowed to depend on the entire noise matrix ${W}$. For example, in spectral methods nguyen2019low, ${\what{P}}$ is a function of the SVD of the entire matrix ${A}$, and

for all $(i, j), (i', j') \in [N] \times [M]$ in general, which implies $\big\{ \what{p}_{i, j} \big \}_{i\in\mathcal{R}_s} ~{\not\!\perp\!\!\!\perp}~ \big\{ \eta_{i, j} \big \}_{i\in\mathcal{R}_s}$, for every $\mathcal{R}_s \subset [N]$. Similarly, in matching methods such as nearest neighbors li2019nearest, ${\what{P}}$ is a function of the matches/neighbors estimated from the entire matrix ${A}$. Dependence structures such as $\what{p}_{i,j} ~{\not\!\perp\!\!\!\perp}~ a_{i,j}$ for any $i, j \in [N] \times [M]$---which is weaker than (ref)---are enough to violate the $\big\{ \what{p}_{i, j} \big \}_{i\in\mathcal{R}_s} \hspace{1mm}{\perp \!\!\! \perp}~ \big\{ \eta_{i, j} \big \}_{i\in\mathcal{R}_s}$ requirement in (ref). Likewise, the requirement $\big\{ \what{\theta}_{i,j}^{(a)} \big \}_{i\in\mathcal{R}_s} \hspace{1mm}{\perp \!\!\! \perp}~ \big\{ \eta_{i, j} \big \}_{i\in\mathcal{R}_s}$ in (ref) can be violated, because ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$ depend respectively on ${{Y}}^{(0),\mathrm{obs}}_{}$ and ${{Y}}^{(1),\mathrm{obs}}_{}$, which themselves depend on the entire matrix ${A}$.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{\texorpdfstring{Cross\allowbreak-\allowbreakFitted\allowbreak-\allowbreakMC}: A meta-cross-fitting algorithm for matrix completion}

We now introduce $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$, a cross-fitting procedure that modifies any $\texttt{MC}$ algorithm to produce $({\what{\Theta}}^{(0)}, {\what{\Theta}}^{(1)}, {\what{P}})$ that satisfy (ref). We employ the following assumption on the noise variables.

For a given block $\mathcal{I}$, (ref) requires the noise in the treatment assignments corresponding to $\mathcal{I}$ to be independent jointly of the noise in the treatment assignments and the potential outcomes corresponding to the remaining three blocks. Likewise, (ref) requires the noise in the treatment assignments corresponding to the remaining three blocks to be independent jointly of the noise in the treatment assignments and the potential outcomes corresponding to $\mathcal{I}$. (ref) leaves unrestricted the dependence of the noise variables across outcomes that belong to the same block.

For notational simplicity, (ref) imposes independence conditions across blocks of outcomes in a partition of $[M]$ into two blocks only. It is important to note, however, that the results in this section hold under more general dependence patterns. In particular, at the cost of additional notational complexity, it is straightforward to extend the result in this section to partitions of outcomes $(\mathcal{C}_0, \mathcal{C}_1, \ldots, \mathcal{C}_m)$ such that for each $k \in \{0, 1, \ldots, m\}$, $s\in\{0,1\}$ and $a\in\{0,1\}$, there exists $k' \in \{0, 1, \ldots, m\} \setminus \{k\}$ with $\{\eta_{i,j}\}_{(i,j) \in \mathcal{R}_s \times \mathcal{C}_k} \hspace{1mm}{\perp \!\!\! \perp}~ \{\eta_{i,j}, \varepsilon_{i,j}^{(a)}\}_{(i,j) \in \mc{R}_{1-s} \times \mathcal{C}_{k'}}$ and $\{\eta_{i,j}\}_{(i,j) \in \mc{R}_{1-s} \times \mathcal{C}_{k'}} \hspace{1mm}{\perp \!\!\! \perp}~ \{\eta_{i,j}, \varepsilon_{i,j}^{(a)}\}_{(i,j) \in \mathcal{R}_s \times \mathcal{C}_k}$. This allows for rather general patterns of dependence across outcomes while preserving independence across specific sets of outcomes (e.g., certain product categories in the retail example of (ref)).

Recall the setup from (ref): Given an observation matrix $S \in\{\mathbb{R} \cup \{\,?\}\}^{N \times M}$, a matrix completion algorithm $\texttt{MC}$ produces an estimate $\what{T} = \texttt{MC}(S) \in \mathbb{R}^{N \times M}$ of a matrix of interest $T$, where $S$ and $T$ are related via (ref). With this background, we now describe the $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ meta-algorithm.

We represent this meta-algorithm succinctly as below:

In summary, $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ produces an estimate $\what{T}$ such that for each block $\mc I \in \mathcal{P}$, the sub-matrix $\what{T}_{\mc I}$ is constructed only using the entries of $S$ corresponding to the remaining three blocks of $\mathcal{P}$. (ref)$(a)$ provides a schematic of the block partition $\mc P$ for $\mathcal{R}_0 = [\left\lfloor N/2 \right \rfloor]$ and $\mathcal{C}_0 = [\left\lfloor M/2 \right \rfloor]$. See (ref)$(b)$ for a visualization of $S \otimes \boldsymbol{1}^{-\mc I}$. The following result, proven in (ref), establishes $({\what{\Theta}}^{(0)},{\what{\Theta}}^{(1)}, {\what{P}})$ generated by $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ satisfy (ref).

A host of $\texttt{MC}$ algorithms are designed to de-noise and impute missing entries of matrices under random patterns of missingness; the most common missingness pattern studied is where each entry has the same probability of being missing, independent of everything else. In contrast, $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ generates patterns where all entries in one block are deterministically missing, as in (ref)$(b)$. A recent strand of research on the interplay between matrix completion methods and causal inference models---specifically, within the synthetic controls framework---has contributed matrix completion algorithms that allow for block missingness athey2021matrix, agarwal2021robustness, bai2021matrix, SI, arkhangelsky2021synthetic, agarwal2023causal, dwivedi2022counterfactual,dwivedi2022doubly. However, it is a challenge to apply known theoretical guarantees for these methods to the setting in this article because of: (i) the use of cross-fitting---which creates blocks where all observations are missing---and (ii) outside of the completely-missing blocks, there can still be missing observations with heterogeneous probabilities of missingness. In the next section, we show how to modify an $\texttt{MC}$ algorithm designed for block missingness patterns so that it can be applied to our setting with cross-fitting and heterogeneous probabilities of missingness outside the folds. For concreteness, we work with the Tall-Wide matrix completion algorithm of bai2021matrix.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{The \texorpdfstring{$\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$} algorithm}

$\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$ is an end-to-end $\texttt{MC}$ algorithm obtained by instantiating the $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ meta-algorithm with the Tall-Wide algorithm of bai2021matrix, which we denote as $\texttt{TW}$. For completeness, we detail the $\texttt{TW}$ algorithm in (ref), and then use it to describe $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$ in (ref). \@startsection{subsubsection}{3}{0mm}{-\baselineskip}{0.05\baselineskip}{\normalfont}{The \texorpdfstring{$\texttt{TW}$} algorithm of \texorpdfstring{bai2021matrix}.} bai2021matrix propose $\texttt{TW}$ to impute missing values in matrices with a set of rows and a set of columns without missing entries. More concretely, for any matrix $S \in \{ \mathbb{R} \cup \{\,?\} \}^{N \times M}$, let $\mathcal{R}_{\mathrm{obs}} \subseteq [N]$ and $\mathcal{C}_{\mathrm{obs}} \subseteq [M]$ denote the set of all rows and all columns, respectively, with all entries observed. Then, all missing entries of $S$ belong to the block $\mathcal{I} = \mathcal{R}_{\mathrm{miss}} \times \mathcal{C}_{\mathrm{miss}}$, where $\mathcal{R}_{\mathrm{miss}} \triangleq [N] \setminus \mathcal{R}_{\mathrm{obs}}$ and $\mathcal{C}_{\mathrm{miss}} \triangleq [M] \setminus \mathcal{C}_{\mathrm{obs}}$.

Given a rank hyper-parameter $r \in [\min\{ | \mathcal{R}_{\mathrm{obs}} |,| \mathcal{C}_{\mathrm{obs}} | \}]$, $\texttt{TW}_r$ produces an estimate of $T$ as follows:

\@startsection{subsubsection}{3}{0mm}{-\baselineskip}{0.05\baselineskip}{\normalfont}{\texorpdfstring{$\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$} algorithm.}

We provide intuition on the key steps of the $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$ algorithm next.

{\bf Computing ${\what{P}}$.} The estimate ${\what{P}}$ comes from applying $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ with $\texttt{TW}$ on ${A}$ and truncating the entries of the resulting matrix to the range $[{\bar\lambda}, 1-{\bar\lambda}]$, in accordance with (ref). The $\texttt{TW}$ sub-routine is directly applicable to ${A}$, because for any block $\mathcal{I} = \mathcal{R}_s \times \mathcal{C}_k \in \mathcal{P}$ the masked matrix ${A} \otimes \boldsymbol{1}^{-\mc I}$ has $[N] \setminus \mathcal{R}_s$ fully observed rows and $[M] \setminus \mathcal{C}_k$ fully observed columns. See (ref)$(a)$ for a visualization of ${A} \otimes \boldsymbol{1}^{-\mc I}$.

{\bf Computing ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$.} The estimates ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$ are constructed by applying $\texttt{Cross}\allowbreak\texttt{-}\allowbreak\texttt{Fitted}\allowbreak\texttt{-}\allowbreak\texttt{MC}$ with $\texttt{TW}$ on ${{Y}}^{(0),\mathrm{full}}_{}$ and ${{Y}}^{(1),\mathrm{full}}_{}$, which do not have missing entries. $\texttt{TW}$ is not directly applicable on ${{Y}}^{(0),\mathrm{obs}}_{}$ and ${{Y}}^{(1),\mathrm{obs}}_{}$, as both matrices may not have any rows and columns that are fully observed. See (ref)$(b)$ and (ref)$(c)$ for visualizations of ${{Y}}^{(0),\mathrm{obs}}_{}\otimes \boldsymbol{1}^{-\mathcal{I}}$ and ${{Y}}^{(1),\mathrm{obs}}_{}\otimes \boldsymbol{1}^{-\mathcal{I}}$, respectively. However, notice that, due to (ref)(ref) and (ref)(ref),

As a result, $\texttt{MC}({{Y}}^{(0),\mathrm{full}}_{})$ and $\texttt{MC}({{Y}}^{(1),\mathrm{full}}_{})$ provide estimates of ${\Theta}^{(0)} \odot (\boldsymbol{1} - {P})$ and ${\Theta}^{(1)} \odot {P}$, respectively---recall the discussion in (ref). To construct ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$, we divide the entries of $\texttt{MC}({{Y}}^{(0),\mathrm{full}}_{})$ and $\texttt{MC}({{Y}}^{(1),\mathrm{full}}_{})$ by the entries of $(\boldsymbol{1} - {\what{P}})$ and ${\what{P}}$, respectively, to adjust for heterogeneous probabilities of missingness jin2021factor,bhattacharya2022matrix,xiong2023large. This inverse probability of treatment weighting adjustment to estimate ${\what{\Theta}}^{(0)}$ and ${\what{\Theta}}^{(1)}$ is distinct and in addition to the augmented IPW procedure that generates $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$ from estimates ${\what{\Theta}}^{(0)}$, ${\what{\Theta}}^{(1)}$ and ${\what{P}}$.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Theoretical guarantees for \texorpdfstring{$\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$}}

To establish theoretical guarantees for $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$, we adopt three assumptions from bai2021matrix. The first assumption imposes a low-rank structure on the matrices ${P}$, ${\Theta}^{(0)}$, and ${\Theta}^{(1)}$, namely that their entries are given by an inner product of latent factors.

(ref) decomposes each of the unobserved confounders (${P}$, ${\Theta}^{(0)}$, and ${\Theta}^{(0)}$) into low-dimensional unit-dependent latent factors ($U$, $U^{(0)}$, and $U^{(1)}$) and measurement-dependent latent factors ($V$, $V^{(0)}$, and $V^{(1)}$). In particular, every unit $i \in [N]$ is associated with three low-dimensional factors: $(i)$ $U_i \in \mathbb{R}^{r_p}$, $(ii)$ $U_i^{(0)} \in \mathbb{R}^{r_{\theta_0}}$, and $(iii)$ $U_i^{(1)} \in \mathbb{R}^{r_{\theta_1}}$. Similarly, every measurement $j \in [M]$ is associated with three factors: $(i)$ $V_j \in \mathbb{R}^{r_p}$, $(ii)$ $V_j^{(0)} \in \mathbb{R}^{r_{\theta_0}}$, and $(iii)$ $V_j^{(1)} \in \mathbb{R}^{r_{\theta_1}}$. Low-rank assumptions are standard in the matrix completion literature.

The second assumption requires that the factors that determine ${P}$, ${\Theta}^{(0)} \odot (\boldsymbol{1} - {P})$, and ${\Theta}^{(1)} \odot {P}$ explain a sufficiently large amount of the variation in the data. This assumption is made on the factors of ${\Theta}^{(0)} \odot (\boldsymbol{1} - {P})$ and ${\Theta}^{(1)} \odot {P}$ instead of ${\Theta}^{(0)}$ and ${\Theta}^{(1)}$ as the $\texttt{TW}$ algorithm is applied on ${{Y}}^{(0),\mathrm{full}}_{}= {Y} \odot (\boldsymbol{1} - {A})$ and ${{Y}}^{(1),\mathrm{full}}_{}= {Y} \odot {A}$, instead of ${{Y}}^{(0),\mathrm{obs}}_{}$ and ${{Y}}^{(1),\mathrm{obs}}_{}$ (see steps (ref) and (ref) of $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$). To determine the factors of ${\Theta}^{(0)} \odot (\boldsymbol{1} - {P})$ and ${\Theta}^{(1)} \odot {P}$, let

where $\boldsymbol{1}_{N}\in \mathbb{R}^N$ and $\boldsymbol{1}_{M} \in \mathbb{R}^M$ are vectors of all $1$'s. Then,

where $\wbar{U}^{(0)} \triangleq \wbar{U} * U^{(0)} \in \mathbb{R}^{N \times r_{\theta_0}(r_p + 1)}$, $\wbar{V}^{(0)} \triangleq \wbar{V} * V^{(0)} \in \mathbb{R}^{M \times r_{\theta_0}(r_p + 1)}$, $\wbar{U}^{(1)} \triangleq U * U^{(1)} \in \mathbb{R}^{N \times r_{\theta_1}r_p}$, and $\wbar{V}^{(1)} \triangleq V * V^{(1)} \in \mathbb{R}^{M \times r_{\theta_1}r_p}$, with the operator $*$ denoting the Khatri-Rao product (see (ref)). We provide details of the derivation of these factors in (ref).

(ref), a classic assumption in the literature on latent factor models, ensures that the factor structure is strong. Specifically, it ensures that each eigenvector of ${P}$, ${\Theta}^{(0)} \odot (\boldsymbol{1} - {P})$, and ${\Theta}^{(1)} \odot {P}$ carries sufficiently large signal.

The third assumption requires a strong factor structure on the sub-matrices of ${P}$, ${\Theta}^{(0)} \odot (\boldsymbol{1} - {P})$, and ${\Theta}^{(1)} \odot {P}$ corresponding to every block $\mathcal{I}$ in the block partition $\mathcal{P}$ from (ref). Further, it also requires that the size $\mathcal{I}$ grows linearly in $N$ and $M$.

The subsequent assumption introduces additional conditions on the noise variables in bai2021matrix than those specified in (ref).

(ref)(ref) and (ref)(ref) requires the noise variables to exhibit only weak dependency across measurements. Still, these assumptions allow the existence of pairs of perfectly correlated outcomes (e.g., $j,j'\in [M]$ such that $a_{i,j}=a_{i,j'}$). (ref)(ref) requires the noise $( {E}^{(a)}, {W} )$ to be jointly independent across units, for every $a \in \{ 0,1 \}$. We are now ready to provide guarantees on the estimates produced by $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$. The proof can be found in (ref).

(ref) implies that the conditions (ref) in (ref) hold whenever $N^{1/2}/M=o(1)$. Then, the DR estimator from (ref) constructed using $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$ estimates ${\what{\Theta}}^{(0)}$, ${\what{\Theta}}^{(1)}$, and ${\what{P}}$ exhibits an asymptotic Gaussian distribution centered at the target causal estimand. Further, (ref) implies that the estimation errors $\mc{E}\big( {\what{P}} \big)$ and $\mc{E}\big( {\what{\Theta}}^{} \big)$ achieve the parametric rate whenever $N/M = O(1)$.

\@startsection{subsection}{2}{0mm}{-\baselineskip}{0.25\baselineskip}{\normalfont\bf}{Application to panel data with staggered adoption} (ref) considered a setting with block independence between noise (formalized in (ref)). (ref) discusses how to extend the proposed doubly-robust framework to a setting of panel data with staggered adoption, where this assumption may not hold. Recall (from (ref)) that in the panel data setting $M$ measurements correspond to $T$ time periods, and $t$ denotes the time index. Then, (ref) considers a setting where a unit remains under control for some period of time, after which it deterministically remains under treatment. In other words, for every unit $i \in [N]$, there exists a time point $t_i \in [T]$ such that $a_{i, t} = 0$ for $t \leq t_i$, and $a_{i, t} = 1$ for $t > t_i$. Such a treatment assignment pattern leads to a heavy dependence in the noise $\{\eta_{i, t}\}_{t \in [T]}$ for every unit $i \in [N]$. (ref) describes an alternative approach to the $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$ algorithm and shows that (ref) still holds for a suitable staggered adoption model.

\@startsection{section}{1}{0mm}{-\baselineskip}{0.5\baselineskip}{\center\normalfont\bf}{Simulations} This section reports simulation results on the performance of the DR estimator of (ref) and the OI and IPW estimators of (ref), respectively.

{\bf Data Generating Process (DGP).} We now briefly describe the DGP for our simulations; (ref) provides details. All simulations set $N=M$. To generate, ${P}$, ${\Theta}^{(0)}$, and ${\Theta}^{(1)}$, we use the latent factor model given in (ref). To introduce unobserved confounding, we set the unit-specific latent factors to be the same across ${P}$, ${\Theta}^{(0)}$, and ${\Theta}^{(1)}$, i.e., $U = U^{(0)} = U^{(1)}$. The entries of $U$ and the measurement-specific latent factors, $V, V^{(0)}, V^{(1)}$ are each sampled independently from a uniform distribution, with hyperparameter $r_p$ equal to the dimension of $U$ and $V$, and hyperparameter $r_p$ equal to the dimension of $U^{(a)}$ and $V^{(a)}$ for $a=0,1$. Further, the entries of the noise matrices ${E}^{(0)}$ and ${E}^{(1)}$ are sampled independently from a normal distribution, and the entries of ${W}$ are sampled independently as in (ref). Then, $y^{(a)}_{i, j}$, $a_{i,j}$, and $y_{i,j}$ are determined from (ref), respectively. The simulation generates ${P}$, ${\Theta}^{(0)}$, and ${\Theta}^{(1)}$ once. Given the fixed values of ${P}$, ${\Theta}^{(0)}$, and ${\Theta}^{(1)}$, the simulation generates $2500$ realizations of $(Y,A)$---that is, only the noise matrices ${E}^{(0)}, {E}^{(1)}, {W}$ are resampled for each of the $2500$ realizations. For each simulation realization, we apply the $\texttt{Cross-}\allowbreak\texttt{Fitted-}\allowbreak\texttt{SVD}$ algorithm with hyper-parameters as in (ref) and ${\bar\lambda} = \lambda = 0.05$ to obtain ${\what{P}}$, ${\what{\Theta}}^{(0)}$, and ${\what{\Theta}}^{(1)}$, and compute $\mathrm{ATE}_{\cdot,j}$ from (ref), and $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{OI}}$, $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{IPW}}$ and $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$ from (ref).

{\bf Results.} (ref) reports simulation results for $N = 1000$, with $r_p = 3$, $ r_{\theta} = 3$ in Panel $(a)$, and $r_p = 5$, $ r_{\theta} = 3$ in Panel $(b)$. (ref) in (ref) reports simulation results for $r_p = 3$, $ r_{\theta} = 5$. In each case, the figure shows a histogram of the distribution of $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}- \mathrm{ATE}_{\cdot,j}$ across $2500$ simulation instances for a fixed $j$, along with the best fitting Gaussian distribution (green curve). The histogram counts are normalized so that the area under the histogram integrates to one. (ref) plots the Gaussian distribution in the result of (ref) (black curve). The dashed blue, red and green lines in (ref) indicate the values of the means of the OI, IPW, and DR error, respectively, across simulation instances. For reference, we place a black solid line at zero. The DR estimator has minimal bias and a close-to-Gaussian distribution. The biases of OI and IPW are non-negligible. In (ref), we compare the biases and the standard deviations of OI, IPW, and DR across many $j$.

Panels $(a)$, $(b)$, and $(c)$ of (ref) report coverage rates over the 2500 simulations for $\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}$-centered nominal $95\%$ confidence intervals with $N = 500$, $N = 1000$, and $N = 1500$, respectively, all with $M=N$ and $r_p = r_{\theta} = 3$. For every $j \in [M]$, panels $(a)$, $(b)$ and $(c)$ show $\what{\mathrm{c}}_j$, the percentage of times $[\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}\pm 1.96 \what{\sigma}_j/\sqrt{N}]$ covers $\mathrm{ATE}_{\cdot,j}$ (in blue), and ${\mathrm{c}}_j$, the percentage of times $[\what{\mathrm{ATE}}{}_{\cdot,j}^{\,\mathrm{DR}}\pm 1.96 {\sigma}_j/\sqrt{N}]$ covers $\mathrm{ATE}_{\cdot,j}$ (in green). Panel $(d)$ shows the means and standard deviations of $\{ \what{\mathrm{c}}_j \}_{j \in [M]}$ and $\{ {\mathrm{c}}_j \}_{j \in [M]}$ for different values of $N$. Confidence intervals based on the large-sample approximation results of (ref) exhibit small size distortion even for fairly small values of $N$.

\@startsection{section}{1}{0mm}{-\baselineskip}{0.5\baselineskip}{\center\normalfont\bf}{Conclusion}

This article introduces a new framework to estimate treatment effects in the presence unobserved confounding. We consider modern data-rich environments, where there are many units, and outcomes of interest per unit. We show it is possible to control for the confounding effects of a set of latent variables when this set is low-dimensional relative to the number of observed treatments and outcomes.

Our proposed estimator is doubly-robust, combining outcome imputation and inverse probability weighting with matrix completion. Analytical tractability of its distribution is gained through a novel cross-fitting procedure for causal matrix completion. We study the properties of the doubly-robust estimator, along with the outcome imputation and inverse probability weighting-based estimators under black-box matrix completion error rates. We show that the decay rate of the error of the doubly-robust estimator dominates those of the outcome imputation and the inverse probability weighting estimators. Moreover, we establish a Gaussian approximation to the distribution of the doubly-robust estimator. Simulation results demonstrate the practical relevance of the formal properties of the doubly-robust estimator.