```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
setwd("~/Github/bigdata/lectures/l16/")
load("../l15/jtpa/jtpa.RData")
library(tidyverse)
```

# Today's Class {.smaller}

1. Review
    - RCTs
    - Individual TEs
2. RCTs: Targeting
3. Observational Methods: IV
    - RD
4. Observational Methods: Diff-in-Diff
    - SCM
5. HW6 Review




# RCTs

![](../l15/dag1context.png)

# RCTs

![](../l15/dag2.png)


# Similarity

Typically the ATE is:

$$ \widehat{ATE} = \bar y_1 - \bar y_0 $$

But we could also define it as the average across individual treatment effects. 

$$\widetilde{ATE} = \frac1n \sum_{i=1}^n \widehat{TE}_i$$

With a lot of rewriting of sums -- we can show that when our individual treatment effects use a constant mean to predict counterfactuals:

$$ \widetilde{ATE} = \widehat{ATE}$$

# Proof in Data

```{r}
ybar0 = mean(jtpa$y[jtpa$offer == 0])
ybar1 = mean(jtpa$y[jtpa$offer == 1])
ybar1-ybar0
jtpa_est = jtpa %>% 
    mutate(y1 = y*offer+(1-offer)*ybar1,
           y0 = y*(1-offer)+offer*ybar0) %>%
    mutate(TE = y1-y0)
mean(jtpa_est$TE)
```


But maybe we can improve on the mean as a prediction. 

# ATE from predictions:
```{r,echo=F}
library(ranger)
jtpa_full = jtpa
jtpa = jtpa %>% select(-Index)
jtpa = jtpa[,-c(3,11:13)] #dropping "train" and other details of the treatment. Probably useful, but interpretation is weird.
```
```{r,message=F}
#Build Treat and control dfs
jtpa_cont = jtpa %>% filter(offer == 0)
jtpa_treat = jtpa %>% filter(offer == 1) 
#Estimate treat and control models
mod_treat = ranger(y~.-offer,data=jtpa_treat)
mod_cont = ranger(y~.-offer,data=jtpa_cont)
#Predict counterfactuals for data from other model
jtpa_cont$confact = predict(mod_treat,data = jtpa_cont)$predictions
jtpa_treat$confact  = predict(mod_cont,data = jtpa_treat)$predictions
#Estimate TEs
jtpa_cont$TE = jtpa_cont$confact - jtpa_cont$y
jtpa_treat$TE = jtpa_treat$y - jtpa_treat$confact
#Recombine
jtpa_est = rbind(jtpa_cont,jtpa_treat)
mean(jtpa_est$TE) #ATE
```

# Individual TEs

```{r,echo=F}
ggplot(jtpa_est,aes(x=TE)) +
    geom_histogram(bins=100) +
    geom_vline(xintercept = mean(jtpa_est$TE),col=2,lty=2)
```


# Individual TEs, Do we care?

It looks like some fraction of individuals lost ~\$50k by engaging in this program. Not to mention the program cost to the government. 

- What if we could target the program to people who benefit?
- In other settings, like marketing, we may wish to target groups for whom the expense of advertising is less than the gain in revenue from those individuals.

$\implies$ Targeting. Can we use the RCT data for targeting?


# Targeting

The goal of using targeting is *to identify subpopulations who benefit more*. 

This will always be about *averages* in some subgroup -- we can't predict outliers. Consequentially, there may be excluded individuals who benefit or suffer from exclusion.

- Do not forget about this. 


# Targeting

It is _also_ at times very questionable, concerning, and even illegal to target _some groups_. Particularly if you are advertising financial products, or developing a government program. 

- Importantly, you can wind up targeting (e.g. race) _even_ if you don't observe relevant variables directly
    - Zip code can be a strong proxy. 
    
- I'm going to drop variables for now -- (and rerun prior code)
```{r}
jtpa = jtpa %>% select(-male,-black,-hispanic)
```
```{r,echo=F, message=F}
#Build Treat and control dfs
jtpa_cont = jtpa %>% filter(offer == 0)
jtpa_treat = jtpa %>% filter(offer == 1) 
#Estimate treat and control models
mod_treat = ranger(y~.-offer,data=jtpa_treat)
mod_cont = ranger(y~.-offer,data=jtpa_cont)
#Predict counterfactuals for data from other model
jtpa_cont$confact = predict(mod_treat,data = jtpa_cont)$predictions
jtpa_treat$confact  = predict(mod_cont,data = jtpa_treat)$predictions
#Estimate TEs
jtpa_cont$TE = jtpa_cont$confact - jtpa_cont$y
jtpa_treat$TE = jtpa_treat$y - jtpa_treat$confact
#Recombine
jtpa_est = rbind(jtpa_cont,jtpa_treat)
```

    

