\section{Introduction}
Unmeasured confounding is a well-known threat to valid causal inference from
observational data. An approach that is sometimes used in practice to assess
residual confounding bias, is to check whether \textit{known null effects} can
be recovered free of bias, by evaluating whether the exposure or treatment of
interest is found to be associated with a so-called negative control outcome
(NCO), upon adjusting for measured confounders \citep{Rosenbaum1989, Lipsitch2010, Shi2020}. 
An observed variable is said to be a valid\textit{ negative control outcome }or\textit{
}more broadly\textit{, }an\textit{ outcome confounding proxy,} to the extent
that it is associated with hidden factors confounding the exposure-outcome
relationship in view, although not directly impacted by the exposure.
Therefore, a negative control outcome which is empirically associated with the
exposure, might suggest the presence of residual confounding. In the event
such an association is present, a natural question is whether the negative
control outcome can be used for bias correction. The most well-established
negative control outcome approach for debiasing observational causal effect
estimates is the difference-in-differences approach (DiD) 
\citep{Snow1856,Snow1936,Angrist2009,Caniglia2020}. 
In fact, DiD may be
viewed as directly leveraging the pre-treatment outcome as a negative control
outcome since it cannot logically be causally impacted by treatment. Identification then follows from an additive equi-confounding assumption that the unmeasured confounder
association with the post-treatment outcome of interest matches that with the pre-treatment
outcome on the additive scale \citep{Sofer2016}. The baseline outcome in DiD is
thus implicitly assumed to be a valid NCO, and equi-confounding is equivalent
to the so-called \textit{parallel trends assumption}, that the average trends
in treatment-free potential outcomes for treatment and untreated units are
parallel. In practice, equi-confounding or equivalently parallel trends may
not be reasonable for a number of reasons, including if the outcome trend is
also impacted by an unmeasured common cause with the treatment. Furthermore,
additive equi-confounding may not be realistic as a broader debiasing method in
non-DiD settings where the NCO is not necessarily a pre-treatment measurement of the outcome of interest, but is instead a post-treatment measurement of a different type of outcome (and might therefore have support on a different scale than the outcome of interest has). 
To address these potential limitations of additive
equi-confounding, \citet{TT2013_COCA} introduced the Control Outcome
Calibration Approach (COCA) as a simple yet formal counterfactual NCO approach to
debias causal effect estimates in observational analyses. At its core, COCA
essentially treats the negative control outcome variable as a proxy
measurement for the treatment-free potential outcome, which therefore is associated with the latter, and which becomes independent of the treatment assignment mechanism, upon conditioning on the treatment-free
counterfactual\ outcome. As the treatment-free potential outcome can be viewed
as an ultimate source of unmeasured confounding, this assumption formalizes
the idea that, as a relevant proxy for the source of residual confounding, the
NCO\ would be made irrelevant for the treatment assignment mechanism if one
were to hypothetically condition on the underlying potential outcome.

For identification and inference for a continuous outcome, the original COCA approach of \citet{TT2013_COCA} 
involves the correct specification of a regression model for the NCO,
conditional on the treatment-free potential outcome and measured confounders, together with a rank-preserving structural model which effectively assumes a constant individual-level treatment
effect. In this paper, we develop a nonparametric COCA identification
framework for the average causal effect for the treated, which equally applies
irrespective of the nature of the primary outcome, whether binary, continuous, or
polytomous. Importantly, as we show, the proposed COCA\ identification framework completely obviates the need
for rank preservation, therefore accommodating arbitrary degrees of effect
heterogeneity across units. Relatedly, an alternative counterfactual approach
named \textit{proximal causal inference} has recently developed in causal
inference literature \citep{Miao2016_arxiv,Miao2018,TT2020_Proximal}, which leverages a pair of
negative treatment and outcome control variables, or more broadly treatment and outcome confounding proxies, to nonparametrically identify
treatment causal effects subject to residual confounding without invoking a
rank-preservation assumption. Importantly, while proximal causal inference
relies on two proxies for causal identification, in contrast, COCA is a
single-proxy control approach which therefore may present practical
advantages. For estimation and inference, we introduce three
strategies to implement COCA\ which improve on prior methods: (i) an
\textit{extended propensity score approach}, (ii) a so-called\textit{ outcome
calibration bridge function approach}, and (iii) \textit{a doubly robust
approach} which carefully combines approaches (i) and (ii) and remains
unbiased for the effect of treatment on the treated, provided that either
approach is also unbiased, without necessarily knowing which method might not
be unbiased. Finally, we illustrate the methods with an application evaluating the causal effect of a Zika outbreak on birth rate in Brazil, and
we conclude with possible extensions to our methods and a brief discussion.

\section{Notation and Brief Review of COCA}

Consider an observational study where, as represented in Figure \ref{fig-1}, one has
observed an outcome variable $Y$, a binary treatment $A$ whose causal effect
on $Y$ is of interest. 
Suppressing measured covariates until a later Section, 
we are concerned that as displayed succinctly in the figure with the bow arc, the
association between $A$ and $Y$ is confounded by hidden factors. 
\begin{figure}[!htp] 
\centering
\scalebox{1}{
\begin{tikzpicture}
\tikzset{line width=1.5pt, outer sep=0pt,
ell/.style={draw,fill=white, inner sep=2pt,
line width=1.5pt},
swig vsplit={gap=5pt,
inner line width right=0.5pt}};
\node[name=A, ell, ellipse] at (0,0){$A$};
\node[name=Y, ell, ellipse] at (4,0) {$Y$};
\begin{scope}[>={Stealth[black]}, every edge/.style={draw=black,very thick}]
    \path [dashed,<->] (A) edge[bend left=90] (Y);
    \path [->] (A) edge (Y);    
\end{scope}
\end{tikzpicture} }
	\caption{A Graphical Illustration of A Simple Causal Model.}
	\label{fig-1}
\end{figure}	

