\documentclass[notitlepage]{article}
\usepackage{Math}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumerate}
\usepackage{bbm}
\usepackage{bm}
\usepackage[all]{xy}
\usepackage{verbatim}
\usepackage{mathtools}
\usepackage{dsfont}

\DeclareMathAlphabet{\mathpzc}{OT1}{pzc}{m}{it}

\makeatletter
\newcommand{\raisemath}[1]{\mathpalette{\raisem@th{#1}}}
\newcommand{\raisem@th}[3]{\raisebox{#1}{$#2#3$}}
\makeatother

\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}}
\def\bH{\mathbf{H}}

\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}}

\def\Xb{\mathbf{X}_{\beta}}
\def\Xp{\mathbf{X}_{\psi}}
\def\beps{\boldsymbol{\epsilon}}
\def\betahat{\widehat \beta}
\def\psihat{\widehat \psi}

\begin{document}

<<setup, cache=FALSE, echo=FALSE>>=
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{Example: Analysis of the Nhanes data}}
\end{center}

In this example, we will explore propensity score based analyses using the publicly available (U.S.) National Health and Nutrition Examination Survey (NHANES). For this, may use the libraries \texttt{NHANES}, \texttt{tableone}, and \texttt{Matching} in \texttt{R}. We will focus our analysis on the question of whether current smoking affects average systolic blood pressure. The variables we will need are:

\begin{itemize}

\item \texttt{BPSysAve} (systolic blood pressure, response),
\item \texttt{SmokeNow} (smoking status: 0--No, 1--Yes),

\item \texttt{Gender},

\item \texttt{Age},

\item \texttt{Race3},

\item \texttt{Education},

\item \texttt{MaritalStatus}, and

\item \texttt{Poverty}

\end{itemize}

where the first two are the outcome and exposure of interest and the remaining are potential confounders. Additionally, we will restrict our attention to adults ($>$ 17 years old) in the second wave of the survey.

<<NH0,comment='+', fig.width=7.25, fig.height=5.5, fig.align='center'>>=
library(NHANES);library(tableone);library(Matching)
NHANES$SmokeNow <- as.numeric(NHANES$SmokeNow)-1
small.nhanes <- na.omit(NHANES[NHANES$SurveyYr=="2011_12" & NHANES$Age > 17,c(3,4,8:11,13,25,61)])
dim(small.nhanes) ## 1377
vars <- c("Gender", "Age", "Race3", "Education", "MaritalStatus", "Poverty")

fit0<-lm(BPSysAve~SmokeNow,data=small.nhanes)
round(coef(summary(fit0)),5)
fit1<-lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=small.nhanes)
round(drop1(fit1,test='F'),5)
@

{\textbf{Checking Confounder balance}
\begin{itemize}
\item In PS--based methods, the goal of the treatment model is to eliminate imbalance in the distribution of covariates between treated and untreated subjects.

\begin{itemize}
    \item Achieving balance on other covariates (particularly strong predictors of treatment) is unhelpful.
\end{itemize}
\item The goal is \textit{not} to build an excellent predictive model for the treatment.

\end{itemize}
}

Common measures of balance include:
\begin{itemize}
    \item Standardized mean difference (SMD) or proportion:
    \[
    \frac{\bar{x}^{1,w}-\bar{x}^{0,w}}{\sqrt{0.5(v^{1,w}+v^{0,w})}}
    \qquad \text{where} \qquad \bar{x}^{z,w} = \frac{1}{n} \sum_{i=1}^n \frac{\Ind_{z}(z_i) x_i}{f_{Z|X}^\Obs(z_i|x_i)},
    \]
        i.e.~the weighted sample mean of variable $X$ among those with treatment value $z$, and similarly $v^{z,w}$ is the weighted variance estimate.
        \begin{itemize}
           \item For all methods of analysis other than IPW, the weights are taken to be 1 for all subjects.
           \item SMD of 0.1 or less typically considered reasonable.
        \end{itemize}
    \item Visual examination of weighted empirical CDFs among the treated and untreated (for binary or categorical treatment).
\end{itemize}

<<NH1,comment='+', fig.width=7.25, fig.height=5.5, fig.align='center'>>=
@
<<NH3,comment='+', fig.width=7.25, fig.height=5.5, fig.align='center'>>=
tabUnmatched <- CreateTableOne(vars = vars, strata = "SmokeNow", data = small.nhanes, test = FALSE)
print(tabUnmatched, smd = TRUE)
Smoke<-small.nhanes$SmokeNow
age0 <- small.nhanes$Age[Smoke==0]
age1 <- small.nhanes$Age[Smoke==1]
ecdf0 <- ecdf(age0)
ecdf1 <- ecdf(age1)
par(mar=c(2,2,2,0))
plot(ecdf0, verticals=TRUE, do.points=FALSE,main='Empirical cdfs for Age',col='blue')
plot(ecdf1, verticals=TRUE, do.points=FALSE, add=TRUE, col='red')
legend(20,0.9,c('Smoke = 0', 'Smoke = 1'),col=c('blue','red'),lty=1)
@
There is clear age imbalance between the two smoking groups.  This implies that if \texttt{Age} is also a possible predictor of the outcome, then it is a possible confounder.