# Simple Targeting

Take an observation, predict outcome under treatment and control, take the difference, compare to some threshold. 

# Simple Targeting

```{r}
pred_treat = predict(mod_treat,data=jtpa_est)$predictions
pred_cont  = predict(mod_cont ,data=jtpa_est)$predictions

jtpa_est$te_pred = pred_treat-pred_cont
```

# Simple Targeting

```{r,echo=F}
ggplot(jtpa_est,aes(x=te_pred)) + geom_histogram(bins=100)
```

# Simple Targeting -- Smoothed Conditional Means

```{r,echo=F,message=F}
ggplot(jtpa_est,aes(x=te_pred,y=TE)) + geom_point(alpha=0.05) + geom_smooth()
```

# Simple Targeting -- Bin Scatter

```{r,echo=F,message=F, warning=F}
nbins = 40
binned_jtpa = jtpa_est %>% mutate(bin = cut(te_pred,breaks=nbins)) %>% group_by(offer,bin) %>% select(bin,y,offer,te_pred) %>% summarize(mu = mean(y),pred = mean(te_pred),se=sqrt(var(y)/n())) %>% pivot_wider(names_from = offer,values_from = c(mu,pred,se)) 
mult = 2 #qnorm(1-0.05/nbins) # this is closer to uniform CIs
binned_jtpa = binned_jtpa %>% mutate(ATE = mu_1-mu_0,midpoint = (pred_1+pred_0)/2,se = sqrt(se_1^2+se_0^2)) %>% mutate(ub = ATE+mult*se,lb = ATE-mult*se) %>% select(-mu_0,-mu_1,-pred_0,-pred_1)
ggplot(binned_jtpa,aes(x=midpoint,y=ATE,ymax=ub,ymin=lb))+
    geom_pointrange(alpha=0.5) + 
    xlab("predicted TE") + ylab("Bin's ATE")+
    geom_abline(slope=1,intercept=0,lty=2,col=2)
```

    
# Simple Targeting -- A threshold

Suppose we picked a cutoff of $0. We can look at which _predicted_ treatment effects are below 0, and which _estimated_ treatment effects are below 0, and think about performance.

```{r}
jtpa_est = jtpa_est %>% mutate(predte0 = te_pred > 0,
                               TE0 = TE > 0)
table(jtpa_est$TE0,jtpa_est$predte0,deparse.level=2)
```


# Simple Targeting  - Quick test 
```{r}
chisq.test(jtpa_est$TE0,jtpa_est$predte0)
```

So these aren't unrelated.


# Simple Targeting -- A threshold

That felt like our classification work. Could we use an ROC curve?

It will take some modification. For starters, our thresholds aren't now restricted to 0,1.

But more importantly, we aren't observing the true outcome of the classification -- even that is a prediction.

- Easier to look at the mean ATE in the two groups that our classification threshold produces. 
- Also important to use OOS predictions. (One step at a time)


# Simple Targeting - Measuring Performance
Continuing with our Threshold of 0 for a moment. Lets calculate mean ATE for each group our predictions and threshold create.
```{r, message = F}
jtpa_TEs = jtpa_est %>% 
    group_by(offer,predte0) %>% 
    summarize(y=mean(y)) %>%
    pivot_wider(names_from = predte0, values_from = y)
jtpa_TEs
```

# Simple Targeting - Measuring Performance

Find the actual "ATE" for each subpopulation:
```{r}
jtpa_TEs[2,2:3]-jtpa_TEs[1,2:3]
```

So for the group the proposed policy would block, the ATE was substantially negative, and for the group it would try to target, the ATE was much more positive than for the entire population.

- Substantial Improvement in our performance. *IF* it holds out of sample. 

# Simple Targeting

But first -- we should consider other thresholds...

Going to plot out possible other thresholds, and the "gain" in performance from using those thresholds. 

"Net Gain" is going to be "benefit-(avoided loss)". For an indiscriminate policy (all eligible), this is the pure ATE. 

For a policy targeting some group and blocking some group (of equal size), this is "ATE(targeted group)-ATE(blocked group)".

- Notice -- if we block a group which would benefit -- this is a loss here. (We are assuming no cost of the policy -- big loss assumptions baked in)

# Simple Targeting - ROC