Throughout, $Y^{a}$ denotes the potential outcome or counterfactual, had possibly contrary to fact, the
exposure been set to $a$ by an external hypothetical intervention.
Furthermore, throughout, we also make the consistency assumption that
$Y^{a}=Y$ if $A=a.$ Hereafter, we aim to make inferences about the causal effect of
treatment on the treated (ETT):
\[
\psi_{0}=E\left(  Y^{a=1}-Y^{a=0}|A=1\right).
\]
Under consistency, $E\left(  Y^{a=1}|A=1\right)$ is identified by $E\left(Y|A=1\right)$; to identify $E\left(  Y^{a=0}|A=1\right)$ requires additional assumptions. Standard methods often resort to the assumption of no unmeasured confounding:%
\[
Y^{a=0}\bot\bot A
\]
a strong assumption we do not make. Instead, we suppose that one has measured
a valid negative control outcome $W,$ which is known a priori to satisfy two
conditions.\ First $W$ must be associated with the mechanism leading to
residual confounding in the sense that it stands for a proxy for the latter;
secondly, $W$ cannot be causally impacted by the treatment of interest. In
this paper, inspired by\ \citet{TT2013_COCA}, we formalized these
conditions as followed:%
\begin{equation}
W^{a}=W;\label{NC.1}%
\end{equation}
where $W^{a\text{ }}$is the potential negative control outcome under an
external intervention that sets $A=a;$ and
\begin{align}
&  W\bot\not \bot Y^{a=0};\label{NC.2}\\
&  W\bot\bot A|Y^{a=0}.\label{NC.3}%
\end{align}
Condition $\left(  \ref{NC.1}\right)  $ encodes the key assumption of a known
null causal effect of the treatment on the NCO\ in potential outcome notation. Condition 2 encodes that W is  relevant for predicting the treatment-free potential outcome of interest.  
Condition 3 states that W is independent of A conditional on treatment-free potential outcome.
These conditions formally encode the assumption that $W$ is a valid proxy for the
treatment-free potential outcome, a source of residual confounding
bias; that is $W$ is only associated with the
treatment mechanism to the extent that it is associated with the confounding mechanism captured by the potential outcome.\ We illustrate these NCO assumptions with the causal graph displayed in figure \ref{fig: Graphical COCA}. The red arrows indicate the deterministic relationship defining the observed outcome in terms of potential outcomes and treatment variable by the consistency assumption. The missing arrows on the graph formally encode the core conditional independence conditions (1)-(3), which can be verified by d-separation.  
\begin{figure}[!htp]
\centering 
\scalebox{1}{
\begin{tikzpicture}
\tikzset{line width=1.5pt, outer sep=0pt,
ell/.style={draw,fill=white, inner sep=2pt,
line width=1.5pt},
swig vsplit={gap=5pt,
inner line width right=0.5pt}};
\node[name=A, ell, ellipse] at (0,0){$\, A\, $};
\node[name=Y0, ell, ellipse] at (4,2) {$\potY{0}{}$}  ;
\node[name=Y1, ell, ellipse] at (1,2) {$\potY{1}{}$}  ;
\node[name=W, ell, ellipse] at (5,0) {$\, W\, $};
\node[name=Y, ell, ellipse] at (2.5,0) {$\, Y\, $};
\begin{scope}[>={Stealth[black]}, every edge/.style={draw=black,very thick}]
        \path [dashed,<->] (Y1) edge[bend left=90] (Y0);
	\path [->] (Y0) edge (A);
        \path [->] (Y1) edge (A);
        \path [->] (Y0) edge (W);    
        \path [-{Stealth[red]}] (Y1) edge[red] (Y); 
        \path [-{Stealth[red]}] (Y0) edge[red] (Y); 
        \path [-{Stealth[red]}] (A) edge[red] (Y); 
\end{scope}
\end{tikzpicture} }
    \caption{\small A Graphical Illustration of the Assumptions for Control Outcome Calibration Approach. Red arrows depict the deterministic relationship between $Y$ and $(\potY{1}{},\potY{0}{},A)$, as established by the consistency assumption $Y=A Y^{a=1}+ (1-A)Y^{a=0}$.}
    \label{fig: Graphical COCA}
\end{figure} 

\ For identification and estimation in the case of a continuous exposure,
\citet{TT2013_COCA} further assumed the rank-preserving structural model:%
\begin{equation}
Y=Y^{a=0}+\psi_{0}A,\label{rank-preserving}%
\end{equation}
which by consistency, implies a constant individual-level causal effect
$\psi_{0}=Y^{a=1}-Y^{a=0}.$ Under this model, he noted that upon defining
$Y(\psi)=Y-\psi A,$ then one can deduce from conditions $\left(
\ref{NC.1}\right)  $-$\left(  \ref{NC.3}\right)  $ that
\[
W\bot\bot A|Y(\psi)
\]
if and only if $\psi=\psi_{0},$ in which case, given equation 4, $Y^{a=0}=Y(\psi_{0})=Y-\psi_{0} A$ which motivates a regression-based
implementation of COCA, that entails searching for the parameter value of
$\psi=\psi_{0}$ such that $E\left(  W|A,Y(\psi)\right)  =E\left(
W|Y(\psi)\right)  .$ A straightforward implementation of the approach uses
linear models whereby for each value of $\psi$ on a sufficiently fine grid$,$ one obtains an
estimate of the regression model $E\left(
W|A,Y(\psi)\right)  =\beta_{1}+\beta_{2}A+\beta_{3}Y(\psi)$  using ordinary least squares (OLS), with estimated
coefficients $\left(  \hat{\beta}_{1}(\psi),\widehat{\beta}_{2}(\psi
),\widehat{\beta}_{3}(\psi)\right)  .$ Then a 95\% confidence interval for
$\psi_{0}$ consists of all values of $\psi$ for which a valid test of the null
hypothesis $\beta_{2}(\psi)=0$ fails to reject at the 0.05 type 1 error level.
\ Such hypothesis test might be performed by verifying whether the interval
$\widehat{\beta}_{2}(\psi)\pm1.96\widehat{\mathrm{SE}}\left(  \widehat{\beta
}_{2}(\psi)\right)  $ covers $0,$ with $\widehat{\mathrm{SE}}\left(
\widehat{\beta}_{2}(\psi)\right)  $ the OLS$\ $estimate of the standard error
of $\widehat{\beta}_{2}(\psi)$. \citet{TT2013_COCA} also describes a
potentially simpler one-shot approach which fits a single regression $E\left(
W|A,Y\right)  =\beta_{1}+\beta_{2}^{\ast}A+\beta_{3}Y$ via OLS where
$\beta_{2}^{\ast}=-\beta_{3}/\psi_{0},$ in which case $\widehat{\psi}%
_{0}=-\widehat{\beta}_{2}^{\ast}/\widehat{\beta}_{3}$ where $\left(
\widehat{\beta}_{2}^{\ast},\widehat{\beta}_{3}\right)  $ are OLS estimates; a
corresponding standard error estimator of $\widehat{\psi}_{0}$ is given in 
Section \ref{sec:supp:COCA_SE} of the Appendix for convenience. 
Though practically convenient, validity of either approach relies on both correct specification of the linear model for
$W$ given $Y^{a=0}$ $,$ and on the rank-preserving structural model, which may
be biologically implausible. In the following, we describe alternative methods
aimed at addressing these limitations.

