day19-Interference.Rmd

---
title: Causal Inference for Models of Interference
date: '`r format(Sys.Date(), "%B %d, %Y")`'
author: ICPSR 2023 Session 1
bibliography:
 - BIB/refs.bib
 - BIB/master.bib
 - BIB/misc.bib
fontsize: 10pt
geometry: margin=1in
graphics: yes
biblio-style: authoryear-comp
output:
  beamer_presentation:
    slide_level: 2
    keep_tex: true
    latex_engine: xelatex
    citation_package: biblatex
    template: icpsr.beamer
    incremental: true
    includes:
        in_header:
           - defs-all.sty
---


<!-- Make this document using library(rmarkdown); render("day12.Rmd") -->


```{r include=FALSE, cache=FALSE}
# Some customization.  You can alter or delete as desired (if you know what you are doing).
# knitr settings to control how R chunks work.
rm(list=ls())

require(knitr)

## This plus size="\\scriptsize" from https://stackoverflow.com/questions/26372138/beamer-presentation-rstudio-change-font-size-for-chunk

knitr::knit_hooks$set(mysize = function(before, options, envir) {
			           if (before){
				      return(options$size)
				   } else { 
					   return("\\normalsize")
				   }
})

knit_hooks$set(plotdefault = function(before, options, envir) {
    if (before) par(mar = c(3, 3, .1, .1),oma=rep(0,4),mgp=c(1.5,.5,0))
})

opts_chunk$set(
  tidy=FALSE,     # display code as typed
  echo=TRUE,
  results='markup',
  strip.white=TRUE,
  fig.path='figs/fig',
  cache=FALSE,
  highlight=TRUE,
  width.cutoff=132,
  size='\\scriptsize',
  out.width='.8\\textwidth',
  fig.retina=FALSE,
  message=FALSE,
  comment=NA,
  mysize=TRUE,
  plotdefault=TRUE)

if(!file.exists('figs')) dir.create('figs')

options(digits=4,
	scipen=8,
	width=132
	)
```

```{r eval=FALSE, include=FALSE, echo=FALSE}
## Run this only once and then not again until we want a new version from github
library('devtools')
library('withr')
with_libpaths('./lib', install_github("markmfredrickson/RItools"), 'pre')
```

```{r echo=FALSE}
library(dplyr)
library(RItools,lib.loc="./lib")
library(coin)
library(dplyr)
library(ggplot2)
library(sandwich)
library(lmtest)
library(estimatr)
library(DeclareDesign)
```

## Today

\begin{enumerate}
\item Agenda: Causal inference when you have Interference.
\item Readings about this for beyond today: \textcite{bowersetal2013}, \textcite{ichino2012deterring}, \textcite{aronow2013estimating}, \cite{liu2014large}
\item Questions arising from the reading or assignments or life?
\item Recap: Sensitivity analysis two ways.
\item Questions and discussion about what should be next for you?
\end{enumerate}

# Remember about approaches to statistical inference for causal quantities

## Design Based Approach 1: Test Models of Potential Outcomes

 1. Make a guess (or model of) about $\tau_i$.
 2. Measure consistency of data with this model given the design.

\centering
  \includegraphics[width=.7\textwidth]{images/cartoonFisherNew1.pdf}

## Design Based Approach 1: Test Models of Potential Outcomes

\centering
  \includegraphics[width=\textwidth]{images/cartoonFisherNew1.pdf}

## Design Based Approach 1: Test Models of Potential Outcomes

\centering
  \includegraphics[width=\textwidth]{images/cartoonFisherNew2.pdf}

## Design Based Approach 1: Test Models of Potential Outcomes

  \centering
  \includegraphics[width=.9\textwidth]{images/cartoonFisher.pdf}

##  Design Based Approach 1: Test Models of Potential Outcomes

```{r}
smdat <- data.frame(Z=c(0,1,0,1),Y=c(16,22,7,4000))

tz_mean_diff <- function(z,y){
	mean(y[z==1]) - mean(y[z==0])
}

tz_mean_rank_diff <- function(z,y){
	ry <- rank(y)
	mean(ry[z==1]) - mean(ry[z==0])
}

newexp <- function(z){
	sample(z)
}
```

##  Design Based Approach 1: Test Models of Potential Outcomes


```{r repexp, cache=TRUE}
set.seed(12345)
rand_dist_md <- replicate(1000,tz_mean_diff(z=newexp(smdat$Z),y=smdat$Y))
rand_dist_rank_md <- replicate(1000,tz_mean_rank_diff(z=newexp(smdat$Z),y=smdat$Y))

```