Built a function, critical elements below
```{r}

threshold_ATES = function(threshold,y,predte,
                          treatment,MC=0,FC=0) {
  n = length(predte)
  targets = predte > threshold
  ATE_treat   = mean(y[treatment == 1 & targets == 1])-
                mean(y[treatment == 0 & targets == 1])
  ATE_untreat = mean(y[treatment == 1 & targets == 0])-
                mean(y[treatment == 0 & targets == 0])
  ntreat = sum(targets)
  nuntreat = n-ntreat
  gain = ntreat*ATE_treat-nuntreat*ATE_untreat
  gain = gain-MC*ntreat-FC 
  gain = gain/n
  c(threshold,ATE_treat,ATE_untreat,ntreat,nuntreat,gain)
}

```


# Gain

```{r,echo=F,warning=F}
nthresholds = 80
thresholds = seq(min(jtpa_est$te_pred),max(jtpa_est$te_pred),length.out=nthresholds)
gains = t(sapply(thresholds,threshold_ATES,
               y=jtpa_est$y,predte=jtpa_est$te_pred,
               treatment=jtpa_est$offer,MC=0,FC=0))

colnames(gains) = c("threshold","ATE_treat","ATE_untreat","ntreat","nuntreat","gain")
gains = as_tibble(gains)
best_thresh = gains$threshold[which.max(gains$gain)]
ggplot(gains,aes(x=threshold,y=gain)) + geom_line() + geom_vline(xintercept = best_thresh,col=2,lty=2)
print(best_thresh)
```

# Costs

Suppose we had some cost for each person treated (call it \$300). And a fixed cost ($1000).

Where is the cutoff?

Changes our "gains" function. 

# New plot


```{r,echo=F, warning=F}
gains = sapply(thresholds,threshold_ATES,
               y=jtpa_est$y,predte=jtpa_est$te_pred,
               treatment=jtpa_est$offer,MC=300,FC=1000)
gains = t(gains) #transpose so the columns and rows are correct

colnames(gains) = c("threshold","ATE_treat","ATE_untreat","ntreat","nuntreat","gain")
gains = as_tibble(gains)
best_thresh = gains$threshold[which.max(gains$gain)]
ggplot(gains,aes(x=threshold,y=gain)) + geom_line() + geom_vline(xintercept = best_thresh,col=2,lty=2)
print(best_thresh)
```

# Interpretation

- First gains plot was "average benefits _to individuals_ of targeting over full treatment" as a function of thresholds
- Second plot was "average benefits to individuals _above costs_ of targeting over full treatment" as a function of thresholds

But maybe we aren't currently implementing anything?

- Then it would be inappropriate to subtract off the ATE for the untreated group. 

# New plot

```{r,echo=F, warning=F}
gains$gain = gains$ATE_treat*gains$ntreat-gains$ntreat*300-1000
best_thresh = gains$threshold[which.max(gains$gain)]
ggplot(gains,aes(x=threshold,y=gain)) + geom_line() + geom_vline(xintercept = best_thresh,col=2,lty=2)
print(best_thresh)
```

# Plot 2

```{r,echo=F}
gains2 = gains %>% select(threshold,ATE_treat,ntreat) %>% pivot_longer(cols = c(ATE_treat,ntreat),names_to = "var",values_to="val")
ggplot(gains2,aes(x=threshold,y=val))+
    geom_line() + facet_grid(var~.,scales="free")
```

# OOS Targeting measures.

**USE CROSS VALIDATION.**

# OOS Targeting

```{r}
OOS_tepreds = function(holdout_ind,data,formula = y~.-offer) {
    holdout = data[holdout_ind,]
    train = data[-holdout_ind,]
    jtpa_cont = train %>% filter(offer == 0)
    jtpa_treat = train %>% filter(offer == 1) 
    #Estimate treat and control models
    mod_treat = ranger(y~.-offer,data=jtpa_treat)
    mod_cont = ranger(y~.-offer,data=jtpa_cont)
    pred_treat = predict(mod_treat,data=holdout)$predictions
    pred_cont  = predict(mod_cont ,data=holdout)$predictions
    te_pred = pred_treat-pred_cont
    cbind(holdout_ind,te_pred)
}
```

# Kfold CV
```{r,echo=F}
set.seed(77261) #34422
```

```{r}
k = 20
fold_ids = sample(rep(1:k,length.out = nrow(jtpa)))
hold_indices = lapply(1:k,function(foldk) which(fold_ids == foldk))
oos_preds = lapply(hold_indices,OOS_tepreds,data=jtpa)
jtpa_est = jtpa
jtpa_est$oos_tepreds = NA
for (i in 1:k){
    jtpa_est$oos_tepreds[oos_preds[[i]][,1]] = oos_preds[[i]][,2]
}
```


# Plots