\section{COCA via Extended Propensity Score Weighting}

In order to introduce the approach, consider the extended propensity score
function (EPS):
\[
\pi\left(  y\right)  =\Pr\left(  A=1|Y^{a=0}=y\right)
\]
which makes explicit the fact that in presence of unmeasured confounding, the
treatment mechanism will generally depend on the treatment-free potential
outcome. Next, we note that were $\pi\left(  y\right)  $ known, identification
of the average treatment-free potential outcome in the treated would then be
empirically identified by the expression
\begin{equation}
E\left(  Y^{a=0}|A=1\right)  =\ \frac{E\left(  \left(  1-A\right)
Y\frac{\pi\left(  Y\right)  }{1-\pi\left(  Y\right)  }\right)  }{E\left(
\left(  1-A\right)  \frac{\pi\left(  Y\right)  }{1-\pi\left(  Y\right)
}\right)  }\label{EPS ID}%
\end{equation}
provided that, as we assume throughout, positivity holds, i.e. $0<\pi(y)<1$ for all $y$; see Appendix \ref{sec:supp:proof:EPSID} for the details. Therefore, the average causal effect of treatment on the treated would be
identified by%
\begin{equation}
E\left( Y^{a=1}-\ Y^{a=0}|A=1\right)  =E\left( Y|A=1\right)
-\frac{E\left(  \left(  1-A\right)  Y\frac{\pi\left(  Y\right)  }{1-\pi\left(
Y\right)  }\right)  }{E\left(  \left(  1-A\right)  \frac{\pi\left(  Y\right)
}{1-\pi\left(  Y\right)  }\right)  }\label{EPS ATT}%
\end{equation}
As $\pi\left(  y\right)  $ is unknown, we next demonstrate how the
NCO\ assumption can be leveraged to identify the latter quantity. Let
$p\left(  w\right)  =\Pr\left(  A=1|W=w\right)$. The proposed approach to
identify $\pi$ is based on the following equality, which we prove in Appendix \ref{sec:supp:proof:EPSeq}:
\begin{result}      \label{result-1}
Under consistency and Assumptions $\left(  \ref{NC.1}%
\right)  -\left(  \ref{NC.3}\right)$, we have that
\begin{align}
\frac{p\left(  w\right)  }{1-p\left(  w\right)  } &  =E\left\{  \frac
{\pi\left(  Y\right)  }{1-\pi\left(  Y\right)  } \, \bigg| \, W=w,A=0\right\}
\label{EPS  equation}\\
&  =\sum_{y}\frac{\pi\left(  y\right)  }{1-\pi\left(  y\right)  }g\left(
y|w,A=0\right)  \nonumber
\end{align}
where $g\left(  y|w,A=0\right)  $ is the law of $Y$ given $Z$ and $A$,
evaluated at $(y,w,a=0)$.
\end{result}
\noindent Result \ref{result-1} provides an expression relating the EPS of
interest to the observed data distribution, as of the three quantities
involved in the expression $\pi,p,g$; two are uniquely determined
by the observed data, mainly $p$ and $g$, which are then related to the
unknown function of interest $\pi$ in Result \ref{result-1}. 
Equation  $\left(  \ref{EPS  equation}\right)  $ is known in mathematics as an
integral equation, more precisely a Fredholm integral equation of the first
kind. We then have the following identification result.

\begin{result}      \label{result-2}
If Assumptions $\left(  \ref{NC.1}\right)  -\left(
\ref{NC.3}\right)  $  hold, and the integral equation $\left(  \ref{EPS  equation}%
\right)$ in Result \ref{result-1} admits a unique solution, then $\pi\left(  Y\right)
$ is nonparametrically identified from the observed data by solving \eqref{EPS equation}; and
 $E\left( Y^{a=0}|A=1\right)  $ is nonparametrically identified by
$\left(  \ref{EPS ID}\right)  .$
\end{result}

Sufficient conditions for the existence, and uniqueness of a solution to
such an equation are well-studied by mathematicians, such conditions were
recently discussed by statisticians in the context of proximal causal
inference \citep{Miao2016_arxiv,Miao2018,TT2020_Proximal},
a recently proposed alternative approach for bias correction leveraging a pair of
negative control variables or proxy variables; a detailed comparison of COCA with proximal causal inference is relegated to Section \ref{sec:Discussion}. Intuitively, the condition that the integral equation admits a unique solution essentially
requires that $W$ is sufficiently relevant for $Y^{a=0}$ in the sense that for
any variation in the latter, there is corresponding variation in the former.
This assumption is akin to the assumption of relevance in the context of
instrumental variable methodology, which states that variation in the instrument should
induce variation in the treatment.
Importantly Result \ref{result-2} applies whether $Y$ is binary, continuous or
polytomous, provided that $W$ is sufficiently relevant for $Y^{a=0}$, to ensure that equation
$\left(  \ref{EPS ID}\right)  $ admits a solution. Also, the result
obviates the need for a rank-preserving structural model, and delivers fully
nonparametric identification of the causal effect of treatment on the
treated.  We now describe simple modeling approaches for estimation and
inference via EPS. 