```{r calc}
obs_md <- tz_mean_diff(z=smdat$Z,y=smdat$Y)
obs_rank_md <- tz_mean_rank_diff(z=smdat$Z,y=smdat$Y)
c(obs_md,obs_rank_md)
table(rand_dist_md)
table(rand_dist_rank_md)
p_md <- mean(rand_dist_md >= obs_md)
p_rank_md <- mean(rand_dist_rank_md >= obs_rank_md)
c(p_md, p_rank_md)
```

##  Design Based Approach 1: Test Models of Potential Outcomes

```{r}
smdat$zF <- factor(smdat$Z)
smdat$rY <- rank(smdat$Y)

md_test_exact <- oneway_test(Y~zF,data=smdat,distribution=exact(),alternative="less")
md_test_sim <- oneway_test(Y~zF,data=smdat,distribution=approximate(nresample=1000),
			   alternative="less")
md_test_asymp<- oneway_test(Y~zF,data=smdat,distribution=asymptotic(),alternative="less")

rank_md_test_exact <- oneway_test(rY~zF,data=smdat,distribution=exact(),alternative="less")
rank_md_test_asymp<- oneway_test(rY~zF,data=smdat,distribution=asymptotic(),alternative="less")
rank_md_test_sim <- oneway_test(rY~zF,data=smdat,distribution=approximate(nresample=1000),
			   alternative="less")
```

##  Design Based Approach 1: Test Models of Potential Outcomes

```{r}
pvalue(md_test_exact)
pvalue(md_test_sim)
pvalue(md_test_asymp)
pvalue(rank_md_test_exact)
pvalue(rank_md_test_sim)
pvalue(rank_md_test_asymp)
```
##  Design Based Approach 1: Test Models of Potential Outcomes

```{r}
s_rank_md_exact <- support(rank_md_test_exact)
d_rank_md_exact <- dperm(rank_md_test_exact,s_rank_md_exact)
d_rank_md_asymp <- dperm(rank_md_test_asymp,seq(-1.6,1.6,.1))

s_md_exact <- support(md_test_exact)
d_md_exact <- dperm(md_test_exact,s_md_exact)
d_md_asymp <- dperm(md_test_asymp,seq(-1.6,1.6,.1))
```

```{r plotexactdists, echo=FALSE}
par(mfrow=c(1,2))
plot(s_rank_md_exact,d_rank_md_exact,type="h",xlim=c(-1.6,1.6),ylim=c(0,.4),lwd=3)
lines(seq(-1.6,1.6,.1),d_rank_md_asymp)
plot(s_md_exact,d_md_exact,type="h",xlim=c(-1.6,1.6),ylim=c(0,.4),lwd=3)
lines(seq(-1.6,1.6,.1),d_md_asymp)
```

## Design Based Approach 2: Estimate Averages of Potential Outcomes

  1. Notice that the observed $Y_i$ are a sample from  the (small, finite) population of $(y_{i,1},y_{i,0})$.
  2. Decide to focus on the average, $\bar{\tau}$, because sample averages, $\hat{\bar{\tau}}$ are unbiased and consistent estimators of population averages.
  3. Estimate $\bar{\tau}$ with the observed difference in means.

\centering
  \includegraphics[width=.5\textwidth]{images/cartoonNeyman.pdf}

## Design Based Approach 2: Estimate Averages of Potential Outcomes

\centering
  \includegraphics[width=.9\textwidth]{images/cartoonNeyman.pdf}


## Design Based Approach 2: Estimate Averages of Potential Outcomes

```{r}

est1 <- difference_in_means(Y~Z,data=smdat)
est1
```

## Design Based Approach 2: Estimate Averages of Potential Outcomes

Is this estimator a good guess about the unobserved causal effect? Let's simulate to learn.

First, make up a simulated research design and underlying potential outcomes:

```{r setupdd}

set.seed(12345)
smdat$x <- runif(4,min=0,max=1)
smdat$e0 <- round(runif(4,min=min(smdat$Y[smdat$Z==0]),max=max(smdat$Y[smdat$Z==0])))
smdat$y0 <- with(smdat,3*sd(e0)*x + e0)
smdat$y1 <- smdat$Y - smdat$y0

thepop <- declare_population(smdat[,c("y0","y1","x")])
theassign <- declare_assignment(m=2)
po_fun <- function(data){
	data$Y_Z_1   <- data$y0
	data$Y_Z_0 <- data$y1
	data
}
thepo <- declare_potential_outcomes(handler=po_fun)
thereveal <- declare_reveal(Y,Z) ## how does assignment reveal potential outcomes
thedesign <- thepop + theassign + thepo + thereveal

oneexp <- draw_data(thedesign)
```

