Presentation_English.Rmd

---
title: "The effect of Tilburg gold on your health"
subtitle: "Presentation Zorginstituut: Combining causal inference and machine learning in practice"
author: "Misja Mikkers & Gertjan Verhoeven"
institute: "Dutch Healthcare Authority (NZa) & Tilburg University"
date: "January 09, 2020"
output:
   beamer_presentation:
    theme: "Goettingen"
    colortheme: "rose"
    fonttheme: "structurebold"
    includes:
          in_header: header.tex
---

```{r setup II, include=FALSE, eval=TRUE}
knitr::opts_chunk$set(eval = TRUE, echo = FALSE, warning = FALSE, message = FALSE)
knitr::opts_chunk$set(fig.width= 6, fig.height = 4) 
```


```{r}
# Packages

library(tidyverse)
library(dagitty)
library(caret)
library(knitr)
library(ranger)
library(stargazer)
library(grf)

```

# Introduction


## Introduction

![](DAGbeach.JPG)\




## Causal inference as missing data

```{r}
X_1 <- c("Male", "Female", "Female", "...")

PO <- as.data.frame(X_1) %>%
  mutate(X_2 =  c("9", "60", "7", "...")) %>%
  mutate(X_3 = c("14", "36", "2", "...")) %>%
  mutate(X_i =  c("1", "0", "1", "...")) %>%
  mutate(I = c("0", "1", "1", "...")) %>%
  mutate(Y0 = c("67", "NA", "NA", "...")) %>%
  mutate(Y1 = c("NA", "113", "54", "..."))



kable(PO)

```



# Causality


## Prediction versus understanding


```{r fig.width= 6, fig.height = 4}
set.seed(123)

N <- 1000

IQ_groep <- rep(c(1:10), times = N/10)

conf <- as.data.frame(IQ_groep) %>%
  mutate(IQ = rnorm(n= N, mean = 75 + 5 * IQ_groep, sd = 2)) %>%
  group_by(IQ_groep) %>%
  mutate(IQ_mean = mean(IQ)) %>%
  ungroup() %>%
  mutate(Aantal_boeken = IQ_mean/3 + rnorm(n= N,  mean = 0), sd = 6) %>%
  mutate(C = 5 * IQ_mean + rnorm(n = N, mean = 0, sd = 10)) %>% # ruwe scores voor Cito
  mutate(Cito = (C - min(C))/max((C - min(C)))*50 + 500) %>% # verander ruwe scores in scores in de relevante range
  mutate(IQ_groep = as.factor(IQ_groep)) %>%
  select(-C)

ggplot(data = conf, aes(x = Aantal_boeken, y = Cito)) +
  geom_point() +
  geom_smooth(method = "lm") + 
  xlab("Number of books possessed by parents") +
  ylab("Test scores of their children") +
  theme_bw()


```





## Prediction versus understanding


```{r}


ggplot(data = conf, aes(x = Aantal_boeken, y = Cito, color = IQ_groep)) +
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("Number of books possessed by parents") +
  ylab("Test scores of their children") +
  labs(color = "IQ Group\n parents") +
  theme_bw()

```





# DAGS



## Causal graph (DAG) building blocks



  

```{r}

Label <- c("X", "Z", "Y", "X", "Z", "Y", "X", "Z", "Y") # Ik gebruik dit niet: geom_text en ggeom_text_repel plaatsen labels niet mooi

DAG <- as.data.frame(Label) %>%
  mutate(X = c(1, 2, 3, 5, 6, 7, 9, 10, 11)) %>%
  mutate(Y = c(1, 3, 1, 1, 3, 1, 1, 3, 1))

ggplot(data = DAG, aes(x = X, y = Y)) +
  geom_point(size = 3, color = "black") +
  geom_segment(x = 2, xend = 1, y = 2.9, yend = 1.1, arrow = arrow(), size = 1, color = "blue") +
  geom_segment(x = 2, xend = 3, y = 2.9, yend = 1.1, arrow = arrow(), size = 1, color = "blue") +
  geom_segment(x = 5, xend = 5.9, y = 1.1, yend = 2.9, arrow = arrow(), size = 1, color = "blue") +
  geom_segment(x = 7, xend = 6.1, y = 1.1, yend = 2.9, arrow = arrow(), size = 1, color = "blue") +
  geom_segment(x = 9, xend = 9.9, y = 1.1, yend = 2.9, arrow = arrow(), size = 1, color = "blue") +
  geom_segment(x = 10.1, xend = 11, y = 2.9, yend = 1.1, arrow = arrow(), size = 1, color = "blue") +
  annotate("text", x = 1, y = 0.8, label = "X") +
  annotate("text", x = 2, y = 3.2, label = "Z") +
  annotate("text", x = 3, y = 0.8, label = "Y") +
  annotate("text", x = 5, y = 0.8, label = "X") +
  annotate("text", x = 6, y = 3.2, label = "Z") +
  annotate("text", x = 7, y = 0.8, label = "Y") +
  annotate("text", x = 9, y = 0.8, label = "X") +
  annotate("text", x = 10, y = 3.2, label = "Z") +
  annotate("text", x = 11, y = 0.8, label = "Y") +
  annotate("text", x = 2, y = 3.5, label = "Fork") +
  annotate("text", x = 6, y = 3.5, label = "Collider") +
  annotate("text", x = 10, y = 3.5, label = "Chain") +
  ylim(0, 5) +
    theme_void()

```