The first approach entails a priori specifying a parametric logistic regression
model for the EPS, say:%
\begin{align*}
& \log\frac{\Pr\left(  A=1|Y^{a=0}=y\right)  }{\Pr\left(  A=0|Y^{a=0}%
=y\right)  }\\
& =\log\frac{\pi\left(  y;\mathbf{\alpha}^{\ast}\right)  }{1-\pi\left(
y;\mathbf{\alpha}^{\ast}\right)  }=\alpha_{0}+\mathbf{\alpha}_{y}%
^{\ast^{\prime}}S_{y}\left(  y\right)  ;\\
\mathbf{\alpha}  & \mathbf{=}\left(  \alpha_{0},\mathbf{\alpha}_{y}^{^{\prime
}}\right)  ^{\prime};
\end{align*}
where $S_{y}\left(  y\right)  $ is a user-specified sufficient statistic for
the log-odds ratio association between $Y^{a=0}$ and $A,$ evaluated at
$y.$\ $\ $ Options for specifying the latter includes both more parsimonious
specifications, e.g. $S_{y}\left(  y\right)  =y$ as a well as more flexible
specifications, say  $S_{y}\left(  y\right)  =\left(  y,y^{2},...,y^{v}\right)
$ or $S_{y}\left(  y\right)  =\left(  I\left(  y\leq y_{1}\right)  ,I\left(
y_{1}<y\leq y_{2}\right)  ,...,I\left(  y_{K}<y\right)  \right)  $ where
$y_{k}$ is the $k^{th}$ user-specified percentile of $Y;$ Alternative
nonparametric smoothing techniques, e.g. Nearest neighbor, kernel smoothing, splines or wavelets might also be considered. While such logistic regression model specifications for the EPS might look familiar, estimation and inference about its unknown parameter
$\mathbf{\alpha}^{\ast}$ via maximum likelihood is complicated by the fact
that $Y^{a=0}$ $\ $is only observed among untreated units with $A=0.$ Instead,
to find an empirical estimate of $\mathbf{\alpha}^{\ast}$ which solves the equation in Result \ref{result-1}, we propose a so-called Generalized Method of Moments (GMM)\ inspired by the
following population moment equation:
\[
E\left\{  \left[  \frac{\left\{  1-A\right\}  }{1-\pi\left(  Y;\mathbf{\alpha
}\right)  }-1\right]  r_{a}\left(  W\right)  \right\}  =0,\text{ for any
user-specified function }r_{a}\text{;}%
\]
which under our conditions, admits a unique solution $\mathbf{\alpha=\alpha
}^{\ast}.$ Note that though the integral equation in Result \ref{result-1} involves the propensity score given $W$ on the left hand-side, GMM obviates the need for estimating the latter, an appealing property. In practice, we might recommend setting $r_{a}\left(  w\right)
=\left(  1,S_{y}^{\prime}\left(  w\right)  \right)  ^{\prime}$ so that the
dimension of $r$ matches that of $\mathbf{\alpha.}$ A key advantage of GMM\ is
that it straightforwardly accommodates settings where multiple NCOs might be
available, say $W_{1},...,W_{k},$ in which case one might specify $r_{a}\left(
w_{1},...,w_{k}\right)  =\left(  1,S_{y}^{\prime}\left(  w_{1}\right)
,S_{y}^{\prime}\left(  w_{2}\right)  ,...,S_{y}^{\prime}\left(  w_{k}\right)
\right)  ^{\prime}$; note that the dimension of $r$ might easily exceed
that of $\mathbf{\alpha;}$ a situation GMM\ is well-suited to handle. Details
about the empirical implementation of GMM\ are relegated to Section \ref{sec:supp:GMM} of the Appendix,
along with corresponding large sample behavior of the resulting EPS\ estimator of the
treatment effect on the treated obtained by substituting the estimated
EPS\ and empirical means in expression $\left(  \ref{EPS ATT}\right)  $.  

\subsection{Estimation via COCA Bridge Function}

In this Section, we introduce an alternative nonparametric identification and estimation approach which does not rely on modeling the EPS, but instead relies on the existence of a so-called outcome
bridge function, which we introduce with the following assumption.

\begin{condition}   \label{condition-1}
For any $y,$ there exist a function (possibly nonlinear), $b(w)$ that satisfies the following equation%
\begin{equation}
y=E\left\{  b\left(  W\right)  |Y=y,A=0\right\}  .\label{Outcome COCA}%
\end{equation}
\end{condition}
 

