\documentclass[notitlepage]{article}
\usepackage{Math598}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumerate}
\usepackage{bbm}
\usepackage{xypic}
\usepackage{amsfonts,amsmath}

\def\E{\Expect}


\lstloadlanguages{R}

\definecolor{keywordcolor}{rgb}{0,0.6,0.6}
\definecolor{delimcolor}{rgb}{0.461,0.039,0.102}
\definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432}

\lstdefinestyle{Rsettings}{
  basicstyle=\ttfamily,
  breaklines=true,
  showstringspaces=false,
  keywords={if, else, function, theFunction, tmp}, % Write as many keywords
  otherkeywords={},
  commentstyle=\itshape\color{Rcommentcolor},
  keywordstyle=\color{keywordcolor},
  moredelim=[s][\color{delimcolor}]{"}{"},
}

\lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}}

\parindent0in

\begin{document}
\begin{center}
{\textsclarge{Causal Contrasts from Outcome Regression}}
\end{center}

In a randomized experimental study, inference concerning the causal effect of binary variable $Z$ on $Y$ can be made by direct comparison of sample averages.  Suppose that $Z \sim Bernoulli(p)$ for $0<p<1$.  Two estimators of the causal contrast $\delta = \E^\Exp[Y(1) - Y(0)]$ can be proposed if the conditional mean model
\[
\E_{Y|X,Z}^\Exp[Y|X=x,Z=z] = \E_{Y|X,Z}^\Obs[Y|X=x,Z=z] = \mu(x,z)
\]
is presumed known.  We have
\begin{equation}\label{eq:tilde}
\widehat{\delta}_{OR} = \frac{1}{n} \sum_{i=1}^n \mu(X_i,1) - \frac{1}{n} \sum_{i=1}^n \mu(X_i,0)
\end{equation}
and
\begin{equation}\label{eq:hat}
\widehat{\delta}_{OR2} = \frac{1}{n \widehat p} \sum_{i=1}^n Z_i \mu(X_i,Z_i) - \frac{1}{n(1-\widehat p)} \sum_{i=1}^n (1-Z_i) \mu(X_i,Z_i)
\end{equation}
where $\widehat p = n^{-1} \sum_{i=1}^n Z_i$.  In the simulation below, we assume that $X \sim Normal(1,1)$ and
\[
Y | X = x, Z = z \sim Normal( 0.5 x + \delta z, 1)
\]
so that $\E^\Exp_{Y|X,Z}[Y|X=x,Z=z] = 0.5 x + \delta z$, and
\[
\E[Y(1) - Y(0)] = \E_X^\Exp \left[ \E^\Exp_{Y|X,Z}[Y|X,Z=1] -  \E_{Y|X,Z}^\Exp[Y|X,Z=0] \right] = \delta .
\]
In this case, $\mu(X,1) - \mu(X,0) = \delta$ which does not depend on $X$, so we can estimate $\delta$ without error using $\delta_{OR}$
<<exp1,comment='+', fig.width=7.5, fig.height=4.5, fig.align='center'>>=
set.seed(23)
nreps<-1000;n<-100;p<-0.5;delta<-2
ests.mat<-matrix(0,nrow=nreps,ncol=2)
for(irep in 1:nreps){
    X<-rnorm(n,1,1)
    Z<-rbinom(n,1,p)
    Y<-rnorm(n,delta*Z+0.5*X,1)
    p.hat<-mean(Z)
    ests.mat[irep,1]<-mean(delta+0.5*X)-mean(0.5*X)
    ests.mat[irep,2]<-sum(Z*(delta*Z+0.5*X))/(n*p.hat)-sum((1-Z)*(delta*Z+0.5*X))/(n*(1-p.hat))
}
apply(ests.mat,2,var)
@