## Fork



```{r}
citodag <- dagitty('dag {
IQ_parents [pos = "1,0"]
Number_of_books [pos = "0,1"]
Test_score_children [pos = "2,1"]

Number_of_books <- IQ_parents -> Test_score_children

}')



plot(citodag)

```


## Confounding / confounder



```{r}
citodag <- dagitty('dag {
IQ_parents [pos = "1,0"]
Number_of_books [pos = "0,1"]
Test_score_children [pos = "2,1"]

Number_of_books <- IQ_parents -> Test_score_children

Number_of_books -> Test_score_children

}')



plot(citodag)

```



## Collider


```{r}
NBAdag <- dagitty('dag {
NBA [pos = "1,0"]
Speed [pos = "0,1"]
Height [pos = "2,1"]

Height -> NBA <- Speed

}')



plot(NBAdag)

```

## Does being tall makes you fast?

```{r}

set.seed(1)

NBA <- data.frame(Lengte = rnorm(n = 500, mean = 1.8, sd = 0.15)) %>%
  mutate(Snelheid = rnorm(n = 500, mean = 10, sd = 3)) %>% # niet gecorreleerd
  mutate(Lengte_stand = Lengte/max(Lengte)) %>%
  mutate(Snelheid_stand = Snelheid/max(Snelheid)) %>%
  mutate(NBA = as.factor(ifelse(Lengte_stand + Snelheid_stand > 1.3,  1 , 0))) # selectie
  
ggplot(data = NBA, aes(x = Lengte, y = Snelheid)) +
  geom_point() +
  geom_smooth(method = "lm", aes(group = 1)) +
  scale_color_manual(values = c("#17408B", "#C9082A")) +
  xlab("Height in meters") +
  ylab("Speed in km per hour") +
  theme_bw()

```


## Does being tall makes you fast?


```{r}

ggplot(data = NBA, aes(x = Lengte, y = Snelheid, color = NBA)) +
  geom_point() +
  geom_smooth(method = "lm", aes(fill = NBA) , alpha = 0.1) +
  scale_color_manual(values = c("#17408B", "#C9082A")) +
  scale_fill_manual(values = c("#17408B", "#C9082A")) +
  xlab("Height in meters") +
  ylab("Speed in km per hour") +
  theme_bw()

```


## Chain / Mediator


```{r}
Voetbaldag <- dagitty('dag {
Team_strength [pos = "1,0"]
Budget [pos = "0,1"]
Points [pos = "2,1"]

Budget -> Team_strength -> Points

}')



plot(Voetbaldag)


```

## Budgets and points

```{r}
set.seed(123)
Budgetten <- rnorm(n = 100, mean = 30, sd = 12)

Vb <- as.data.frame(Budgetten) %>%
  mutate(Teamsterkte = (3 * Budgetten)/(3*mean(Budgetten)) * 100 + rnorm(n = 100, mean = 0, sd = 10)) %>%
  mutate(Hulppunten = 2 * Teamsterkte + rnorm(n= 100, mean = 0, sd = 15)) %>%
  mutate(Punten = round((Hulppunten - min(Hulppunten))/max((Hulppunten - min(Hulppunten)))*63 + 23, 0 )) %>% # punten in de relevante range
  select(- Hulppunten) %>%
  mutate(Teamsterkte_groep = as.factor(ntile(Teamsterkte, 10)))
  

ggplot(data = Vb, aes(x = Budgetten, y = Punten)) +
  geom_point() + 
  geom_smooth(method = "lm") +
  xlab("Budgets") +
  ylab("Points") +
  theme_bw()

```


## Budgets and Points