## Design Based Approach 2: Estimate Averages of Potential Outcomes

Next specify our estimator and compare to a couple of others:

```{r}
theestimand <- declare_estimand(ATE = mean(Y_Z_1 - Y_Z_0))
theest1 <- declare_estimator(Y~Z, estimand=theestimand, model=difference_in_means,
			     label="est1: diff means1")
theest2 <- declare_estimator(Y~Z, estimand=theestimand, model=lm_robust,
			     label="est2: diff means2")
theest3 <- declare_estimator(Y~Z+x, estimand=theestimand, model=lm_robust,
			     label="est3: covadj diff means")

thedesign_plus_est <- thedesign + theestimand + theest1 + theest2 + theest3
theest1(oneexp)
theest2(oneexp)
theest3(oneexp)

```

## Design Based Approach 2: Estimate Averages of Potential Outcomes

Finally, repeat the design, applying our estimators and comparing to our known potential outcomes.

```{r diagnose1, warnings=FALSE, cache=TRUE}
set.seed(12345)

thediagnosands <- declare_diagnosands(
     bias = mean(estimate - estimand),
     rmse = sqrt(mean((estimate - estimand) ^ 2)),
     power = mean(p.value < .25),
     coverage = mean(estimand <= conf.high & estimand >= conf.low),
     mean_estimate = mean(estimate),
     sd_estimate = sd(estimate),
     mean_se = mean(std.error),
     mean_estimand = mean(estimand)
     )

diagnosis <- diagnose_design(thedesign_plus_est,sims=1000,bootstrap_sims=0,
			     diagnosands = thediagnosands)
```


## Design Based Approach 2: Estimate Averages of Potential Outcomes

```{r}
kable(reshape_diagnosis(diagnosis)[,-c(1:2,4)])
```

```{r simmethod, eval=FALSE, echo=FALSE}
thedesign_sims<- simulate_design(thedesign_plus_est,sims=1000)
res <- thedesign_sims %>% group_by(estimator_label) %>% summarize(bias=mean(estimate-estimand))
```



## Model Based Approach 1: Predict Potential Outcomes}

  \smallskip
  \centering
  \includegraphics[width=.9\textwidth]{images/cartoonBayes.pdf}

# What about interference? How define, estimate and test causal effects.

## Statistical inference with interference?

  \centering
  \includegraphics[width=.3\textwidth]{complete-graph.pdf}

  \only<1>{

    \includegraphics[width=.95\textwidth]{interference-example.pdf}

    On estimation see (Sobel, Aronow, Samii, Hudgens, Ogburn,
    VanderWeele, Toulis, Kao, Coppock, Sicar, Raudenbush, Hong, \ldots). What is the function of potential
    outcomes that we can estimate using observed data? WWFisherD? WWNeymanD?
  }

  \only<2>{
    \includegraphics[width=.95\textwidth]{interference-example-2.pdf}

    Introducing the \textbf{uniformity trial} $\equiv \by_{i,0000}$
    (Rosenbaum, 2007).
  }



# A Simple Design-based Estimation Approach

## Voter Registration in Ghana 2008

\begin{columns}
  \begin{column}{.4\textwidth}
    \centering
    \includegraphics[height=.75\textheight]{Ghana2008BigNetwork.png}
  \end{column}
  \begin{column}{.6\textwidth}
    \begin{block}{}
      \begin{itemize}
        \item Presidential and parliamentary elections in December 2008.
        \item 13 day voter registration exercise in August 2008.
        \item Estimated 800,000 people newly eligible to vote, but \textbf{2 million} new voters registered.
        \item Term-limited president, election expected to be very close.  Decided by less than 50,000 votes out of more than 9 million votes cast.
      \end{itemize}
    \end{block}

  \end{column}
\end{columns}

## Voter Registration in Ghana 2008