```{r,echo=F,warning=F}
gains = t(sapply(thresholds,threshold_ATES,
               y=jtpa_est$y,predte=jtpa_est$oos_tepreds,
               treatment=jtpa_est$offer,MC=0,FC=0))
colnames(gains) = c("threshold","ATE_treat","ATE_untreat","ntreat","nuntreat","gain")
gains = as_tibble(gains)
best_thresh = gains$threshold[which.max(gains$gain)]
ggplot(gains,aes(x=threshold,y=gain)) + geom_line() + geom_vline(xintercept = best_thresh,col=2,lty=2)
print(best_thresh)
```

# Targeting

- Powerful tool in many domains
    - Optimal Policy Design
    - Marketing
    - Medical treatment choice
    - Many more

- Care and Caution needed
    - Discrimination issues
    - OOS/overfitting issues
        - We may think we are finding good targeting things -- and really be failing
        - Our standard tools for examining and improving predictions will help here.
        
# RCT Wrapup

RCTs can have problems arise in several ways. 

- Placebo effects: I know I recieved fancy treatment
- Spillover effects: My friend got tutoring and I benefit from that
- and others.
    
They can be great -- but they aren't flawless.

# Observational Causal Inference

Frequently we care about causal issues in contexts where we can't run an experiment.

There are still tools for trying to assess a causal effect. They will rely on _something_ being **"as if"** random. 

# Instrumental Variables

The simplest example here is "instrumental variables". 

Recall, when we use an RCT, we assume some random number generator influences X without having any way to effect Y, and without being affected by Y. 

- Thus we can use the RNG output to predict Y _only_ because of the causal effect of X

# RCT DAG

![](../l15/dag2.png)

# Instrumental Variables

IV takes the same idea. What if I had some other variable -- not a random number -- which affects X, but can't affect Y or be affected by Y. Then I could predict Y as a function of that random variable -- and any relationship would be only because of X.

# IV DAG

![](dag_iv.png)


# IV - Example

Many Examples out there.