```{r}
ggplot(data = Vb, aes(x = Budgetten, y = Punten, color = Teamsterkte_groep)) +
  geom_point() + 
  geom_smooth(method = "lm") +
  xlab("Budgets") +
  ylab("Points") +
  labs(color = "Team strength\nGroup") +
  theme_bw()

```


## Example of a more complex DAG



```{r}


thesismodel10 <- dagitty('dag {
  GSB [pos="0,5"]
  FTE [pos="1,7"]
  DR [pos="2,7"]
  GB [pos="1,5"]
  BB [pos="1,3"]
  PRB [pos="3,5"]
  LSB [pos="4,5"]
  GM [pos="5,3"]
  PRM [pos="6,3"]
  LSM [pos="7,3"]
  FPh [pos="3,1"]
  FPc [pos="7,1"]
  SC [pos="7,7"]
  UR [pos="7,14"]  
  AF [pos="7,6"]
  BP [pos="7,15"]
  BM [pos="5,2"]
  SS [pos="7,8"]
  ST [pos="7,9"]
  AS [pos="7,10"]
  SM [pos="7,11"]
  CM [pos="7,13"]
  PQ [pos="9,5"]
  
  GSB-> GB -> PRB -> LSB -> GM -> PRM -> LSM -> PQ
  BB -> PRB
  BM -> PRM
  LSB -> LSM
  FTE -> GB
  FPh -> GB
  FPh -> FPc
  LSB -> SS
  LSB -> SM
  LSB -> ST
  LSB -> AS
  LSB -> SC
  LSB -> AF
  DR -> GB
  DR -> PRB
  LSB -> PQ
  SS -> PQ
  SM -> PQ
  ST -> PQ
  AS -> PQ
  SC -> PQ
  CM -> PQ
  UR -> PQ
  AF -> PQ
  BP -> PQ
  FPc -> PQ
  CM -> SM
}')

plot(thesismodel10)

```

# Machine learning and causality


## Machine learning and causality


```{r}


g <- dagitty('dag {
    Health_care_cost [pos="1,0"]
    Treatment [pos="0,1"]
    Proportional_shortfall [pos="2,1"]
    V_i [pos="3,0"]
   
    
V_i -> Proportional_shortfall -> Health_care_cost <- Treatment    

    

}')
plot(g)

```


## Approach

The variable $Proportional\_Shortfall$ is based on:

\begin{multline*}
\begin{split}
Proportional\_Shortfall = abs(scale(V1^3 + 2 * V2 + \\ 3 *  V3^2 + 4 * V4 + 
5 * V5 \\ + 6 * V6*V7)) + \epsilon
\end{split}
\end{multline*}

\footnotesize{

1. Fit Random Forest model on the data
2. Determine the average treatment effect with generalized random forests (grf)

We will fit 2 models for each of these steps:

a.  An analysis with all variables ("the wrong model")
b.  An analysis with all variables, except the collider $Health\_care\_cost$ ("the right model")
}




## Summary statistics

\captionsetup[table]{labelformat=empty}

```{r, results = 'asis'}

set.seed(123)

X <- matrix(sample.int(1000, size = 1000 * 8 , replace = TRUE), 
            nrow = 1000, ncol = 8)


d1 <- as.data.frame(X) %>%
  mutate(Treatment = sample(c(0, 1), 1000, replace = T)) %>%
  mutate(Ziektelast = abs(scale(V1^3 + 2 * V2 + 3 * V3^2 + 4 * V4 + 5 * V5 + 6 * V6*V7) + rnorm(n = 1000, mean = 0, sd = 1))) %>%
  mutate(Zorgkosten = 3 * Treatment + 5 * Ziektelast  + rnorm(n = 1000, mean = 0, sd = 1)) %>%
  select(Ziektelast, Treatment, Zorgkosten, everything())




stargazer(d1, type= "latex", header = FALSE, font.size = "tiny", covariate.labels = c("Proportional Shortfall", "Treatment", "Health care cost"), omit.summary.stat = c("p25", "p75"))

```





## Predictions Random Forest