\begin{columns}
  \begin{column}{.4\textwidth}
    \centering
    \only<1>{\includegraphics[width=\textwidth]{Ghana2008NetworkSlice-Annotated.png}}
    \only<2>{\includegraphics[width=1.8\textwidth]{Ghana2008NetworkSlice-Annotated.png}}
    \only<3>{\includegraphics[width=1.8\textwidth]{ichinoMap.pdf}}
  \end{column}
  \begin{column}{.5\textwidth}
    \begin{block}{}
      \begin{itemize}
        \item Coalition of Domestic Election Observers (CODEO) organize registration observers.  Registration day was generally \textbf{not} routinely monitored.
        \item Design: 4 regions (non-random); within-region, 13 blocks by 2004 parliamentary results; 1 of 3 constituencies in each block receives observers (random).
        \item Randomly assign observers to approximately 25\% of election
          polling stations (ELAs) in selected constituency (77 of 868).
        \item Party agents seen approaching treated ELAs in buses, and then
          driving away toward control ELAs.
      \end{itemize}
    \end{block}
  \end{column}
\end{columns}


## Assessing the sharp null hypothesis of no effects.

\begin{figure}[th]
  \centering
  \includegraphics[width=.8\textwidth]{outcome-desc.pdf}
  %  \caption{Left:  The number of registered voters was lower in 2008 than in 2004 in treated ELAs (centered within experimental block).  Right: QQ-plot of change in the number of registered voters in control and treated ELAs. The distributions differ largely at the tails.}
  \label{fig:observedoutcome}
\end{figure}

Q: What is the probability of seeing
as large an observed difference between the treated and control groups, if the
observers had no effect at all --- recalling that no effect means no interference as well as no other effect?

A:$p=0.018$ (using a mean-difference test-statistic).


## Approaches for going beyond the sharp-null of no effects


**Estimation:**

   - Use design to isolate units
   - Or weight average differences by model of propagation / spillover \parencite{aronow2013estimating,toulis2013estimation}

**Testing:**

   - Assess implications of models of network-propagation effects \parencite{bowersetal2013}.
   - Invert hypothesis tests comparing levels/ranks of treatment outcomes to the uniformity trial \parencite{rosenbaum2007a}.


## Estimation restricting interference by design

Imagine that $Z_i \in \{U,C,T\}$ where  $T$ is treatment
(election observers), $C$ is control with possible spillover and and $U$ is
"uniformity trial" or control with no possible spillover. 

\medskip


Thus, if you have isolated units and randomization (such that *all units have
positive probability of $Z_i \in \{U,C,T\}$*) we have $y_{i,T}$, $y_{i,C}$, and
$y_{i,U}$ for each unit.^[The two-level design
\parencite{sinclair2012detecting}. See also \textcite[Chap 8]{gerbergreen2012}
or generalized saturation design \parencite{baird2014designing}.
\textcite{liu2014large} for some nice theory.]

\pause
\medskip

And you can define and estimate $\bar{\tau}_{\text{spillover}}=\bar{y}_{C} -
\bar{y}_{U}$ or $\bar{\tau}_{\text{Direct Effect}} = \bar{y}_{T} - \bar{y}_U$
etc..

## Estimation restricting interference by design

\igrphx{ichinoMap.pdf}

## Estimation restricting interference by design

How would we use this data to estimate direct, indirect, or spillover effects?

```{r loadghanastuff, echo=FALSE, cache=FALSE}
library(sp)
library(maptools)
load("~/Documents/PROJECTS/GhanaModels/data/ELAs.rda")
load("~/Documents/PROJECTS/GhanaModels/data/network.rda")
ELAs.df <- as.data.frame(ELAs)
outcome<-ELAs$outcome
treatment<-ELAs$tela
options(scipen=10)
```

```{r}
table(Z=ELAs.df$tela,TrtRegion=ELAs.df$NSF_Const_registration_Treat)
tmpdat <- group_by(ELAs.df,block,ZLv2=tela,ZLv1=NSF_Const_registration_Treat) %>% 
  summarise(barYb=round(mean(reg2008ELA - reg2004ELA),5),
            nb=n(), nTb=sum(tela),
            barY08=mean(reg2008ELA),
            barY04=mean(reg2004ELA))
tmpdat
```

# The Model of Effects Approach

[See the ColumbiaStats Talk]


## Summary

  - The sharp null implies no interference. So no need for an *assumption* of no interference.
  - Careful design can allow estimation when the number of potential outcomes is small-ish per unit.
  - Models of effects can specify flexible theoretical models of propagation over networks and randomization can justify statistical inference and causal inference.

## References