Estimation_of_Causal_Effects.Rnw

\RequirePackage[l2tabu, orthodox]{nag} % warn about outdated packages
\documentclass[11pt,leqno]{article}
\usepackage{microtype} %
\usepackage{setspace}
\onehalfspacing
\usepackage{xcolor, color, ucs}     % http://ctan.org/pkg/xcolor
\usepackage{natbib}
\usepackage{booktabs}          % package for thick lines in tables
\usepackage{amsfonts,amsthm,amsmath,amssymb}          % AMS Stuff
\usepackage[linewidth=1pt]{mdframed}
\usepackage{mdframed}
\usepackage{empheq}            % To use left brace on {align} environment
\usepackage{graphicx}          % Insert .pdf, .eps or .png
\usepackage{enumitem}          % http://ctan.org/pkg/enumitem
\usepackage[mathscr]{euscript}          % Font for right expectation sign
\usepackage{tabularx}          % Get scale boxes for tables
\usepackage{float}             % Force floats around
\usepackage{afterpage}% http://ctan.org/pkg/afterpage
\usepackage[T1]{fontenc}
\usepackage{rotating}          % Rotate long tables horizontally
\usepackage{bbm}                % for bold betas
\usepackage{csquotes}           % \enquote{} and \textquote[][]{} environments
\usepackage{subfig}
\usepackage{titling}            % modify maketitle in latex
% \usepackage{mathtools}          % multlined environment with size option
\usepackage{verbatim}
\usepackage{geometry}
\usepackage{bigfoot}
\usepackage[format=hang,
            font={small},
            labelfont=bf,
            textfont=rm]{caption}

\geometry{verbose,margin=2cm,nomarginpar}
\setcounter{secnumdepth}{2}
\setcounter{tocdepth}{2}

\usepackage{url}
\usepackage{relsize}            % \mathlarger{} environment
\usepackage[unicode=true,
            pdfusetitle,
            bookmarks=true,
            bookmarksnumbered=true,
            bookmarksopen=true,
            bookmarksopenlevel=2,
            breaklinks=false,
            pdfborder={0 0 1},
            backref=page,
            colorlinks=true,
            hyperfootnotes=true,
            hypertexnames=false,
            pdfstartview={XYZ null null 1},
            citecolor=blue!70!black,
            linkcolor=red!70!black,
            urlcolor=green!70!black]{hyperref}
\usepackage{hypernat}

\usepackage{multirow}
\usepackage[noabbrev]{cleveref} % Should be loaded after \usepackage{hyperref}

\parskip=12pt
\parindent=0pt
\delimitershortfall=-1pt
\interfootnotelinepenalty=100000

\makeatletter
\def\thm@space@setup{\thm@preskip=0pt
\thm@postskip=0pt}
\makeatother