```{r}
fullrun <- 0
set.seed(123)
tr_control <- trainControl(method = "repeatedcv",
                            number = 10,
                            repeats = 5)
 
if(fullrun){
 rf_fit_all <- train(Ziektelast ~ . ,
                 data= d1,
                 method = "ranger",
                 importance= "permutation",
                 trControl = tr_control)
 saveRDS(rf_fit_all, "rf_fit_all.rds")
 
 rf_fit <- train(Ziektelast ~ .-Zorgkosten ,
                 data= d1,
                 method = "ranger",
                 importance= "permutation",
                 trControl = tr_control)
 saveRDS(rf_fit, "rf_fit.rds")
} else {
 rf_fit <- readRDS("rf_fit.rds")
 rf_fit_all <- readRDS("rf_fit_all.rds")
}

dresamps <- as.data.frame(rf_fit_all$resample$Rsquared) %>%
  rename(`Wrong model`= `rf_fit_all$resample$Rsquared`) %>%
  mutate(`Right model` = rf_fit$resample$Rsquared) %>%
  gather(Model, Rsquared)

ggplot(data = dresamps, aes(x= reorder(Model, Rsquared), y = Rsquared*100)) + 
  geom_boxplot(width = 0.1) +
  ylab("R squared (%)") + 
  xlab(" ") +
  #ylim(50, 100) +
  theme_bw() + 
  coord_flip()

```

## Average treatment effect


```{r}
fullrun <- 0
# We schatten het interventie-effect met grf (causal forest)

## Verkeerde model

Treatment <- d1$Treatment
Ziektelast = d1$Ziektelast
X = as.matrix(d1[,-(1:2)])

if(fullrun){
  Y.forest <- regression_forest(X, Ziektelast)
  Y.hat <- predict(Y.forest)$predictions
  T.forest <- regression_forest(X, Treatment)
  W.hat <- predict(T.forest)$predictions
  
  cf.raw <- causal_forest(X, Ziektelast, Treatment,
                         Y.hat = Y.hat, W.hat = W.hat)
  varimp <- variable_importance(cf.raw)
  selected.idx <- which(varimp > mean(varimp))
  
  cf <- causal_forest(X[,selected.idx], Ziektelast, Treatment,
                     Y.hat = Y.hat, W.hat = W.hat,
                     tune.parameters = TRUE)
  saveRDS(cf, "cf.rds")
} else {cf <- readRDS("cf.rds")}



ATE <- average_treatment_effect(cf)


Effect_0 <- as.data.frame(ATE[1]) %>%
  rename(Effect = `ATE[1]`) %>%
  mutate(Upper = Effect + 1.96 * ATE[2]) %>%
  mutate(Lower = Effect - 1.96 * ATE[2]) %>%
  mutate(Model = as.factor("Wrong model"))

```


```{r}
fullrun <- 0

## Goede model

Treatment <- d1$Treatment
Ziektelast = d1$Ziektelast
X <- as.matrix(d1[,-(1:3)])

if(fullrun){
  Y.forest <- regression_forest(X, Ziektelast)
  Y.hat <- predict(Y.forest)$predictions
  T.forest <- regression_forest(X, Treatment)
  W.hat <- predict(T.forest)$predictions
  
  cf.raw_1 <- causal_forest(X, Ziektelast, Treatment,
                         Y.hat = Y.hat, W.hat = W.hat)
  varimp_1 = variable_importance(cf.raw_1)
  selected.idx_1 = which(varimp_1 > mean(varimp_1))
  
  cf_1 = causal_forest(X[,selected.idx_1], Ziektelast, Treatment,
                     Y.hat = Y.hat, W.hat = W.hat,
                     tune.parameters = TRUE)
  saveRDS(cf_1, "cf_1.rds")
} else{ cf_1 <- readRDS("cf_1.rds")
}

ATE_1 <- average_treatment_effect(cf_1)



Effect_1 <- as.data.frame(ATE_1[1]) %>%
  rename(Effect = `ATE_1[1]`) %>%
  mutate(Upper = Effect + 1.96 * ATE_1[2]) %>%
  mutate(Lower = Effect - 1.96 * ATE_1[2]) %>%
  mutate(Model = as.factor("Right model"))



```


```{r}
## Combine and make graph

d4 <- rbind(Effect_0, Effect_1)

ggplot(data = d4, aes(x = Model, y = Effect)) +
    geom_hline(yintercept = 0, col = "aliceblue", size = 2) +
  geom_point() +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1) + 
  xlab("") +
  ylab("Average Treatment Effect") +
  ylim(-1,1) +
  theme_bw() +
  coord_flip()


```


# Conclusion

## Conclusion

- Causal models are necessary for inference
- It is tempting to use all variables in a machine learning model
- However, this could lead to misleading conclusions

# Blog and code

## Blog and code (in Dutch)

https://misjamikkers.github.io/post/causaliteit-en-machine-learning/

https://github.com/misjamikkers/Meetup_Informatieberaad