\medskip

\textbf{Building the propensity score:}  We attempt to build and assess a propensity score model using the available covariates:

<<NH4,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
ps.mod <- glm(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,
                data=small.nhanes,family="binomial")
ps.lr <- predict(ps.mod,type="response")
summary(ps.lr)
ps0<-ps.lr[Smoke==0]
ps1<-ps.lr[Smoke==1]
quints <- c(0,quantile(ps.lr,seq(.2,1,.2)))
rbind(table(cut(ps.lr[Smoke==0],quints)),
      table(cut(ps.lr[Smoke==1],quints)))
par(mar=c(2,2,2,0))
hist(ps0, col=rgb(0,0,1,0.5), breaks=seq(0,1,by=0.05), ylim=c(0,150),
        main="Propensity Score overlap", xlab="PS")
hist(ps1, col=rgb(1,0,0,0.5), breaks=seq(0,1,by=0.05), add=T);box()
legend(0,0.9,c('Smoke = 0', 'Smoke = 1'),col=c(rgb(0,0,1,0.5),rgb(1,0,0,0.5)),lty=1)
boxplot(ps0,ps1,ylab="PS",xlab="Treatment Group",names=c('SmokeNow==0','SmokeNow==1'),
            col=c(rgb(0,0,1,0.5),rgb(1,0,0,0.5)));abline(h=quints[2:5],col="red")
@
We can therefore proceed to check for balance knowing we have sufficient numbers of smokers and non-smokers in each quintile
to ensure the stratum-specific estimates are not too unstable.

<<NH5,comment='+', fig.width=7.00, fig.height=7.0, fig.align='center'>>=
Pcat<-as.numeric(cut(ps.lr,quints,include.lowest=T))
par(mar=c(2,3,4,0),mfrow=c(3,2),oma=c(0,0,2,0))
for(k in 1:5){
    age0 <- small.nhanes$Age[Smoke==0 & Pcat==k]
    age1 <- small.nhanes$Age[Smoke==1 & Pcat==k]
    ecdf0 <- ecdf(age0)
    ecdf1 <- ecdf(age1)
    plot(ecdf0, verticals=TRUE, do.points=FALSE,
        main=substitute(paste('Quintile ',k),list(k=k)),col='blue')
    plot(ecdf1, verticals=TRUE, do.points=FALSE, add=TRUE, col='red')
}
plot(age0,type='n',ylim=range(0,1),axes=FALSE)
title("ECDFs for Age by PS quintile",outer = TRUE)
legend(30,0.75,c('Smoke = 0', 'Smoke = 1'),col=c('blue','red'),lty=1)
@

Balance does not appear to have been achieved: we have SMDs $>$ 0.1 for at least three quintiles for all variables, and the empirical CDFs of age do not overlap in several quintiles.

<<NH6,comment='+', fig.width=7.25, fig.height=8.5, fig.align='center'>>=
W<-(Smoke==0)/(1-ps.lr)+(Smoke==1)/ps.lr
smd.mat<-ExtractSmd(tabUnmatched)
for(k in 1:5){
    nhanesQ<-small.nhanes[Pcat == k,]
    tabQs <- CreateTableOne(vars = vars, strata = "SmokeNow", data = nhanesQ, test = FALSE)
    smd.mat<-cbind(smd.mat,ExtractSmd(tabQs))
}
colnames(smd.mat)<-c('Original','Q1','Q2','Q3','Q4','Q5')
round(smd.mat,4)
@
We may try to find evidence in smaller PS strata, such as those formed by deciles of the propensity score.

<<NH7,comment='+', fig.width=7.25, fig.height=8.5, fig.align='center'>>=
dec <- c(0,quantile(ps.lr,seq(0.1,1,.1)))
Pcat10 <- cut(ps.lr,dec,labels=1:10)
SMD.10.table <- ExtractSmd(tabUnmatched)
for(k in 1:10) {
   tabPSdec <- CreateTableOne(vars = vars, strata = "SmokeNow",
            data = small.nhanes[Pcat10==k,], test = FALSE)
   SMD.10.table <- cbind(SMD.10.table,ExtractSmd(tabPSdec))
}
colnames(SMD.10.table)<-c('Original','Q1','Q2','Q3','Q4','Q5','Q6','Q7','Q8','Q9','Q10')
round(SMD.10.table,4)
@
This does not assist with providing evidence of balance.