Intuitively, upon noting that under assumptions\ $\left(  \ref{NC.1}\right)
-\left(  \ref{NC.3}\right)  ,$ Condition \ref{condition-1} can equivalently be stated in terms
of potential outcomes%
\begin{equation}
Y^{a=0}=E\left\{  b\left(  W\right)  |Y^{a=0}\right\}
\label{Counterfactual bridge}%
\end{equation}
which essentially formalizes the idea that $W$ is a sufficiently relevant
proxy for the potential outcome $Y^{a=0}$ if there exist a (possibly unknown)
transformation of $W$ whose conditional expectation given $Y^{a=0},$ recovers
$Y^{a=0}.$ As $b(w)$ provides a bridge between the observed data equation
$\left(  \ref{Outcome COCA}\right)  ,$ and its potential outcome counterpart
$\left(  \ref{Counterfactual bridge}\right)  ,$ we aptly refer to $b\left(
W\right)  $ as a \textit{COCA confounding bridge function}. Note that classical measurement error is a special case of the equation in the display above in which case $b$ is the identity map and $W=Y^{a=0}+e$ where $e$ is an independent mean zero error.  The condition can therefore be viewed as nonparametric generalization of classical measurement error which allows $W$ and $Y$ to be of arbitrary nature and does not assume the error to be unbiased on the additive scale.  We further  illustrate the assumption in the case of binary $W$ and $Y$.
As shown in Appendix \ref{sec:supp:proof:Binary}, 
in which case, the following $b(W)$ satisfies \eqref{Outcome COCA} provided that $\Pr\left(  W=1|Y=1,A=0\right)  -\Pr\left(  W=1|Y=0,A=0\right)  \not =0$:
\begin{align*}
    b(W)
    =
    \frac{ -(1-W) \cdot \Pr\left(  W=1|Y=0,A=0\right) + W \cdot  \Pr\left(  W=0|Y=0,A=0\right)  }{\Pr\left(
W=1|Y=1,A=0\right)  -\Pr\left(  W=1|Y=0,A=0\right)  }
\end{align*}
This latter condition is equivalently
expressed under\ $\left(  \ref{NC.1}\right)  -\left(  \ref{NC.3}\right)  $ as%
\[
\Pr\left(  W=1|Y^{a=0}=1\right)  \not =\Pr\left(  W=1|Y^{a=0}=0\right)
\]
encoding the requirement that $W$ cannot be independent of $Y^{a=0}.$ Beyond
the binary case, for more general outcome types,
Condition \ref{condition-1} likewise formally defines a Fredholm integral equation of the
first kind, for which sufficient conditions for existence of a solution are
well characterized in functional analysis textbooks; we again
refer the reader to \citep{Miao2016_arxiv}.   We are now ready to state our result, which we prove in Appendix \ref{sec:supp:proof:OutcomeBridge}:
\begin{result}  \label{result-3}
    Suppose that Assumptions $\left(  \ref{NC.1}\right)  -\left(
\ref{NC.3}\right)  $  and Condition \ref{condition-1} hold, then 
\begin{equation}
E\left( Y^{a=0}|A=1\right)  =E\left(  b^{\ast}\left(  W\right)
|A=1\right)  \label{Outcome bridge id}%
\end{equation}
where $b^{\ast}$ satisfies $\left(  \ref{Outcome COCA}%
\right)  .$
\end{result}
 

In the binary example discussed above where $b^{*}(w)$ was uniquely identified, we have that
\begin{align*}
E\left(  Y^{a=0}|A=1\right)    & =E\left(  b^{\ast}\left(  W\right)
|A=1\right)  \\
& =\frac{\Pr\left(  W=1|A=1\right)  -\Pr\left(  W=1|Y=0,A=0\right)  }%
{\Pr\left(  W=1|Y=1,A=0\right)  -\Pr\left(  W=1|Y=0,A=0\right)  }%
\end{align*}


We briefly highlight a key feature of the above result reflected in its proof,
which is that $b^{\ast}\left(  W\right)  $ need not be uniquely identified by
equation  $\left(  \ref{Outcome COCA}\right)  ,$ and that any such solution
leads to a unique value for $E\left(  Y^{a=0}|A=1\right)  .$ Interestingly, the identifying formula in the display above was also obtained by \cite{TT2013_COCA} in the binary case, although he did not emphasize the key role of the bridge function as a general framework for identification beyond the binary case. 

For estimation, we primarily consider the special case where $b^{\ast}\left(
W\right)  $ is in fact uniquely identified by equation $\left(
\ref{Outcome COCA}\right)  .$ In this case, a practical approach might entail
specifying a flexible model, say  $b\left(  W;\eta\right)  =\eta^{\prime}%
m_{w}\left(  W\right)  $ for user-specified function $m_{w}\left(  W\right)
,$ assuming that $b^{\ast}\left(  W\right)  =\eta^{\ast\prime}m_{w}\left(
W\right)  $ for an unknown unique value $\eta^{\ast}.$ Similar to
$S_{y}\left(  y\right)  $ in the previous Section, $m_{w}\left(  W\right)  $
can likewise account for nonlinearities by specifying specific basis
functions, such as polynomials, splines, or wavelets. We then note that a
standard least-squares regression $Y$ on $W$ would generally fail to recover
$\eta^{\ast}$ even if the model is correctly specified. In other words, under
our assumptions it will generally be the case that $b^{\ast}\left(  W\right)
\neq E\left(Y|A=0,W\right)  ,$ as this would in fact imply that  $E\left(
E\left(  Y|A=0,W\right)  |A=0,Y=y\right)  =y$ which does not follow from our assumptions.
A correct approach in fact is based on the following consequence of our
assumptions:%
\[
E\left\{  \left(  1-A\right)  \left(  \eta^{\prime}m_{w}\left(  W\right)
-Y\right)  r_{W}\left(  Y\right)  \right\}  =0,\text{ for any user-specified
function }r_{W}\text{;}%
\]
if and only if $\eta=\eta^{\ast}.$ 
In practice, analogous to the prior
Section, it might be convenient to take $r_{W}\left(  Y\right)  =m_{w}\left(
Y\right)  ,$ which yields
\[
\eta^{\ast}=E\left(  m_{w}\left(  W\right)  m_{W}^{\prime}\left(  Y\right)
|A=0\right)  ^{-1}E\left(  Ym_{W}\left(  Y\right)  |A=0\right)
\]
assuming $E\left(  m_{w}\left(  W\right)  m_{W}^{\prime}\left(  Y\right)
|A=0\right)  $ is invertible.\ Details for empirical implementation and
inference based on this approach is relegated to Appendix \ref{sec:supp:GMM}, 
together with inference about the ETT obtained upon substituting unknown population
parameters with their empirical analogs in equation $\left(
\ref{Outcome bridge id}\right)  .$ 

\subsection{ Doubly Robust COCA}

