\documentclass[notitlepage]{article} \usepackage{Math} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{bm} \usepackage{bbold} \usepackage{verbatim} \usepackage{amsfonts,amsmath,amssymb} \newtheorem{alg}{Algorithm} \newtheorem{ex}{Example} \newtheorem{note}{Note} \newenvironment{proof}{\par\noindent{\normalfont{\bf{Proof }}}} {\hfill \small$\blacksquare$\\[2mm]} \setlength{\parindent}{0in} \def\E{\mathbb{E}} \def\Var{\text{Var}} \def\Cov{\text{Cov}} \def\Corr{\text{Corr}} \def\bW{\mathbf{W}} \newcommand{\bs}[1]{\bm{#1}} \newcommand{\Rho}{\mathrm{P}} \newcommand{\cif}{\text{ if }} \newcommand{\Ra}{\Longrightarrow} \newcommand{\ra}{\longrightarrow} \def\bphi{\boldsymbol \phi} \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}} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) @ \begin{center} {\textsclarge{AIPW Estimation via Regression}} \end{center} For binary treatment $Z$, the doubly robust AIPW estimator of the average treatment effect (ATE) $\delta$ \[ \delta = \mu(1) - \mu(0) = \E[Y(1)-Y(0)] \] is \[ \widehat \delta = \frac{1}{n} \sum_{i=1}^n \left\{ \mu(X_i,1) - \mu(X_i,0) \right\} + \frac{1}{n} \sum_{i=1}^n \frac{Z_i}{e(X_i)} (Y_i - \mu(X_i,Z_i)) - \frac{1}{n} \sum_{i=1}^n \frac{1-Z_i}{1-e(X_i)} (Y_i - \mu(X_i,Z_i)) \] where $\mu(x,z)$ is the modelled conditional expectation of $Y$ given $Z=z$ and $X=x$. This is the difference of the two estimators of the average potential outcomes (APOs) for treatments 1 and 0 ie for $z=0,1$ \begin{equation} \widetilde \mu (z) = \frac{1}{n} \sum_{i=1}^n \mu(X_i,z) + \frac{1}{n} \sum_{i=1}^n \left( z \frac{Z_i}{e(X_i)} + (1-z) \frac{1-Z_i}{1-e(X_i)} \right) (Y_i - \mu(X_i,Z_i)) \label{eq:DR1}. \end{equation} For example, for $z=1$, the corresponding estimating equation is \begin{equation} \sum_{i=1}^n \left[ \mu(X_i,1) + \frac{Z_i}{e(X_i)} (Y_i - \mu(X_i,Z_i)) - \mu (1) \right] = 0\label{eq:EE1}. \end{equation} \medskip If the true conditional expectation model is $\mu(x,z)$ and the true propensity model is $e(x)$, but that models $m(x,z)$ and $g(x)$ are used, in their place then \begin{align*} \E[\widetilde \mu (z)] & = \E[m(X_i,z)] + \E \left[ \left( z \frac{Z_i}{g(X_i)} + (1-z) \frac{1-Z_i}{1-g(X_i)} \right) (Y_i - m(X_i,Z_i)) \right] \\[6pt] & = \E[m(X_i,z)] + \E \left[ z \frac{Z_i}{g(X_i)} (Y_i - m(X_i,Z_i)) \right] + \E \left[ (1-z) \frac{1-Z_i}{1-g(X_i)} (Y_i - m(X_i,Z_i)) \right] \\[6pt] & = \E[m(X_i,z)] + \E \left[ z \frac{e(X_i)}{g(X_i)} (Y_i - m(X_i,1)) \right] + \E \left[ (1-z) \frac{1-e(X_i)}{1-g(X_i)} (Y_i - m(X_i,0)) \right] \\[6pt] & = \E[m(X_i,z)] + \E \left[ z \frac{e(X_i)}{g(X_i)} (\mu(X_i,1) - m(X_i,1)) \right] + \E \left[ (1-z) \frac{1-e(X_i)}{1-g(X_i)} (\mu(X_i,0) - m(X_i,0)) \right]. \end{align*} \begin{itemize} \item If $m(x,z) \equiv \mu(x,z)$, then the latter two terms are zero, and \[ \E[\widetilde \mu (z)] = \E[\mu(X_i,z)]. \] \item If $g(x) \equiv e(x)$, \[ \E[\widetilde \mu (z)] = \E[m(X_i,z)] + \E \left[ z (\mu(X_i,1) - m(X_i,1)) \right] + \E \left[ (1-z) (\mu(X_i,0) - m(X_i,0)) \right] = \E[\widetilde\mu(X_i,z)]. \] \end{itemize} Hence the estimator is doubly robust. \medskip \textbf{Regression Approach:} The \textit{augmented outcome regression} with \[ \E[Y|Z,X] = \mu_{}(Z,X) + \phi_1 \frac{Z}{e(X)} + \phi_0 \frac{1-Z}{1-e(X)} = \mu_A(Z,X) \] say, could be used to estimate the quantity $\delta$ using OLS. This is easily seen as the two estimating equations for the parameters $\phi_1$ and $\phi_0$ are \begin{align*} \sum_{i=1}^n \frac{Z_i}{e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i) - \phi_1 \frac{Z_i}{e(X_i)} - \phi_0 \frac{1-Z_i}{1-e(X_i)} \right) & = 0 \\[6pt] \sum_{i=1}^n \frac{1-Z_i}{1-e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i) - \phi_1 \frac{Z_i}{e(X_i)} - \phi_0 \frac{1-Z_i}{1-e(X_i)} \right) & = 0. \end{align*} or equivalently \begin{align} \sum_{i=1}^n \frac{Z_i}{e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i) - \phi_1 \frac{Z_i}{e(X_i)} \right) & = 0 \label{eq:DR1-1}\\[6pt] \sum_{i=1}^n \frac{1-Z_i}{1-e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i) - \phi_0 \frac{1-Z_i}{1-e(X_i)} \right) & = 0 \label{eq:DR1-0}. \end{align} Thus \begin{equation} \widehat \phi_1 = \dfrac{\sum\limits_{i=1}^n \dfrac{Z_i}{e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i)\right) }{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} \right)^2 } \qquad \qquad \widehat \phi_0 = \dfrac{\sum\limits_{i=1}^n \dfrac{1-Z_i}{1-e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i)\right) }{\sum\limits_{i=1}^n \left(\dfrac{1-Z_i}{1-e(X_i)} \right)^2 } \label{eq:phihats}. \end{equation} The corresponding estimators are derived from the fitted values from this model. For example \begin{align} \widetilde \mu(1) & = \frac{1}{n} \sum_{i=1}^n \mu_{}(X_i,1) + \widehat \phi_1 \frac{1}{n} \sum\limits_{i=1}^n \dfrac{1}{e(X_i)} \nonumber\\[6pt] & = \frac{1}{n} \sum_{i=1}^n \mu_{}(X_i,1) + \dfrac{\sum\limits_{i=1}^n \dfrac{Z_i}{e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i)\right) }{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} \right)^2 } \frac{1}{n} \sum\limits_{i=1}^n \dfrac{1}{e(X_i)} \nonumber\\[6pt] & = \frac{1}{n} \sum_{i=1}^n \mu_{}(X_i,1) + \dfrac{\sum\limits_{i=1}^n \dfrac{1}{e(X_i)}}{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} \right)^2 } \frac{1}{n} \sum\limits_{i=1}^n \dfrac{Z_i}{e(X_i)} \left( Y_i - \mu_{}(X_i,Z_i)\right) \label{eq:DR2} \end{align} Now, \[ \frac{1}{n} \sum\limits_{i=1}^n \dfrac{1}{e(X_i)} \CiP \E \left[ \dfrac{1}{e(X)} \right] \qquad \qquad \frac{1}{n} \sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} \right)^2 \CiP \E \left[ \left( \dfrac{Z}{e(X)} \right)^2 \right] = \E \left[ \dfrac{1}{e(X)} \right] \] so the two estimators in \eqref{eq:DR1} and \eqref{eq:DR2} are asymptotically identical. The doubly robust estimating equation being solved here is the equivalent to \eqref{eq:EE1}, that is \begin{equation} \sum_{i=1}^n \left[ \widetilde \mu_A(X_i,1) + \frac{Z_i}{e(X_i)} (Y_i - \widetilde \mu_A(X_i,Z_i)) - \mu (1) \right] = 0\label{eq:EE2} \end{equation} where \[ \widetilde \mu_A(X_i,1) = \mu(X_i,1) + \widehat \phi_1 \frac{1}{e(X_i)}. \] By \eqref{eq:DR1-1}, we have that \[ \sum_{i=1}^n \frac{Z_i}{e(X_i)} (Y_i - \widetilde \mu_A(X_i,Z_i)) = \sum_{i=1}^n \left( Y_i - \mu_{}(X_i,Z_i) - \widehat \phi_1 \frac{Z_i}{e(X_i)} \right) \frac{Z_i}{e(X_i)} = 0. \] The double robustness property relates to the correct specification of $\mu_A(x,z)$ as the conditional expectation of $Y$ given $Z=z$ and $X=x$ (or of $e(x)$ as the propensity model). However, if in fact $\mu(x,z)$ is the correct conditional expectation, then from \eqref{eq:phihats} it is clear that $\widehat \phi_1 \CiP 0$ as $n \lra \infty$. \pagebreak \textbf{Simulation study:} In the following study we examine the performance of six estimators: \begin{enumerate}[1.] \item $\mu(x,z)$ and $e(x)$ correctly specified, augmented outcome regression estimator; \item $\mu(x,z)$ and $e(x)$ correctly specified, direct AIPW estimator; \item $e(x)$ correctly specified, $\mu(x,z)$ incorrectly specified, augmented outcome regression estimator; \item $e(x)$ correctly specified, $\mu(x,z)$ incorrectly specified, direct AIPW estimator; \item $\mu(x,z)$ correctly specified, $e(x)$ incorrectly specified, augmented outcome regression estimator; \item $\mu(x,z)$ correctly specified, $e(x)$ incorrectly specified, direct AIPW estimator; \end{enumerate} First we compute the true value of the ATE by simulation using a large sample: <>= library(MASS) set.seed(23987) n<-1000000 D<-diag(c(0.25,0.5,0.75)) Mu<-c(1,-2,-1) Sigma<-D %*% (matrix(c(1,0.9,-0.1,0.9,1,-0.2,-0.1,-0.2,1),3,3) %*% D) X<-mvrnorm(n,mu=Mu,Sigma) theta<-2.0 al<-c(6.0,-0.2,0.7,2) Xal<-cbind(1,X[,1],X[,2],X[,1]*X[,2]) expit<-function(x){return(1/(1+exp(-x)))} ps.true<-expit(Xal %*% al) Z<-rbinom(n,1,ps.true) be<-c(0,theta,-2.0,1.2,0.6) Xbe<-cbind(1,Z,X[,1],X[,2],X[,1]*Z) mu.true<- Xbe %*% be Y<-rnorm(n)*5+mu.true X0<-cbind(1,0,X[,1],X[,2],X[,1]*0) X1<-cbind(1,1,X[,1],X[,2],X[,1]*1) mu0<-X0 %*% be mu1<-X1 %*% be ATE.true<-mean(mu1)-mean(mu0) ATE.true X1<-X[,1]; X2<-X[,2]; X3<-X[,3] @ In this simulation, the true ATE value is \Sexpr{round(ATE.true,3)}. <>= Yoff<-Y-mu.true C1<-Z/ps.true phi1.hat<-sum(C1*Yoff)/sum(C1^2) mu1.hat<-mean(mu1)+phi1.hat*mean(1/ps.true) C0<-(1-Z)/(1-ps.true) phi0.hat<-sum(C0*Yoff)/sum(C0^2) mu0.hat<-mean(mu0)+phi0.hat*mean(1/(1-ps.true)) ATE.hat<-mu1.hat-mu0.hat Yoff<-Y C1<-Z/ps.true phi1.hat<-sum(C1*Yoff)/sum(C1^2) mu1.hat<-phi1.hat*mean(1/ps.true) C0<-(1-Z)/(1-ps.true) phi0.hat<-sum(C0*Yoff)/sum(C0^2) mu0.hat<-phi0.hat*mean(1/(1-ps.true)) ATE.hat.mis<-mu1.hat-mu0.hat c(ATE.hat,ATE.hat.mis) @ The AOR and the AIPW estimators exhibit apparent consistent estimation; however, they are not numerically identical. <>= nreps<-10000 res.mat<-matrix(0,nrow=nreps,ncol=6) nsamp<-1000 for(irep in 1:nreps){ isub<-sample(1:n,size=nsamp,rep=T) Yoff<-Y[isub]-mu.true[isub] W1.sub<-Z[isub]/ps.true[isub] W0.sub<-(1-Z[isub])/(1-ps.true[isub]) mu.sub<-mu.true[isub] Z.sub<-Z[isub] mu1.sub<-mu1[isub] mu0.sub<-mu0[isub] #AIPW estimator, correct specification ATE.direct<-mean(mu1.sub)-mean(mu0.sub)+mean(W1.sub*Yoff)-mean(W0.sub*Yoff) #AOR estimator, correct specification phi1.hat<-sum(W1.sub*Yoff)/sum(W1.sub^2) mu1.hat<-mean(mu1.sub)+phi1.hat*mean(1/ps.true[isub]) phi0.hat<-sum(W0.sub*Yoff)/sum(W0.sub^2) mu0.hat<-mean(mu0.sub)+phi0.hat*mean(1/(1-ps.true[isub])) ATE.hat<-mu1.hat-mu0.hat #Incorrect mean model Yoff<-Y[isub] #AIPW estimator, incorrect mean specification ATE.direct.mis<-mean(W1.sub*Yoff)-mean(W0.sub*Yoff) #AOR estimator, incorrect mean specification phi1.hat<-sum(W1.sub*Yoff)/sum(W1.sub^2) mu1.hat<-phi1.hat*mean(1/ps.true[isub]) phi0.hat<-sum(W0.sub*Yoff)/sum(W0.sub^2) mu0.hat<-phi0.hat*mean(1/(1-ps.true[isub])) ATE.hat.mis<-mu1.hat-mu0.hat #Incorrect PS model Yoff<-Y[isub]-mu.true[isub] ps.mis<-fitted(glm(Z[isub]~X1[isub],family=binomial)) W1.sub<-Z[isub]/ps.mis W0.sub<-(1-Z[isub])/(1-ps.mis) #AIPW estimator incorrect ps specification ATE.direct.ps.mis<-mean(mu1.sub)-mean(mu0.sub)+mean(W1.sub*Yoff)-mean(W0.sub*Yoff) #AOR estimator, incorrect ps specification phi1.hat<-sum(W1.sub*Yoff)/sum(W1.sub^2) mu1.hat<-mean(mu1.sub)+phi1.hat*mean(1/ps.mis) phi0.hat<-sum(W0.sub*Yoff)/sum(W0.sub^2) mu0.hat<-mean(mu0.sub)+phi0.hat*mean(1/(1-ps.mis)) ATE.ps.mis<-mu1.hat-mu0.hat res.mat[irep,]<-c(ATE.direct,ATE.hat,ATE.direct.mis,ATE.hat.mis,ATE.direct.ps.mis,ATE.ps.mis) } @ <>= br<-seq(0.0,5.0,by=0.1);par(mar=c(2,2,2,0));res<-res.mat[,1];res<-res[res>0 & res < 5] hist(res,ylim=range(0,2000),breaks=br,main='AIPW correct specification',col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) br<-seq(0.0,5.0,by=0.1);par(mar=c(2,2,2,0));res<-res.mat[,2];res<-res[res>0 & res < 5] hist(res,ylim=range(0,2000),breaks=br,main='AOR correct specification',col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) br<-seq(0.0,5.0,by=0.1);par(mar=c(2,2,2,0));res<-res.mat[,3];res<-res[res>0 & res < 5] hist(res,ylim=range(0,2000),breaks=br, main=expression(paste('AIPW ',mu(x,z),' incorrect')),col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) br<-seq(0.0,5.0,by=0.1);par(mar=c(2,2,2,0));res<-res.mat[,4];res<-res[res>0 & res < 5] hist(res,ylim=range(0,2000),breaks=br, main=expression(paste('AOR ',mu(x,z),' incorrect')),col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) br<-seq(0.0,5.0,by=0.1);par(mar=c(2,2,2,0));res<-res.mat[,5];res<-res[res>0 & res < 5] hist(res,ylim=range(0,2000),breaks=br, main=expression(paste('AIPW ',pi(x),' incorrect')),col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) br<-seq(0.0,5.0,by=0.1);par(mar=c(2,2,2,0));res<-res.mat[,6];res<-res[res>0 & res < 5] hist(res,ylim=range(0,2000),breaks=br, main=expression(paste('AOR ',pi(x),' incorrect')),col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) <>= par(mar=c(2,2,0,0)) nvec<-c("AIPW","AOR","AIPW-mu","AOR-mu","AIPW-ps","AOR-ps") boxplot(res.mat,names=nvec,pch=19,cex=0.5);abline(h=ATE.true,col='red',lwd=1) (apply(res.mat,2,mean)-ATE.true)*sqrt(nsamp) #Bias apply(res.mat,2,var)*nsamp #Variance apply((res.mat-ATE.true)^2,2,mean)*nsamp #MSE @ The second Bang and Robins approach proposes taking $\phi_0 = -\phi_1$, and using the AOR model \[ \E[Y|Z,X] = \mu_{}(Z,X) + \phi \left( \frac{Z}{e(X)} - \frac{1-Z}{1-e(X)} \right) \] to produce the ATE estimator \[ \frac{1}{n} \sum_{i=1}^n \left\{ \mu(X_i,1) - \mu(X_i,0) \right\} + \widehat \phi \sum_{i=1}^n \left( \frac{1}{e(X_i)} + \frac{1}{1-e(X_i)} \right). \] By elementary calculations we have that \[ \widehat \phi = \dfrac{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right) \left( Y_i - \mu(X_i,Z_i)\right) }{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right)^2 } \] so that \begin{align}\label{eq:AOR2} \widehat \delta & = \frac{1}{n} \sum_{i=1}^n \left\{ \mu(X_i,1) - \mu(X_i,0) \right\} + \frac{1}{n} \sum_{i=1}^n \left( \frac{1}{e(X_i)} + \frac{1}{1-e(X_i)} \right) \dfrac{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right) \left( Y_i - \mu(X_i,Z_i)\right) }{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right)^2 } \nonumber\\[6pt] & = \frac{1}{n} \sum_{i=1}^n \left\{ \mu(X_i,1) - \mu(X_i,0) \right\} + \dfrac{\sum\limits_{i=1}^n \left( \dfrac{1}{e(X_i)} + \dfrac{1}{1-e(X_i)} \right) }{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right)^2 }\frac{1}{n} \sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right) \left( Y_i - \mu(X_i,Z_i)\right) . \end{align} As $n \lra \infty$ we have that \[ \dfrac{\sum\limits_{i=1}^n \left( \dfrac{1}{e(X_i)} + \dfrac{1}{1-e(X_i)} \right) }{\sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right)^2}\CiP \dfrac{\E\left[\dfrac{1}{e(X)} + \dfrac{1}{1-e(X)} \right]}{\E\left[\left(\dfrac{Z}{e(X)}\right)^2 + \left(\dfrac{1-Z}{1-e(X)}\right)^2 \right]} = 1 \] so that $\widehat \delta$ is asymptotically identical to \[ \frac{1}{n} \sum_{i=1}^n \left\{ \mu(X_i,1) - \mu(X_i,0) \right\} + \frac{1}{n} \sum\limits_{i=1}^n \left(\dfrac{Z_i}{e(X_i)} - \dfrac{1-Z_i}{1-e(X_i)} \right) \left( Y_i - \mu(X_i,Z_i)\right) \] that is, the AIPW estimator. <>= nreps<-10000 res.mat2<-matrix(0,nrow=nreps,ncol=3) phi.ests<-matrix(0,nrow=nreps,ncol=3) nsamp<-1000 for(irep in 1:nreps){ isub<-sample(1:n,size=nsamp,rep=T) Yoff<-Y[isub]-mu.true[isub] C.sub<-(Z[isub]/ps.true[isub]-(1-Z[isub])/(1-ps.true[isub])) mu.sub<-mu.true[isub] Z.sub<-Z[isub] mu1.sub<-mu1[isub] mu0.sub<-mu0[isub] #AOR estimator, correct specification phi.hat<-sum(C.sub*Yoff)/sum(C.sub^2) mu1.hat<-mean(mu1.sub)+phi.hat*mean(1/ps.true[isub]) mu0.hat<-mean(mu0.sub)+phi.hat*mean(1/(1-ps.true[isub])) ATE.hat<-mu1.hat-mu0.hat phi.ests[irep,1]<-phi.hat #Incorrect mean model Yoff<-Y[isub] #AOR estimator, incorrect mean specification phi.hat<-sum(C.sub*Yoff)/sum(C.sub^2) mu1.hat<-phi.hat*mean(1/ps.true[isub]) mu0.hat<--phi.hat*mean(1/(1-ps.true[isub])) ATE.hat.mis<-mu1.hat-mu0.hat phi.ests[irep,2]<-phi.hat #Incorrect PS model Yoff<-Y[isub]-mu.true[isub] ps.mis<-fitted(glm(Z[isub]~X1[isub],family=binomial)) C.sub<-Z[isub]/ps.mis-(1-Z[isub])/(1-ps.mis) #AOR estimator, incorrect ps specification phi.hat<-sum(C.sub*Yoff)/sum(C.sub^2) mu1.hat<-mean(mu1.sub)+phi.hat*mean(1/ps.mis) mu0.hat<-mean(mu0.sub)-phi.hat*mean(1/(1-ps.mis)) ATE.ps.mis<-mu1.hat-mu0.hat phi.ests[irep,3]<-phi.hat res.mat2[irep,]<-c(ATE.hat,ATE.hat.mis,ATE.ps.mis) } @ <>= br<-seq(1.5,3.5,by=0.05);par(mar=c(3,2,2,0)) res<-res.mat2[,1];res<-res[res>1.5 & res < 3.5] hist(res,ylim=range(0,2000),breaks=br, main='AOR correct specification',col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) res<-res.mat2[,2];res<-res[res>1.5 & res < 3.5] hist(res,ylim=range(0,2000),breaks=br, main=expression(paste('AOR ',mu(x,z),' incorrect')),col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) res<-res.mat2[,3];res<-res[res>1.5 & res < 3.5] hist(res,ylim=range(0,2000),breaks=br, main=expression(paste('AOR ',pi(x),' incorrect')),col='gray',xlab='ATE') box();abline(v=ATE.true,col='red',lwd=1) @ <>= par(mar=c(2,2,0,0)) nvec<-c("Correct",expression(paste(mu(x,z),' incorrect')),expression(paste(e(x),' incorrect'))) boxplot(res.mat2,names=nvec,pch=19,cex=0.5);abline(h=ATE.true,col='red',lwd=1) apply(res.mat2,2,var)*nsamp apply((res.mat2-ATE.true)^2,2,mean)*nsamp @ <>= par(mar=c(2,2,2,0)) nvec<-c("Correct",expression(paste(mu(x,z),' incorrect')),expression(paste(e(x),' incorrect'))) boxplot(phi.ests,names=nvec,pch=19,cex=0.5);abline(h=ATE.true,col='red',lwd=1) abline(h=0,col='red',lty=2) title(expression(paste('Distribution of ',hat(phi),' under the correct and incorrect specifications'))) @ \end{document}