\bigskip

\textbf{Matching:}  We now attempt matching using the propensity score.  The function \texttt{MatchBalance} from the \texttt{Matching} library provides many more details than \texttt{CreateTableOne}, including:
\begin{itemize}
  \item mean, median, and maximum difference in empirical CDF plots,
  \item mean, median, and maximum difference in empirical QQ plots,
  \item Kolmogorov-Smirnov statistics,
  \item ratio of variances,
  \item p-value for t-test.
  \item Note that SMDs reported are multiplied by 100.
\end{itemize}

<<NH8,comment='+', fig.width=7.25, fig.height=8.5, fig.align='center'>>=
small.nhanes$ps.lr<-ps.lr
ps.lr.match <- Match(Tr=small.nhanes$SmokeNow,X=small.nhanes$ps.lr,estimand="ATE",ties=FALSE)
matched.samp <- small.nhanes[c(ps.lr.match$index.control,ps.lr.match$index.treated),]
table(table(c(ps.lr.match$index.control, ps.lr.match$index.treated)))
tabMatched <- CreateTableOne(vars = vars, strata = "SmokeNow",data = matched.samp, test = FALSE)
MatchBalance(SmokeNow~Gender+Age+Education+MaritalStatus+Poverty,data=small.nhanes,match.out=ps.lr.match)
@

<<NH9,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
age0 <- matched.samp$Age[matched.samp$SmokeNow==0]
age1 <- matched.samp$Age[matched.samp$SmokeNow==1]
ecdf0 <- ecdf(age0)
ecdf1 <- ecdf(age1)
par(mar=c(2,2,2,0))
plot(ecdf0, verticals=TRUE, do.points=FALSE,main='Empirical cdfs for Age in matched sample',col='blue')
plot(ecdf1, verticals=TRUE, do.points=FALSE, add=TRUE, col='red')
legend(20,0.9,c('Smoke = 0', 'Smoke = 1'),col=c('blue','red'),lty=1)
@

\textbf{Inverse Weighting:}  We can construct an inverse weighted (IPW) sample and examine it using the \texttt{survey} and \texttt{Hmisc} packages.
<<NH10,comment='+', fig.width=7.25, fig.height=4.25, fig.align='center'>>=
library(survey);library(Hmisc);library(htmlwidgets)
ps.lr.weight <- small.nhanes$SmokeNow/ps.lr + (1-small.nhanes$SmokeNow)/(1-ps.lr)
nhanes.IPW.lr <- svydesign(ids=~0, data=small.nhanes, weights=ps.lr.weight)
tabIPW <- svyCreateTableOne(vars = vars, strata = "SmokeNow",data = nhanes.IPW.lr, test = FALSE)
print(tabIPW, smd = TRUE)
temp0 <- Ecdf(small.nhanes$Age[Smoke==0],weights=ps.lr.weight[Smoke==0],pl=FALSE)
temp1 <- Ecdf(small.nhanes$Age[Smoke==1],weights=ps.lr.weight[Smoke==1],pl=FALSE)
par(mar=c(2,2,2,0))
plot(temp0$x,temp0$y,ylab="ECDF(Age)",xlab="Age",
        main='Empirical cdfs for Age in weighted sample',col='blue',type="s",lwd=1)
lines(temp1$x,temp1$y,col="red",lwd=1,type='s')
smd.mat<-cbind(smd.mat,ExtractSmd(tabMatched))
smd.mat<-cbind(smd.mat,ExtractSmd(tabIPW))
colnames(smd.mat)<-c('Original','Q1','Q2','Q3','Q4','Q5',"Match",'IPW')
round(smd.mat,4)
@
Here the matched and inverse weighted samples exhibit good balance.


\pagebreak

\textbf{Summary: }
\begin{itemize}
\item Creating or restoring confounder balance is essential to estimating a causal effect.
\item It can be hard to assess overlap or achieve balance in high dimensions.
\item The propensity score, a scalar summary of confounding variables, simplifies this task.
\item However:
\begin{itemize}
\item fitting a model for treatment does not guarantee balance,
\item fitting a model that predicts treatment with a high degree of precision can be unhelpful.
\end{itemize}
\end{itemize}

\bigskip

\textbf{Estimating the ATE:} We proceed now to estimating the average treatment effect (ATE), using:
\begin{itemize}
\item outcome regression,
\item PS stratification,
\item PS matching,
\item PS regression,
\item IPW.
\end{itemize}
We will use the PS estimated via logistic regression, as this provided the best balance.