We have thus far described two separate COCA estimators: (i) an EPS approach, which relies on correct specification of a model
for the exposure mechanism, and (ii) an outcome bridge COCA estimator, which
instead relies on correct specification of a model for a bridge function which
if correctly specified, should be conditionally unbiased for the
treatment-free counterfactual outcome. Thus, misspecification of either model
will likely result in biased inferences about the ETT. As in most practical
situations one cannot be sure that either model is correctly specified, the
most one can probably hope for is to obtain unbiased inferences about the causal effect of interest
if either model is correct, without a priori knowledge of which model if any is
incorrect. An estimator with such a property is said to be doubly robust
\citep{Scharfstein1999, Lunceford2004, Bang2005}. 
Such doubly robust estimators have previously been proposed under
unconfoundedness conditions 
\citep{Bang2005, Tsiatis2006, TTRR2010}, 
and more recently for proximal causal
inference leveraging a pair of proxies \citep{Cui2020}. 
Here, we propose a doubly
robust COCA estimator of the effect of treatment on the treated, which
involves obtaining a doubly robust estimator of the average exposure-free
counterfactual outcome in the treated parameter $E\left(  Y^{a=0}   |A=1\right), $ which we now give in the following result; see Appendix \ref{sec:supp:proof:DR} for the proof:
\begin{result} \label{result-4}
    Under Assumptions $\left(  \ref{NC.1}\right)  -\left(
\ref{NC.3}\right)  ,$ we have that %
\[
E\left(  Y^{a=0}  |A=1\right)  =E\left\{  Ab\left(
W;\eta\right)  +\frac{\left(  1-A\right)  \pi\left(  Y;\mathbf{\alpha}\right)
}{  1-\pi\left(  Y;\mathbf{\alpha}\right)     }\left\{  Y-b\left(
W;\eta\right)  \right\}  \right\}  / \Pr \left(  A=1\right)
\]
if either $E\left(  b\left(  W;\eta\right)  |A=0,Y=y\right)  =y$ or
$\pi\left(  y;\mathbf{\alpha}\right)  =\Pr\left(  A=1|Y^{a=0}=y\right)  .$
\end{result}
Result \ref{result-4} motivates a doubly robust estimator which we obtain by substituting
unknown quantities $\pi\left(  Y;\mathbf{\alpha}\right)  $ and $b\left(
W;\eta\right)  $ with corresponding estimators obtained using methods
described in previous sections. Details are relegated to Appendix \ref{sec:supp:GMM}, together
with large sample behavior of the resulting estimator, and corresponding
inference. 





\section{Data Application: Zika Virus Outbreak in Brazil}            \label{sec:Data}

The Zika virus, which can be transmitted from a pregnant woman to her fetus, can cause serious brain abnormalities, including microcephaly (i.e., an abnormally small head) \citep{Zika2016}. Brazil is one of the countries hardest hit by the Zika virus. In particular, the outbreak in 2015 resulted in  over 200,000 cases in Brazil by 2016 \citep{Zika2018}. As a result, many prior works \citep{Zika2018_3, Zika2018_2, Zika2022, TTPR2023} asked whether the Zika virus outbreak caused a drop in birth rates. 

We re-analyze the dataset analyzed in \citet{Zika2022} and \citet{TTPR2023}. In the dataset, we focused on municipalities in two states of Brazil, Pernambuco and Rio Grande do Sul, which are northeastern and southernmost states. Out of the 1248 cases of microcephaly that occurred in Brazil by November 28, 2015, 51.8\% (646 cases) were reported in Pernambuco (PE), while there were no reported cases in Rio Grande do Sul (RS) \citepalias{Zika2015}, which shows that PE was severely impacted by the Zika virus outbreak, while RS was minimally affected. Based on their epidemiologic histories, we define the two states as treated and control groups, respectively. We use the post-epidemic municipality-level birth rate in 2016 as the outcome $Y$, where the birth rate is defined as the total number of live human births per 1,000 persons. We use the pre-epidemic municipality-level birth rate in 2014 as the outcome proxy (i.e., NCO), denoted $W$. To be a valid proxy, the birth rate in 2014 must satisfy the following three conditions: (i) birth rate in 2014 cannot be causally impacted by the Zika virus epidemic which occurred in 2015, (ii) birth rate in 2014 is   correlated with what the birth rate in 2016 would have been had there not been a Zika virus epidemic, and (iii) birth rate in 2014 would become irrelevant and therefore independent of a municipality's Zika epidemic status, upon conditioning on its Zika virus epidemic-free potential birth rate in 2016.
The first two conditions are uncontroversial, while the third condition largely relies on the extent to which pre-epidemic birth rates can accurately be viewed as a proxy for the counterfactual birth rate had the pandemic not occurred, and as such would not further be predictive of whether the municipality experienced a high rate of Zika virus incidence, conditional on the regions epidemic free counterfactual birth rate in 2016. Although one might consider this last assumption reasonable, ultimately, it is empirically untestable without making an alternative assumption. Nevertheless, in the next Section we describe a straightforward sensitivity analysis to evaluate the extent to which violation of the assumption might impact inference. We focus the analysis on 673 municipalities with complete data on the outcome, outcome proxy, and treatment, where 185 and 488 municipalities belong to PE and RS, respectively. The crude estimand $E(\obsY{1} \, \big| \, A=1) - E(\obsY{1} \, \big| \, A=0)$ is estimated to be equal to $3.384$, suggesting that municipalities in the PE region (with higher incidence of Zika virus) experienced a higher birth rate than RS regions in 2016 during the Zika virus outbreak. An immediate concern is that this crude association between $A$ and $Y$ might be subject to significant confounding bias, leading us to conduct two separate analyses geared at addressing residual confounding bias; the proposed COCA methods, which we compare with a more standard difference-in-differences analysis. Thus, we estimate the average treatment effect on the treated using the three COCA estimators. In particular, because both $Y$ and $W$ are continuous, we specify the EPS and the outcome bridge function by specifying flexible models for discrete versions of $Y$ and $W$, respectively, where we vary the number of discrete bins as 5 and 10, leading to a form of histogram-like nonparametric estimator, respectively. As previously mentioned, this approach yields an increasingly flexible model specification for both the outcome bridge function and the EPS functions. In addition, implementation of COCA with either EPS or the bridge function might benefit from a regularization via $\ell_2$-penalization (i.e., ridge penalty) for numerical stability. In fact, EPS estimation in the application appeared to benefit from such a small penalty, while bridge function estimation did not appear to require one; see Appendix \ref{sec:supp:data} for details.

