The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
81,964 characters
Covariate-Adjusted Deep Causal Learning for Heterogeneous Panel Data Models
\newcounter{appendixcounter}
\refstepcounter{appendixcounter}
\label{append:staggered-design}
\refstepcounter{appendixcounter}
\label{section:proofs}
\refstepcounter{appendixcounter}
\label{append:simulation}
\refstepcounter{appendixcounter}
\label{append: covid19}
\title{Covariate-Adjusted Deep Causal Learning for Heterogeneous Panel Data Models}
\author{
Guanhao Zhou, Yuefeng Han, and Xiufan Yu \\ University of Notre Dame
}
\date{}
\maketitle{}
\pagestyle{plain}
\setstretch{1.2}
\begin{abstract}
This paper studies the task of estimating heterogeneous treatment effects in causal panel data models, in the presence of covariate effects.
We propose a novel \textbf{Co}variate-Adjusted \textbf{DE}ep C\textbf{A}usal \textbf{L}earning (\textbf{CoDEAL}) for panel data models, that employs flexible model structures and powerful neural network architectures to cohesively deal with the underlying heterogeneity and nonlinearity of both panel units and covariate effects.
The proposed CoDEAL integrates nonlinear covariate effect components (parameterized by a feed-forward neural network) with nonlinear factor structures (modeled by a multi-output autoencoder) to form a heterogeneous causal panel model.
The nonlinear covariate component offers a flexible framework for capturing the complex influences of covariates on outcomes.
The nonlinear factor analysis enables CoDEAL to effectively capture both cross-sectional and temporal dependencies inherent in the data panel.
This latent structural information is subsequently integrated into a customized matrix completion algorithm, thereby facilitating more accurate imputation of missing counterfactual outcomes.
Moreover, the use of a multi-output autoencoder explicitly accounts for heterogeneity across units and enhances the model interpretability of the latent factors.
We establish theoretical guarantees on the convergence of the estimated counterfactuals, and demonstrate the compelling performance of the proposed method using extensive simulation studies and a real data application.
\end{abstract}
\noindent Keywords: causal panel data models, counterfactual estimation, heterogeneous treatment effects, matrix completion, missing not at random, multi-output autoencoders, nonlinear factor models
\newpage
\setstretch{1.5}
\section{Introduction}
Causal inference in panel data settings has attracted growing interest in recent years, with broad applications across various fields such as economics \citep{clarke2024synthetic}, life science
\citep{helske2024estimating}, political science \citep{imai2021use}, and social science \citep{imbens2024causal}.
Causal panel data analysis is inherently a missing data problem due to the \emph{fundamental problem of causal inference}: for each unit and time period, we could only observe either the treated or the untreated outcome, but never both. The panel data structure brings unique challenges to causal analysis. Different units may get exposed to the treatment at various time periods, forming a block-wise non-random missing pattern \citep{athey2022design}. The presence of both temporal and cross-sectional dependencies could bias the estimation if not properly addressed \citep{sun2021estimating}, and unit-specific heterogeneity can confound causal interpretation \citep{millimet2023fixed}, especially if additional covariates are present.
To estimate causal effects, one natural approach is to first impute the missing counterfactual, and then examine the difference between the imputed counterfactuals and the actual observed outcomes. \citet{athey2021matrix} formulated the imputation of missing counterfactuals as a matrix completion task with nuclear norm minimization.
It is important to note that standard matrix completion approaches are not directly suitable for causal imputation \citep{choi2024matrix}. Most conventional matrix completion methods operate under the assumption that missing entries are \emph{missing at random}, which enables estimation of missing values without explicitly modeling the missingness mechanism.
However, in causal panel data settings where different units receive treatment at different times, missing entries correspond to counterfactual outcomes that are systematically unobserved due to treatment assignments, invalidating the missing-at-random assumption.
As a result, applying standard matrix completion blindly without properly accounting for this structured, non-random missingness can lead to biased estimation \citep{agarwal2023causal}.
In recent years, a growing body of work has focused on developing causal matrix completion methods that handle missing-not-at-random panel designs.
Besides the optimization-based causal matrix completion \citep{athey2021matrix,choi2024matrix}, another line of research leverages factor models \citep{bai2003,bai2009panel} to capture the latent low-rank structure of observed entries.
Most factor-model-based causal matrix completion methods in the literature \citep{xu2017generalized,bai2021matrix,agarwal2023causal,xiong2023large,yan2024entrywise} posit a linear factor structure. Linear factor analysis helps capture cross-sectional and temporal dependencies to some extent \citep{han2024cp,han2024tensor,yu2024power}, but is often limited in capturing complex nonlinear relationships, motivating the generalization to nonlinear factor models \citep{yalcin2001nonlinear}.
Another key aspect in panel data analysis is to account for unit-specific heterogeneity \citep{millimet2023fixed,semenova2023inference}. Most relevant works in the literature still focus on the overall average treatment effects across units \citep{athey2021matrix,athey2025identification}, as opposed to the unit-specific heterogeneous treatment effects. Pooled estimates may mask important variations across units. Units in panel data often differ in characteristics and respond to the treatment differently \citep{wager2018estimation,nandy2023generalized}, which may bias the estimates of treatment effects, motivating a critical need to estimate unit-specific heterogeneous treatment effects and adjust for the covariate effects if possible.
In this paper, we develop \textbf{Co}variate-Adjusted \textbf{DE}ep C\textbf{A}usal \textbf{L}earning (\textbf{CoDEAL}) for panel data models, a deep-learning-based causal learning method for estimating heterogeneous treatment effects. CoDEAL captures covariate effects with a deep neural network regression, and models the covariate-adjusted outcomes via a nonlinear factor structure.
Our contributions can be summarized as follows.
\begin{itemize}[leftmargin = *]
\item CoDEAL provides \textbf{a causal matrix completion method} that addresses the non-random missingness mechanism in causal panel data analysis.
By leveraging nonlinear factor analysis, CoDEAL \textbf{exploits both cross-sectional and temporal dependencies}.
This latent structural information is subsequently integrated into the matrix completion procedure, thereby facilitating more accurate imputation of missing counterfactual outcomes.
\item CoDEAL employs flexible model structures and powerful neural network architectures to cohesively deal with the underlying \textbf{heterogeneity} and \textbf{nonlinearity} of \textbf{both panel units and covariate effects}.
Algorithmically, we begin with a feed-forward deep neural network (DNN) to flexibly remove nonlinear covariate influences, followed by a multi-output autoencoder (AE) to recover nonlinear latent factors and explicitly take the unit-specific heterogeneity into account.
\item CoDEAL offers \textbf{a unifying framework} for heterogeneous causal panel models in the presence of covariates. By combining nonlinear covariate adjustment with deep latent factor models, CoDEAL generalizes beyond traditional linear factor models and unifies a class of factor-model-based causal matrix completion approaches, e.g., \cite{xu2017generalized,bai2021matrix,agarwal2023causal,xiong2023large,yan2024entrywise}. Furthermore, with different choices of factor models and neural network architectures, the proposed framework can be further generalized to accommodate various data types such as spatial and tensor data.
\end{itemize}
\noindent \textbf{Related Work.} \textbf{(a) An overview of causal learning methods in panel data models.} One traditional approach to estimate treatment effects in causal panel models is \emph{the difference-in-differences} \citep{imai2021use,athey2022design,wing2024designing} that compares the average changes in outcomes over time between the treated and untreated units under a \emph{parallel trend} assumption.
Alternatively, one can estimate the causal effects by first imputing the missing counterfactual outcomes,
and then examining the difference between the imputed counterfactuals and the actual observed outcomes.
Following this strand, two approaches have emerged in recent causal panel literature: \emph{the uncounfoundedness-based approach} (also known as horizontal regression) \citep{rosenbaum1983central,imbens2015causal} and \emph{the synthetic control method} (also known as vertical regression) \citep{abadie2003economic,abadie2010synthetic,abadie2015comparative,doudchenko2016balancing,abadie2021using}.
The close connections between the difference-in-differences, unconfoundedness-based methods, synthetic control approaches, and causal matrix completion estimators were revealed by \citet{athey2021matrix}. Despite their distinct appearances, \citet{athey2021matrix} showed that the linear versions of all four estimators can be characterized as solutions to the same optimization problem with the exact same objective function, but subject to different restrictions on parameters of that optimization.
We refer readers to \cite{arkhangelsky2024causal} for a comprehensive review of past and recent progress on causal panel models.
\noindent \textbf{(b) Causal matrix completion and deep matrix completion.} A surge of interest in the recent past has led to new research on causal matrix completion methods capable of addressing missing-not-at-random panel data, e.g., optimization-based methods \citep{athey2021matrix,choi2024matrix} and factor-model-based methods \citep{xu2017generalized,bai2021matrix,agarwal2023causal,xiong2023large,yan2024entrywise}.
There has been a handful of deep-learning-based methods for matrix completion, e.g., \cite{fan2018matrix,radhakrishnan2022simple,xiu2024deep,fan2024neuron}. However, most of them assume missing-at-random, and therefore, not directly suitable for causal imputation in panel data models.
\vspace{-2ex}
\section{Methodology}
\vspace{-2ex}
\subsection{Problem Setup}
\noindent \textbf{Heterogeneous Treatment Effects \footnote{Here, the ``heterogeneous treatment effect'' refers to the heterogeneity across units.
}
in Causal Panel Data Models. }
We begin by setting up the causal learning problem in panel data. Consider a panel data setting with $N$ units over $T$ periods.
Let $[n]$ denote the set $\{1, 2, \ldots, n\}$.
For each unit $i \in [N]$ and time point $t \in [T]$, there are two potential outcomes, $Y_{it}(1)$ and $Y_{it}(0)$, representing the outcomes the $i$-th unit would experience at time $t$ under treatment and control, respectively.
The treatment assignment is captured by a binary indicator variable $W_{it} \in \{0, 1\}$, where $W_{it} = 1$ indicates that the unit $i$ receives the treatment at time $t$ and $W_{it} = 0$ otherwise. Suppose we also have access to some unit-specific covariates $\mathrm{\bf X}_i \in \mathbb{R}^P$ that could potentially affect the observations $\{Y_{it}, t \in [T]\}$. We denote the matrix forms of the observed outcomes, treatment indicators, and covariates by $\mathrm{\bf Y} = (Y_{it}) \in \mathbb{R}^{N \times T}$, $\boldsymbol{\Omega} = (W_{it}) \in \mathbb{R}^{N \times T}$, and $\mathrm{\bf X} = (\mathrm{\bf X}_1, \ldots, \mathrm{\bf X}_N)^\top \in \mathbb{R}^{N \times P}$, respectively. Let $\mathrm{\bf Y}(1) = (Y_{it}(1))\in \mathbb{R}^{N \times T}$ and $\mathrm{\bf Y}(0) = (Y_{it}(0))\in \mathbb{R}^{N \times T}$ be the matrices of potential outcomes of the treated and untreated, respectively.
Our aim\footnote{Aligning with relevant works in the literature \citep{athey2021matrix,athey2025identification},
we focus on the unit-specific average treatment effects on the treated (ATT) as our primary estimand of interest to demonstrate our proposed model and algorithm. The proposed method can be also be used to many other estimands,
such as the overall (non-unit-specific) ATT $\mathbb{E}_{(i,t)} [Y_{it}(1) - Y_{it}(0) \mid W_{it}=1]$, the unit-specific average treatment effects (ATE) $ \mathbb{E}_{t} [Y_{it}(1) - Y_{it}(0)]$, and the overall (non-unit-specific) ATE $\mathbb{E}_{(i,t)} [Y_{it}(1) - Y_{it}(0)]$. } in this work is to estimate the unit-specific (potentially heterogeneous) average treatment effects on the treated (ATT), $\tau_i^\star = \mathbb{E}_{t} [Y_{it}(1) - Y_{it}(0) \mid W_{it} = 1]$ for each treated unit $i$.
\bigskip
\noindent \textbf{Staggered Adoption Design. }
This work focuses on staggered adoption design \citep{athey2022design}, where treatments can be rolled out at different times across various units and are irreversible once initiated.
We assume there is at least one never-treated unit in the panel.
Details about staggered design are in Appendix \ref{append:staggered-design}.
\bigskip
\noindent \textbf{Causal Assumptions. } To ensure the validity of the estimators' interpretation, we adopt the following assumptions throughout the analysis. (i) (Stable Unit Treatment Value Assumption (SUTVA).) The potential outcomes for a given unit at a particular time depend solely on the treatment administered to that unit at that time, and are unaffected by treatments assigned to other units or at other time points. (ii) (Static Treatment Effects.) The treatment effects are non-dynamic, i.e., do not change over time.
\bigskip
\noindent \textbf{Challenges and Limitations of Existing Methods. } The \emph{fundamental problem of causal inference} implies that $Y_{it}(0)$ and $Y_{it}(1)$ cannot be observed simultaneously. The panel data forms two partially observed data matrices $\mathrm{\bf Y}(0)$ and $\mathrm{\bf Y}(1)$, where the missing pattern is determined according to the treatment assignments. To estimate the treatment effects, one natural idea is to estimate the unobserved counterfactuals, in other words, to impute the missing potential outcomes. Since our focus is the unit-specific ATT and $\{Y_{it}(1): W_{it}=1\}$ is observed, it suffices to estimate $\{Y_{it}(0): W_{it}=1\}$ to impute missing elements in $\mathrm{\bf Y}(0)$. Let $\{\widehat{Y}_{it}(0): W_{it}=1\}$ denote imputed values.
The unit-specific ATT can then be estimated by $\widehat\tau_i = \sum_{t: W_{it} = 1} [Y_{it}(1)-\widehat{Y}_{it}(0)]/\sum_{t} W_{it}$ for each treated unit $i$.
There are three main limitations of existing methods that this work aims to overcome in the task of modeling unit-specific ATT and imputing missing potential outcomes.
First, the vast majority of matrix completion methods rely on the assumption of \emph{missing at random}, an assumption that is typically violated in staggered adoption designs, thereby leading to biased estimates.
Second, existing matrix-completion-based causal panel models mostly focus on the overall (non-unit-specific) ATT, without adequately accounting for unit-level heterogeneity.
Third, many existing methods often lack the capacity to disentangle the causal impact of treatment from confounding covariate effects and to model the underlying complex nonlinear dependencies, limiting both the precision of counterfactual predictions and the interpretability of causal analysis.
\subsection{Covariate-Adjusted Deep Causal Learning (\textbf{CoDEAL}) for Panel Data Models}
\label{subsec:method-CoDEAL}
To tackle the above challenges, we propose the CoDEAL method that (i) leverages the staggered adoption structure to fully exploit both cross-sectional and temporal dependencies for more accurate imputation of missing potential outcomes, (ii) models the unit-specific heterogeneity explicitly to accommodate variation in treatment effects across units, and (iii) involves flexible model structures and employs powerful neural network architectures to efficiently capture complex nonlinear relationships.
For each unit $i \in [N]$ at time $t \in [T]$, we model the observed outcome by
\vspace{-1ex}
\begin{equation}
Y_{it}= \phi_i^\star(\mathrm{\bf F}^\star_{t}) + g^\star_{t}(\mathrm{\bf X}_i) + \tau_i^\star \cdot \mathds{1}\{W_{it}=1\} + \varepsilon_{it},
\label{eq:CoDEAL-Model}
\vspace{-1ex}
\end{equation}
where $\mathrm{\bf F}^\star_t \in \mathbb{R}^K$ is a $K$-dimension vector of latent factors, $\phi^\star_{i}(\cdot): \mathbb{R}^K \rightarrow \mathbb{R}$ is a potentially nonlinear factor loading function, $g^\star_t(\cdot): \mathbb{R}^P \rightarrow \mathbb{R}$ is a potentially nonlinear function that captures how covariates influence the outcome, and $ \mathds{1}\{\cdot\}$ is an indicator function, and $\varepsilon_{it}$ is the idiosyncratic noise. CoDEAL employs a nonlinear factor model \citep{yalcin2001nonlinear} to characterize an underlying low-rank structure of the covariate-adjusted outcomes. By letting $\phi^\star_{i}(\cdot)$ take different forms, CoDEAL unifies a class of factor-model-based causal matrix completion approaches, e.g., \citep{xu2017generalized,bai2021matrix,agarwal2023causal,xiong2023large,yan2024entrywise}.
As shown by \eqref{eq:CoDEAL-Model},
CoDEAL accommodates heterogeneous but non-dynamic treatment effects. The covariate effects can be heterogeneous across units and time-varying over periods.
Here, we allow $N$ and $T$ to be large and diverging, while $P$ and $K$ are fixed and small.
For clarity, we first present the learning procedures using the simplest form of a $2\times 2$ staggered adoption (also known as the four-block design) in Section \ref{subsec:CoDeal-Four-Block}, and generalize the algorithm to the more general staggered adoption design in Section \ref{subsec:CoDeal-Staggered-Design}. More details about the four-block and staggered adoption design as well as their graphical illustrations are in Appendix \ref{append:staggered-design}. A graphical illustration of CoDEAL in four-block design is provided in Figure \ref{fig:method}.
The complete CoDEAL method in general staggered adoption design is summarized in Algorithm \ref{alg:staggered}.
\subsubsection{CoDEAL in Four-Block Design} \label{subsec:CoDeal-Four-Block}
We begin with the four-block design, the simplest form of the staggered adoption and often regarded as foundational building blocks for more general staggered adoption designs. The four-block design is a $2 \times 2$ staggered adoption design, where a subset of $N_1$ units never receive the treatment while the remaining $N_2 = N-N_1$ units are exposed to the irreversible treatment starting at time $t=T_1+1$. Without loss of generality, we rearrange the panel such that the first $N_1$ units ($i=1,\ldots, N_1$) are untreated. We call the time periods $t\in [T_1]$ as pre-treatment periods and $t\in [T]\backslash [T_1]$ as post-treatment periods. Let $T_2 = T-T_1$ denote the time duration of the post-treatment periods. The four-block design naturally decomposes the panel into four blocks:
\begin{equation}
\mathrm{\bf Y} =
\begin{bmatrix}
\mathrm{\bf Y}_\mathcal{A} & \mathrm{\bf Y}_\mathcal{B} \\
\mathrm{\bf Y}_\mathcal{C} & \mathrm{\bf Y}_\mathcal{D}
\end{bmatrix}
=
\begin{bmatrix}
\mathrm{\bf Y}_{\mathcal{A}}(0) & \mathrm{\bf Y}_{\mathcal{B}}(0) \\
\mathrm{\bf Y}_{\mathcal{C}}(0) & \mathrm{\bf Y}_{\mathcal{D}}(1)
\end{bmatrix}
\quad \text{and } \quad
\mathrm{\bf Y}(0) =
\begin{bmatrix}
\mathrm{\bf Y}_{\mathcal{A}}(0) & \mathrm{\bf Y}_{\mathcal{B}}(0) \\
\mathrm{\bf Y}_{\mathcal{C}}(0) & ?
\end{bmatrix}
\label{eq:four-block}
\end{equation}
where $\mathrm{\bf Y}_\mathcal{A} \in \mathbb{R}^{N_1\times T_1}$, $\mathrm{\bf Y}_\mathcal{B} \in \mathbb{R}^{N_1\times T_2}$, $\mathrm{\bf Y}_\mathcal{C} \in \mathbb{R}^{N_2\times T_1}$, and $\mathrm{\bf Y}_\mathcal{D} \in \mathbb{R}^{N_2\times T_2}$.
By construction, $\mathrm{\bf Y}_\mathcal{A}$, $\mathrm{\bf Y}_\mathcal{B}$, $\mathrm{\bf Y}_\mathcal{C}$ are observed outcomes of the untreated with $Y_{it} = Y_{it}(0)$ for $(i,t)\in \mathcal{A} \cup \mathcal{B} \cup \mathcal{C}$, while $\mathrm{\bf Y}_\mathcal{D}$ are observed outcomes of the treated with $Y_{it} = Y_{it}(1)$ for $(i,t)\in \mathcal{D}$, leaving a missing block in $\mathrm{\bf Y}(0)$. To estimate $\tau_i^\star$, it suffices to impute the missing values of $\mathrm{\bf Y}_{\mathcal{D}}(0)$. Denote the imputed block as $\widehat{\mathrm{\bf Y}}_{\mathcal{D}}(0)$. The unit-specific ATT can be then estimated by
$\widehat{\tau}_i = \frac{1}{T_2} \sum_{t: (i,t)\in\mathcal{D}} [Y_{it}(1) -\widehat{Y}_{it}(0)] \text{ for } i \in [N]\backslash [N_1]$.
\begin{figure}[tb]
\centering
\includegraphics[width=\textwidth]{Figures-Finalized/codeal5.png}
\caption{A graphical illustration of the proposed CoDEAL in a four-block design. Gray and orange refer to the untreated and treated blocks, respectively. Inputs are observed outcomes \(\mathbf{Y}\in\mathbb{R}^{N\times T}\), covariates \(\mathrm{\bf X}\in\mathbb{R}^{N\times P}\), indicator matrix \(\boldsymbol{\Omega}\in\mathbb{R}^{N\times T}\). Output is the estimated unit-specific ATT $\widehat{\tau}_i$.}
\label{fig:method}
\end{figure}
\bigskip
\noindent \textbf{Modeling Covariate Effects. } To account for the potentially nonlinear, heterogeneous, and time-varying effects brought by covariates, we consider a fully connected DNN for $g_t^\star(\cdot)$.
We construct our model using a fully connected DNN with ReLU activation,
due to its strong empirical performance. We refer to this architecture as a \textit{deep ReLU network}. Let $L\in\mathbb N$ denote the network depth and let $\mathrm{\bf d} = (d_1, \ldots, d_{L+1}) \in \mathbb{N}^{L+1}$ specify the layer widths. A deep ReLU network is a function mapping $\mathbb{R}^{d_0}$ to $\mathbb{R}^{d_{L+1}}$ and takes the form
\begin{equation}\label{eq:relu}
h(x) = \mathcal{L}_{L+1} \circ \bar{\sigma}_{L} \circ \mathcal{L}_L \circ \bar{\sigma}_{L-1} \circ \cdots \circ \mathcal{L}_2 \circ \bar{\sigma}_{1} \circ \mathcal{L}_1(x),
\end{equation}
where each $\mathcal{L}_\ell(z) = M_\ell z + b_\ell$ is an affine transformation with weight matrix $M_\ell \in \mathbb{R}^{d_\ell \times d_{\ell-1}}$ and bias vector $b_\ell \in \mathbb{R}^{d_\ell}$, and $\bar{\sigma}_\ell : \mathbb{R}^{d_\ell} \to \mathbb{R}^{d_\ell}$ applies the ReLU activation function entrywise. For simplicity, we refer to both $M_\ell$ and $b_\ell$ as the weights of the deep ReLU network.
\begin{definition}[Deep ReLU network class]
For any $L \in \mathbb{N}$, $\mathrm{\bf d} \in \mathbb{N}^{L+1}$, and constants $B, C>0$, we define the class of deep ReLU networks with depth $L$, width parameter $\mathrm{\bf d}$, and weights $\mathrm{\bf M}=\{M_1,...,M_{L+1},b_1,...,b_{L+1}\}$ bounded by $C$ as
\[
\mathcal{G}_{n_0}^{n_{L+1}}(L,\mathrm{\bf d}, B, C) = \left\{ \tilde{h}(x) = \text{sgn}(h(x))(|h(x)| \wedge B) : h \text{ of form } \eqref{eq:relu}, \, \|M_\ell\|_{\infty} \leq C, \, \|b_\ell\|_{\infty} \leq C \right\} .
\]
\end{definition}
The covariate effects can be estimated via
\begin{equation}
\widehat{g}_t(\cdot) = \underset{g_t \in \mathcal{G}_{P}^{1}(L,\mathrm{\bf d}, B, C) }{\arg\min} \sum_{i: W_{it}=0} (Y_{it} - g_t(\mathrm{\bf X}_i) )^2 ,
\label{eq:g-hat}
\end{equation}
for certain functional class $\mathcal{G}_{P}^{1}(L,\mathrm{\bf d}, B, C)$ specified in Appendix~\ref{section:proofs}.
We define the covariate-adjusted observed outcomes for unit $i \in [N]$ at time $t \in [T]$ as $\widetilde{Y}_{it} = Y_{it} - \widehat{g}_t(\mathrm{\bf X}_i).$
\bigskip
\noindent \textbf{Estimating Counterfactuals by Matrix Completion via a Multi-Output AE.} With covariate-adjusted outcomes, it remains to model the nonlinear factor structure $\phi^\star_i(\mathrm{\bf F}_t^\star)$ in \eqref{eq:CoDEAL-Model}. Inspired by \cite{xiu2024deep}, we consider tackling $\phi^\star_i(\mathrm{\bf F}_t^\star)$ using a multi-output AE.
The architecture of the multi-output autoencoder (AE) consists of a single DNN encoder and multiple independent DNN decoders. The encoder $\rho(\cdot) \in \mathcal{G}^{K_1}_{N}(L_1, \mathrm{\bf d}_1, B,C)$ maps the input to a $K_1$-dimensional
bottleneck representation. For each of the $N$ output components, a distinct DNN decoder is used; specifically, the decoder for the $i$-th output is given by $\phi_i(\cdot) \in \mathcal{G}^{1}_{K_1}(L_2, \mathrm{\bf d}_2, B, C)$ for $i \in [N]$. Consequently, the $i$-th output of the multi-output AE for an input $x$ is given by $\phi_i \circ \rho(x)$. We formally define the function class of the multi-output AE as:
\begin{align}\label{eq:ae-class}
\mathcal{G}_{\text{AE}}^{K_1} := \left\{ (\phi_1, \ldots, \phi_N) \circ \rho : \rho \in \mathcal{G}^{K_1}_{N}(L_1, \mathrm{\bf d}_1, B, C), \, \phi_i \in \mathcal{G}^{1}_{K_1}(L_2, \mathrm{\bf d}_2, B, C), \, i \in [N] \right\},
\end{align}
where $K_1\ge K$. Further details on this function class are provided in Appendix~\ref{section:proofs}.
With the multi-output AE, the encoder $\rho^\star(\cdot)$ effectively approximates the latent factors, facilitating the factor-model-based matrix completion. The unit-specific decoders $\phi_i^\star(\cdot)$, $i \in [N]$, are better equipped to capture the heterogeneity across units, in contrast to a traditional single-output AE which limits $\phi_1^\star = \ldots =\phi_N^\star$.
DNNs involved in the multi-output AE are estimated by
\begin{equation}
(\widehat{\phi}_1,\dots,\widehat{\phi}_N)\circ\widehat{\rho}(\cdot)=\underset{(\phi_1,\dots,\phi_N)\circ{\rho}(\cdot)\in \mathcal{G}_{\text{AE}}^{K_1} }{\arg \min}\sum_{(i,t):W_{it}=0}\Bigr(\phi_i(\rho(\widetilde{\mathrm{\bf Y}}_{\cdot t}))-\widetilde{Y}_{it}\Bigr)^2 .
\label{eq:multi-output-ae}
\end{equation}
The counterfactuals $\mathrm{\bf Y}_{\mathcal{D}}(0)$ can then be imputed by $\widehat{Y}_{it}(0) = \widehat{\phi}_i (\widehat\rho(\widetilde{\mathrm{\bf Y}}_{\cdot t})) + \widehat{g}_t(\mathrm{\bf X}_i)$ for any $(i,t) \in \mathcal{D}$.
The complete steps of imputing counterfactuals in the four-block design are summarized as a standalone algorithm in Algorithm \ref{alg:four-block} to facilitate further generalizations to staggered design.
\begin{algorithm}[h]
\caption{Estimating Counterfactuals by Matrix Completion in the Four-Block Design}
\label{alg:four-block}
\begin{algorithmic}[1]
\STATE \textbf{Input:} Observed outcomes \(\mathrm{\bf Y}\), covariates $\mathrm{\bf X}$, the estimated covariate effects $\widehat{g}_t(\cdot)$ from \eqref{eq:g-hat}, treatment indicators $\boldsymbol{\Omega}$ with a four-block design in \eqref{eq:four-block}, latent factor dimension $K$
\STATE \textbf{Output:} Imputed outcome matrix $\widehat{\mathrm{\bf Y}}_{\mathcal{D}}(0)$.
\STATE Obtain the covariate-adjusted observed outcomes by $\widetilde{Y}_{it} = Y_{it} - \widehat{g}_t(\mathrm{\bf X}_i)$ for $i\in[N]$, $t\in[T]$.
\STATE Fit a multi-output autoencoder using the entries in the untreated part $\{\widetilde{Y}_{it}(0): W_{it}=0\}$ and obtain $(\widehat{\phi}_1, \ldots, \widehat{\phi}_N, \widehat{\rho})$ following \eqref{eq:multi-output-ae}
\STATE Estimate the missing elements
by $\widehat{Y}_{it}(0) = \widehat{\phi}_i (\widehat\rho(\widetilde{\mathrm{\bf Y}}_{\cdot t})) + \widehat{g}_t(\mathrm{\bf X}_i)$ for any $(i,t) \in \mathcal{D}$.
\RETURN $\widehat{\mathrm{\bf Y}}_{\mathcal{D}}(0)$.
\end{algorithmic}
\end{algorithm}
\begin{algorithm}[h]
\caption{\textbf{CoDEAL}: Staggered Adoption Design}
\label{alg:staggered}
\begin{algorithmic}[1]
\STATE \textbf{Input:} Observed outcomes \(\mathbf{Y}\), covariates \(\mathrm{\bf X}\), treatment indicators \(\boldsymbol{\Omega}\), latent factor dimension \(K\).
\STATE \textbf{Output:} Estimated unit-specific ATT $\widehat{\tau}_i$ for each treated unit $i$.
\STATE Obtain the estimated covariate effect $\widehat{g}_t(\cdot)$ by \eqref{eq:g-hat}
\STATE Initialize $\widehat{\mathrm{\bf Y}}(0)\in \mathbb{R}^{N\times T}$ by setting $\widehat{Y}_{it}(0) = Y_{it}$ for the untreated $\{(i,t): W_{it} = 0\}$.
\STATE Extract block partitions \(\{N_\xi \}_{1\le \xi \le r}\) and \(\{T_\eta\}_{1\le \eta\le r}\) from \(\boldsymbol{\Omega}\).
\FOR{\(\xi_0 = 2, \dots, r\)}
\FOR{\(\eta_0 = r+2-\xi_0, \dots, r\)}
\STATE Construct the four-block submatrix \( \mathrm{\bf Y}^{(\xi_0,\eta_0)}\) and \(\boldsymbol{\Omega}^{(\xi_0,\eta_0)}\) according to Equation~\eqref{eq:staggered_design}.
\STATE Call Algorithm~\ref{alg:four-block} with inputs \( \mathrm{\bf Y}^{(\xi_0,\eta_0)}\), $\mathrm{\bf X}^{(\xi_0)}$, $\widehat{g}_t(\cdot)$, \(\boldsymbol{\Omega}^{(\xi_0,\eta_0)}\), and \(K\) to obtain the imputed block $\widehat{\mathrm{\bf Y}}_{\mathcal{D}_{(\xi_0,\eta_0)}}(0)$, and extract the estimated block \(\widehat{\mathrm{\bf Y}}_{(\xi_0,\eta_0)}(0)\) from $\widehat{\mathrm{\bf Y}}_{\mathcal{D}_{(\xi_0,\eta_0)}}(0)$.
\STATE Update the matrix \(\widehat{\mathrm{\bf Y}}(0)\) with the imputed block \(\widehat{\mathrm{\bf Y}}_{(\xi_0,\eta_0)(0)}\).
\ENDFOR
\ENDFOR
\STATE Compute unit-specific estimated ATT as $\widehat\tau_i = \sum_{t: W_{it} = 1} [Y_{it}(1)-\widehat{Y}_{it}(0)]/\sum_{t} W_{it}$.
\RETURN Estimated unit-specific ATT $\widehat{\tau}_i$.
\end{algorithmic}
\end{algorithm}
\subsubsection{CoDEAL in Staggered Adoption Design} \label{subsec:CoDeal-Staggered-Design}
The extension to staggered adoption design is a natural generalization of the four-block design.
Under the staggered adoption design, we rearrange the panels by sorting the $N$ units based on the time at which they began receiving the treatment in a descending order (with never-treated units at the top, followed by the latest adopters, and then progressively earlier adopters towards the bottom). This sorting procedure yields a block structure in the treatment indicator matrix $\boldsymbol{\Omega}$. By grouping the units based on their adoption timing and segmenting the overall time periods into relevant intervals, we partition the indicator matrix $\boldsymbol{\Omega}$ (as well as the outcome matrix $\mathrm{\bf Y}$) into submatrices $\boldsymbol{\Omega}_{(\xi, \eta)} \in \mathbb{R}^{N_{\xi} \times T_{\eta}}$ and $\mathrm{\bf Y}_{(\xi, \eta)} \in \mathbb{R}^{N_{\xi} \times T_{\eta}}$ for $\xi=1,\ldots, r$ and $\eta =1,\ldots, r$, where $\sum_{\xi=1}^{r} N_{\xi} =N $ and $\sum_{\eta=1}^{r} T_{\eta} = T $.
Examples of block partitions in the staggered design are included in Appendix \ref{append:staggered-design}.
Let $\mathcal{I}_{tr} = \{(\xi,\eta): \boldsymbol{\Omega}_{(\xi, \eta)} = 1\}$ denote the index set of the treated blocks. By the construction of these partitions, the indices in $\mathcal{I}_{tr}$ satisfy $\xi+\eta > r+1$.
To estimate the unit-specific ATT, it suffices to estimate the counterfactuals of $\mathrm{\bf Y}_{(\xi,\eta)}$ for all $(\xi,\eta) \in \mathcal{I}_{tr}$. The imputation of $\mathrm{\bf Y}_{(\xi,\eta)}(0)$ can be reduced to a four-block design problem \citep{yan2024entrywise}, therefore the method proposed in the previous Section \ref{subsec:CoDeal-Four-Block} can be directly applied. Specifically, to estimate $\mathrm{\bf Y}_{(\xi_0, \eta_0)}(0)$ for some $(\xi_0, \eta_0) \in \mathcal{I}_{tr}$ , we construct a four-block data matrix $ \mathrm{\bf Y}^{(\xi_0,\eta_0)}$ as follows
\begin{equation}
\mathrm{\bf Y}^{(\xi_0,\eta_0)} =
\begin{bmatrix}
\mathrm{\bf Y}_{\mathcal{A}_{(\xi_0,\eta_0)}} & \mathrm{\bf Y}_{\mathcal{B}_{(\xi_0,\eta_0)}} \\
\mathrm{\bf Y}_{\mathcal{C}_{(\xi_0,\eta_0)}} & \mathrm{\bf Y}_{\mathcal{D}_{(\xi_0,\eta_0)}}
\end{bmatrix}
=
\begin{bmatrix}
\mathrm{\bf Y}_{\mathcal{A}_{(\xi_0,\eta_0)}}(0) & \mathrm{\bf Y}_{\mathcal{B}_{(\xi_0,\eta_0)}}(0) \\
\mathrm{\bf Y}_{\mathcal{C}_{(\xi_0,\eta_0)}}(0) & \mathrm{\bf Y}_{\mathcal{D}_{(\xi_0,\eta_0)}}(1)
\end{bmatrix}
\label{eq:staggered_design}
\end{equation}
where
$\mathcal{A}_{(\xi_0,\eta_0)} = \{(\xi,\eta): 1 \leq \xi \leq k_1, 1 \leq \eta \leq k_2 \}$, $\mathcal{B}_{(\xi_0,\eta_0)} = \{(\xi,\eta):1\leq \xi \leq k_1, k_2+1 \leq \eta \leq \eta_0 \}$, $\mathcal{C}_{(\xi_0,\eta_0)} = \{(\xi,\eta): k_1+1 \leq \xi \leq \xi_0, 1 \leq j \leq k_2 \}$, and $\mathcal{D}_{(\xi_0,\eta_0)} = \{(\xi,\eta): k_1+1 \leq \xi \leq \xi_0, k_2+1 \leq j \leq \eta_0 \}$
with \( k_1 = r + 1 - \eta_0 \) and \( k_2 = r + 1 - \xi_0 \).
The associated treatment indicator matrix $\boldsymbol{\Omega}^{(\xi_0,\eta_0)}$ is defined analogously and the associated covariate matrix $\mathrm{\bf X}^{(\xi_0)}$ can be obtained by retaining the units included in $\mathrm{\bf Y}^{(\xi_0,\eta_0)}$.
The missing potential outcomes $\mathrm{\bf Y}_{\mathcal{D}_{(\xi_0,\eta_0)}}(0)$ can be directly estimated using Algorithm \ref{alg:four-block}.
By construction, $\mathrm{\bf Y}_{\mathcal{D}_{(\xi_0,\eta_0)}} \supseteq \mathrm{\bf Y}_{(\xi_0,\eta_0)}$. In the case when $\mathrm{\bf Y}_{\mathcal{D}_{(\xi_0,\eta_0)}} \supset \mathrm{\bf Y}_{(\xi_0,\eta_0)}$, we treat all the blocks of $\mathrm{\bf Y}_{\mathcal{D}_{(\xi_0,\eta_0)}}$ as missing and impute their values, but only extract $\widehat{\mathrm{\bf Y}}_{(\xi_0,\eta_0)}(0)$ from the imputed $\widehat{\mathrm{\bf Y}}_{\mathcal{D}_{(\xi_0,\eta_0)}}(0)$; see Figure \ref{fig:eg-construction-of-4-block-design} for examples.
With the four-block construction in \eqref{eq:staggered_design} and the CoDEAL method in four-block design introduced in Section \ref{subsec:CoDeal-Four-Block}, we obtain the imputed values of $\mathrm{\bf Y}_{(\xi, \eta)}(0)$ for all $(\xi,\eta) \in \mathcal{I}_{tr}$. The unit-specific ATT can then be estimated by $\widehat\tau_i = \sum_{t: W_{it} = 1} [Y_{it}(1)-\widehat{Y}_{it}(0)]/\sum_{t} W_{it}$ for each treated unit $i$. The complete algorithm of CoDEAL in the staggered adoption design is summarized as Algorithm \ref{alg:staggered}.
\begin{figure}[bt]
\centering
\resizebox{0.47\linewidth}{!}{
\begin{tabular}{|c|c|c|c||c||c|}\hline
& $T_1$ & $T_2$ & $T_3$ & $T_4$ & $T_5$ \\ \hline \hline
$N_1$ & \cellcolor{green!20}$\mathrm{\bf Y}_{(1,1)}$ & \cellcolor{green!20}$\mathrm{\bf Y}_{(1,2)}$& \cellcolor{green!20}$\mathrm{\bf Y}_{(1,3)}$& \cellcolor{red!15}$\mathrm{\bf Y}_{(1,4)}$ & $\mathrm{\bf Y}_{(1,5)}$\\
\hline
$N_2$ & \cellcolor{green!20}$\mathrm{\bf Y}_{(2,1)}$& \cellcolor{green!20}$\mathrm{\bf Y}_{(2,2)}$& \cellcolor{green!20}$\mathrm{\bf Y}_{(2,3)}$& \cellcolor{red!15}$\mathrm{\bf Y}_{(2,4)}$& \\
\hline \hline
$N_3$ & \cellcolor{blue!15}$\mathrm{\bf Y}_{(3,1)}$ & \cellcolor{blue!15}$\mathrm{\bf Y}_{(3,2)}$& \cellcolor{blue!15}$\mathrm{\bf Y}_{(3,3)}$& \cellcolor{orange!30}?& \\
\hline \hline
$N_4$ &$\mathrm{\bf Y}_{(4,1)}$& $\mathrm{\bf Y}_{(4,2)}$&& & \\
\hline
$N_5$ & $\mathrm{\bf Y}_{(5,1)}$& & & & \\
\hline
\end{tabular}
}
\quad
\resizebox{0.47\linewidth}{!}{
\begin{tabular}{|c||c||c|c||c|c|} \hline
& $T_1$ & $T_2$ & $T_3$ & $T_4$ & $T_5$ \\ \hline \hline
$N_1$ & \cellcolor{green!20}$\mathrm{\bf Y}_{(1,1)}$ & \cellcolor{red!15}$\mathrm{\bf Y}_{(1,2)}$& \cellcolor{red!15}$\mathrm{\bf Y}_{(1,3)}$ & $\mathrm{\bf Y}_{(1,4)}$& $\mathrm{\bf Y}_{(1,5)}$\\ \hline
$N_2$ & \cellcolor{green!20}$\mathrm{\bf Y}_{(2,1)}$ & \cellcolor{red!15}$\mathrm{\bf Y}_{(1,2)}$& \cellcolor{red!15}$\mathrm{\bf Y}_{(2,3)}$ & $\mathrm{\bf Y}_{(2,4)}$& \\ \hline
$N_3$ & \cellcolor{green!20}$\mathrm{\bf Y}_{(3,1)}$& \cellcolor{red!15}$\mathrm{\bf Y}_{(3,2)}$& \cellcolor{red!15}$\mathrm{\bf Y}_{(3,3)}$& & \\ \hline \hline
$N_4$ & \cellcolor{blue!15}$\mathrm{\bf Y}_{(4,1)}$ & \cellcolor{orange!30}$\mathrm{\bf Y}_{(4,2)}$& \cellcolor{orange!30}& & \\ \hline
$N_5$ & \cellcolor{blue!15}$\mathrm{\bf Y}_{(5,1)}$& \cellcolor{orange!30}& \cellcolor{orange!30}?& & \\ \hline
\end{tabular}
}
\caption{Examples of constructing a four-block submatrix \(\mathrm{\bf Y}^{(\xi_0,\eta_0)}\) to estimate \(\mathrm{\bf Y}_{(\xi_0,\eta_0)}(0)\). Here, \( r=5 \). The question mark denotes the block of interest, with $(\xi_0,\eta_0)$ = (3,4) (Left) and (5,3) (Right). The associated four blocks are designed by green (\(\mathrm{\bf Y}_{\mathcal{A}_{(\xi_0,\eta_0)}}\)), red (\(\mathrm{\bf Y}_{\mathcal{B}_{(\xi_0,\eta_0)}}\)), blue (\(\mathrm{\bf Y}_{\mathcal{C}_{(\xi_0,\eta_0)}} \)), and orange (\(\mathrm{\bf Y}_{\mathcal{D}_{(\xi_0,\eta_0)}}\)) as in \eqref{eq:staggered_design}.
}
\label{fig:eg-construction-of-4-block-design}
\end{figure}
\subsection{Theoretical Properties} \label{subsec:theoretical-results}
\begin{assumption}\label{asmp:density}
Suppose $\mathrm{\bf X}_i$ is supported on $[0,1]^P$ and has a density function $f$ satisfying
$\sup_{x \in [0,1]^P} f(x) =: f_{\max} < \infty,$
where $f_{\max}$ is independent of dimension $P$. We assume there exists a constant $B > 0$ such that $\mathbb{E} F_t^\star=0$ and $\|F_t^\star\|_\infty \leq B$ holds almost surely. We further assume that $\varepsilon_{it}$ is sub-Gaussian with sub-Gaussian norm $\sigma_\epsilon^2$, and that $\varepsilon_{it}, F_t^\star,\mathrm{\bf X}_i$ are mutually independent.
\end{assumption}
\begin{assumption}\label{asmp:smooth}
The functions $\{g_t^\star\}$ are assumed to be $(\beta, C)$-H\"older smooth, that is, each function is $\lfloor \beta \rfloor$ times continuously differentiable, and its $\lfloor \beta \rfloor$th derivative is H\"older continuous with exponent $\beta - \lfloor \beta \rfloor$ and constant $C$. Here, $\lfloor \beta \rfloor$ denotes the greatest integer less than $\beta$.
\end{assumption}
\begin{assumption}[Pervasiveness] \label{asmp:factor}
There exist a matrix $M^\star \in \mathbb{R}^{K \times N}$ and a function $\rho^\star$, defined on the image of the mapping $M^\star \phi^\star$, denoted by $M^\star \phi^\star([-B,B]^K)$. Assume that $\rho^\star $ is also $(\beta, C)$-H\"older smooth, and $\|M^\star\|_\infty \lesssim n^{-1}$, $\|M^\star\|_0 \asymp n$ for some diverging integer $n>0$. We further assume that
{\small
\begin{equation}\label{eq:factor}
\sup_{x \in [-B,B]^K} \| \rho^\star(M^\star \phi^\star(x)) - x \|^2 \lesssim n^{-1}.
\end{equation}
}
\end{assumption}
Assumption \ref{asmp:density} is a standard assumption in the literature of nonparametric regression. Note that assuming $\mathrm{\bf X}_i$ has compact support is equivalent to assuming that the support is $[0,1]^P$ via centering and scaling.
The smoothness condition is standard in the nonparametric regression literature (cf. \cite{tsybakov2003introduction}, Chapter 1), as it governs the complexity of the underlying function class. A formal definition is provided in Appendix \ref{section:proofs}.
Assumption \ref{asmp:factor} generalizes the classical pervasiveness condition (with $n = N$) from linear factor models \citep{bai2003}. It effectively enables approximate recovery of the latent factors from the observed input: when $x$ denotes the latent factors, condition \eqref{eq:factor} ensures the existence of a reconstruction map $\rho^\star \circ M^\star$ that approximately inverts the encoding $\phi^\star(x)$ and recovers $x$ from the noiseless data.
To illustrate, consider the linear factor model setup in \cite{bai2003}, where $\phi^\star(x) = A x$ with bounded loadings satisfying $\|A\|_\infty \lesssim 1$ and $\|A\|_0 \asymp N$ indicating that each factor affects nearly all observed variables. Assumption~\ref{asmp:factor} is satisfied by setting $M^\star = N^{-1} A^\top$ and $\rho^\star(x) = N (A^\top A)^{-1} x$.
The following theorem establishes the convergence of the estimated counterfactuals under a four-block design. The result can be readily extended to general staggered adoption designs.
\begin{theorem}\label{thm1}
Suppose that Assumptions \ref{asmp:density}-\ref{asmp:factor} hold, and let $K_1\ge K$, with $K_1$ defined in \eqref{eq:ae-class}. Assume that $\log(N+T)=o(n)$, $N_1\lesssim N_2$ and $T_1\asymp T_2$. Then, under a suitable choice of parameters for the DNN function class (specified in Appendix~\ref{section:proofs}), we have
\begin{align*}
& \frac{1}{T_2 N_2} \sum_{t=T_1+1}^{T}\sum_{i=N_1+1}^N \mathbb{E} \left( \widehat Y_{it}(0) - \phi_i^\star(\mathrm{\bf F}_{t}^\star) - g_t^\star(\mathrm{\bf X}_i) \right)^2 \nonumber \\
=& O_{\mathbb{P}} \left( \left( T^{-\frac{2\beta}{2\beta + K}} + T^{-1} n + N^{-1}K_1 + n^{-1} \right) \log^4(N T) +N^{-\frac{2\beta}{2\beta + k}} \log^5 N \right).
\end{align*}
\end{theorem}
Additional theoretical properties, technical conditions and proofs are detailed in Appendix~\ref{section:proofs}.
\vspace{-2ex}
\section{Simulation Studies}
\label{sec:simulation}
In this section, we use extensive simulation studies to demonstrate the finite-sample performance of the proposed CoDEAL\footnote{Our code for implementing CoDEAL is available in the supplementary material.}.
\vspace{-1.5ex}
\subsection{Benchmark Methods} \label{append:simu-benchmark}
We consider five benchmark methods, including a single-output AE (\textbf{Single AE}) (detailed below), neuron-enhanced AE (\textbf{AEMC-NE}) \citep{fan2024neuron}, matrix completion with nuclear norm minimization (\textbf{MC-NNM}) \citep{athey2021matrix}, vertical regression (\textbf{Vert-Reg}) \citep{doudchenko2016balancing}, and difference-in-differences (\textbf{DiD}) \citep{imai2021use}.
Additionally, we evaluate the algorithm's performance in terms of the covariate effect removal. To this end, we compare the performance of CoDEAL using different methods to remove the covariate effects, specifically, DNN-based removal (as illustrated in Section \ref{subsec:method-CoDEAL}), and linear-regression(LR)-based removal (detailed below).
We evaluate the algorithm's performance using two evaluation metrics:
the mean absolute error
$MAE=(\sum_{(i, t): W_{it} = 1}|Y_{it}(1)-\widehat{{Y}}_{it}(0)-\tau_{i}^\star|)/(\sum_{(i,t)} W_{it})$
and the mean squared error $MSE=(\sum_{(i,t): W_{it} = 1}(Y_{it}(1)-\widehat{{Y}}_{it}(0)-\tau_{i}^\star)^2)/(\sum_{(i,t)} W_{it})$.
\bigskip
\noindent {\textbf{Matrix Completion by Single-Output AE (as a benchmark method).}}
Our proposed CoDEAL utilizes a multi-output autoencoder (as illustrated in Figure \ref{fig:method}) to recover nonlinear latent factors while explicitly modeling the unit-specific heterogeneity through multiple decoders.
The multi‐output autoencoder uses a shared encoder network \(\rho(\cdot)\) to map each column \(\widetilde{\mathbf{Y}}_{\!\cdot t}\) into a \(K\)-dimensional latent space, followed by \(N\) separate decoder networks \(\{\phi_1,\dots,\phi_N\}\) to reconstruct each unit’s value.
To demonstrate CoDEAL's ability to account for the unit-specific heterogeneity, we include an additional benchmark that replaces the multi-output autoencoder with a single-output autoencoder in the CoDEAL algorithm. This benchmark is referred to as the Single AE in later discussions.
A single-output autoencoder uses the same encoder \(\rho(\cdot)\) but only one decoder \(\phi(\cdot)\) (i.e. $\phi=\phi_1=\dots=\phi_N$) for all units. By design, the common decoder makes the single-output AE inherently unable to model the unit-specific heterogeneity effectively. For implementation, DNNs in the single-output AE are estimated by
\begin{equation*}
\widehat{\phi}\circ\widehat{\rho}(\cdot)=\underset{\phi\circ{\rho}(\cdot)\in \mathcal{G}_{\text{AE}}^{K_1} }{\arg \min}\sum_{(i,t):W_{it}=0}\Bigr(\phi(\rho(\widetilde{\mathrm{\bf Y}}_{\cdot t}))-\widetilde{Y}_{it}\Bigr)^2 .
\end{equation*}
\bigskip
\noindent {\textbf{Covariate Effect Removal by Linear Regression (as a benchmark method).}}
Our proposed CoDEAL uses a DNN to remove the potentially nonlinear, heterogeneous, and time-varying effects brought by covariates. To evaluate the algorithm's performance in terms of the covariate effect removal, we include a linear-regression(LR)-based removal as a benchmark. The LR-based removal assumes that covariates have a linear impact on the outcome.
With covariates $\mathbf{X}_i\in\mathbb{R}^{P}$, $i\in [N]$, we denote the effect of covariates on the outcome by a coefficient matrix $\mathrm{\bf H} = [\mathrm{\bf H}_1, \ldots, \mathrm{\bf H}_T] \in\mathbb{R}^{P\times T}$, and we estimate $\mathrm{\bf H}_t \in \mathbb{R}^P$ separately for each period.
For each pre-treatment time $t\in\{1,\dots,T_1\}$, we estimate the covariate effect vector $\widehat{\mathrm{\bf H}}_t\in\mathbb{R}^{P}$ using OLS on all units:
$\widehat{\mathrm{\bf H}}_t = \arg\min_{\mathrm{\bf H}_t} \sum_{i=1}^{N}\left(Y_{it}-\mathrm{\bf X}_i^\top\mathrm{\bf H}_t\right)^2.$
For each post-treatment period $t\in\{T_1+1,\dots,T\}$, we use only control units to estimate the covariate effect:
$\widehat{\mathrm{\bf H}}_t = \arg\min_{\mathrm{\bf H}_t} \sum_{i: W_{it}=0}\left(Y_{it}-\mathrm{\bf X}_i^\top\mathrm{\bf H}_t\right)^2.$
With the estimated $\widehat{\mathrm{\bf H}}_t, t\in [T]$, we obtain the LR-based covariate-adjusted
observed outcomes by
$\widetilde{Y}_{it} = Y_{it} - \mathrm{\bf X}_i^\top\widehat{\mathrm{\bf H}}_t$ for each $i \in [N]$ and $t \in [T]$.
\subsection{Data Generating Process and Simulation Configurations } \label{append:simulation-DGP}
We perform extensive simulation studies under various linear and nonlinear generating mechanisms, with different setups of unit numbers $N$, time periods $T$, matrix ranks $K$, and treatment assignments. Specifically, we consider the following configurations 1-4:
\begin{description}
\item[] 1. Four-Block Design: $(N,T,N_1,T_1,P,K)=(100,200,50,100,3,4)$
\item[] 2. Four-Block Design: $(N,T,N_1,T_1,P,K)=(200,120,100,60,5,3)$
\item[] 3. Staggered Adoption Design: $(N,T,K)=(100,120,4)$, $r=\{5,10\}$
\item[] 4. Staggered Adoption Design: $(N,T,K)=(200,120,3)$, $r=\{5,10\}$.
\end{description}
For the generation of factor effects and covariate effects, we first generate covariates
$\mathrm{\bf X} = (X_{it}) \in \mathbb{R}^{N\times P}$, factors $\mathrm{\bf F} = (F_{it}) \in \mathbb{R}^{T\times K}$, factor loadings $\boldsymbol{\Lambda} = (\Lambda_{it})\in \mathbb{R}^{N\times K}$ with $X_{it},F_{it}, \Lambda_{it}\overset{\text{i.i.d.}}{\sim}N(0,1)$. With these building blocks, we consider various linear and nonlinear generating mechanisms to generate factor effects $\phi_i^\star(\mathrm{\bf F}^\star_{t})$ and covariate effects $g^\star_t(\mathbf{X}_i)$. Specifically, we experiment with 4 variants for generating factor effects and 6 variants for generating covariate effects, as listed in Table \ref{tab:data_gen}. The unit-specific causal effect $\tau_i^\star \overset{i.i.d}{\sim}N(12,5) $ is designed to be added to treated units in treated periods. The observed outcome matrix $\mathrm{\bf Y} = (Y_{it}) \in \mathbb{R}^{N \times T}$ is then obtained by
$Y_{it}= Y^\star_{it} + \tau^\star_i \cdot \mathcal{I}\{W_{it}=1\} + \varepsilon_{it}$
where $W_{it}$ are determined based on the treatment assignments according to the simulation configuration, and the idiosyncratic error $\varepsilon_{it}$ follows $\varepsilon_{it}\overset{i.i.d}{\sim}\mathcal{N}(0,0.5^2)$.
\begin{table}[h]
\centering
\caption{Data Generation Processes for Latent Factor Effects and Covariate Effects}
\vspace{-1ex}
\resizebox{\textwidth}{!}{
\begin{tabular}{ll|ll}
\toprule
\multicolumn{2}{c|}{\textbf{Factor Effect Generation} ($\phi_i^\star(\mathrm{\bf F}^\star_{t})$)} & \multicolumn{2}{c}{\textbf{Covariate Effect Generation} ($g^\star_t(\mathbf{X}_i)$)} \\
\midrule
\textbf{Linear} & $0.5C_1\boldsymbol{\Lambda}_i^\top \mathbf{F}_t$, $C_1\overset{\text{i.i.d.}}{\sim}N(0,1)$& \textbf{Matrix Linear} & $\mathrm{\bf X}_i^\top\mathrm{\bf U}_t$, $\mathrm{\bf U} = (\mathrm{\bf U}_1, \ldots, \mathrm{\bf U}_T)\overset{\text{i.i.d.}}{\sim}N_{P\times T}(1,1)$\\
\textbf{Sine} & $2C_1\sin(\boldsymbol{\Lambda}_i^\top \mathbf{F}_t)$ & \textbf{Vector Linear} & $\mathrm{\bf X}_i^\top \mathrm{\bf U}$, $\mathrm{\bf U} \overset{\text{i.i.d.}}{\sim}N_P(1,1)$\\
\textbf{Polynomial} & $0.2C_1\boldsymbol{\Lambda}_i^\top \mathbf{F}_t + 0.2C_2\mathbf{F}_t^\top \mathbf{F}_t$, $C_2\overset{\text{i.i.d.}}{\sim}N(0,1)$& \textbf{tanh} & $[\tanh(|\mathbf{X}_i^\top \boldsymbol{w}_t|)]^{1/2}+b_t$, $\boldsymbol{w}_t\sim N_P(0,1)$\\
\textbf{ReLU MLP} & $\mathrm{\bf R}_i^{(2)}\text{ReLU}(\mathrm{\bf R}_i^{(1)}\mathbf{F}_t+\mathbf{b}_i^{(1)})+b_i^{(2)}$ & \textbf{poly} & $|\mathbf{X}_i^\top \boldsymbol{w}_t|^{1/2}+b_t$ \\
& $b_i^{(2)}\overset{\text{i.i.d.}}{\sim}N(0,0.5^2)$, $\mathrm{\bf R}_i^{(1)}\overset{\text{i.i.d.}}{\sim}N_{h_1\times r}(0,0.5^2)$& \textbf{log} & $\log(|\mathbf{X}_i^\top \boldsymbol{w}_t+b_t|)$, $b_{t}\overset{\text{i.i.d.}}{\sim}N(0,1)$ \\
& $\mathbf{b}_i^{(1)},\mathrm{\bf R}_i^{(2)}\overset{\text{i.i.d.}}{\sim}N_{h_1}(0,0.5^2)$, $h_1=10$& \textbf{ReLU} & $\text{ReLU}(\mathrm{\bf X}_i\mathrm{\bf R}_1^{(t)}+\mathbf{b}_1^{(t)})\mathrm{\bf R}_2^{(t)}+b_2^{(t)}, b_2^{(t)} \overset{\text{i.i.d.}}{\sim}N(0,1)$\\
& & &$\mathrm{\bf R}_1^{(t)}\overset{\text{i.i.d.}}{\sim}N_{P \times h_2}(0,\frac{2}{P}), \mathbf b_1^{(t)} \overset{\text{i.i.d.}}{\sim}N_{h_2}(0,1)$\\
& & &$\mathrm{\bf R}_2^{(t)} \overset{\text{i.i.d.}}{\sim}N_{h_2}(0,\frac{2}{h_2})$, $h_2=32$\\
\bottomrule
\end{tabular}}
\label{tab:data_gen}
Note: The vector linear covariate effect generation setting mimics a situation in which covariate effects are heterogeneous across panel units but remain time-invariant over time. The matrix linear covariate effect generation setting, however, allows for heterogeneous and time-varying covariate effects.
\end{table}
\subsection{Evaluation on the Performance of CoDEAL in the Removal of Covariate Effects}
\label{append:simu-covariate-removal}
To evaluate the algorithm’s performance in terms of the covariate effect removal, we compare the performance of CoDEAL using different methods, specifically, DNN (as illustrated in Section \ref{subsec:method-CoDEAL}) and LR (detailed in Section \ref{append:simu-benchmark}), to remove the covariate effects.
Six scenarios of covariate effects are considered as detailed in Table \ref{tab:data_gen}. We also include an oracle benchmark, when no covariate effects are present and no covariate-adjusted operations are needed.
\begin{table}[p]
\small
\centering
\caption{Comparisons of imputation accuracy (measured by MAE and MSE) with covariate effects removed by \textbf{LR} and \textbf{DNN} under various covariate generations in Config.1. Results are averaged over 50 replications. Numbers in parentheses are standard errors.}
\resizebox{\textwidth}{!}{
\begin{tabular}{llccccccll}
\toprule
\multirow{2}{*}{\textbf{Covariate}} & \multirow{2}{*}{\textbf{Method}} & \multicolumn{2}{c}{\textbf{LR}, Linear Factor}&\multicolumn{2}{c}{\textbf{NN}, Linear Factor}& \multicolumn{2}{c}{\textbf{LR}, Sine Factor}& \multicolumn{2}{c}{\textbf{NN}, Sine Factor}\\
& &MAE& MSE&MAE& MSE&MAE& MSE& \text{ }\text{ }\text{ }\text{ }MAE&\text{ }\text{ }\text{ }\text{ }\text{ }MSE\\
\midrule
\multirow{6}{*}{\textbf{No covariate}}
& CoDEAL & \textbf{0.560(0.009)} &0.578(0.028)&\textbf{0.560(0.009)} &0.578(0.028)&\textbf{0.933(0.015)} & \textbf{1.722(0.066)} & \textbf{0.933(0.015)} &\textbf{1.722(0.066)} \\
& Single AE & 0.724(0.009)&1.081(0.034)&0.724(0.009)&1.081(0.034)&1.071(0.013)& 2.131(0.056)& 1.071(0.013)&2.131(0.056)\\
& AEMC-NE & 0.734(0.009)&1.142(0.039)&0.734(0.009)&1.142(0.039)&1.090(0.013)& 2.172(0.056)& 1.090(0.013)&2.172(0.056)\\
& DiD & 0.804(0.010)&1.364(0.043)&0.804(0.010)&1.364(0.043)&1.141(0.013)& 2.343(0.059)& 1.141(0.013)&2.343(0.059)\\
& Vert-Reg & 0.601(0.002)&\textbf{0.574(0.004)} &0.601(0.002)&\textbf{0.574(0.004)} &1.024(0.012)& 2.228(0.069)& 1.024(0.012)&2.228(0.069)\\
& MC-NNM & 0.580(0.010)&0.593(0.030)&0.580(0.010)&0.593(0.030)&1.093(0.013)& 2.147(0.054)& 1.093(0.013)&2.147(0.054)\\
\midrule
\multirow{6}{*}{\textbf{Matrix Linear}}
& CoDEAL & \textbf{0.585(0.013)}&\textbf{0.615(0.024)}&\textbf{0.646(0.009)} &\textbf{0.722(0.023)} &\textbf{0.983(0.011)}& \textbf{1.889(0.046)}& \textbf{1.079(0.012)}&\textbf{2.135(0.054)} \\
& Single AE & 0.751(0.008)&1.131(0.032)&0.821(0.010)&1.320(0.039) &1.123(0.014)& 2.283(0.060)& 1.226(0.013)&2.626(0.062)\\
& AEMC-NE & 0.779(0.010)&1.244(0.041)&0.838(0.009)&1.394(0.041) &1.145(0.013)& 2.336(0.057)& 1.209(0.013)&2.564(0.060) \\
& DiD & 0.837(0.010)&1.444(0.045)&0.892(0.010)&1.587(0.045) &1.187(0.013)& 2.490(0.060)& 1.244(0.013)&2.696(0.063)\\
& Vert-Reg & 0.686(0.006)&0.766(0.015)&0.757(0.006)&0.957(0.019) &1.103(0.012)& 2.498(0.066)& 1.169(0.012)&2.647(0.064)\\
& MC-NNM & 0.604(0.009)&0.664(0.025)&0.699(0.006)&0.895(0.020) &1.140(0.013)& 2.293(0.055)& 1.201(0.013)&2.510(0.059)\\
\midrule
\multirow{6}{*}{\textbf{Vector Linear}}
& CoDEAL & \textbf{0.577(0.010)}&\textbf{0.592(0.026)}&\textbf{0.652(0.006)}&\textbf{0.715(0.017)}&\textbf{0.964(0.016)}& \textbf{1.775(0.060)}& \textbf{1.053(0.011)}&\textbf{2.035(0.045)} \\
& Single AE & 0.746(0.007)&1.106(0.026)&0.841(0.008)&1.373(0.030) &1.106(0.012)& 2.226(0.054)& 1.205(0.012)&2.549(0.053)\\
& AEMC-NE & 0.778(0.010)&1.205(0.032)&0.823(0.008)&1.317(0.029)&1.135(0.013)& 2.298(0.053)& 1.191(0.012)&2.497(0.053)\\
& DiD & 0.825(0.008)&1.373(0.031)&0.872(0.009)&1.489(0.033)&1.173(0.012)& 2.437(0.053)& 1.224(0.012)&2.620(0.056) \\
& Vert-Reg & 0.669(0.005)&0.727(0.012)&0.733(0.006)&0.887(0.017) &1.085(0.011)& 2.413(0.057)& 1.145(0.011)&2.520(0.053) \\
& MC-NNM & 0.591(0.008)&0.617(0.020)&0.678(0.007)&0.826(0.020) &1.106(0.017)& 2.219(0.099)& 1.175(0.012)&2.412(0.053)\\
\midrule
\multirow{6}{*}{\textbf{tanh}}
& CoDEAL & 0.778(0.008)&1.168(0.029)&\textbf{0.701(0.008)}&\textbf{0.871(0.031)}&\textbf{1.142(0.013)}& 2.366(0.059)& \textbf{1.080(0.013)}&\textbf{2.175(0.059)} \\
& Single AE & 0.889(0.008)&1.554(0.035)&0.851(0.008)&1.383(0.033) &1.299(0.014)& 2.910(0.067)& 1.257(0.014)&2.748(0.066) \\
& AEMC-NE & 0.988(0.008)&1.822(0.036)&0.862(0.009)&1.445(0.036) &1.305(0.013)& 2.925(0.064)& 1.239(0.014)&2.677(0.065) \\
& DiD & 0.855(0.008)&1.473(0.035)&0.916(0.009)&1.642(0.037) &1.199(0.013)& 2.536(0.061)& 1.274(0.014)&2.817(0.067) \\
& Vert-Reg & \textbf{0.721(0.005)} &\textbf{0.844(0.015)} &0.798(0.007)&1.075(0.024) &1.146(0.012)& 2.660(0.071)& 1.212(0.013)&2.787(0.072) \\
& MC-NNM & 0.837(0.007)&1.131(0.021)&0.735(0.005)&1.007(0.017)&1.153(0.013)& \textbf{2.339(0.056)} & 1.234(0.013)&2.636(0.063) \\
\midrule
\multirow{6}{*}{\textbf{log}}
& CoDEAL & \textbf{1.020(0.006)} &\textbf{1.873(0.026)} &\textbf{0.905(0.010)} &\textbf{1.517(0.030} &\textbf{1.345(0.012)}& \textbf{3.233(0.065)} & \textbf{1.253(0.013)}&\textbf{2.849(0.064)} \\
& Single AE & 1.119(0.006)&2.266(0.031)&1.013(0.008)&1.943(0.033) &1.427(0.012)& 3.502(0.065)& 1.386(0.013)&3.332(0.068) \\
& AEMC-NE & 1.137(0.007)&2.386(0.036)&1.024(0.008)&2.007(0.036) &1.421(0.012)& 3.483(0.064)& 1.372(0.014)&3.281(0.069) \\
& DiD & 1.127(0.007)&2.371(0.037)&1.069(0.008)&2.182(0.038) &1.408(0.013)& 3.437(0.066)& 1.401(0.014)&3.402(0.071) \\
& Vert-Reg & 1.100(0.006)&2.087(0.024)&1.055(0.007)&1.940(0.028) &1.485(0.012)& 4.054(0.077)& 1.439(0.012)&3.737(0.074) \\
& MC-NNM & 1.044(0.006)&1.998(0.028)&0.939(0.005)&1.649(0.020) &1.373(0.012)& 3.266(0.061)& 1.366(0.013)&3.233(0.067)\\
\midrule
\multirow{6}{*}{\textbf{poly}}
& CoDEAL & 0.833(0.007)&1.290(0.028)&\textbf{0.715(0.007)} &\textbf{0.913(0.035)} &\textbf{1.182(0.013)}& \textbf{2.507(0.059)} & \textbf{1.024(0.013)}&\textbf{2.086(0.062)} \\
& Single AE & 0.946(0.008)&1.704(0.036)&0.854(0.008)&1.395(0.033) &1.345(0.014)& 3.088(0.068)& 1.261(0.014)&2.760(0.066) \\
& AEMC-NE & 1.060(0.008)&2.038(0.037)&0.868(0.008)&1.461(0.036) &1.348(0.012)& 3.088(0.063)& 1.245(0.014)&2.701(0.065) \\
& DiD & 0.906(0.008)&1.596(0.036)&0.921(0.009)&1.653(0.037) &1.235(0.013)& 2.659(0.062)& 1.279(0.014)&2.832(0.067) \\
& Vert-Reg & \textbf{0.756(0.006)} &\textbf{0.934(0.017)} &0.800(0.007)&1.078(0.023) &1.281(0.012)& 2.983(0.072)& 1.216(0.013 &2.790(0.070) \\
& MC-NNM & 0.847(0.006)&1.305(0.032)&0.747(0.005)&1.014(0.017) & 1.290(0.013)& 2.763(0.057)& 1.238(0.013)&2.652(0.063) \\
\midrule
\multirow{6}{*}{\textbf{ReLU}}
& CoDEAL & 0.946(0.009)&1.702(0.042)&\textbf{0.717(0.008)} &\textbf{0.916(0.025)} &$1.269 (0.013)$& $2.817 (0.064)$& \textbf{1.121(0.012)} &\textbf{2.269(0.054)} \\
& Single AE & 1.020(0.010)&1.971(0.046)&0.858(0.010)&1.407(0.039) &$1.329 (0.013)$& $3.044 (0.064)$& 1.222(0.014)&2.617(0.062) \\
& AEMC-NE & 1.288(0.011)&2.982(0.053)&0.868(0.010)&1.471(0.043) &$1.539 (0.014)$& $4.010 (0.075)$& 1.235(0.013)&2.660(0.062) \\
& DiD & 0.960(0.009)&1.777(0.046)&0.922(0.010)&1.670(0.047) &$1.275 (0.013)$& $2.819 (0.062)$& 1.269(0.014)&2.791(0.064) \\
& Vert-Reg & \textbf{0.779(0.006)} &\textbf{1.013(0.020} &0.782(0.006)&1.027(0.021) &\textbf{1.204(0.012)}& $2.846 (0.069)$& 1.199(0.012)&2.753(0.067) \\
& MC-NNM & 0.967(0.009)&1.893(0.035)&0.764(0.009)&1.073(0.024) &1.235(0.012)& \textbf{2.639(0.057)} & 1.227(0.013)&2.608(0.060)\\
\bottomrule
\label{tab:lrnn1}
\end{tabular}
}
\end{table}
\begin{table}[p]
\small
\centering
\caption{Comparisons of imputation accuracy (measured by MAE and MSE) with covariate effects removed by \textbf{LR} and \textbf{DNN} under various covariate generations in Config.2. Results are averaged over 50 replications. Numbers in parentheses are standard errors.}
\resizebox{\textwidth}{!}{
\begin{tabular}{llccccccll}
\toprule
\multirow{2}{*}{\textbf{Covariate}} & \multirow{2}{*}{\textbf{Method}} & \multicolumn{2}{c}{\textbf{LR}, Linear Factor}&\multicolumn{2}{c}{\textbf{NN}, Linear Factor}& \multicolumn{2}{c}{\textbf{LR}, Sine Factor}& \multicolumn{2}{c}{\textbf{NN}, Sine Factor}\\
& &MAE& MSE&MAE& MSE&MAE& MSE& \text{ }\text{ }\text{ }\text{ }MAE&\text{ }\text{ }\text{ }\text{ }\text{ }MSE\\
\midrule
\multirow{6}{*}{\textbf{No covariate}}
& CoDEAL & \textbf{0.507(0.004)} &\textbf{0.425(0.010)} &\textbf{0.507(0.004)} &\textbf{0.425(0.010)} &\textbf{0.779(0.007)} & \textbf{1.173(0.027)} & \textbf{0.779(0.007)} &\textbf{1.173(0.027)} \\
& Single AE & 0.622(0.006)&0.761(0.021)&0.622(0.006)&0.761(0.021)&0.963(0.010)& 1.772(0.038)& 0.963(0.010)&1.772(0.038)\\
& AEMC-NE & 0.647(0.006)&0.850(0.023)&0.647(0.006)&0.850(0.023)&1.004(0.009)& 1.884(0.036)& 1.004(0.009)&1.884(0.036)\\
& DiD & 0.713(0.007)&1.043(0.028)&0.713(0.007)&1.043(0.028)&1.085(0.010)& 2.163(0.041)& 1.085(0.010)&2.163(0.041)\\
& Vert-Reg & 0.651(0.002)&0.671(0.004)&0.651(0.002)&0.671(0.004)&0.865(0.007)& 1.422(0.034)& 0.865(0.007)&1.422(0.034)\\
& MC-NNM & 0.535(0.005)&0.495(0.013)&0.535(0.005)&0.495(0.013)&1.020(0.009)& 1.856(0.034)& 1.020(0.009)&1.856(0.034)\\
\midrule
\multirow{6}{*}{\textbf{Matrix Linear}}
& CoDEAL & \textbf{0.510(0.005)} &\textbf{0.436(0.013)} &\textbf{0.609(0.004)} &\textbf{0.622(0.010)} &\textbf{0.781(0.007)} & \textbf{1.177(0.025)} & \textbf{0.896(0.008)}&\textbf{1.456(0.029)} \\
& Single AE & 0.645(0.007)&0.807(0.022)&0.777(0.007)&1.147(0.024) &1.007(0.009)& 1.885(0.037)& 1.146(0.010)&2.330(0.042) \\
& AEMC-NE & 0.682(0.007)&0.919(0.024)&0.759(0.007)&1.093(0.025) &1.051(0.009)& 2.018(0.037)& 1.127(0.009)&2.259(0.040) \\
& DiD & 0.736(0.008)&1.093(0.029)&0.810(0.007)&1.262(0.029) &1.122(0.010)& 2.275(0.041)& 1.191(0.009)&2.503(0.042) \\
& Vert-Reg & 0.680(0.003)&0.733(0.006)&0.761(0.004)&0.936(0.011) &0.917(0.007)& 1.550(0.031)& 1.110(0.007)&2.024(0.027) \\
& MC-NNM & 0.537(0.005)&0.498(0.013)&0.650(0.005)&0.736(0.015) &1.025(0.009)& 1.889(0.033)& 1.106(0.009)&2.147(0.037) \\
\midrule
\multirow{6}{*}{\textbf{Vector Linear}}
& CoDEAL & \textbf{0.508(0.005)} &\textbf{0.428(0.010)} &\textbf{0.589(0.004)}&\textbf{0.574(0.010)}&\textbf{0.784(0.007)} & \textbf{1.184(0.026)} & \textbf{0.883(0.007)}&\textbf{1.426(0.025)}\\
& Single AE & 0.659(0.006)&0.843(0.018)&0.772(0.007)&1.135(0.024)&1.010(0.008)& 1.890(0.031)& 1.147(0.009)&2.333(0.036)\\
& AEMC-NE & 0.700(0.009)&0.952(0.024)&0.761(0.008)&1.099(0.026)&1.071(0.010)& 2.073(0.036)& 1.130(0.009)&2.265(0.034)\\
& DiD & 0.737(0.006)&1.093(0.023)&0.800(0.007)&1.234(0.025) &1.125(0.009)& 2.279(0.036)& 1.186(0.009)&2.478(0.036)\\
& Vert-Reg & 0.675(0.003)&0.724(0.007)&0.742(0.004)&0.884(0.010)&0.915(0.007)& 1.555(0.031)& 0.995(0.007)&1.768(0.030)\\
& MC-NNM & 0.536(0.005)&0.499(0.012)&0.640(0.006)&0.714(0.017) &1.019(0.010)& 1.870(0.036)& 1.099(0.008)&2.124(0.030) \\
\midrule
\multirow{6}{*}{\textbf{tanh}}
& CoDEAL & 0.700(0.007)&0.920(0.024)&\textbf{0.644(0.005)} &\textbf{0.704(0.015)} &1.012(0.008)& 1.913(0.034)& \textbf{0.904(0.008)}&\textbf{1.492(0.030)} \\
& Single AE & 0.880(0.009)&1.442(0.032)&0.795(0.007)&1.190(0.027) &1.227(0.011)& 2.638(0.047)& 1.180(0.010)&2.446(0.043) \\
& AEMC-NE & 0.902(0.008)&1.489(0.028)&0.771(0.007)&1.120(0.026) &1.243(0.009)& 2.681(0.041)& 1.160(0.010)&2.365(0.042) \\
& DiD & 0.758(0.007)&1.139(0.028)&0.828(0.008)&1.309(0.030) &1.135(0.010)& 2.313(0.041)& 1.224(0.010)&2.613(0.044) \\
& Vert-Reg & 0.731(0.003)&0.849(0.008)&0.770(0.005)&0.968(0.015) &\textbf{0.968(0.007)} & \textbf{1.698(0.032)} & 1.048(0.007)&1.921(0.029) \\
& MC-NNM & \textbf{0.570(0.005)} &\textbf{0.555(0.013)} &0.697(0.006)&0.859(0.019) &1.136(0.009)& 2.317(0.036)& 1.145(0.009)&2.281(0.040) \\
\midrule
\multirow{6}{*}{\textbf{log}}
& CoDEAL & 1.030(0.005)&1.867(0.017)&\textbf{0.878(0.004)} &\textbf{1.385(0.013)} &\textbf{1.236(0.008)} & \textbf{2.683(0.038)} & \textbf{1.107(0.007)} &\textbf{2.264(0.030)} \\
& Single AE & 1.085(0.006)&2.125(0.026)&0.985(0.006)&1.808(0.026) &1.367(0.008)& 3.241(0.040)& 1.338(0.009)&3.116(0.042) \\
& AEMC-NE & 1.078(0.006)&2.115(0.028)&0.968(0.007)&1.751(0.027) &1.360(0.008)& 3.210(0.040)& 1.319(0.009)&3.034(0.042) \\
& DiD & 1.079(0.007)&2.135(0.030)&1.013(0.007)&1.919(0.030) &1.382(0.009)& 3.312(0.043)& 1.374(0.009)&3.268(0.045) \\
& Vert-Reg & 1.275(0.005)&2.722(0.020)&1.108(0.005)&2.079(0.017) &1.516(0.008)& 3.956(0.049)& 1.388(0.008)&3.301(0.039) \\
& MC-NNM & \textbf{0.999(0.005)} &\textbf{1.795(0.021)} &0.907(0.005)&1.503(0.019) &1.306(0.008 & 2.957(0.038)& 1.305(0.009)&2.951(0.040) \\
\midrule
\multirow{6}{*}{\textbf{poly}}
& CoDEAL & 0.798(0.006)&1.128(0.022)&\textbf{0.656(0.005)} &\textbf{0.725(0.013)} &\textbf{1.082(0.009)}& \textbf{2.105(0.038)} & \textbf{0.962(0.008)}&\textbf{1.640(0.030)} \\
& Single AE & 0.965(0.007)&1.672(0.029)&0.808(0.007)&1.223(0.027) &1.298(0.009)& 2.895(0.044)& 1.187(0.009)&2.463(0.041) \\
& AEMC-NE & 0.999(0.007)&1.756(0.027)&0.789(0.007)&1.162(0.025) &1.313(0.010 & 2.934(0.044)& 1.173(0.009)&2.412(0.041) \\
& DiD & 0.848(0.007)&1.337(0.029)&0.841(0.008)&1.339(0.030) &1.197(0.009)& 2.510(0.041)& 1.234(0.010)&2.650(0.044) \\
& Vert-Reg & 0.809(0.004)&1.045(0.010)&0.791(0.005)&1.017(0.015) &1.159(0.007)& 2.288(0.033)& 1.072(0.008)&2.001(0.032) \\
& MC-NNM & \textbf{0.710(0.005)} &\textbf{0.853(0.014)} &0.714(0.006)&0.896(0.018) & 1.204(0.009)& 2.321(0.036)& 1.158(0.009)&2.324(0.039) \\
\midrule
\multirow{6}{*}{\textbf{ReLU}}
& CoDEAL & 0.771(0.007)&1.205(0.027)&\textbf{0.612(0.005)} &\textbf{0.629(0.012)} &$1.269 (0.013)$ & $2.817 (0.064)$ & \textbf{0.902(0.008)} &\textbf{1.478(0.029)} \\
& Single AE & 1.078(0.011)&2.110(0.043)&0.764(0.006)&1.069(0.020) &$1.329 (0.013)$ & $3.044 (0.064)$ & $1.121 (0.010)$ &$2.232 (0.040)$\\
& AEMC-NE & 1.234(0.011)&2.692(0.048)&0.784(0.007)&1.148(0.024) &$1.539 (0.014)$ & $4.010 (0.075)$ & $1.156 (0.009)$ &$2.355 (0.040)$\\
& DiD & 0.868(0.007)&1.407(0.029)&0.840(0.007)&1.334(0.029) &$1.275 (0.013)$ & $2.819 (0.062)$ & $1.220 (0.010)$ &$2.603 (0.042)$\\
& Vert-Reg & 0.808(0.004)&1.051(0.011)&0.783(0.004)&0.992(0.011) &\textbf{1.204(0.012)} & $2.846 (0.069)$ & $1.049 (0.007)$ &$1.938 (0.030)$\\
& MC-NNM & \textbf{0.766(0.006)} &\textbf{1.034(0.019)} &0.651(0.005)&0.737(0.016) &$1.235 (0.012)$ & \textbf{2.639(0.057)} & $1.107 (0.009)$ &$2.149 (0.036)$\\
\bottomrule
\label{tab:lrnn2}
\end{tabular}
}
\end{table}
\vspace{-2ex}
Table \ref{tab:lrnn1} and Table \ref{tab:lrnn2} report the comparisons of MAE and MSE using the two covariate effect removal methods in Config.1 and Config.2, respectively. Under each configuration, we consider both linear factor effects and nonlinear Sine factor effects, with all six types of covariate effects.
As shown in Tables \ref{tab:lrnn1} and \ref{tab:lrnn2}, LR-based removal and DNN-based removal behave very differently depending on the underlying true types of covariate effects. In cases where the covariate effect is strictly linear, LR often outperforms NN by a small margin,
indicating that a simple linear adjustment suffices for purely linear effects.
In contrast, whenever the covariate effects enter nonlinearly (e.g.\ tanh, log, polynomial, or ReLU), the proposed DNN covariate removal achieves substantially lower errors, with MAE reduced by 10–30 \% and MSE by similar amounts relative to LR.
These results confirm that LR is adequate for linear covariate effects but that the DNN approach is essential to capture complex, nonlinear covariate–outcome relationships.
\vspace{-2ex}
\subsection{Evaluation on the Performance of CoDEAL in Causal Matrix Completion}
\label{append:simu-CoDEAL-performance}
\vspace{-1ex}
Table \ref{tab:simu-CoDeal} reports the comparisons of prediction accuracy between the proposed CoDEAL and five benchmark methods (detailed in Section \ref{append:simu-benchmark}) in Config.1-4. Under each configuration, we consider all four types of factor effects (listed in Table \ref{tab:data_gen}), with linear covariate effects and nonlinear tanh covariate effects.
\begin{onehalfspace}
\begin{table}[p]
\centering
\caption{\small Comparisons of imputation accuracy (measured by MAE and MSE) under various factor and covariate generations in Config.1-4. Results are averaged over 50 replications. Numbers in parentheses are standard errors.}
\resizebox{\textwidth}{!}{
\begin{tabular}{llccccccll}
\toprule
\multirow{2}{*}{\textbf{Factor}}& \multirow{2}{*}{\textbf{Method}} & \multicolumn{2}{c}{Config. 1, linear covariate}& \multicolumn{2}{c}{Config. 1, nonlinear covariate}& \multicolumn{2}{c}{Config. 2, linear covariate}& \multicolumn{2}{c}{Config. 2, nonlinear covariate}\\
& &MAE& MSE&MAE& MSE &MAE& MSE&\text{ }\text{ }\text{ }\text{ } MAE&\text{ }\text{ }\text{ }\text{ }MSE \\
\midrule
\multirow{6}{*}{\textbf{Linear}}
& CoDEAL & \textbf{0.646(0.009)} &\textbf{0.722(0.023)} &\textbf{0.701(0.008)}&\textbf{0.871(0.031)}&\textbf{0.609(0.004)} & \textbf{0.622(0.010)} & \textbf{0.644(0.005)} &\textbf{0.704(0.015)} \\
& Single AE & 0.821(0.010)&1.320(0.039)&0.851(0.008)&1.383(0.033)&0.777(0.007)& 1.147(0.024)& 0.795(0.007)&1.190(0.027)\\
& AEMC-NE & 0.838(0.009)&1.394(0.041)&0.862(0.009)&1.445(0.036)&0.759(0.007)& 1.093(0.025)& 0.771(0.007)&1.120(0.026)\\
& DiD & 0.892(0.010)&1.587(0.045)&0.916(0.009)&1.642(0.037)&0.810(0.007)& 1.262(0.029)& 0.828(0.008)&1.309(0.030)\\
& Vert-Reg & 0.757(0.006)&0.957(0.019)&0.798(0.007)&1.075(0.024)&0.761(0.004)& 0.936(0.011)& 0.770(0.005)&0.968(0.015)\\
& MC-NNM & 0.699(0.006)&0.895(0.020)&0.735(0.005)&1.007(0.017)&0.650(0.005)& 0.736(0.015)& 0.697(0.006)&0.859(0.019)\\
\midrule
\multirow{6}{*}{\textbf{Sine}}& CoDEAL & \textbf{1.078(0.012)} &\textbf{2.136(0.055)} &\textbf{1.117(0.011)}&\textbf{ 2.238(0.050)}&\textbf{0.894(0.007)}& \textbf{1.454(0.028)} & \textbf{0.945(0.008)}&\textbf{1.593(0.030)}\\
& Single AE & 1.223(0.013)&2.618(0.061)&1.250(0.012)&2.705(0.053)&1.142(0.009)& 2.312(0.039)& 1.185(0.009)&2.458(0.039)\\
& AEMC-NE & 1.209(0.013)&2.566(0.061)&1.235(0.012)&2.645(0.053)&1.128(0.009)& 2.256(0.039)& 1.167(0.010)&2.389(0.041)\\
& DiD & 1.244(0.013)&2.697(0.063)&1.272(0.012)&2.785(0.055)&1.191(0.009)& 2.498(0.042)& 1.234(0.009)&2.650(0.041)\\
& Vert-Reg & 1.169(0.012)&2.633(0.063)&1.208(0.011)&2.736(0.059)&1.012(0.007)& 1.829(0.031)& 1.054(0.008)& 1.935(0.034)\\
& MC-NNM & 1.201(0.013)&2.511(0.059)&1.231(0.012)&2.610(0.053)&1.105(0.009)& 2.142(0.036)& 1.153(0.010)&2.305(0.040)\\
\midrule
\multirow{6}{*}{\textbf{Polynomial}}
& CoDEAL & \textbf{0.908(0.013)}&\textbf{2.250(0.100)}&\textbf{0.952(0.015)}&\textbf{2.477(0.126)}&\textbf{0.772(0.010)}& \textbf{1.399(0.076)}& \textbf{0.807(0.011)}&\textbf{1.539(0.070)} \\
& Single AE & 0.962(0.013)&2.727(0.116)&0.997(0.016)&3.014(0.164)&0.820(0.010)& 1.759(0.087)& 0.854(0.011)&1.964(0.084)\\
& AEMC-NE & 0.997(0.014)&3.198(0.142)&1.038(0.017)&3.500(0.192)&0.856(0.011)& 2.132(0.110)& 0.887(0.012)&2.339(0.111)\\
& DiD & 0.995(0.013)&2.708(0.111)&1.025(0.016)&2.939(0.144)&0.854(0.010)& 1.860(0.083)& 0.886(0.011)&2.065(0.092)\\
& Vert-Reg & 0.953(0.013)&1.906(0.076)&1.019(0.016)&2.327(0.113)&0.855(0.007)& 1.289(0.032)& 0.884(0.010)&1.456(0.042)\\
& MC-NNM & 0.951(0.013)&2.419(0.101)&0.989(0.015)&2.627(0.130)& 0.819(0.009)& 1.654(0.07)& 0.851(0.011)&1.839(0.078)\\
\midrule
\multirow{6}{*}{\textbf{ReLU MLP}}
& CoDEAL & \textbf{0.797(0.009)}&\textbf{1.097(0.027)}&\textbf{0.838(0.010)}& \textbf{1.217(0.032)} &\textbf{0.698(0.005)}& \textbf{0.820(0.014)}& \textbf{0.743(0.006)}&\textbf{0.929(0.015)}\\
& Single AE & 0.989(0.009)&1.712(0.034)&1.026(0.011)& 1.845(0.044)&0.891(0.007)& 1.375(0.022)& 0.930(0.008)&1.502(0.028)\\
& AEMC-NE & 1.045(0.010)&1.919(0.038)&1.088(0.011)& 2.084(0.048)&0.959(0.007)& 1.605(0.026)& 0.982(0.008)&1.679(0.029)\\
& DiD & 0.930(0.008)&1.499(0.028)&0.963(0.010)& 1.616(0.036)&0.842(0.005)& 1.218(0.017)& 0.870(0.007)&1.303(0.021)\\
& Vert-Reg & 0.888(0.010)&1.315(0.035)&0.958(0.013)& 1.544(0.046)&0.847(0.005)& 1.166(0.014)& 0.894(0.006)&1.302(0.020)\\
& MC-NNM & 0.843(0.007)&1.227(0.025)&0.892(0.009)& 1.453(0.031)&0.767(0.005)& 1.001(0.013)& 0.806(0.006)&1.110(0.019)\\
\toprule
\multirow{2}{*}{\textbf{Factor}}& \multirow{2}{*}{\textbf{Method}}& \multicolumn{2}{c}{Config. 3, $r=5$}& \multicolumn{2}{c}{Config. 3, $r=10$} & \multicolumn{2}{c}{Config. 4, $r=5$}& \multicolumn{2}{c}{Config. 4, $r=10$} \\
& & MAE & MSE & MAE & MSE& MAE & MSE& MAE & MSE \\
\midrule
\multirow{6}{*}{\textbf{Linear}}
& CoDEAL &0.614(0.007) &0.721(0.023) & 0.625(0.005)& 0.754(0.018)&\textbf{0.532(0.003)}&\textbf{0.506(0.009)}&\textbf{0.560(0.004)}&\textbf{0.591(0.012)} \\
& Single AE &0.789(0.011) &1.263(0.044) & 0.795(0.008)& 1.270(0.033)&0.693(0.005) &0.948(0.017) &0.714(0.006) &1.009(0.023) \\
& AEMC-NE &0.731(0.010) &1.125(0.042) & 0.730(0.007)& 1.110(0.029)&0.647(0.005) &0.846(0.018) &0.660(0.005) &0.886(0.021) \\
& DiD &0.801(0.011) &1.341(0.047) & 0.802(0.008)& 1.326(0.033)&0.711(0.005) &1.031(0.021) &0.726(0.006) &1.081(0.025) \\
& Vert-Reg & 0.940(0.005)&6.880(2.378)& 0.865(0.023)&8.238(2.625) &0.691(0.003) &0.888(0.012) &0.750(0.016) &3.071(0.914) \\
& MC-NNM &\textbf{0.599(0.007)} &\textbf{0.685(0.023)} & \textbf{0.616(0.006)}& \textbf{0.733(0.019)}&0.537(0.003) &0.520(0.009) &0.565(0.004) &0.596(0.014) \\
\hline
\multirow{6}{*}{\textbf{Sine}}
& CoDEAL &\textbf{0.973(0.010)} &\textbf{1.798(0.041)} &\textbf{0.996(0.010)}&\textbf{1.907(0.039)} &\textbf{0.850(0.007)} &\textbf{1.389(0.024)} &\textbf{0.869(0.008)} &\textbf{1.463(0.028)} \\
& Single AE &1.087(0.011) &2.138(0.047) &1.122(0.011) &2.288(0.047) &1.057(0.008) &2.050(0.030) &1.061(0.008) &2.082(0.035) \\
& AEMC-NE &1.055(0.012) &2.029(0.047) &1.078(0.011) &2.128(0.044) &1.011(0.008) &1.890(0.029) &1.008(0.009) &1.897(0.034) \\
& DiD &1.110(0.012) &2.209(0.049) &1.146(0.012) &2.361(0.048) &1.091(0.008) &2.172(0.031) &1.088(0.009) &2.189(0.037) \\
& Vert-Reg & 1.549(0.090)&8.954(2.462)& 1.524(0.123)&15.689(4.057) &0.974(0.009) &2.153(0.059) &1.047(0.013) &4.384(0.544) \\
& MC-NNM &1.066(0.011) &2.033(0.046) &1.097(0.011) &2.164(0.043) &0.998(0.008) &1.819(0.031) &1.003(0.008) &1.858(0.032) \\
\hline
\multirow{6}{*}{\textbf{Polynomial}}
& CoDEAL &\textbf{0.785(0.011)} &\textbf{2.146(0.112)} &\textbf{0.811(0.010)} &\textbf{2.379(0.161)} &\textbf{0.662(0.007)} &\textbf{1.311(0.073)} &\textbf{0.699(0.007)} &\textbf{1.530(0.064)} \\
& Single AE &0.895(0.015) &2.423(0.117) &0.926(0.016) &2.710(0.172) &0.739(0.011) &1.582(0.086) &0.785(0.010) &1.798(0.071) \\
& AEMC-NE &0.860(0.014) &2.921(0.160) &0.897(0.013) &3.273(0.220) &0.720(0.009) &1.889(0.108) &0.764(0.009) &2.166(0.093) \\
& DiD &0.887(0.013) &2.422(0.116) &0.917(0.014) &2.693(0.165) &0.739(0.009) &1.613(0.082) &0.773(0.008) &1.801(0.068) \\
& Vert-Reg & 1.155(0.035)&6.771(1.168)& 1.131(0.032)&5.530(0.666) &0.788(0.005) &1.246(0.028) &0.881(0.014) &3.254(0.526) \\
& MC-NNM &0.839(0.012) &2.117(0.104) &0.865(0.012) &2.308(0.129) &0.703(0.008) &1.397(0.071) &0.735(0.007) &1.569(0.058) \\
\hline
\multirow{6}{*}{\textbf{ReLU MLP}}
& CoDEAL &\textbf{0.632(0.004)} &\textbf{0.886(0.012)} &\textbf{0.733(0.004)} &\textbf{0.953(0.013)} &\textbf{0.633(0.004)} &\textbf{0.698(0.011)}&\textbf{0.646(0.004)} &\textbf{0.730(0.011)} \\
& Single AE &0.876(0.005) &1.338(0.017) &0.913(0.006) &1.448(0.019) &0.797(0.005) &1.108(0.015) &0.809(0.005) &1.142(0.015) \\
& AEMC-NE &0.934(0.006) &1.574(0.021) &0.968(0.007) &1.677(0.024) &0.860(0.006) &1.330(0.023) &0.878(0.006) &1.381(0.022) \\
& DiD &0.817(0.004) &1.163(0.013) &0.842(0.004) &1.238(0.014) &0.742(0.004) &0.962(0.013) &0.748(0.005) &0.976(0.014) \\
& Vert-Reg & 1.045(0.031)&5.157(0.846)& 0.948(0.026)&10.387(6.292) &0.755(0.003) &1.060(0.015) &1.070(0.253) &8.399(2.377)\\
& MC-NNM &0.721(0.003) &0.897(0.009) &0.750(0.003) &0.976(0.010) &0.656(0.004) &0.741(0.011) &0.662(0.004) &0.756(0.011) \\
\bottomrule
\label{tab:simu-CoDeal}
\end{tabular}
}
\end{table}
\end{onehalfspace}
\bigskip
\noindent \textbf{CoDEAL in Four-Block Design.} As shown in Config.1 and Config.2 in Table \ref{tab:simu-CoDeal}, CoDEAL consistently \textbf{achieves the lowest MAE and MSE across nearly all settings} of factor effects and covariate effects, and achieves substantial improvements in MAE and MSE relative to all competing methods. Under the simplest linear factor model, CoDEAL reduces MAE by approximately 5–10\% and MSE by about 10–20\% compared to MC-NNM, the second-best performing method, and significantly outperforms other baseline methods by even larger margins. For scenarios involving nonlinear factor models, CoDEAL continues to show superior performance, typically surpassing MC-NNM by approximately 5–20\% reduction in MAE and 10–20\% reduction in MSE.
Overall, these findings clearly demonstrate that CoDEAL’s integration of DNN-based covariate adjustment and deep autoencoder-based latent factor extraction is robust and highly effective across diverse linear and nonlinear data-generating settings.
\bigskip
\noindent \textbf{CoDEAL in Staggered Adoption Design.} As Config.3 and Config.4 of Table \ref{tab:simu-CoDeal} demonstrate, CoDEAL achieves substantial percentage gains in MAE and MSE relative to all competing methods across staggered adoption designs.
Under \emph{Linear} factors, CoDEAL CoDEAL \textbf{{ties MC-NNM within 3\% while reduces MAE by 16–35\% over other benchmarks}}. In staggered design with nonlinear factors, CoDEAL delivers evident advantages. Especially when under \emph{ReLU}-based factors, CoDEAL delivers its largest margins, with \textbf{{MAE reductions of 12–40\% and MSE reductions of 15–41\% versus all competing approaches}}.
Finally, as $r$ increases from 5 to 10, the missingness pattern grows more complex and extensive, so all methods show some drop in accuracy. However, CoDEAL’s errors only suffer a slight increase, whereas simpler methods like Vert-Reg show a dramatic rise in both MAE and MSE when missingness increases. Above results reemphasize CoDEAL’s consistent superiority and robustness in staggered‐treatment settings.
\bigskip
\noindent \textbf{In summary}, CoDEAL achieves performance comparable to MC-NNM in linear settings, while exhibiting substantial improvements over all considered methods in the presence of nonlinear effects.
By integrating DNN-based covariate adjustments with deep multi-output-AE–based nonlinear factor analysis, CoDEAL is robustly effective across a wide range of linear and nonlinear scenarios.
\section{Real Data: Oxford COVID-19 Government Response Tracker (OxCGRT)} \label{sec:realdata}
\vspace{-1ex}
In this section, we illustrate our methodology using data from the OxCGRT \citep{Hale2021}. The data is publicly available at \href{https://github.com/OxCGRT}{https://github.com/OxCGRT}. OxCGRT systematically collects panel data on the timing and intensity of 24 different government interventions, such as school closures, travel restrictions, and vaccination mandates, across over 185 countries from January 2020 through 2022.
We examine the impact of two specific policy interventions: (1) the \textbf{implementation of mandatory vaccination}, and (2) the \textbf{termination of internal travel restrictions} (results in Appendix \ref{append: covid19}). Our analysis covers their impacts on COVID-19 confirmed cases and deaths across 64 states and territories in the United States and Canada from April 16, 2021, to December 11, 2021. This time frame aligns with the period following the detection of the Delta variant and precedes the rapid emergence of the Omicron variant in North America. In our analysis, we use four OxCGRT indices (overall government response, containment and health, stringency, and economic support), averaged over the time frame, as covariates, resulting in data dimensions of $N=64$, $T=240$, and $P=4$.
The policy is adopted in different states at various times, forming a staggered adoption design (visualized in Figure \ref{fig:MV_mask}).
Figure \ref{fig:large_MV} plots the comparisons of the total confirmed cases and deaths across policy-in-effect states between observed results (with policy executed) and estimated counterfactuals (if no policy).
As shown in Figure \ref{fig:large_MV}, if the mandatory vaccination policy was not implemented, there would be an increase in both confirmed cases and deaths, with a notably larger impact on mortality. This observation aligns with recent research highlighting a greater reduction in mortality than in infection rates following vaccination \citep{hernandez2023impact,xu2023systematic,mesleEtAl2024estimated,wu2023long}. Although the visual differences in the log-transformed graphs appear modest, converting these back to their original scales reveals substantial impacts: on December 11, 2021, the actual cumulative counts of confirmed cases and deaths for states with mandatory vaccination policies were ${26.57\text{ }million, 0.39\text{ }million}$, respectively. In contrast, our estimates suggest that without these policies, the counts would have risen to ${29.81\text{ }million(\uparrow 12.2\%), 0.52\text{ }million(\uparrow 34.6\%)}$.
\begin{figure}[p]
\centering
\includegraphics[width=0.8\textwidth]{Figures-Finalized/missing_MA.png}
\captionsetup{justification=centering}
\caption{Indicator matrix of the implementation of mandatory vaccination policy. \\ Blocks with darker color refer to policy executed.}
\label{fig:MV_mask}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[width=1\textwidth]{Figures-Finalized/MA1.png}
\caption{Comparison of the total confirmed cases (Left) and deaths (Right) across policy-in-effect states between observed results with the mandatory vaccination policy executed and estimated counterfactuals if no policy.
Values are reported in ${\log_{10}}$ scale, and plots are visualized with 14-day moving averages.}
\label{fig:large_MV}
\end{figure}
\section{Discussion of Limitations and Future Extensions} \label{sec:discussion}
In this paper, we propose CoDEAL, a deep-learning-based causal matrix completion method and a unifying framework for estimating heterogeneous causal effects in panel data models.
One limitation of CoDEAL is that it is only designed for binary treatment settings. Causal panel data with multiple treatments is an open research area with limited existing work in the literature \citep{abadie2021using}, yet this topic is gradually gaining attention due to its growing relevance in practical applications \citep{agarwal2020synthetic,squires2022causal}.
To extend CoDEAL to multiple treatment settings, one possible direction is to adopt a similar strategy to that mentioned by \citep{agarwal2020synthetic}, in which the treatment type is incorporated as an additional data dimension in the analysis. Together with the original panel data structures, the inclusion of treatment dimension forms a three-dimensional tensor. Methodologically, \textbf{\emph{with different choices of factor models and advanced neural network architectures, the proposed CoDEAL framework can be further generalized to accommodate various data types with more complex data structures}}, such as tensors, images, and networks. For example, CoDEAL can be naturally extended to three-dimensional or higher-order tensors by involving tensor factorization methodologies \citep{han2022tensor, chen2024estimation, zhou2025factor,chen2025diffusion}, providing a more general approach compared to matrix factor analysis \citep{bai2021matrix, luo2022inverse,yu2022nonparametric,yu2024dynamic}. That being said, to extend to multiple treatments, how to handle the increased level of missingness brought by the treatment assignments in a higher-dimensional tensor remains challenging and requires careful modeling strategies. We leave the further investigation on multiple treatment settings as future work.
{
\bibliographystyle{plainnat}
\bibliography{ref}
}