\documentclass{beamer} \usepackage{verbatim} \usecolortheme{seagull} \usepackage{pstricks,pst-node,pst-tree,amscd,amsfonts,verbatim,geometry,bbm} \usepackage{fancyvrb,fancybox} \usepackage{color} \usepackage{xcolor} \usepackage{colortbl} \usepackage{graphicx} %\RequirePackage{etex} %\usepackage{etex} %\usepackage{xypic} \usepackage[all,color]{xy} \usepackage{natbib} \usepackage{tkz-graph} \usepackage{tikz} \usetikzlibrary{arrows,automata} \usetikzlibrary{graphs,arrows.meta} \usetikzlibrary{positioning} \usepackage{caption} \captionsetup[figure]{labelfont={it},textfont={color=black,it},labelformat={default},labelsep=period,name={}} \usepackage[latin1]{inputenc} \usepackage{undertilde} \usepackage{bm} \usepackage{amssymb} \usepackage{amsmath} \usepackage{lastpage} \usepackage{rotating} \usepackage{multirow} \usepackage{graphicx} \usepackage{bbm} \usepackage{bbold} \usepackage{upgreek} \usepackage{fancybox} \usepackage{fancyhdr} \usepackage{amsfonts} \usepackage{dsfont} %\usepackage{arydshln} \usepackage{listings} \usepackage{courier} \lstloadlanguages{R} %\lstset{basicstyle=\tiny,basicstyle=\ttfamily, numbers=none, % literate={~} {$\sim$}{2}} %\usepackage{mathptmx} %\usepackage{bookman} %\usepackage[T1]{fontenc} %\usepackage{newpxtext} %\usepackage[T1]{fontenc} %\usepackage[urw-garamond]{mathdesign} %\usepackage{garamondx} %\usepackage[sc]{mathpazo} %\linespread{1.05} % Palladio needs more leading (space between lines) %\usepackage[T1]{fontenc} %\usepackage{arev} %\usepackage[T1]{fontenc} \usepackage[scaled=.85]{DejaVuSerifCondensed} \usepackage[T1]{fontenc} \usefonttheme{serif} \DeclareFontFamily{OT1}{pzc}{} \DeclareFontShape{OT1}{pzc}{m}{it}{<-> s * [1.30] pzcmi7t}{} \DeclareMathAlphabet{\mathpzc}{OT1}{pzc}{m}{it} %\theoremstyle{myNote} %\newtheorem{myNote}{Note} \usepackage{mathtools} \usepackage{bookmark} \usepackage{amsthm} \usepackage{ulem} \usepackage{soul} \usepackage{hanging} \setbeamertemplate{footnote}{\hangpara{2em}{1}\makebox[2em][l]{\insertfootnotemark}\footnotesize\insertfootnotetext\par} \addtobeamertemplate{footnote}{}{\vspace{5ex}} \usepackage{calc} \newsavebox\CBox \newcommand\hcancel[2][0.25pt]{% \ifmmode\sbox\CBox{$#2$}\else\sbox\CBox{#2}\fi% \makebox[0pt][l]{\usebox\CBox}% \rule[0.5\ht\CBox-#1/2]{\wd\CBox}{#1}} \newcommand{\msout}[1]{\text{\sout{\ensuremath{#1}}}} \makeatletter \newcommand*{\shifttext}[2]{% \settowidth{\@tempdima}{#2}% \makebox[\@tempdima]{\hspace*{#1}#2}% } \makeatother \newcommand\oldmyzeta[1]{{\textcolor{red}{\protect\mathpzc{z}_{#1}} \shifttext{-10pt}{\raisebox{0.1ex}{\tiny{\sim}}}}} \newcommand\firstmyzeta[1]{\textcolor{black}{{\protect\mathpzc{z}_{#1}}}} \newcommand\myzeta[1]{{\textcolor{mcgillred}{{\protect\textsf{z}_{#1}}}}} \newcommand\tildemyzeta[1]{\textcolor{mcgillred}{{\protect\widetilde{\textsf{z}}_{#1}}}} \theoremstyle{plain} \newtheorem{thm}{Theorem} % reset theorem numbering for each chapter \theoremstyle{definition} \newtheorem{defn}[thm]{Definition} % definition numbers are dependent on theorem numbers \newtheorem{Note}{Note} % same for example numbers \setbeamertemplate{theorems}[numbered] \def\U{\mathcal{U}} \def\u{\upsilon} \def\Ind{\mathds{1}} \def\transpose{\top} \def\Xvec{{\skew0\bm{X}}} \def\Xveca{{\skew0\bm{X}}_1} \def\Xvecb{{\skew0\bm{X}}_2} \def\Yvec{{\skew0\bm{Y}}} \def\xvec{{\skew0\bm{x}}} \def\xveca{{\skew0\bm{x}}_1} \def\xvecb{{\skew0\bm{x}}_2} \def\mvec{{\skew0\bm{m}}} \def\muvec{{\skew0\bm{\mu}}} \def\cvec{{\skew0\mathbf{c}}} \def\Xbar{\overline{X}} \def\xbar{\overline{x}} \def\ybar{\overline{y}} \def\yhat{\widehat{y}} \def\betahat{\widehat{\beta}} \def\betavec{\utilde{\beta}} \def\betahatvec{\skew3\widehat{\utilde{\beta}}} \def\sigmahat{\widehat{\sigma}} \def\ehat{\widehat{e}} \def\zhat{\widehat{z}} \def\nhat{\widehat{n}} \def\xbarA{\overline{x}_A} \def\xbarB{\overline{x}_B} \def\calR{\mathcal{R}} \def\calS{\mathcal{S}} \def\calE{{\mathcal{S}_{e}}} \def\calB{\mathcal{B}} \def\calU{\mathcal{U}} \def\calI{\mathcal{I}} \def\B{\mathcal{B}} \def\Qbb{\mathbb{Q}} \def\Pbb{\mathbb{P}} \def\E{\mathbbm{E}} \def\Ehat{\widehat{\mathbbm{E}}} \def\bhatint{\widehat{\beta}_0} \def\bhatsl{\widehat{\beta}_1} \def\d{\:\textrm{d}} \def\ps{{\textsf{e}}} \def\bs{{\textsf{b}}} \def\bX{\bm{X}} \def\Ident{\bm{I}} \newcommand\xbarBi[1]{\overline{x^{\text{\tiny{(B)}}}_{#1}}} \newcommand\xbarib[1]{\overline{x^{\text{\tiny{(B)}}}_{#1}}} \def\Fstat{\frac{(SSE_R - SSE_C)/(k-g)}{SSE_C/(n-k-1)}} \def\approxsim{\sim \! \! \! : \;} \def\hsp{\vspace{0.1 in}} \def\fsp{\vspace{0.2 in}} %\def\dfrac#1#2{{\displaystyle {#1 #2}}}% \def\dint{\mathop{\displaystyle \int}}% \def\dsum{\mathop{\displaystyle \sum}}% \newcommand\betaZ{$\beta_{1}$} \newcommand\betaA[1]{$\beta_{#1}^{(A)}$} \newcommand\betaB[1]{$\beta_{#1}^{(B)}$} \newcommand\gammaAB[1]{$\gamma_{#1}^{(AB)}$} \def\SST{\mathit{SST}} \def\MST{\mathit{MST}} \def\SSE{\mathit{SSE}} \def\MSE{\mathit{MSE}} \def\SSB{\mathit{SSB}} \def\MSB{\mathit{MSB}} \def\SSI{\mathit{SSI}} \def\MSI{\mathit{MSI}} \def\SS{\mathit{SS}} \def\Ztilde{{\widetilde Z}} \def\ztilde{{\widetilde z}} \def\bmu{\boldsymbol\mu} \newcounter{Lecture} \def\calA{\mathcal{A}} \renewcommand{\thelecture}{\arabic{lecture}} \newcommand\xbari[1]{\overline{x}_{#1}} \newcommand\subsectnum[1]{\textbf{\arabic{section}.\arabic{subsection}} \textbf{\uppercase{#1}}} \newcommand\textsclarge[1]{\fontshape{sc}\fontsize{14}{15} \selectfont{#1}} \renewcommand\emph[1]{\textcolor{mcgillred}{\textit{#1}}} \def\Exp{\mathcal{E}} \def\Ex{\mathcal{E}} \def\Exzero{\mathcal{E}_0} \def\Obs{\mathcal{O}} \def\Ob{\mathcal{O}} \usepackage{tikz} \newcommand*\circled[1]{\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {#1};}} \def\WideSpacing{\itemsep=2.0ex\topsep=0.5ex\partopsep=0.0ex\parskip=0.0ex\parsep=0.0ex} \makeatletter \let\orig@Enumerate =\enumerate \newenvironment{enumerateWide}{\orig@Enumerate\WideSpacing}{\endlist} \makeatother \definecolor{midgray}{rgb}{0.75,0.75,0.75} \definecolor{lightgray}{rgb}{0.85,0.85,0.85} \definecolor{darkgray}{rgb}{0.654,0.724,0.795} \definecolor{lightdarkgray}{rgb}{0.727,0.805,0.883} \definecolor{purpleblue}{rgb}{0.50,0.50,1.00} \definecolor{lightblue}{rgb}{0.90,0.90,1.00} \definecolor{kugreen}{RGB}{50,93,61} \definecolor{kugreenlys}{RGB}{132,158,139} \definecolor{kugreenlyslys}{RGB}{173,190,177} \definecolor{kugreenlyslyslys}{RGB}{214,223,216} \definecolor{mcgillred}{rgb}{.926,0.105,0.184} \definecolor{mcgillreddark}{rgb}{0.617,0.035,0.094} \definecolor{mcgillredlight}{rgb}{0.926,0.153,0.241} \definecolor{verylightred}{rgb}{1,0.95,0.95} \definecolor{darkgreen}{rgb}{0.008, 0.294, 0.188} \definecolor{lightgreen}{rgb}{0.75,1.0,0.75} \definecolor{verylightgreen}{rgb}{0.90,1.0,0.90} \definecolor{darkblue}{rgb}{0.00, 0.00, 0.444} \definecolor{lightblue}{rgb}{0.933,0.933,1.00} \definecolor{psyc}{rgb}{0.87, 0.0, 1.0} \definecolor{mayablue}{rgb}{0.45, 0.76, 0.98} \definecolor{paleblue}{rgb}{0.69, 0.93, 0.93} \definecolor{spirodiscoball}{rgb}{0.06, 0.75, 0.99} \definecolor{lightcyan}{rgb}{0.88, 1.0, 1.0} \definecolor{frenchblue}{rgb}{0.0, 0.45, 0.73} \definecolor{bubbles}{rgb}{0.91, 1.0, 1.0} \definecolor{azurewm}{rgb}{0.94, 1.0, 1.0} \definecolor{airforceblue}{rgb}{0.36, 0.54, 0.66} \definecolor{darkmidnightblue}{rgb}{0.0, 0.2, 0.4} \definecolor{darkpastelblue}{rgb}{0.47, 0.62, 0.8} \definecolor{bleudefrance}{rgb}{0.19, 0.55, 0.91} \definecolor{babyblue}{rgb}{0.54, 0.81, 0.94} \setbeamercolor{background}{bg=red!15!white} \setbeamercovered{transparent} \setbeamercolor{title in sidebar}{fg=lightgray} \setbeamercolor{sidebar}{bg=lightgray,fg=purpleblue} \setbeamercolor{author in sidebar}{fg=lightgray} \setbeamercolor{section in sidebar}{fg=white,bg=mcgillred} \setbeamercolor{title}{fg=white,bg=mcgillred} \setbeamercolor{author}{use=structure,fg=black,bg=darkgray} \setbeamercolor{part name}{fg=black,bg=blue!30!white} \setbeamercolor{institute}{use=structure,fg=black,bg=lightdarkgray} \setbeamercolor{frametitle}{fg=white,bg=mcgillred} \setbeamercolor{subsection in sidebar}{fg=purpleblue} \setbeamercolor{block title}{use=structure,fg=black,bg=blue!30!white} \setbeamercolor{block body}{use=structure,fg=black,bg=blue!20!white} \setbeamertemplate{navigation symbols}{} \setbeamercovered{transparent} \mode \setbeamertemplate{part page} { \begin{centering} %\begin{beamercolorbox}[sep=16pt,center]{author} % \usebeamerfont{part title}\partname~\insertpartnumber \par %\end{beamercolorbox} \vskip1em\par \begin{beamercolorbox}[sep=16pt,center]{institute} \usebeamerfont{part title} \insertpart\par \end{beamercolorbox} \end{centering} } \makeatletter \AtBeginPart{% \beamer@tocsectionnumber=0\relax \setcounter{section}{0} \numberwithin{section}{part} {% start a group to keep the template change local \setbeamertemplate{background canvas}{% \color{white}\rule{\paperwidth}{\paperheight}% }% \setbeamercolor{frametitle}{fg=white,bg=white} \setbeamertemplate{footline}{ \leavevmode% \hbox{\begin{beamercolorbox}[wd=.5\paperwidth,ht=2.5ex,dp=1.125ex,leftskip=.3cm plus1fill,rightskip=.3cm]{author in head/foot}% \usebeamerfont{author in head/foot} \end{beamercolorbox}% \begin{beamercolorbox}[wd=.5\paperwidth,ht=2.5ex,dp=1.125ex,leftskip=.3cm,rightskip=.3cm plus1fil]{title in head/foot}% \usebeamerfont{title in head/foot}\hfill\insertpagenumber \end{beamercolorbox}}% \vskip0pt% } %\addtocounter{framenumber}{-1} \begin{frame}[noframenumbering,c]{\partpage}\end{frame}% }% end group } \makeatother \setbeamercolor{author in head/foot}{fg=white, bg=mcgillred} \setbeamercolor{title in head/foot}{fg=white, bg=mcgillred} \setbeamertemplate{footline} {% %\rule{\paperwidth}{0.025cm} \leavevmode% \hbox{\begin{beamercolorbox}[left, wd=.5\paperwidth,ht=2.5ex,dp=1.125ex,leftskip=.3cm]{title in head/foot}% \usebeamerfont{author in head/foot}\insertpartnumber.\insertsectionnumber: \insertsection \end{beamercolorbox}% \begin{beamercolorbox}[wd=.5\paperwidth,ht=2.5ex,dp=1.125ex,leftskip=.3cm,rightskip=.3cm plus1fil]{title in head/foot}% \usebeamerfont{title in head/foot}\hfill\insertpagenumber \end{beamercolorbox}}% \vskip0pt% } %\usecolortheme[named=kugreen]{structure} %\useinnertheme{circles} %\setbeamercovered{transparent} %\setbeamertemplate{blocks}[rounded][shadow=true] \setbeamertemplate{blocks}[default] \setbeamerfont{itemize/enumerate subbody}{size=\footnotesize} \setbeamertemplate{itemize items}[circle] \newenvironment<>{Ex}[1]{% \begin{actionenv}#2% \def\insertblocktitle{Example: #1}% \par% \mode{% \setbeamercolor{block title}{use=structure,fg=darkmidnightblue,bg=babyblue} \setbeamercolor{block body}{use=structure,fg=black,bg=bubbles} \setbeamertemplate{itemize item}{\scriptsize\raise1.25pt\hbox{\donotcoloroutermaths$\bullet$}} }% \usebeamertemplate{block begin}\justifying} {\par\usebeamertemplate{block end}\end{actionenv}} \newenvironment<>{Thm}[1]{% \begin{actionenv}#2% \def\insertblocktitle{Theorem: #1}% \par% \mode{% \setbeamercolor{block title}{use=structure,fg=white,bg=darkgreen} \setbeamercolor{block body}{use=structure,fg=black,bg=lightgreen} \setbeamercolor{itemize item}{fg=orange!20!black} \setbeamertemplate{itemize item}[circle] }% \usebeamertemplate{block begin}\justifying} {\par\usebeamertemplate{block end}\end{actionenv}} \newenvironment<>{Nt}[1]{% \begin{actionenv}#2% \def\insertblocktitle{Note}% \par% \mode{% \setbeamercolor{block title}{use=structure,fg=white,bg=darkblue} \setbeamercolor{block body}{use=structure,fg=black,bg=azurewm} \setbeamercolor{itemize item}{fg=orange!20!black} \setbeamertemplate{itemize item}[circle] }% \usebeamertemplate{block begin}\justifying} {\par\usebeamertemplate{block end}\end{actionenv}} \newcounter{myNote} \newcommand{\myNote}{\arabic{myNote}} \setcounter{myNote}{1} \renewcommand*{\thefootnote}{[\arabic{footnote}]} \renewenvironment<>{Note}{% \begin{actionenv}% \def\insertblocktitle{Note \myNote.}% \par% \mode{% \setbeamercolor{block title}{use=structure,fg=white,bg=darkblue} \setbeamercolor{block body}{use=structure,fg=black,bg=lightblue} \setbeamercolor{itemize item}{fg=orange!20!black} \setbeamertemplate{itemize item}[circle] }% \usebeamertemplate{block begin}\addtocounter{myNote}{1}\justifying} {\par\usebeamertemplate{block end}\end{actionenv}} \date{ } % % The following info should normally be given in you main file: % \usepackage{tikz} \newcommand{\topline}{% \tikz[remember picture,overlay] {% \draw[midgray] ([yshift=-1cm]current page.north west) -- ([yshift=-1cm,xshift=\paperwidth]current page.north west);}} \setbeamertemplate{frametitle}{% \nointerlineskip \begin{beamercolorbox}[sep=0.3cm,ht=1.8em,wd=\paperwidth,rightskip=0cm]{frametitle}% \usebeamerfont{frametitle}\usebeamercolor[fg]{frametitle} \vbox{}\vskip-2ex% \strut\insertframetitle\strut \vskip-1.2ex% \end{beamercolorbox}% } \newcommand\mypart[1]{#1} %\setbeamertemplate{frametitle continuation}[from %second][(\lowercase{\insertcontinuationcountroman})] \setbeamertemplate{frametitle continuation}[from second][] \setbeamertemplate{section in toc}{\insertpartnumber.\inserttocsectionnumber\hspace*{1em}\inserttocsection} \title{Propensity Score Methods, Models and Adjustment} \author{Dr David A. Stephens} \institute{ Department of Mathematics \& Statistics\\ McGill University\\ Montreal, QC, Canada. \vspace{0.1 in} \texttt{david.stephens@mcgill.ca}\\ \texttt{\texttt{www.math.mcgill.ca/dstephens/SISCER2020/}} } \titlegraphic{\includegraphics[width=1.93cm,height=2cm]{McGillLogo}} \setbeamertemplate{itemize item}{\scriptsize\raise1.25pt\hbox{\donotcoloroutermaths$\bullet$}} \setbeamertemplate{itemize subitem}{\tiny\raise1.25pt\hbox{\donotcoloroutermaths$\blacktriangleright$}} \usepackage{etoolbox} %\usepackage{setspace} %\addtocontents{toc}{\protect\setstretch{0.9}} \usepackage{ragged2e} \apptocmd{\frame}{}{\justifying}{} % Allow optional arguments after frame. \usepackage[bookmarks=true]{hyperref} \lstset{basicstyle=\ttfamily\tiny, numbers=left, numberstyle=\tiny, stepnumber=1, numbersep=5pt} \begin{document} \part{NHANES Example} \setbeamercolor{author in head/foot}{fg=white, bg=mcgillred} \setbeamercolor{title in head/foot}{fg=white, bg=mcgillred} { \setbeamertemplate{footline} {% \leavevmode% \hbox{\begin{beamercolorbox}[wd=.5\paperwidth,ht=2.5ex,dp=1.125ex,leftskip=.3cm plus1fill,rightskip=.3cm]{author in head/foot}% \usebeamerfont{author in head/foot} \end{beamercolorbox}% \begin{beamercolorbox}[wd=.5\paperwidth,ht=2.5ex,dp=1.125ex,leftskip=.3cm,rightskip=.3cm plus1fil]{title in head/foot}% \usebeamerfont{title in head/foot}\hfill\insertpagenumber \end{beamercolorbox}}% \vskip0pt% } <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE, size = "tiny") options(scipen=5) options(repos=c(CRAN="https://cloud.r-project.org/")) @ %\frame{\titlepage} \begin{frame}[fragile,allowframebreaks]\frametitle{NHANES example} We study the publicly available (U.S.) National Health and Nutrition Examination Survey (NHANES). \begin{itemize} \item \texttt{BPSysAve} (systolic blood pressure, \emph{outcome}), \item \texttt{SmokeNow} (smoking status: 0--No, 1--Yes, \emph{exposure}), \item \texttt{Gender}, \item \texttt{Age}, \item \texttt{Race3}, \item \texttt{Education}, \item \texttt{MaritalStatus}, and \item \texttt{Poverty} \end{itemize} We restrict our attention to 1377 adults ($>$ 17 years old) in the second wave of the survey. \framebreak First the R packages are loaded <>= #Libraries and dependencies muse be pre-loaded in R suppressPackageStartupMessages({ library(NHANES) #Contains the NHANES data library(tableone) #For comparing the treated and untreated library(Matching) #For comparison based on matching library(survey) #For handling reweighted samples library(splines) #For fitting a flexible PS model library(Hmisc) #Some plotting functions library(htmlwidgets) #Some plotting functions }) @ \framebreak <>= #Selecting the subset of data for analysis small.nhanes <- na.omit(NHANES[NHANES$SurveyYr=="2011_12" & NHANES$Age > 17,c(3,4,8:11,13,25,61)]) small.nhanes$SmokeNow <- as.numeric(small.nhanes$SmokeNow)-1 #dim(small.nhanes) ## 1377 vars <- c("Gender", "Age", "Race3", "Education", "MaritalStatus", "Poverty") @ \framebreak Regression: <>= fit0<-lm(BPSysAve~SmokeNow,data=small.nhanes) round(coef(summary(fit0)),5) fmod<-formula(BPSysAve~SmokeNow+Gender+Age+Race3+ Education+MaritalStatus+HHIncome+Poverty) fit1<-lm(fmod,data=small.nhanes) round(coef(summary(fit1))[1:2,],5) @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{Checking confounder balance} Typical ways of assessing balance include inspection of \begin{itemize} \item \emph{Standardized mean difference} (SMD) or proportion: \[ \frac{\bar{x}^{1,w}-\bar{x}^{0,w}}{\sqrt{(v^{1,w}+v^{0,w})/2}} \qquad \text{where}, \] where, if $w_i$ is a case weight, the quantities \[ \bar{x}^{z,w} \quad \text{and} \quad v^{z,w} \] are the weighted sample mean and sample variance of variable $X$ among those with treatment value $z$; SMD $< 0.1$ indicates reasonable balance. \item Empirical CDFs in the treated and untreated groups. \end{itemize} \framebreak There is clear Age imbalance between the two groups. <>= tabUnmatched<-CreateTableOne(vars=vars[2],strata="SmokeNow", data=small.nhanes,test=FALSE) print(tabUnmatched, smd = TRUE) @ \texttt{CreateTableOne} can be used for all predictors <>= tabUnmatched<-CreateTableOne(vars=vars,strata="SmokeNow", data=small.nhanes,test=FALSE) @ \framebreak <>= 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) @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{Building the propensity score} We attempt to build and assess a propensity score model using the available covariates <>= #Fit the PS using logistic regression 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") @ \framebreak Investigating the PS quintiles: <>= ps0<-ps.lr[Smoke==0];ps1<-ps.lr[Smoke==1] quints <- c(0,quantile(ps.lr,seq(.2,1,.2))) ps.tab<-rbind(table(cut(ps.lr[Smoke==0],quints)), table(cut(ps.lr[Smoke==1],quints))) rownames(ps.tab)<-c('Smoke=0','Smoke=1') ps.tab @ The table indicates that there is good representation of the exposed and unexposed groups in the five quintile-based strata. \framebreak We look for overlap of the propensity score histograms. <>= par(mar=c(4,2,2,0)) colvec<-c(rgb(0,0,1,0.5),rgb(1,0,0,0.5)) 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.6,150,c('Smoke = 0', 'Smoke = 1'),col=colvec,pch=15) boxplot(ps0,ps1,ylab="PS",xlab="Treatment Group",ylim=range(0,1), names=c('Smoke=0','Smoke=1'),col=colvec); abline(h=quints[2:5],col="red",lty=2) @ \framebreak There seems to be good overlap. \medskip 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. \framebreak We can proceed to check for balance within PS quintiles <>= Pcat<-as.numeric(cut(ps.lr,quints,include.lowest=T)) 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('All','Q1','Q2','Q3','Q4','Q5') round(smd.mat,4) @ \framebreak 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. \framebreak <>= par(mar=c(2,3,2,0),mfrow=c(2,3),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) @ \framebreak <>= 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('All','Q1','Q2','Q3','Q4','Q5','Q6','Q7','Q8','Q9','Q10') @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{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. \end{itemize} \framebreak Create the PS-matched sample: <>= small.nhanes$ps.lr<-ps.lr ps.lr.match <- Match(Tr=small.nhanes$SmokeNow,X=small.nhanes$ps.lr, estimand="ATE",ties=FALSE) length(ps.lr.match$index.dropped) @ \medskip The matched sample can be checked for balance: <>= MatchBalance(SmokeNow~Gender+Age+Education+MaritalStatus+Poverty, data=small.nhanes,match.out=ps.lr.match) @ (results not shown). \framebreak <>= matched.samp <- small.nhanes[c(ps.lr.match$index.control, ps.lr.match$index.treated),] tabMatched<-CreateTableOne(vars = vars, strata = "SmokeNow", data = matched.samp, test = FALSE) @ \framebreak <>= 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) @ \framebreak Matching seems to have matched the age profiles for the smoking and non-smoking groups if we compare the SMDs: <>= smd.mat<-cbind(smd.mat,ExtractSmd(tabMatched)) colnames(smd.mat)<-c('All','Q1','Q2','Q3','Q4','Q5',"Match") round(smd.mat,4) @ \bigskip \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{Inverse Weighting} We may study an \emph{inverse weighted} sample using the \texttt{survey} and \texttt{Hmisc} packages. <>= 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) @ \framebreak <>= 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(tabIPW)) colnames(smd.mat)<-c('All','Q1','Q2','Q3','Q4','Q5','Match','IPW') round(smd.mat,4) @ Matching and weighting samples yield good balance. \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{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} \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{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 \& AIPW. \end{itemize} \framebreak Note that the IPW and AIPW estimators can suffer from the fact that the inverse weights based on the propensity score can get large if the propensity score is near zero or one. \medskip For these estimators, \emph{truncation} can be used, where propensity score values more extreme than a given threshold (0.001 or 0.999 say) are set equal to these thresholds. \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{Outcome regression} \textbf{Outcome regression, unadjusted:} <>= coef(summary(lm(BPSysAve~SmokeNow,data=small.nhanes)))[2,] @ \framebreak \textbf{Outcome regression no interaction:} We use the model with other predictors included additively: <>= fmod<-formula(BPSysAve~SmokeNow+Gender+Age+Race3+Education+ MaritalStatus+HHIncome+Poverty) coef(lm(fmod,data=small.nhanes))[2] @ <>= nhanes.allsmoke <- small.nhanes nhanes.allsmoke$SmokeNow <- 1 nhanes.nosmoke <- small.nhanes nhanes.nosmoke$SmokeNow <- 0 mod1.lm <- lm(fmod,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 #Note: same as for previous answer regression @ The naive conditional effect estimate is more than 3 times greater than its confounder-adjusted counterpart. \framebreak \textbf{Outcome regression with interactions:} We now add interaction terms <>= fmodi<-formula(BPSysAve~SmokeNow+Gender+Age+Race3+Education+ MaritalStatus+HHIncome+Poverty+ SmokeNow:HHIncome+SmokeNow:Gender+SmokeNow:Age) mod1.lmX <- lm(fmodi,data=small.nhanes) @ \framebreak <>= round(coef(summary(mod1.lmX))[-c(3:30),],3) @ \framebreak Here the coefficients alone \emph{do not} provide the ATE estimate, and we need to use the \texttt{predict} method to compute the ATE estimate as \[ \frac{1}{n} \sum_{i=1}^n m(x_i,1;\widehat \beta) - \frac{1}{n} \sum_{i=1}^n m(x_i,0;\widehat \beta) \] \framebreak <>= APO.lmX.1 <- mean(predict(mod1.lmX,nhanes.allsmoke)) APO.lmX.0 <- mean(predict(mod1.lmX,nhanes.nosmoke)) APO.lmX.1 - APO.lmX.0 @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{PS stratification} We first fit the propensity score model based on logistic regression and then use \emph{stratification} based on PS quintiles. <>= 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") 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]) } sum(ATE.strat*p.strat) @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{PS matching} <>= 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),] mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 1]) - mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 0]) @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{PS regression} We now attempt three propensity score regression analyses based on increasingly complex specifications \[ \begin{array}{lll} \text{PSR 1} & Y & = \ \beta_0 + \beta_1 e(X;\widehat \alpha) + \psi_1 Z + \epsilon \\[6pt] \text{PSR 2} & Y & = \ \beta_0 + \beta_1 e(X;\widehat \alpha) + \beta_2 \{e(X;\widehat \alpha) \}^2 + \\ && \qquad \qquad \psi_1 Z + \psi_2 Z e(X;\widehat \alpha) + \psi_3 Z \{e(X;\widehat \alpha)\} + \epsilon \\[6pt] \text{PSR 3} & Y & = \ s(e(X;\widehat \alpha)) + \psi_1 Z +\psi_2 Z s(e(X;\widehat \alpha)) + \epsilon \end{array} \] where $s(.)$ is a flexible spline function. \framebreak PSR 1: <>= 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 @ \framebreak PSR 2: <>= psmod2<-formula(BPSysAve~(1+SmokeNow)*(ps.lr+I(ps.lr^2))) mod1.PSlm2 <- lm(psmod2,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 @ \framebreak PSR 3: <>= psmod3<-formula(BPSysAve~ (1+SmokeNow)*bs(ps.lr,df=4,Boundary.knots = range(0,1))) mod1.PSlm3 <- lm(psmod3,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 @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{Inverse Weighting} The IPW estimator can be computed using a direct calculation based on the weight \[ w_i = \frac{Z_i}{e(X_i;\widehat \alpha)} + \frac{(1-Z_i)}{1-e(X_i;\widehat \alpha)} \] computed after the PS model has been estimated. \framebreak \textcolor{red}{IPW0}: Ordinary weight form \[ \frac{1}{n} \sum_{i=1}^n w_i Z_i Y_i - \frac{1}{n} \sum_{i=1}^n w_i (1-Z_i) Y_i \] \medskip \textcolor{red}{IPW1}: Standardized weight form \begin{equation*} \frac{1}{\sum\limits_{i=1}^n w_i} \left(\sum_{i=1}^n w_i Z_i Y_i - \sum_{i=1}^n w_i (1-Z_i) Y_i \right) \end{equation*} obtained via weighted least squares. \framebreak <>= ps.lr.weight <- Smoke/ps.lr + (1-Smoke)/(1-ps.lr) mean(Smoke*Y*ps.lr.weight) -mean((1-Smoke)*Y*ps.lr.weight) #IPW0 coef(lm(Y ~ Smoke, weights = ps.lr.weight))[2] #IPW1 @ \framebreak The augmented IPW (AIPW) estimator adds a proposed mean model for the APOs, that is \begin{align*} \widehat \mu(1) & = \dfrac{1}{n} \sum_{i=1}^n \dfrac{Z_i (Y_i-m(X_i,1;\widehat \beta))}{e(X_i;\widehat \alpha)} + \dfrac{1}{n} \sum_{i=1}^n m(X_i,1;\widehat \beta) \\[6pt] \widehat \mu(0) & = \dfrac{1}{n} \sum_{i=1}^n \dfrac{(1-Z_i) (Y_i-m(X_i,0;\widehat \beta))}{1-e(X_i;\widehat \alpha)} + \dfrac{1}{n} \sum_{i=1}^n m(X_i,0;\widehat \beta) \end{align*} or computed via an augmented outcome regression with the two additional predictors \[ \dfrac{Z_i}{e(X_i;\widehat \alpha)} \qquad \dfrac{(1-Z_i)}{1-e(X_i;\widehat \alpha)}. \] \framebreak Here we use the outcome model with interactions fitted previously: <>= fmodi<-formula(BPSysAve~SmokeNow+Gender+Age+Race3+Education+ MaritalStatus+HHIncome+Poverty+ SmokeNow:HHIncome+SmokeNow:Gender+SmokeNow:Age) mod1 <- lm(fmodi,data=small.nhanes) m1<-predict(mod1,nhanes.allsmoke) m0<-predict(mod1,nhanes.nosmoke) ps.trunc<-pmin(pmax(0.001,ps.lr),0.999) APO.1<-mean(Smoke*(Y-m1)/ps.trunc)+mean(m1) APO.0<-mean((1-Smoke)*(Y-m0)/(1-ps.trunc))+mean(m0) APO.1-APO.0 @ Here we use the \emph{augmented outcome regression} (AOR) approach \framebreak <>= small.nhanes$ps <- ps.lr fmodia<-formula(BPSysAve~SmokeNow+Gender+Age+Race3+Education+ MaritalStatus+HHIncome+Poverty+ SmokeNow:HHIncome+SmokeNow:Gender+SmokeNow:Age+ I(SmokeNow/ps)+I((1-SmokeNow)/(1-ps))) mod1.lmX <- lm(fmodia,data=small.nhanes) nhanes.allsmoke <- nhanes.nosmoke <- small.nhanes; nhanes.allsmoke$SmokeNow <- 1 nhanes.nosmoke$SmokeNow <- 0 mean(predict(mod1.lmX,nhanes.allsmoke))- mean(predict(mod1.lmX,nhanes.nosmoke)) @ \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{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} \framebreak \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; however, that the bootstrap is \textit{not} in general valid for matching. \end{itemize} \end{frame} \begin{frame}[fragile,allowframebreaks]\frametitle{Bootstrap analysis} 1000 bootstrap samples are used to approximate the sampling distribution of the various estimators. <>= nboot<-1000 set.seed(45) ests.mat<-matrix(0,nrow=nboot,ncol=10) n<-nrow(small.nhanes) library(splines) for(ib in 1:nboot){ ind<-sample(1:n,size=n,rep=T) boot.nhanes<-small.nhanes[ind,] boot.allsmoke <- boot.nosmoke <- boot.nhanes boot.allsmoke$SmokeNow <- 1 boot.nosmoke$SmokeNow <- 0 ests.mat[ib,1]<-coef(lm(BPSysAve~SmokeNow,data=boot.nhanes))[2] ests.mat[ib,2]<-coef(lm(fmod,data=boot.nhanes))[2] mod1.lmX <- lm(fmodi,data=boot.nhanes) APO.lmX.1 <- mean(predict(mod1.lmX,boot.allsmoke)) APO.lmX.0 <- mean(predict(mod1.lmX,boot.nosmoke)) ests.mat[ib,3]<-APO.lmX.1 - APO.lmX.0 ps.mod <- glm(fmods,data=boot.nhanes,family="binomial") ps.lr <- predict(ps.mod,type="response") ps0<-ps.lr[boot.nhanes$SmokeNow==0];ps1<-ps.lr[boot.nhanes$SmokeNow==1] quints <- c(0,quantile(ps.lr,seq(.2,1,.2))) Y<-boot.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[boot.nhanes$SmokeNow == 1 & ps.lr.quints==j])- mean(Y[boot.nhanes$SmokeNow == 0 & ps.lr.quints==j]) } ests.mat[ib,4]<-sum(ATE.strat*p.strat) boot.nhanes$ps.lr<-ps.lr ps.lr.match <- Match(Tr=boot.nhanes$SmokeNow,X=boot.nhanes$ps.lr,estimand="ATE",ties=FALSE) matched.samp <- boot.nhanes[c(ps.lr.match$index.control,ps.lr.match$index.treated),] ests.mat[ib,5]<-mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 1]) - mean(matched.samp$BPSysAve[matched.samp$SmokeNow == 0]) mod1.PSlm1 <- lm(BPSysAve~SmokeNow+ps.lr,data=boot.nhanes) APO.PSlm1.1 <- mean(predict(mod1.PSlm1,boot.allsmoke)) APO.PSlm1.0 <- mean(predict(mod1.PSlm1,boot.nosmoke)) ests.mat[ib,6]<-APO.PSlm1.1 - APO.PSlm1.0 mod1.PSlm2 <- lm(BPSysAve~(1+SmokeNow)*(ps.lr+I(ps.lr^2)),data=boot.nhanes) APO.PSlm2.1 <- mean(predict(mod1.PSlm2,boot.allsmoke)) APO.PSlm2.0 <- mean(predict(mod1.PSlm2,boot.nosmoke)) ests.mat[ib,7]<-APO.PSlm2.1 - APO.PSlm2.0 mod1.PSlm3 <- lm(BPSysAve~(1+SmokeNow)*bs(ps.lr,df=4,Boundary.knots = range(0,1)), data=boot.nhanes) APO.PSlm3.1 <- mean(predict(mod1.PSlm3,boot.allsmoke)) APO.PSlm3.0 <- mean(predict(mod1.PSlm3,boot.nosmoke)) ests.mat[ib,8]<-APO.PSlm3.1 - APO.PSlm3.0 ps.lr.weight <- boot.nhanes$SmokeNow/ps.lr + (1-boot.nhanes$SmokeNow)/(1-ps.lr) ests.mat[ib,9]<-coef(lm(Y ~ boot.nhanes$SmokeNow, weights = ps.lr.weight))[2] boot.allsmoke <- boot.nosmoke <- boot.nhanes boot.allsmoke$SmokeNow <- 1 boot.nosmoke$SmokeNow <- 0 mod1 <- lm(fmodi,data=boot.nhanes) m1<-predict(mod1,boot.allsmoke) m0<-predict(mod1,boot.nosmoke) Smoke<-boot.nhanes$SmokeNow ps.trunc<-pmin(pmax(0.001,ps.lr),0.999) APO.1<-mean(Smoke*(Y-m1)/ps.trunc)+mean(m1) APO.0<-mean((1-Smoke)*(Y-m0)/(1-ps.trunc))+mean(m0) ests.mat[ib,10]<-APO.1-APO.0 } @ \framebreak <>= par(mar=c(3,2,1,0)) lvec<-c('Reg.','Reg. 2','Reg. 3','Strat.','Match.','PSR 1','PSR 2','PSR 3','IPW','AIPW') boxplot(ests.mat,names=lvec,pch=19,cex=0.5,ylim=range(-10,5),cex.axis=0.75) @ \framebreak The boxplot demonstrates largely similar estimates from quite different analysis methods. The estimates and standard errors are as follows. <>= ATE<-apply(ests.mat,2,mean) ATE.se<-sqrt(apply(ests.mat,2,var)) ATE.res<-cbind(ATE,ATE.se) rownames(ATE.res)<-lvec colnames(ATE.res)<-c('Est.','s.e.') ATE.res @ \end{frame} \end{document}