Table \ref{Tab-1} summarizes corresponding results. We find that the six estimates (two per COCA estimator) vary between $-1.695$ and $-2.602$, meaning between $1.695$ and $2.602$ birth per 1,000 persons were reduced in PE due to the Zika virus outbreak, an empirical finding better aligned with the scientific hypothesis that Zika may likely adversely impact exposed populations birthrates. Compared to the crude estimate of 3.384, the negative effect estimates indeed provide compelling evidence of potential confounding. We also obtain an estimate using the difference-in-difference estimator under a standard parallel trends assumption (e.g., \citet{CardKrueger1994, Angrist2009}), which yields a considerably smaller effect estimate of -1.191; noting that the DiD estimator requires the assumption of equi-confounding of the $A-Y$ association in the pre- and post-periods, while our proposed estimator does not (but instead requires conditions (i)-(iii) outlined above). Regardless of the estimator, all 7 estimates appear to be consistent with the anticipated adverse causal impact of the Zika Virus epidemic, although the doubly robust COCA estimator produced a more efficient estimator than the other COCA estimators. Consequently, we conclude that based on inferences aimed at accounting for confounding (DiD and COCA), the Zika virus outbreak likely led to a decline in the birthrate of affected regions in Brazil, which agrees with similar findings in the literature \citep{Zika2018_3,Zika2018_2,Zika2022,TTPR2023}.


\begin{table}[!htp]
\renewcommand{\arraystretch}{1.2} \centering
		\setlength{\tabcolsep}{7pt}
\begin{tabular}{ccccc}
\hline
\multicolumn{5}{|c|}{COCA}                                                                                                                                                                                                           \\ \hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Number of bins}} & \multicolumn{1}{c|}{\multirow{2}{*}{Statistic}} & \multicolumn{3}{c|}{Estimator Type}                                                                                        \\ \cline{3-5} 
\multicolumn{1}{|c|}{}                                & \multicolumn{1}{c|}{}                           & \multicolumn{1}{c|}{EPS}                 & \multicolumn{1}{c|}{ORB}               & \multicolumn{1}{c|}{DR}                \\ \hline
  \multicolumn{1}{|c|}{\multirow{2}{*}{5}}  &    \multicolumn{1}{c|}{Estimate}  &            \multicolumn{1}{c|}{-1.695}  &            \multicolumn{1}{c|}{-2.470}  &            \multicolumn{1}{c|}{-1.858} \\ \cline{2-5}
                     \multicolumn{1}{|c|}{}  &  \multicolumn{1}{c|}{95$\%$ CI}  &  \multicolumn{1}{c|}{(-2.585, -0.804)}  &  \multicolumn{1}{c|}{(-3.846, -1.094)}  &  \multicolumn{1}{c|}{(-2.864, -0.851)}      \\ \hline
 \multicolumn{1}{|c|}{\multirow{2}{*}{10}}  &    \multicolumn{1}{c|}{Estimate}  &            \multicolumn{1}{c|}{-1.941}  &            \multicolumn{1}{c|}{-2.602}  &            \multicolumn{1}{c|}{-1.817} \\ \cline{2-5}
                     \multicolumn{1}{|c|}{}  &  \multicolumn{1}{c|}{95$\%$ CI}  &  \multicolumn{1}{c|}{(-2.757, -1.126)}  &  \multicolumn{1}{c|}{(-4.618, -0.586)}  &  \multicolumn{1}{c|}{(-2.946, -0.688)}      \\ \hline
                                                      &                                                 &                                          &                                        &                                        \\[-0.4cm] \hline
\multicolumn{5}{|c|}{Difference-in-Differences under Parallel Trends}   \\ \hline
\multicolumn{2}{|c|}{Estimate}                                                                         & \multicolumn{3}{c|}{-1.191}                                                                                                \\ \hline
\multicolumn{2}{|c|}{95$\%$ CI}                                                                         & \multicolumn{3}{c|}{(-1.507,  -0.876)}                                                                                     \\ \hline
&                                                 &                                          &                                        &                                        \\[-0.4cm] \hline
\multicolumn{5}{|c|}{Crude Estimate $E(Y_1|A=1) - E(Y_1|A=0)$}   \\ \hline
\multicolumn{2}{|c|}{Estimate}                                                                         & \multicolumn{3}{c|}{3.384}                                                                                                \\ \hline
\multicolumn{2}{|c|}{95$\%$ CI}                                                                         & \multicolumn{3}{c|}{(2.953, 3.816)}                                                                                     \\ \hline
\end{tabular}
\caption{Summary of Data Analysis. EPS, ORB, DR stand for the extended propensity score, outcome regression bridge, doubly robust estimators of the average treatment effect on the treated.}
\label{Tab-1}

\end{table}

\bigskip