\pagebreak
Suppose now that
\[
\E_{Y|X,Z}^\Exp[Y|X=x,Z=z] = 0.5 x + \delta z + 2 x z
\]
so that
\[
\E_{Y|X,Z}^\Exp[Y|X=x,Z=1] -  \E_{Y|X,Z}^\Exp[Y|X=x,Z=0] = \delta + 2 x  .
\]
and thus
\[
\E^\Exp[Y(1) - Y(0)] = \E^\Exp[\delta + 2 X] = \delta + 2 \E^\Exp[X] = \delta + 2
\]
<<exp2,comment='+', fig.width=7.5, fig.height=4.5, fig.align='center'>>=
set.seed(23)
nreps<-1000;n<-100;p<-0.5;delta<-2
ests.mat<-matrix(0,nrow=nreps,ncol=2)
for(irep in 1:nreps){
    X<-rnorm(n,1,1)
    Z<-rbinom(n,1,p)
    Y<-rnorm(n,delta*Z+0.5*X + 2*Z*X,1)
    p.hat<-mean(Z)
    ests.mat[irep,1]<-mean(delta+0.5*X+2*X)-mean(0.5*X)
    ests.mat[irep,2]<-sum(Z*(delta*Z+0.5*X+ 2*X*Z))/(n*p.hat)-
                      sum((1-Z)*(delta*Z+0.5*X+2*X*Z))/(n*(1-p.hat))
}
apply(ests.mat,2,var)
par(mar=c(4,2,1,1))
boxplot(ests.mat,ylim=range(2.5,5));abline(h=delta+2,lty=2,col='red')
@

\pagebreak
Suppose we \textbf{mis-specify} the modelled mean using
\[
m(x,z) = x - 2 x^2 + \delta z + 2 x z
\]
that is, where the dependence on $z$ is correct, but it is not the true conditional mean.  We can still unbiasedly estimate $\delta$ using the same estimator:
<<exp3,comment='+', fig.width=7.5, fig.height=4.5, fig.align='center'>>=
set.seed(23)
nreps<-1000;n<-100;p<-0.5;delta<-2
ests.mat<-matrix(0,nrow=nreps,ncol=2)
for(irep in 1:nreps){
    X<-rnorm(n,1,1)
    Z<-rbinom(n,1,p)
    Y<-rnorm(n,delta*Z+0.5*X + 2*Z*X,1)
    p.hat<-mean(Z)
    ests.mat[irep,1]<-mean(delta+X-X^2+2*X)-mean(X-X^2)
    ests.mat[irep,2]<-sum(Z*(delta*Z+X-X^2+2*Z*X))/(n*p.hat)-
                      sum((1-Z)*(delta*Z+X-X^2+2*Z*X))/(n*(1-p.hat))
}
apply(ests.mat,2,var)
par(mar=c(4,2,1,1))
boxplot(ests.mat,ylim=range(2.5,5));abline(h=delta+2,lty=2,col='red')
@
\pagebreak
However, if we \textbf{mis-specify} the modelled mean using
\[
m(x,z) = x - 2 x^2 + \delta z + x z
\]
that is, where the dependence on $z$ is incorrect, and it is not the true conditional mean, we can no longer unbiasedly estimate $\delta$ using the same estimator:
<<exp4,comment='+', fig.width=7.5, fig.height=4.5, fig.align='center'>>=
set.seed(23)
nreps<-1000;n<-100;p<-0.5;delta<-2
ests.mat<-matrix(0,nrow=nreps,ncol=2)
for(irep in 1:nreps){
    X<-rnorm(n,1,1)
    Z<-rbinom(n,1,p)
    Y<-rnorm(n,delta*Z+0.5*X + 2*Z*X,1)
    p.hat<-mean(Z)
    ests.mat[irep,1]<-mean(delta+X-X^2+X)-mean(X-X^2)
    ests.mat[irep,2]<-sum(Z*(delta*Z+X-X^2+Z*X))/(n*p.hat)-
                      sum((1-Z)*(delta*Z+X-X^2+Z*X))/(n*(1-p.hat))
}
apply(ests.mat,2,var)
par(mar=c(4,2,1,1))
boxplot(ests.mat,ylim=range(2.5,5));abline(h=delta+2,lty=2,col='red')
@

\end{document}