One classic example is testing the relationship between [lead and crime](https://www.sciencedirect.com/science/article/abs/pii/S0014498316300109) (or other things like test scores)

Essentially, we can look at a city's proximity to lead foundries to predict their usage of lead pipes. Then we can try to use the proximity to lead foundries to predict crime. 

"In the absence of other plausible relationships between lead foundry proximity and crime, any effect we find comes from the lead".

# IV - Estimation

Simplest method ("2SLS") runs a linear regression between Z and X, getting $\beta_X$, then another regression between Z and Y, getting $\beta_Y$. 

The estimated effect of X on Y is $\beta_Y/\beta_X$. 

# IV - Problems

- Weak IV: If $\beta_X$ is near 0, our estimates for the effect head towards $\infty$ or $-\infty$. This is a bigger problem than it may seem. 
    - Z needs to be relevant to X
- Exclusion: If Z has other channels for affecting Y, then we can't interpret this as the causal effect of X on Y.
    - E.g. if lead foundries disproportionately went out of business, causing unemployed workers to form mobs that committed crimes. 

# IV - Other Uses

Implicitly, that 2SLS procedure of dividing one coefficient by another is how we estimated "average treatment effect on the treated" on Tuesday. 

> Divide the ATE by the fraction recieving treatment


# Regression Discontinuity

Regression discontinuity is similar to IV.

It takes advantage of some policy threshold, on one side of which individuals receive treatment, and on one side of which they don't.

# RD Example

Classic examples involve looking at students.

- E.g. Some schools have a GPA cutoff for being put on academic probation. We may want to know the effect of 'academic probation' on    
    - dropout rates
    - future GPA
    - Etc

# RD

![](d0.png)

# RD

![](d6.png)

# RD - Estimation

Typically you use nonparametric tools (local linear regressions above), with carefully chosen bandwidths (choose with CV if you like).

That will let you estimate $E[Y|X=threshold,T=1]$ and for $T=0$ (i.e. the intercept on each side). Then we take the difference.

# RD - Notes

RD gives you a Conditional Average Treatment Effect -- the RD estimate of the treatment effect is the estimated _average_ treatment effect _for individuals at the threshold_. 

# RD - Problems

RD assumes individuals are fundamentally similar on either side of the threshold.

- But if individuals are aware of the threshold, and can control their position relative to it, then they may not be fundamentally similar
    - Individuals on right side may have chosen to be there
    - And on the left side, not chosen to be on the right side
    - Selection problems


# Difference-in-Differences

Suppose we observe a state implement some policy.

- E.g. California imposes a large tax (\$0.25) on cigarettes. 

And we want to know the effect of the cigarette tax on consumption. 

What could we do?

# Difference-in-Differences

Recall the cigarette data. Per capita cigarette consumption in 38 states over 30 years. California policy implemented in 1989. 

```{r, message=F}
df = read_csv("https://codowd.com/bigdata/lectures/l13/cigarettes.csv")
df
```

# Quick Plot

```{r,echo=F}
df = df %>% pivot_longer(!year,names_to = "state",values_to = "consumption")
```
```{r,echo=F}
df = df %>% mutate(treated = state == "CA")
ggplot(df,aes(x=year,y=consumption,group=state)) + geom_line() +
    geom_line(data=df%>%filter(state=="CA"),col=2,lwd=2)+
    geom_vline(xintercept=1989,lty=2,col="blue")
```

# Idea: State Differences

We could look at the difference between California's consumption and consumption in other states after the policy takes effect. 

- Suppose California was randomly selected for the policy -- this is how we would run an RCT analysis. 
- But suppose California already consumed a different amount -- we would be biased by the scale of pre-existing differences.

# Idea: Time Differences

We could look at the difference between California's consumption before and after the policy takes effect.

- This would avoid concerns about differences between states. 
- But if there was some kind of time trend, we would be adding a bias from the time trend.

# Quick Plot

```{r,echo=F}
ggplot(df,aes(x=year,y=consumption,group=state)) + geom_line() +
    geom_line(data=df%>%filter(state=="CA"),col=2,lwd=2)+
    geom_vline(xintercept=1989,lty=2,col="blue")
```

# Difference-in-Difference

Diff-in-diff says: "why not both?"

- We can look at the difference between California and other states, and see how it changes over time. 
    - Or we could look at each state's changes over time, and see how California compares to other states.
        - These are equivalent. 
        
# Diff-in-Diff: Year means in treatment and control groups

```{r, message=F,echo=F }
dfdd = df %>% group_by(year,treated) %>% summarize(consumption = mean(consumption))
ggplot(dfdd,aes(x=year,y=consumption,col=treated))+geom_line()+
    geom_vline(xintercept=1989,col="blue",lty=2)
```

# Diff-in-Diff: Difference between year means

```{r,echo=F}
dfdd %>% pivot_wider(names_from = treated,values_from = consumption) %>% mutate(Difference = `TRUE`-`FALSE`) %>% ggplot(aes(x=year,y=Difference))+geom_line()+geom_vline(xintercept=1989,col="blue",lty=2)
```

# Diff-in-Diff

The difference looks to have gotten bigger - so there was probably an effect. The tax probably reduced consumption.

# Diff-in-Diff: Notes

- Relies on the 'parallel trends' assumption.
    - Both treatment and control groups have time trends that are similar.
        - In that case, we can estimate the treatment effect.
- Usually relies on naive means in each of the treatment and control groups.
    - I.e. equal weighted means
    
# Synthetic Controls

Synthetic controls builds on Diff-in-Diff by saying "can we improve those weights?"

- This should sound familiar.

We may be able to find weights for the control group that make it better at predicting the treatment group.

- This will let us relax the 'parallel trends' assumption somewhat
- It will also improve the statistical efficiency. 

# Synthetic Controls

I'm not going to go through code here. SCM is mostly beyond the scope of this class. 

- But its useful for you _because the basic idea is familiar_
    - "Can we improve those weights for out of sample prediction"

# Observational Causal Inference Wrapup

RCTs, IV, RD all rely on "exogenous variation" coming from something. Variation which is unrelated to the outcome of interest.

- RCTs: Varation from a random number generator
- IV: Variation from some instrument
- RD: Variation from randomness around some threshold

Diff-in-diff and SCM rely on finding the two sources of bias, and trying to eliminate them.

- Bias 1: our treated unit is fundamentally different
- Bias 2: out treated time periods are fundamentally different

These are different approaches with different strengths and weaknesses. 

# Causal Inference Wrapup

There is no magic bullet.

We have questions that desperately need answers. 

We have data. That data _may_ be able to help. But it may only mislead us. 

- Critical to make decisions and implement policies which _fail gracefully_
    - We want things that work wonderfully if we're right about the world
    - But which aren't tremendously destructive if we're wrong.
    
# HW 6 Review

If we have time I'll review HW6 in a minute.

# Wrap up

Homework 7 will be posted this afternoon, due next Wednesday. It is about causal inference.  

Prediction competition 3 submissions must be in by midnight tomorrow. 

Next Tuesday we will look Neural Nets, and next Thursday I'll try to review the whole of this course in 90 minutes. 

See you next week!