\makeatletter
% align all math after the command
\newcommand{\mathleft}{\@fleqntrue\@mathmargin\parindent}
\newcommand{\mathcenter}{\@fleqnfalse}
% tilde with text over it
\newcommand{\distas}[1]{\mathbin{\overset{#1}{\kern\z@\sim}}}%
\newsavebox{\mybox}\newsavebox{\mysim}
\newcommand{\distras}[1]{%
  \savebox{\mybox}{\hbox{\kern3pt$\scriptstyle#1$\kern3pt}}%
  \savebox{\mysim}{\hbox{$\sim$}}%
  \mathbin{\overset{#1}{\kern\z@\resizebox{\wd\mybox}{\ht\mysim}{$\sim$}}}%
}
\makeatother

\newtheoremstyle{newstyle}
{12pt} %Aboveskip
{12pt} %Below skip
{\itshape} %Body font e.g.\mdseries,\bfseries,\scshape,\itshape
{} %Indent
{\bfseries} %Head font e.g.\bfseries,\scshape,\itshape
{.} %Punctuation afer theorem header
{ } %Space after theorem header
{} %Heading

\theoremstyle{newstyle}
\newtheorem{thm}{Theorem}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{lem}{Lemma}
\newtheorem{cor}{Corollary}
\newcommand*\diff{\mathop{}\!\mathrm{d}}
\newcommand*\Diff[1]{\mathop{}\!\mathrm{d^#1}}
\newcommand*{\QEDA}{\hfill\ensuremath{\blacksquare}}%
\newcommand*{\QEDB}{\hfill\ensuremath{\square}}%
\DeclareMathOperator{\E}{\mathbb{E}}
\DeclareMathOperator{\Var}{\rm{Var}}
\DeclareMathOperator{\Cov}{\rm{Cov}}
% \DeclareMathOperator{\Pr}{\rm{Pr}}

% COLORS FOR GRAPHICS (3-class Set1)
\definecolor{Blue}{RGB}{55,126,184}
\definecolor{Red}{RGB}{228,26,28}
\definecolor{Green}{RGB}{77,175,74}

% COLORS FOR EQUATIONS (3-class Dark2)
\definecolor{eqgreen}{RGB}{27,158,119}
\definecolor{eqblue}{RGB}{117,112,179}
\definecolor{eqred}{RGB}{217,95,2}


\title{Estimation of Causal Effects}
\author{Jake Bowers, Ben Hansen \& Tom Leavitt}
\date{\today}

\begin{document}

\maketitle

\tableofcontents

<< setup, include = FALSE, eval = TRUE >>=

if (!require("pacman")) { install.packages("pacman") }

pacman::p_load(
  plyr, # data wrangling
  dplyr, # data manipulation
  magrittr, # pipes
  tidyr, # gather and other data reshaping
  haven, # read_dta, _sas, _sav files
  broom, # tidy() and glance() 
  randomizr, # performs different random assignments
  lmtest,
  ggplot2,
  knitr,
  xtable,
  devtools, # sourcing and working with github repos
  optmatch,
  RItools,
  brglm
)

opts_chunk$set(tidy = TRUE,
               echo = TRUE,
               results = 'markup',
               strip.white = TRUE,
               fig.path = 'figs/fig',
               cache = FALSE,
               highlight = TRUE,
               width.cutoff = 100,
               size = 'footnotesize',
               out.width = '1\\textwidth',
               message = FALSE,
               warning = FALSE,
               comment = NA)

options(width = 100,
        digits = 3)

@

\newpage

\section{Estimation of Mean Unit-Level Causal Effects}

<< >>=

## Load data
news_df <- read.csv("http://jakebowers.org/PS531Data/news.df.csv")

## Create potential outcomes
news_df %<>% rename(y = r) %>%
  mutate(yc = ifelse(test = z == 0,
                     yes = y,
                     no = NA),
         yt = ifelse(test = z == 1,
                     yes = y,
                     no = NA))

kable(news_df[, c(1, 3:4, 8:9, 7)])

@

\vspace{5mm}
\begin{mdframed}
\textbf{Question for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Now let's imagine the \textit{true} unit-level treatment effect is some specific positive natural number that differs for almost every unit. What do we mean by ``\textit{true} unit-level treatment effect'' as opposed to a unit-level treatment effect that we hypothesize and subsequently assume to be true in order to assess evidence in favor of or against that hypothesis?
\end{itemize}
\end{mdframed}

<< >>=

news_df %<>% mutate(true_tau = c(6, 4, 19, 3, 9, 9, 13, 15),
                    yc = ifelse(test = is.na(yc),
                                yes = yt - true_tau,
                                no = yc),
                    yt = ifelse(test = is.na(yt),
                                yes = yc + true_tau,
                                no = yt))

## Other ways to get the true ATE if you knew the true individual level effects.
## true_ate<-with(news_df,mean(yt) - mean(yc))
## true_ate<-news_df %$% {mean(yt) - mean(yc)}

true_ate <- news_df %$% mean(true_tau)

true_ate

kable(news_df[, c(1, 3:4, 8:9, 7)])

@

\vspace{5mm}
\begin{mdframed}
\textbf{Question for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item I claim that the difference-in-means estimator and the stratified difference-in-means estimator are both good estimators. Name some criteria that people use to assess the quality of estimators . . .
\end{itemize}
\end{mdframed}

In the code chunk below, we create a function that calculates $\widehat{ATE}$ on the outcomes we observe under different treatment assignment permutations.

<< >>=

treatment_permutations <- function(z, yc, yt, s){

	if(missing(s)){
		## If this is not a block-randomized experiment then do the following

		## Permute treatment assignment
		## Student Question: What is the difference between Z and z? 
	  ## Why is one uppercase and the other lower?
		Z = sample(z)

		Y = Z * yt + (1 - Z) * yc

		## Calculate unstratified test statistic
		ate_hat_unstrat = coef(lm(Y ~ Z))[["Z"]]

		return(ate_hat_unstrat)

	}

	else { # If the experiment is block randomized

		## Permute treatment assignment WITHIN blocks
		Z = unsplit(lapply(split(x = z, f = s), sample), s)

		Y = Z * yt + (1 - Z) * yc

		## Calculate stratified test-statistic
		ate_hat_strat = coef(lm(Y ~ Z + s))[["Z"]]

		return(ate_hat_strat)
	}

}

@

\vspace{5mm}
\begin{mdframed}
\textbf{Question for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Describe in words what the function above is doing.
\end{itemize}
\end{mdframed}

Now we are repeatedly estimating the ATE under 1000 treatment assignment permutations in order to assess unbiasedness. 

\vspace{5mm}
\begin{mdframed}
\textbf{Question for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Which two quantities would we like to compare in order to assess whether or not the difference-in-means estimator is unbiased?
\end{itemize}
\end{mdframed}

<< >>=

## Set seed for simulations
set.seed(1:5)

## If we use blocks
block_randomization_distribution <- data.frame(ate = replicate(10^3,
							       treatment_permutations(z = news_df$z,
										      yc = news_df$yc,
										      yt = news_df$yt,
										      s = news_df$s)))

colMeans(block_randomization_distribution)

## If we don't use blocks
randomization_distribution <- data.frame(ate = replicate(10^3,
							 treatment_permutations(z = news_df$z,
										yc = news_df$yc,
										yt = news_df$yt)))

colMeans(randomization_distribution)

block_randomization_distribution$type <- 'Stratified'

randomization_distribution$type <- 'Unstratified'

randomization_dists <- as.data.frame(rbind(block_randomization_distribution, randomization_distribution))

ggplot(randomization_dists, aes(x = ate, fill = type)) +
	geom_histogram(alpha = 0.5, aes(), position = 'identity') +
	xlab("Estimated ATEs") +
	ylab("Count") +
	ggtitle("Randomization Distributions") +
	scale_fill_discrete(name = "Random \nAssignment")

## If you want to add the observed ate hat
# +  geom_vline(xintercept = obs_ate_hat, colour = "black", linetype = "longdash")

@

Notice that the randomization distributions look a little strange; that is, both randomization distributions are very different from the normal distribution to which we are accustomed. This is important. Nothing so far has invoked any claim about the shape of the randomization distribution, such as its normality if the finite population central limit theorem obtains \citep{hajek1960}. What matters (in terms of unbiasedness) is whether the mean of the randomization distribution---i.e., the mean of all the estimated
difference in means---is equal to the true average treatment effect.\footnote{Notice that we do not have simulated results that are \emph{exactly}	equal to the truth. Why is that? What might we do to move them closer to the truth? If I ran this simulation with one seed, and you ran it with another seed, by how much would we expect our two averages to	differ?}

\vspace{5mm}
\begin{mdframed}
\textbf{Questions for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Is the variance of the randomization distribution tighter when treatment is assigned within blocks or when treatment is \textit{not} assigned within blocks? And why?
\item Our \textit{single} $\widehat{ATE}$ from the experiment is our ``best guess'' about the true ATE. But what do we mean by ``best guess''? And can one's ``best guess'' be misleading?
\end{itemize}
\end{mdframed}

<< >>=

## Let's calculate out single observed test statistic that accounts for blocks
obs_block_ate_hat <- coef(lm(y ~ z + s, data = news_df))[["z"]]

## This is our SINGLE estimated ATE that we actually observe
obs_block_ate_hat

## This is the expected value of the estimated ATE
mean(block_randomization_distribution$ate)

## This is the true ATE (not estimated) that we are trying to estimate
true_ate

@


\vspace{5mm}
\begin{mdframed}
\textbf{Question for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item How would you relate this \texttt{[R]} exercise to the analytic proof we did in class?
\end{itemize}
\end{mdframed}

\section{Binary Outcomes}

Now let's pretend that we have an experiment with a binary outcome:

\citet{freedman2008b} famously argues that ``randomization does not justify logistic regression." Let's recall that the coefficient of a logistic regression model is the difference in potential log odds.

<< >>=

experiment <- data.frame(unit = seq(from = 1, to = 8, by = 1))

experiment %<>% mutate(Z = complete_ra(N = length(unit))) %>%
	arrange(desc(Z))

experiment %<>% mutate(y1 = c(1, 1, 0, 1, 1, 0, 0, 1),
		       y0 = c(1, 1, 0, 1, 0, 0, 1, 0),
		       Y = Z * y1 + (1 - Z) * y0)

N <- length(experiment$Z)

Z <- experiment$Z

m <- sum(experiment$Z == 1)

y1 <- experiment$y1

y0 <- experiment$y0

Y <- experiment$Y

true_ate_binary <- mean(y1) - mean(y0)

true_ate_binary

## Log odds
true_logit <- log(mean(y1)/(1 - mean(y1))) - log(mean(y0)/(1 - mean(y0)))

true_logit

est_ate <- mean(Y[Z == 1]) - mean(Y[Z == 0])

est_logit <- log(mean(Y[Z == 1])/(1 - mean(Y[Z == 1]))) -
	log(mean(Y[Z == 0])/(1 - mean(Y[Z == 0])))

all.equal(est_logit, coef(glm(Y ~ Z, family = binomial(link = "logit")))[["Z"]])

new_experiment <- function(z, y0, y1){

	Z = sample(z)

	Y = Z * y1 + (1 - Z) * y0

	est_ate = coef(lm(Y ~ Z))[["Z"]]

	est_logit = coef(glm(Y ~ Z,
                  family = binomial(link = "logit")))[["Z"]]

  return(c(est_ate, est_logit))
}

set.seed(1:5)

randomization_dists <- data.frame(t(replicate(10^3,
                                              new_experiment(z = Z,
                                                             y0 = y0,
                                                             y1 = y1)))) %>%
  rename(est_ates = X1,
         est_logits = X2)

colMeans(randomization_dists)

@

\vspace{5mm}
\begin{mdframed}
\textbf{Questions for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Is the difference-in-means estimator still unbiased (give or take some simulation error) when potential outcomes are binary? Is the coefficient from a logistic regression model unbiased?
\item To which quantities (defined as \texttt{[R]} objects above) are we comparing the means of the two respective randomization distributions in order to assess unbiasedness?
\item How does the randomization distribution above differ from the \textit{null} randomization distribution?
\end{itemize}
\end{mdframed}

\vspace{5mm}
\begin{mdframed}
\textbf{Exercise for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item In a code chunk below, plot both randomization distributions (of the estimated ATEs and estimated logit coefficients) in \texttt{[R]}:
\end{itemize}
\end{mdframed}

<< >>=






@

\section{Application of Unbiased Inference to Matched Designs}

\subsection{Make Your Own Matched Design}

\vspace{5mm}
\begin{mdframed}
\textbf{Exercise for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Please create a matched design with a matched design with the data from the
\citet{cerdaetal2012} paper. We made one of our own, that we call, `fm1,' but
you should make your own.
\end{itemize}
\end{mdframed}

<<>>=

## This our data with the matched sets that you DON'T have access to
load("matched_data.Rdata")

@

<< eval = FALSE >>=

## To load the meddat data and then make your own matched design
load(url("http://jakebowers.org/Matching/meddat.rda"))

@

\subsection{Perform Outcome Analysis After Matching}

Now let's think about how we can draw upon the principles above in order to actually do outcome analysis once we have a matched design.

We want to perform outcome analysis (in this case to estimate average treatment effects) as if we have a \textit{block} randomized experiment in which the matched sets are the experiment's blocks. We therefore want to estimate an average treatment effect \textit{within} each block. But then we need some scheme to weight each block-specific $\widehat{ATE}$ to yield one overall $\widehat{ATE}$. There are three different weighting possibilities below.

\subsubsection{Harmonic Mean Weighting}

Let's first estimate the average treatment effect \textit{within} each block:

<< >>=

blocks <- unique(meddat$fm1)

ATE_hat_by_block <- sapply(blocks, function(x) {
    coef(lm(HomRate08 ~ nhTrt, data = meddat, subset = fm1 == x))[["nhTrt"]]
})

ATE_hat_by_block

@

\vspace{5mm}
\begin{mdframed}
\textbf{Question for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Describe what the output above is giving us . . .
\end{itemize}
\end{mdframed}

Now let's calculate an overall $\widehat{ATE}$ that uses harmonic mean weighting. Harmonic mean weighting weights each block by $\frac{ \frac{2}{n_{st}^{-1} + n_{sc}^{-1}} }{\sum \limits_{s = 1}^{S} \frac{2}{n_{st}^{-1} + n_{sc}^{-1}} }$. 

\vspace{5mm}
\begin{mdframed}
\textbf{Question and Exercise for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Why is there a $2$ in the numerator? 
\item Pick a block with more than $2$ units and calculate its ``effective sample size'' by hand using the appropriate formula---i.e., $\frac{2}{n_{ts}^{-1} + n_{cs}^{-1}}$.
\end{itemize}
\end{mdframed}

Now let's estimate an overall ATE using harmonic mean weights:

<< >>=

num_treated <- meddat %$% {sapply(split(nhTrt, fm1), sum)}

num_control <- meddat %$% {sapply(split(1 - nhTrt, fm1), sum)}

harm_mean <- 2/(1/num_treated + 1/num_control)

obs_overall_ate_hat_hm <- sum(ATE_hat_by_block * (harm_mean/sum(harm_mean)))

@

What are some other weighting schemes one could use? Let's try two others.

\subsubsection{Block Size Weighting}

\vspace{5mm}
\begin{mdframed}
\textbf{Exercise for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Write a mathematical expression for the weights attached to each block-specific $\widehat{ATE}$ in which each block-specific $\widehat{ATE}$ is weighted by the total number of units in each block.
\end{itemize}
\end{mdframed}

Now let's implement this estimation strategy in \texttt{[R]}:

<<>>=

num_units_per_block <- sapply(blocks, function(x){
  length(meddat$nhTrt[meddat$fm1 == x])
})

obs_overall_ate_hat_bs <- sum(ATE_hat_by_block * (num_units_per_block/nrow(meddat)))

@

\subsubsection{ETT Weighting}

Let's try one final weighting scheme (although there are surely many other possibilities.) Let's weight each block-specific $\widehat{ATE}$ by the number of \textit{treated} units within that block. 

That is, the weights associated with each block-specific $\widehat{ATE}$ are: $W_{ETT} = \frac{n_{st}}{\sum \limits_{s = 1}^S n_{st}}$.

<<>>=

num_treated_units_per_block <- sapply(blocks, function(x){
  sum(meddat$nhTrt[meddat$fm1 == x])
})

obs_overall_ate_hat_ett <- sum(ATE_hat_by_block * (num_treated_units_per_block/sum(meddat$nhTrt)))

@

\vspace{5mm}
\begin{mdframed}
\textbf{Questions for Students:}
\vspace{-5mm}
\begin{itemize}\itemsep1pt
\item Now, compare all the overall estimates of the ATE from the $3$ different weighting schemes. Are they all similar?
\item Notice that thus far we have not talked about standard errors, confidence intervals, or hypothesis testing. We have focused on only \textit{estimating} causal effects. How would you describe difference between \textit{estimation} and \textit{testing}?
\end{itemize}
\end{mdframed}

\section{Agenda for Tomorrow}

Tomorrow we will discuss standard errors and \textit{testing} hypotheses about causal effects.

\section{Recommended Readings}

\subsection{Readings on Weighting Schemes}

\subsubsection{Harmonic Mean Weighting}

\begin{itemize}

\item \citet[Section 9.3.3]{hansen2011}

\item \citet[Section 2.1]{hansenbowers2008}

\end{itemize}

\subsubsection{Block-Size Weights}

\begin{itemize}

\item \citet[Section 3.6.1]{gerbergreen2012}

\end{itemize}

\subsubsection{ATT Weighting}

\begin{itemize}

\item
	\url{http://egap.org/methods-guides/10-types-treatment-effect-you-should-know-about}

\end{itemize}

\subsection{Readings on Standard Errors for Tomorrow's Class}

\begin{itemize}

\item \citet[Chapter 3]{gerbergreen2012}

\item \citet[Chapter 6]{dunning2012}

\item \citet[Chapter 33, footnote 11]{fpp2007}

\item \citet{lin2013}

\item \citet{aronowetal2014}

\end{itemize}

\newpage

\bibliographystyle{chicago}
\begin{singlespace}
\bibliography{refs}   % name your BibTeX data base
\end{singlespace}

\newpage

\end{document}