\section{Discussion and Possible Extensions}        \label{sec:Discussion}
We have described a COCA nonparametric identification framework, therefore extending previous results of \citet{TT2013_COCA} 
to a more general setting accommodating outcomes of arbitrary nature and obviating the need for an assumption of constant treatment effects, i.e. rank preservation. We have proposed three estimation strategies including a doubly robust method which has appealing robustness properties. 
While we have ignored the role of measured confounders (say denoted by $X$), when available such covariates can readily be included in our analyses by specifying their functional form via either $S_y(W,X)$ and $r_a(W,X)$ in EPS approach or their analogs in the confounding bridge-based approach or respectively in the doubly robust approach. 
Interestingly, the COCA central identifying assumption, that conditioning on the treatment-free counterfactual would in principle shield the treatment assignment from any association with the negative control outcome is isomorphic to an analogous assumption 
in the missing data literature where an outcome might be missing not at random; however, a fully observed so-called \emph{shadow variable} (the missing data analog of an NCO) reasonably assumed to be conditionally independent of the missing data process given the value of the potentially missing outcome. For example, \cite{Zahner1992} considered a study of the children’s mental health evaluated through their
teachers’ assessments in Connecticut. However, the data for the teachers’ assessments are subject to nonignorable missingness. As a proxy of the teacher’s assessment, a separate parent
report is available for all children in this study. The parent report is likely to be correlated with
the teacher’s assessment, but is unlikely to be related to the teacher’s response rate given the
teacher’s assessment and fully observed covariates. Hence, the parental assessment is regarded
as a shadow variable for the teacher's assessment in this study.  The literature on shadow variables is fast-growing \citep{DH2010, Kott2014, Wang2014, MLTTG2015, MTT2016, LMTT2021},
the methods developed in this paper have close parallels to shadow variable counterparts in this literature. This connection to shadow variables is particularly salient in the bridge function approach in the context of multiple NCOs, in which case, the outcome bridge function might not be uniquely defined, thus complicating inference; fortunately the methods developed by \citep{LMTT2021} for the analogous shadow variable setting directly apply to our setting.  Although the more general framework of shadow variables is not directly discussed in the current work, we refer the interested reader to this latter work for further details. It is worth noting that the doubly robust estimator proposed in this paper appears to be completely knew, and different from those of \citet{MLTTG2015, MTT2016, LMTT2021}  and therefore may also be of use in  shadow variable applications. Likewise, the doubly robust estimators proposed by the latter works can equally be applied to current the COCA setting as an alternative approach.  

Additionally, as mentioned in the previous Section, the key assumption that conditioning on the treatment-free potential outcome, would in principle make the NCO or outcome proxy irrelevant to treatment mechanism is ultimately untestable, and may in certain settings not hold exactly. In fact, this would be the case if the NCO were in fact explicitly used in assigning the treatment in which case the assumption might be violated. In order to address such eventuality, the analyst might consider several candidate proxies/NCOs when available, and may even perform an over-identification test, by inspecting the extent to which the estimated causal effect depends on the choice of proxy. Alternatively, a sensitivity analysis might also be performed to evaluate the potential impact of a hypothesize departure from the assumption. To perform such sensitivity analysis, one might consider the slightly modified EPS model obtained by including the NCO into the EPS to encode a departure from the assumption, say:
\begin{align*}
& \log\frac{\Pr\left(  A=1|Y^{a=0}=y, W=w\right)  }{\Pr\left(  A=0|Y^{a=0}%
=y, W=w\right)  } 
=\alpha_{0}+\mathbf{\alpha}_{y}%
^{\ast^{\prime}}S_{y}\left(  y\right) +\mathbf{\alpha}_{w} S_{w}\left(  W\right);
\end{align*}
 where $S_{w}$ is specified by the user and  $\mathbf{\alpha}_{w}$ is a pre-specified sensitivity parameter encoding a hypothetical violation of a valid NCO assumption in the direction of  $S_{w}$. For instance, one might let  $S_{w}(W)=W$ and vary $\mathbf{\alpha}_{w}$ over a grid of values in the neighborhood of zero (which corresponds to the valid NCO assumption). For each hypothetical value of $\mathbf{\alpha}_{w}$, one would then re-estimate the EPS using methods described in the paper with the sensitivity function specified as an offset. Such a sensitivity analysis is implemented in context of the Zika virus Application for illustration and reported  with corresponding discussion in Appendix \ref{sec:supp:sensitivity}.
 
 Finally, as previously mentioned in the Introduction, COCA offers an alternative approach to proximal causal inference \citep{Miao2016_arxiv,TT2020_Proximal}, for debiasing observational estimates of the effect of treatment on the treated by leveraging negative control or valid confounding proxies. A key difference highlighted earlier between these two frameworks is that COCA relies on a single valid NCO which directly proxies the treatment-free potential outcome, while proximal causal inference requires both valid negative control outcome and negative control treatment variables that proxy an underlying unmeasured confounder. Importantly, COCA takes advantage of the fact that the treatment-free potential outcome is observed in the untreated, while in proximal causal inference, the unmeasured confounders for which proxies are available, are themselves never observed, arguably a more challenging identification task. Despite the practical advantage of needing one rather than two proxies, it is important to note that though COCA identifies the effect of treatment on the treated, it fails to nonparametrically identify the population average treatment effect (ATE), identified by proximal causal inference, without an additional assumption. In contrast, proximal causal inference provides nonparametric identification of both causal parameters and thus can be interpreted as providing richer identification opportunities. A key reason for this difference in the scope of identification is the fact that in the current paper, we have emphasized an interpretation of the NCO as a proxy for the treatment-free potential outcome but not for the potential outcome under treatment, in the sense that under our conditions, it must be that the treatment-free potential outcome does not only shield the treatment from the NCO, but also shields the potential outcome under treatment from the latter. Two potential strategies to recover COCA identification of the population ATE, might be either (i) evoke a rank-preservation assumption which if appropriate would imply that the ETT and the ATE are equal (this is the assumption made in \cite{TT2013_COCA}); or (ii) identify a second proxy negative control outcome which is a valid proxy for the counterfactual outcome under treatment. The second condition would be needed if $Y^{a=1}$ can also be viewed as a hidden confounder.  By a symmetry argument, one can show that (ii) in fact would provide identification of the average counterfactual outcome under treatment for the untreated. A weighted average of both counterfactual means would then provide identification of the ATE. Details are not provided, but can easily be deduced from the presentation. 






\section{Acknowledgment}

The authors would like to thank Jamie Robins, Thomas Richardson and Ilya Shpitser for helpful discussions,



\bibliographystyle{apa}