\begin{itemize}

\item \textbf{Linear regression:} 
<<NH11,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
coef(lm(BPSysAve~SmokeNow,data=small.nhanes))[2]
coef(lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,data=small.nhanes))[2]
@
The naive conditional effect estimate is more than 3 times greater than its confounder-adjusted counterpart.

\item \textbf{Outcome regression no interaction:}
<<NH12,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
nhanes.allsmoke <- small.nhanes
nhanes.allsmoke$SmokeNow <- 1
nhanes.nosmoke <- small.nhanes
nhanes.nosmoke$SmokeNow <- 0
mod1.lm <- lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty,
      data=small.nhanes)
APO.lm.1 <- mean(predict(mod1.lm,nhanes.allsmoke))
APO.lm.0 <- mean(predict(mod1.lm,nhanes.nosmoke))
APO.lm.1 - APO.lm.0
@
Conditional and marginal effect are the same in a linear model with no interaction.

\item \textbf{Outcome regression with interactions:}
<<NH13,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
mod1.lmX <- lm(BPSysAve~SmokeNow+Gender+Age+Race3+Education+
     MaritalStatus+HHIncome+Poverty+SmokeNow:HHIncome+SmokeNow:Gender+SmokeNow:Age,data=small.nhanes)
APO.lmX.1 <- mean(predict(mod1.lmX,nhanes.allsmoke))
APO.lmX.0 <- mean(predict(mod1.lmX,nhanes.nosmoke))
APO.lmX.1 - APO.lmX.0
@

\item \textbf{PS stratification}
<<NH14,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
Y<-small.nhanes$BPSysAve
ps.lr.quints <- cut(ps.lr,quints,labels=1:5)
p.strat <- table(ps.lr.quints)/length(ps.lr.quints)
ATE.strat <- rep(NA,5)
for(j in 1:5) {
    ATE.strat[j] <- mean(Y[Smoke == 1 & ps.lr.quints==j])-mean(Y[Smoke == 0 & ps.lr.quints==j])
}
ATE.strat
sum(ATE.strat*p.strat)
@



\item \textbf{PS matching}
<<NH15,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
ps.lr.match <- Match(Tr=small.nhanes$SmokeNow,
      X=small.nhanes$ps.lr,estimand="ATE",ties=FALSE)
matched.samp <- small.nhanes[c(ps.lr.match$index.control,ps.lr.match$index.treated),]

dim(matched.samp)

mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 1]) -
    mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 0])
@

\item \textbf{PS regression}
<<NH16,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
library(splines)
mod1.PSlm1 <- lm(BPSysAve~SmokeNow+ps.lr,data=small.nhanes)
APO.PSlm1.1 <- mean(predict(mod1.PSlm1,nhanes.allsmoke))
APO.PSlm1.0 <- mean(predict(mod1.PSlm1,nhanes.nosmoke))
APO.PSlm1.1 - APO.PSlm1.0

mod1.PSlm2 <- lm(BPSysAve~SmokeNow+ps.lr+I(ps.lr^2),
     data=small.nhanes)
APO.PSlm2.1 <- mean(predict(mod1.PSlm2,nhanes.allsmoke))
APO.PSlm2.0 <- mean(predict(mod1.PSlm2,nhanes.nosmoke))
APO.PSlm2.1 - APO.PSlm2.0

mod1.PSlm3 <- lm(BPSysAve~SmokeNow+bs(ps.lr,df=4),
     data=small.nhanes)
APO.PSlm3.1 <- mean(predict(mod1.PSlm3,nhanes.allsmoke))
APO.PSlm3.0 <- mean(predict(mod1.PSlm3,nhanes.nosmoke))
APO.PSlm3.1 - APO.PSlm3.0
@


\item \textbf{IPW}
<<NH17,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
ps.lr.weight <- Smoke/ps.lr + (1-Smoke)/(1-ps.lr)
mean(Smoke*Y*ps.lr.weight) -mean((1-Smoke)*Y*ps.lr.weight)
coef(lm(Y ~ Smoke, weights = ps.lr.weight))[2]
@

\end{itemize}



{\textbf{Additional considerations:}
\begin{itemize}
\item All of the PS approaches considered rely on substitution estimators.
\begin{itemize}
\item In PS regression, we plug in an estimated PS as a covariate.
\item In IPW, we plug in estimated weights.
\end{itemize}
\item We need to account for this when estimating standard errors and/or confidence intervals.
\item Analytically derived asymptotic variances can be used, but are not provided in many standard software packages.
\item The easiest approach is to bootstrap.
\item Note, however, that the bootstrap is \textit{not} valid for matching.
\end{itemize}
}

\end{document}