\documentclass[notitlepage,11pt]{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=FALSE, autodep=TRUE)
options(scipen=6)
options(repos=c(CRAN="https://cloud.r-project.org/"))
@
\begin{center}
{\textsclarge{Matching analysis of the Nhanes data}}
\end{center}

The variables are:

\begin{itemize}

\item \texttt{BPSysAve} (systolic blood pressure, response),
\item \texttt{SmokeNow} (smoking status: 0--No, 1--Yes),

\item \texttt{Gender}, \texttt{Age}, \texttt{Race3}, \texttt{Education}, \texttt{MaritalStatus}, and \texttt{Poverty}

\end{itemize}

where the first two are the outcome and exposure of interest and the remaining are potential confounders. We restrict our attention to 1377 adults ($>$ 17 years old) in the second wave of the survey.

<<NH0,comment='+', fig.width=7.25, fig.height=5.5, fig.align='center'>>=
suppressPackageStartupMessages({
library(NHANES); library(tableone); library(Matching)   #Must be pre-loaded in R
library(survey); library(Hmisc); library(htmlwidgets)   #Must be pre-loaded in R
})
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)])
vars <- c("Gender", "Age", "Race3", "Education", "MaritalStatus", "Poverty")
@


<<NH1,comment='+', fig.width=6.25, fig.height=3.9, fig.align='center'>>=
tabUnmatched<-CreateTableOne(vars=vars,strata="SmokeNow",data=small.nhanes,test=FALSE)
print(tabUnmatched, smd = TRUE)
Smoke<-small.nhanes$SmokeNow
@

\bigskip

\textbf{Propensity score matching:}  
<<NH4,comment='+', fig.width=6.25, fig.height=4.25, fig.align='center'>>=
#Build the PS model
fmods<-formula(SmokeNow~Gender+Age+Race3+Education+MaritalStatus+HHIncome+Poverty)
ps.mod <- glm(fmods,data=small.nhanes,family="binomial")
ps.lr <- predict(ps.mod,type="response")
small.nhanes$ps.lr<-ps.lr

#Matching for the ATE
set.seed(1234)
ps.lr.match <- Match(Tr=small.nhanes$SmokeNow,X=small.nhanes$ps.lr,
                        estimand="ATE",ties=FALSE)
summary(ps.lr.match)

match.details<-cbind(ps.lr.match$index.control,
      small.nhanes$SmokeNow[ps.lr.match$index.control],
      ps.lr.match$index.treated,small.nhanes$SmokeNow[ps.lr.match$index.treated],
      small.nhanes$ps.lr[ps.lr.match$index.control],
      small.nhanes$ps.lr[ps.lr.match$index.treated])
cols<-c('Control ID','Control Z','Treated ID','Treated Z','Control PS','Treated PS')
colnames(match.details)<-cols

match.details[1:18,]

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)
print(tabMatched)

###################################################################################
#Matching for the ATT
set.seed(1234)
ps.lr.match <- Match(Tr=small.nhanes$SmokeNow,X=small.nhanes$ps.lr,
                        estimand="ATT",ties=FALSE)
summary(ps.lr.match)

match.details<-cbind(ps.lr.match$index.control,
      small.nhanes$SmokeNow[ps.lr.match$index.control],
      ps.lr.match$index.treated,small.nhanes$SmokeNow[ps.lr.match$index.treated],
      small.nhanes$ps.lr[ps.lr.match$index.control],
      small.nhanes$ps.lr[ps.lr.match$index.treated])
colnames(match.details)<-cols

match.details[1:18,]

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)
print(tabMatched)
@

<<NH9a,comment='+', fig.width=6.25, fig.height=4.25, fig.align='center',echo=FALSE>>=
smd.mat<-ExtractSmd(tabMatched)
colnames(smd.mat)<-c("Match")
@


Estimation for the ATE:
<<NH15,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
ps.lr.match.ate <- Match(Y=small.nhanes$BPSysAve,Tr=small.nhanes$SmokeNow,
                            X=small.nhanes$ps.lr,estimand="ATE",ties=FALSE)
summary(ps.lr.match.ate)
@

Estimation for the ATT:
<<NH16,comment='+', fig.width=7.25, fig.height=4.5, fig.align='center'>>=
ps.lr.match.att <- Match(Y=small.nhanes$BPSysAve,Tr=small.nhanes$SmokeNow,
                            X=small.nhanes$ps.lr,estimand="ATT",ties=FALSE)
summary(ps.lr.match.att)
@
\end{document}
