overview_causal_inference.Rmd

---
title: "A Brief Overview of Causal Inference"
subtitle: ""
author: 'Todd R. Johnson, PhD <br> Professor<br><br>'
institute: "The University of Texas School of Biomedical Informatics at Houston" 
date: "`r Sys.Date()`<br><span class=footnote>Use arrow keys to navigate, press h for help</span>"
output:
  xaringan::moon_reader:
    css: ["default",  "my-fonts.css", "my-theme.css"]
    lib_dir: libs
    self-contained: false
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
      ratio: '16:9'
      beforeInit: "https://platform.twitter.com/widgets.js"
      navigation:
        scroll: false # disable slide transitions by scrolling
    includes:
      in_header: header.html  
        
---

```{r setup, include=FALSE}
options(htmltools.dir.version = FALSE)
library(dagitty)
library(ggdag)
library(tidyverse)
library(ggpmisc)
#library(hexbin)
library(plotly)
library(widgetframe)
library(DT)
library(kableExtra)
library(pcalg)
library(Rgraphviz)
library(rcausal)
library(DOT)
library(DiagrammeR)
library(corrplot)
library(lavaan)
```

# Notes

Use the arrow keys to navigate. Press `h` for help. Press `p` to see speaker notes. These notes often elaborate on slide content and include additional material and references.

My primary objective for creating this is to help me understand the basics of causal infererence. I hope to use it to help others as well, especially machine learning researchers who tend to make predictive associative models, rather than causal models.

This presentation is hosted on Github Pages at: https://tjohnson250.github.io/overview_causal_infererence/overview_causal_infererence.html#1

All of these slides, including model outcomes, are generated in R using R Markdown and the [xaringan](https://github.com/yihui/xaringan) presentation package. The source code to generate the presentation is at: https://github.com/tjohnson250/overview_causal_infererence

A shorter version, suitable for a one hour seminar, may be checked out from the `short` branch: https://github.com/tjohnson250/overview_causal_inference/tree/short

If you find errors or have suggestions for improvement, open an issue on this project's Github repo.

---
# Pearl's Three Layer Causal Hierarchy.fn[1]

```{r, echo=FALSE}
text_tbl <- data.frame(
Level = c("1. Association<br><&nbsp>\\(\\mathbf{P(y|x)}\\)", "2. Intervention<br>\\(\\mathbf{P(y|do(x), z)}\\)", "3. Counterfactuals<br>\\(\\mathbf{P(y_x|x',y')}\\)"),
"Typical Activity" = c(
"Seeing,<br>Observing",
"Doing,<br>Intervening",
"Imagining,<br>Retrospecting"
),
"Typical Questions" = c("How would seeing X change my belief in Y?", "What if I do X?", "Was it X that caused Y?<br>What if I had acted differently?"),
Example = c("What is my chance of having prostate cancer if my PSA test is 3 ng/ml?", "Will taking aspirin daily reduce my risk of stroke?", "If I had stopped smoking 10 years earlier, could I have prevented lung cancer?")
)
ladder <- text_tbl %>% knitr::kable(format="html", booktabs = FALSE, escape=FALSE, booktabs = TRUE) %>%
kable_styling(full_width = F) %>%
column_spec(1, bold = T, color = "#1c5253") %>%
column_spec(2, width = "8em") %>%
column_spec(3, width = "8em") 
ladder

```

.footnote[[1] Pearl, 2019]

???
I have modified this with what I think are more illustrative questions. This is based on the paper below, but also appears in other publications by Pearl.

Pearl J. The seven tools of causal inference, with reflections on machine learning. Communications of the ACM. 2019 Feb 21;62(3):54-60.

---
class: middle

### Part I: Causation

### Part II: Identifying and Estimating Causal Effects from Observational Data

### Part III: Causal Model Evaluation

### Part IV: Causal Discovery from Observational Data

### Part V: Counterfactual Reasoning

---
class: center, middle

# Part I: Causation

---
# Part I Learning Objectives: Causation
- Define and identify interventions, causal effects, and average causal effects
- Given a generative model, draw the corresponding causal DAG
- Given a causal DAG describe the DAG that would result from an intervention on one of the variables
- Given data on two interventions, calculate the average causal effect

---
# Example: Age, Exercise and Cholesterol
Suppose that the following model accurately describes how age, exercise, and cholesterol are related

.pull-left[
```{r codeCAEObs, tidy=FALSE}
n <- 10000
a <- rnorm(n, mean = 50, sd = 10)
e <- 0.3*a + rnorm(n)
{{c <- 0.5*a + -1*e + rnorm(n, mean = 100, sd = 5)}}
```
```{r, echo=FALSE, fig.height=4}
dobs <- tibble(E = e, A = a, C=c)
gobs <- dagitty('dag{
  A [pos="1,1"]
  E [pos="0,0"]
  C [pos="2,0"]
  A -> C A->E  E->C }')
pgobs <- ggdag(gobs, text_size=11) + theme_dag_gray()
pgobs
```
]

.pull-right[
From the model, the causal effect of exercise on cholesterol is $-1e$:

  - For every unit increase in excercise, cholesterol decreases by $1$
  
Causal DAG (Directed Acyclic Graph)
  - Graphical representation of the data generating model
  - Links represent *possible* direct causal effects
  - Missing links indicate the strong assumption of no direct causal effect
  - Must include all common causes of any pair of variables already in the DAG
]

---
## The Observational Model Produces Observational Data
```{r, echo=FALSE, out.width="80%", comment=NA}
DT::datatable(head(dobs, n=8), fillContainer = FALSE, class = "nowrap hover", options = list(pageLength=8, dom = "t")) %>% formatRound(c("E", "A", "C"), 2)
```
---
.pull-left[
### The observational data show that exercise increases cholesterol!
```{r, echo=FALSE, out.width="70%"}
formula <- y ~ x
pobs <- ggplot(dobs, aes(x = E, y = C)) + geom_point() + geom_smooth(method = "lm", formula = formula)
pobs
fitobs = lm(C ~ E, dobs)
Cobs <- coef(fitobs)
```


$c = `r round(Cobs[1], digits=2)`+`r round(Cobs[2], digits=2)`e \space `r ifelse(summary(fitobs)$coefficients[2,4] < .001, "(p<.001)", "")`$
]

--

.pull-right[
### But in the model that generated the data exercise decreases cholesterol!
```{r, eval=FALSE}
n <- 10000
a <- rnorm(n, mean = 50, sd = 10)
e <- 0.3*a + rnorm(n)
{{c <- 0.5*a + -1*e + rnorm(n, mean = 100, sd = 5)}}
```
```{r, echo=FALSE, fig.height=3}
pgobs
```
Why does the data show the opposite association?
]

---

# Definitions

- **Causal Infererence** is the process of inferring the **causal effect** of one random variable, $X$, on a second random variable. $Y$
- The **Causal Effect** of $X$ on $Y$ is the expected value of $Y$ when we **intervene** to directly set the value of $X$
- **Intervening** on $X$ means to force $X$ to take a specific value, independent of any of the variables that normally influence the value of $X$
    - Suppose that the decision to treat is normally influenced only by gender, *intervening* means ignoring gender when deciding on treatment, such as by randomly assigning whether to treat or not  
- This is why Randomized Controlled Trials (RCT's) are considered the gold standard for identifying causal effects


---
## Intervening Changes the Generative Model
Intervene by randomly setting $E$ from $0 - 20$
.pull-left[
### Observational Model
```{r, eval=FALSE}
n <- 10000
a <- rnorm(n, mean = 50, sd = 10)
{{e <- 0.3*a + rnorm(n)}}
c <- 0.5*a + -1*e + rnorm(n, mean = 100, sd = 5)
```
```{r, echo=FALSE, fig.height=3}
pgobs
```
]

.pull-right[
### Interventional Model
```{r codeCAEint}
n <- 10000
a <- rnorm(n, mean = 50, sd = 10)
{{e <- sample(0:20, n, replace = TRUE)}}
c <- 0.5*a + -1*e + rnorm(n, mean = 100, sd = 5)
```
```{r, echo=FALSE, fig.height=3}
dint <- tibble(E = e, A = a, C=c)
gint <- dagitty('dag{
  A [pos="1,1"]
  E [pos="0,0"]
  C [pos="2,0"]
  A -> C E->C }')
pgint <- ggdag(gint, text_size=11) + theme_dag_gray()
pgint
```
]
???
Intervening on $E$, means replacing the expression for $e$ from the observational model with a specific value, here one randomly assigned between $0$ and $20$. This changes the generative model and also changes the causal DAG by removing the link $A \rightarrow E$, since in the intervential model, $e$'s value no longer directly depends on $a$. 

---
## Intervening Changes The Distribution

.pull-left[
### Observational Distribution
```{r, echo=FALSE, out.width="70%"}
formula <- y ~ x
pobs <- ggplot(dobs, aes(x = E, y = C)) + geom_point() + geom_smooth(method = "lm", formula = formula)
pobs
#p <- ggplotly(p)
#frameWidget(p, width='350', height='350')
fitobs = lm(C ~ E, dobs)
Cobs <- coef(fitobs)
```

$c = `r round(Cobs[1], digits=2)`+`r round(Cobs[2], digits=2)`e \space `r ifelse(summary(fitobs)$coefficients[2,4] < .001, "(p<.001)", "")`$
]

--

.pull-right[
### Interventional Distribution
```{r, echo=FALSE, out.width="70%"}
formula <- y ~ x
pint <- ggplot(dint, aes(x = E, y = C)) + geom_point() + geom_smooth(method = "lm", formula = formula) +
   stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), sep = "~~~~")),
                formula = formula, parse = TRUE)       
pint
fitint = lm(C ~ E, dint)
Cint <- coef(fitint)
```

$c = `r round(Cint[1], digits=2)`+`r round(Cint[2], digits=2)`e \space `r ifelse(summary(fitint)$coefficients[2,4] < .001, "(p<.001)", "")`$
]

???
The coefficient for $e$ on the observational distribution is way off (indicating that exercising increases cholesterol), but the coefficient for the interventional distribution is close to the true value. Both coefficents are statistically significant at the p = .001 level.

---
# Difference Between Seeing and Doing
.pull-left[
### Observational Association
$E(C|e)$: The expected value of Cholesterol given that we observe a specific level of exercise.

$E(C|e)$ is confounded, because age $A$ affects both the amount of exercise $E$ a person gets and cholesterol levels $C$.

In the observational data
$$E(C|e) \neq E(C|do(e))$$
]
.pull-right[
### Interventional Association (Causal Effect)
$E(C|do(e))$: The expected value of Cholesterol given that we intervene and set exercise to a specific value, independent of any of the variables that normally influence the value of $E$.

$do()$ is called the do-operator.

In the interventional data
$$E(C|e) = E(C|do(e))$$
]

---
### Average Causal Effect (ACE) for a Population
Contrast of the mean values of $Y$ (outcome) given two specific interventions on $X$ (treatment)

- Dichotomous (binary) treatment and outcomes: 

    $$P(Y=1|do(X = 1)-P(Y=1|do(X = 0)))$$
    
    $$P(Survival|do(Treat))-P(Survival|do(Don't Treat))$$

- Continuous treatment and outcomes:  

    $$E(Y|do(X=x)) - E(Y|do(X=x^{\prime})$$

First equation is a special case, since for dichotomous outcomes 
  $$E(Y|do(X=x)) = P(Y=1|do(X = x))$$
???
ACE is also called the Average Treatment Effect (ATE)

Contrast ACE with the individual causal effect: the difference between counterfactual outcomes for a single patient--the outcome if treated vs. the outcome if not treated. Lack of an average causal effect does not imply lack of an individual causal effect. However, in most cases it is impossible to estimate individual causal effects from observational data because we only ever observe the effect of one level of the treatment for an individual. 

A crossover randomized experiment in which we give different treatments over time to the same person can sometimes measure individual causal effects when treatment and outcome are short term and do not affect the next treatment period. This is commonly used in cognitive science  and psychology in within-subject designs. For example, we might measure a subject's ability to monitor different numbers of patients using a telemonitoring system to determine the effect of patient load on time to respond to events, error rate, etc.

See Hernan and Robbins, Causal Inference, Chapter 1 and Fine Point 2.1 for more details in individual vs. population causal effects. https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/

---
## Estimating ACE of Exercise on Cholesterol Using the Interventional Dataset
Step 1: Estimate functional relationship between C and E using linear regression on the **interventional dataset**:

```{r, echo=FALSE}
fit = lm(C ~ E, dint)
new <- data.frame(E = c(1, 2))
est = predict(fit, new)
```

$$c = `r round(Cint[1], digits=2)`+`r round(Cint[2], digits=2)`e \space `r ifelse(summary(fit)$coefficients[2,4] < .001, "(p<.001)", "")`$$

$$E(C|do(E=e))=`r round(Cint[1], digits=2)`+`r round(Cint[2], digits=2)`e \space `r ifelse(summary(fit)$coefficients[2,4] < .001, "(p<.001)", "")`$$

Step 2: Use regression equation to calculate ACE

$$
\begin{align}
  ACE &= E(Y|do(E=1))-E(Y|do(E=2)) \\
      &= (`r round(Cint[1], digits=2)`+`r round(Cint[2], digits=2)`*1) - (`r round(Cint[1], digits=2)`+`r round(Cint[2], digits=2)`*2) \\
      &= `r round(est[2], 2)` - `r round(est[1], 2)` \\
      &= `r round(est[2]-est[1], 2)`
\end{align}
$$
Equivalently: The regression coefficent on $e$ directly tells us how much $C$ will change for each unit increase in $E$.

---
class: center, middle, Large

We say that $\mathbf{X}$ **causes** $\mathbf{Y}$ in a population when the **average causal effect** of $X$ on $Y$ is non-zero.

---
class: center, middle

## Part II: Identifying and Estimating Causal Effects from
## Observational Data

---
# Part 2 Learning Objectives: Identifying and Estimating Causal Effects from
## Observational Data

- Distinguish between total and direct causal effects with respect to a Causal DAG
- Define adjustment and show how to adjust for a variable when doing linear regression
- Given a causal DAG identify the causal and non-causal paths between two variables
- Apply the backdoor criterion to identify which variables must be adjusted for in order to identify direct and total causal effects
- Write a regression expression to estimate the causal effect of $X$ on $Y$ given a set of variables that meet the backdoor criterion
- Given the results of a linear regression that satisfies the backdoor criterion for the causal effect of $X$ on $Y$ give the causal effect estimate of $X$ on $Y$

---
## With Observational Data...
.pull-left[
#### using E to predict C gives an incorrect estimate,
```{r, echo=FALSE, out.width="70%"}
#formula <- y ~ x
#p <- ggplot(dobs, aes(x = E, y = C)) + geom_point() + geom_smooth(method = "lm", formula = formula)
pobs
#p <- ggplotly(p)
#frameWidget(p, width='350', height='350')
fit = lm(C ~ E, dobs)
Cobs <- coef(fit)
```

$c = `r round(Cobs[1], digits=2)`+`r round(Cobs[2], digits=2)`e \space `r ifelse(summary(fit)$coefficients[2,4] < .001, "(p<.001)", "")`$
]

--

.pull-right[
#### but controlling for A gives the correct estimate
```{r, echo=FALSE, out.width="70%") #, fig.width=8, fig.height=5}
#dobsstrat <- dobs %>% mutate(abin=cut(a, 10))
#ggplot(d, aes(E, C)) + geom_point(aes(color=abin))
ggplot(dobs, aes(E, C)) + geom_point() + geom_smooth(method="lm", aes(group = cut(A, 10), color=cut(A, 10)), se = FALSE)
fitobsstrat = lm(C ~ E + A, dobs)
Cintadjusted <- coef(fitobsstrat)
# Save coefficients for C = ci + ce*E + ca*A
ca <- round(Cintadjusted[3], 2)
ce <- round(Cintadjusted[2], 2)
ci <- round(Cintadjusted[1], 2) # intercept
```

$c = `r round(Cintadjusted[1], 2)`+`r round(Cintadjusted[2], 2)`e + `r round(Cintadjusted[3], 2)`a \space `r ifelse(summary(fitobsstrat)$coefficients[2,4] < .001, "(p<.001)", "")`$
]
???
The graph on the right shows the observational data with a regression line for each of several age groups. As you can see, most of the lines show the appropriate negative correlation between $C$ and $E$. Lines near the left and right are incorrect due to sparse data. Stratifying by age (here age group) is one method for controlling or adjusting for a confounder, such as $A$.

The regression equation also shows a very close estimate of the correct causal effect of $E$ on $C$. Regressing $C$ on both $E$ and $A$ controls or adjusts for $A$, returning the close estimate of the true causal effect of $E$ on $C$.

---
## Adjustment
Stratifying by age $(A)$ or including $A$ as a covariate in the regression $(C = {\beta_1}+{\beta_2}A + {\beta_3}E)$ is called **adjusting**, **conditioning**, or **controlling for A**

In R code, regression using the observational dataset (`dobs`):
.pull-left[
Unadjusted for $A$
```{r}
lm(C ~ E, dobs)
```
]
.pull-right[
Adjusted for $A$
```{r}
lm(C ~ E + A, dobs)
```
]
---
### What is the Causal Effect of Age on Cholesterol?
Using the observational dataset, regressing $C$ on $A$ and $E$ gave us: $c = `r round(Cintadjusted[1], 2)`+`r round(Cintadjusted[2], 2)`e + `r round(Cintadjusted[3], 2)`a$
.pull-left[
```{r, tidy=FALSE, eval=FALSE}
n <- 10000
a <- rnorm(n, mean = 50, sd = 10)
e <- 0.3*a + rnorm(n)
{{c <- 0.5*a + -1*e + rnorm(n, mean = 100, sd = 5)}}
```
```{r, echo=FALSE, fig.height=4}
pgobs
```
]
.pull-right[
From the model, the **direct causal effect** of $A$ on $C$ is $0.5a$

  - Regressing $C$ on $A$ and $E$ closely estimates this effect.

  - However $A$ also affects $C$ through $E$

The **total causal effect** of $A$ on $C$ includes the effect of $A$ through both paths,  which from the model is:

$$
\begin{align}
  c &= .5a + -1\mathbf{e} \\
    &= 0.5a + -1\mathbf{(0.3a)} \\
      &= 0.5a\mathbf{- 0.3a} \\
      &= 0.2a \\
\end{align}
$$

]


???
To get the total causal effect of A on C, we must regress C on A alone (`lm(C ~ A, dobs)`), controlling for E in the regression (`lm(C ~ A + E, dobs)`) gives us the direct effect of A on C.

---
### A Single Analysis is Usually Insuffient to Measure All Effects
.pull-left[
```{r, echo=FALSE, out.width="70%"}
formula <- y ~ x
p <- ggplot(dobs, aes(x = A, y = C)) + geom_point() + geom_smooth(method = "lm", formula = formula)
p

fitCAobs = lm(C ~ A, dobs)
CAobs <- coef(fitCAobs)
```

`lm(C ~ A, dobs)`

$c = `r round(CAobs[1], digits=2)`+`r round(CAobs[2], digits=2)`a \space `r ifelse(summary(fitCAobs)$coefficients[2,4] < .001, "(p<.001)", "")`$
]
.pull-right[
This simple 3 variable dataset requires two different regression analyses to estimate the causal effects of $A$ on $C$ and $E$ on $C$
  
- Total effect of $E$, Direct effect of $A$:

    `lm(C ~ A + E, dobs)`

    $c = `r round(Cintadjusted[1], 2)`+`r round(Cintadjusted[2], 2)`e + `r round(Cintadjusted[3], 2)`a \space `r ifelse(summary(fitobsstrat)$coefficients[2,4] < .001, "(p<.001)", "")`$
-  Total effect of $A$:
    `lm(C ~ A, dobs)`

    $c = `r round(CAobs[1], digits=2)`+`r round(CAobs[2], digits=2)`a \space `r ifelse(summary(fitCAobs)$coefficients[2,4] < .001, "(p<.001)", "")`$
  
Simply throwing all covariates into a single regression is insufficient... but this is what most machine learning models do
]
???
To get the direct effect of $A$ and $E$ (thereby estimating all path coefficients) we would have to regress $E$ on $A$. Later slides show that we don't need to adjust for $C$ because it is a collider. Adjusting for $C$ would give an incorrect estimate of the effect of $A$ on $E$.

---
class: center, middle

.Large[How do we know when and what to adjust for?]

---
## Paths in Causal DAGs
.pull-left[
Observational Causal DAG
```{r echo=FALSE, out.width="70%"}
pgobs
```
]
.pull-right[
A **path** is a sequence of non-intersecting adjacent edges
  - Non-intersecting means a path cannot cross a node more than once
  - Direction of the arrows doesn't matter

Examples:
  - $A \rightarrow E$
  - $A \rightarrow E \rightarrow C$
  - $E \leftarrow A \rightarrow C$
  
All associations are transmitted along paths,
but not all paths transmit association!
]

???
Source: https://www.ssc.wisc.edu/~felwert/causality/wp-content/uploads/2013/06/2-elwert_dags.pdf
See here for how to cite sources using refmanager package: https://github.com/yihui/xaringan/wiki/Bibliography-and-citations

---
class: center, middle
.Large[Two Kinds of Paths: **Causal** and **Non-Causal**]

---
## Causal Paths
.pull-left[
Observational Causal DAG
```{r echo=FALSE, out.width="70%"}
pgobs
```
]
.pull-right[
**Causal Path:** A path in which all arrows point away from the treatment to the outcome
 - Causal path for Treatment $E$ and Outcome $C$
    * $E \rightarrow C$
  
 - Causal paths for Treatment $A$ and Outcome $C$
    * $A \rightarrow E \rightarrow C$
    * $A \rightarrow C$
    
The **total causal effect** of a treatment on an outcome consists of association transmitted along all causal paths connecting them
]

---
## Non-Causal Paths
.pull-left[
Observational Causal DAG
```{r echo=FALSE, out.width="70%"}
pgobs
```
]
.pull-right[
**Non-Causal Path:** A path from the treatment to the outcome in which at least one arrow points back to the treatment
 - Non-Causal path for Treatment E and Outcome C
    * $E \leftarrow A \rightarrow C$
  
 - Non-Causal paths for Treatment $A$ and Outcome C
    * None
    
Non-causal paths are potential spurious sources of association between treatment and outcome
]
???
Technically $A$ cannot be a treatment, because we cannot externally manipulate age. Body Mass Index, Obesity, and Cholesterol levels are similar kinds of variables, because none are directly manipulable. This has led to a bit of feud between different groups of causal inference researchers on the question of whether it even makes sense to estimate the causal effect of non-manipulable causes. For some background on this, see the papers below, but also check out Pearl and Hernan's twitter feeds for a look at their different perspectives. Much of this centers on the need to test well-defined interventions. Obesity can only be manipulated indirectly, so why not just define a precise intervention and measure the causal effect of that intervention on some outcome of interest? Pearl argues that we still want to know the effect of obesity, because if we find that indirect changes in obesity has no effect on say life expectancy, then we might as well stop worrying about obesity, at least with respect to life expectancy.

Pearl J. Does obesity shorten life? Or is it the soda? On non-manipulable causes. Journal of Causal Inference. 2018 Sep 25;6(2).

Hernán MA, Taubman SL. Does obesity shorten life? The importance of well-defined interventions to answer causal questions. International journal of obesity. 2008 Aug 11;32(S3):S8.

Hernán MA. Does water kill? A call for less casual causal inferences. Annals of epidemiology. 2016 Oct 1;26(10):674-80.

---
### Causal and Non-Causal Paths for Treatment E and Outcome C
.pull-left[
Observational Causal DAG
```{r echo=FALSE, out.width="70%"}
pgobs
```
]
.pull-right[
Causal path:
  - $E \rightarrow C$
  
Non-Causal path:
  - $E \leftarrow A \rightarrow C$

    
The non-causal path is a source of confounding because it transmits spurious association between $E$ and $C$

Adjusting for $A$ in the observational data blocks this path, allowing us to use the observational data to estimate the causal effect of $E$ on $C$
]

---
## Identifiability and Estimation
A causal effect is **identifiable** if the controlled (post-intervention) distribution can be estimated from data drawn from the observational (pre-intervention) distribution.

- The **Backdoor Criterion** is one set of constraints for determining whether an effect is identifiable:

    - The causal effect of $T$ on $Y$ is **identifiable** if
we can adjust for a set of variables that

      - blocks all non-causal paths between $T$ and $Y$
      - without blocking any causal paths between $T$ and $Y$
 
**Estimation** is the act of quantifying an identifiable causal effect from finite data

- Uses statistical estimation techniques

---
### Blocking and Unblocking Paths
All Causal DAGs are a combination of these three patterns:

.pull-left[
```{r, echo=FALSE, out.width="70%", fig.height=2}
glinear <- dagitty('dag{
  A [pos="0,0"]
  E [pos="1,0"]
  C [pos="2,0"]
  A -> E E->C}')
ggdag(glinear, text_size=11) + theme_dag_gray()
```
```{r, echo=FALSE, out.width="70%", fig.height=3}
gdiverging <- dagitty('dag{
  E [pos="1,2"]
  A [pos="0,1"]
  C [pos="1,0"]
  A -> E A->C}')
ggdag(gdiverging, text_size=11) + theme_dag_gray()
```
```{r, echo=FALSE, out.width="70%", fig.height=3}
gconverging <- dagitty('dag{
  A [pos="0,2"]
  E [pos="0,0"]
  C [pos="1,1"]
  A -> C E->C}')
ggdag(gconverging, text_size=11) + theme_dag_gray()
```
]
.pull-right[
**Linear Path (or Chain)**
- Adjusting for $E$ (a mediator) *blocks* the path between $A$ and $C$

**Common Cause**
- Adjusting for $A$ *blocks* the path between $E$ and $C$

$$\space$$

**Common Effect**
- $C$ is called a collider
- Unadjusted colliders *block* the path
- Adjusting for $C$ (or a descendent of $C$)  **unblocks or opens** the path
]
???
In the linear chain, adjusting for a child of the mediator $E$ will block part of the causal path from $A$ to $C$. A child of a mediator is sometimes called a descending proxy. Here $F$ is a descending proxy:

```{r, echo=FALSE, out.width="70%", fig.height=3}
gdescproxy <- dagitty('dag{
  A [pos="0,0"]
  E [pos="1,0"]
  F [pos="1, -1"]
  C [pos="2,0"]
  A -> E E->C E->F}')
ggdag(gdescproxy, text_size=11) + theme_dag_gray()
```

If $F$ was a parent of $E$ adjusting for $F$ would not affect the path from $A$ to $C$.

---
## Identifying the Total and Direct Effects of $A$ on $C$
.pull-left[
Observational Causal DAG
```{r echo=FALSE, out.width="70%"}
pgobs
```
]
.pull-right[
**Total Effect of $A$ on $C$**

- No adjustments needed
    - All paths from $A$ to $C$ are causal
        - $A \rightarrow C$
        - $A \rightarrow E \rightarrow C$
    - `lm(C ~ A, dobs)`

**Direct Effect of $A$ on $C$**

- Adjusting for $E$ blocks the path $A \rightarrow E \rightarrow C$ leaving only the direct path $A \rightarrow C$
    - `lm(C ~ A + E, dobs)`
]
```{r, echo=FALSE}
AEadjustedfit = lm(E ~ A, dobs)
AEadjusted <- coef(AEadjustedfit)
# Save coefficients for E = ei + ea*A
ea <- round(AEadjusted[2], 2)
ei <- round(Cintadjusted[1], 2) # intercept
```
---
## Identifying the Total Effect of $A$ on $E$
.pull-left[
Observational Causal DAG
```{r echo=FALSE, out.width="70%"}
pgobs
```
]
.pull-right[
**Total (and Direct) Effect of $A$ on $E$**

- No adjustments needed
    - Causal path
        - $A \rightarrow E$
    - Non-causal path, blocked by collider $C$
        - $A \rightarrow C \leftarrow  E$
        - Adjusting for $C$ would distort the effect estimate
    - `lm(E ~ A, dobs)`
]

???
Show the effects and how adjusting distorts the effect
---
## Backdoor Criterion (Pearl 1995)
A set of variables $\{Z\}$ (possibly empty) satisfies the **backdoor criterion** relative to an ordered pair of variables $\{X, Y\}$ in a DAG if:

  1. no node in $\{Z\}$ is a descendent of $X$, and
  
  2. $\{Z\}$ blocks every path between $X$ and $Y$ that contain an arrow into $X$ (the "backdoor paths")
  
If $\{Z\}$ meets the backdoor criterion, the total causal effect of $X$ on $Y$ is non-parametrically identifiable given $\{Z\}$, such that:

$$ P(Y|do(X)) = \sum_{Z} P(Y|X,Z)P(Z)$$

The backdoor criterion recognizes that all paths that descend from $X$ are either causal or blocked
---

## Multiple Adjustment Sets
.pull-left[
```{r, echo=FALSE, out.width="70%"}
g1 <- dagitty( "dag {
    Y <- X <- Z1 <- V -> Z2 -> Y
    Z1 <- W1 -> W2 -> Z2
    X <- W1 -> Y
    X <- W2 -> Y
}")
pg1 <- ggdag(g1, text_size=11) + theme_dag_gray()
pg1
```
]
.pull-right[
Adjustment sets for $X \rightarrow Y$:

```{r, echo=FALSE}
print( adjustmentSets(g1, 'X', 'Y'))
```

Adjusting for any one of these sets of variables eliminates confouding for the causal effect of $X$ on $Y$

For a linear model this means regressing $Y$ on $X$, including as covariates only one set of these variables 

In R, any of these statements will estimate the total effect of $X$ on $Y$ from observational data `obsdata`:

```{r, eval=FALSE}
lm(obsdata, Y ~ X + W1 + W2 + Z2)
lm(obsdata, Y ~ X + V + W1 + W2)
lm(obsdata, Y ~ X + W1 + W2 + Z1)
```
]

---
# Adjustment Sets for $X \rightarrow Y$

```{r, echo=FALSE, fig.width=15}
ggdag_adjustment_set(g1, exposure='X', outcome='Y') + theme_dag_gray()
```
???
Bold edges are unadjusted, light edges are adjusted

The backdoor criterion implicitly prevents adjusting for children of mediators. For instance, suppose we want to determine the causal effect of $X$ on $Y$ in the model:

$$X \rightarrow Z \rightarrow Y$$
$$Z \rightarrow Z_2$$
where $Z$, a mediator, causes $Z_2$. It is clear that adjusting for $Z$ would block the path between $X$ and $Y$, but adjusting for $Z_2$ would partially block the path.

---
exclude: FALSE
# Frontdoor Adjustment Criterion
.pull-left[
U unmeasured: Cannot block backdoor path from $X \rightarrow Y$

```{r, echo=FALSE, out.width="70%"}
gfront <- dagitty( 'dag {
    X [pos="0,0"]
    Z [pos="1,0"]
    Y [pos="2,0"]
    U [pos="1,1"]
    X -> Z -> Y
    U -> X
    U -> Y
}')
pgfront <- ggdag(gfront, text_size=11) + theme_dag_gray()
pgfront
```
]
.pull-right[
$X \rightarrow Y$ is still identifiable using a two step process:

Step 1. Identify $X \rightarrow Z$:

  * Backdoor path blocked by collider $Y$
  
$$X \leftarrow U \rightarrow Y \leftarrow Z$$
  * In R: `lm(obsdata, Z ~ X)`
  
Step 2. Identify $Z \rightarrow Y$:

  * Backdoor path blocked by adjusting for $X$

$$Z \leftarrow X \leftarrow U \rightarrow Y$$
  * In R: `lm(obsdata, Y ~ Z + X)`
]
???
Since $U$ is unmeasured, we cannot adjust for it. However, in this case, the frontdoor adjustment criterion allows us to identify the causal effect of $X$ on $Y$ in two steps.

---
exclude: false
## Frontdoor Adjustment Criterion (Pearl 1995)
A set of variables $\{Z\}$ satisfies the **frontdoor criterion** relative to an ordered pair of variables $\{X, Y\}$ in a DAG if:

  1. $\{Z\}$ intercepts all directed paths from $X$ to $Y$,
  
  2. there is no backdoor path between $X$ and $Z$,
  
  3. every backdoor path between $Z$ and $Y$ is blocked by $X$
  
If $\{Z\}$ meets the frontdoor criterion, the causal effect of $X$ on $Y$ is identifiable given $\{Z\}$, such that:

$$P(y|do(x)) = \sum_{z} P(z|x)\sum_{x'}P(y|x', z)P(x')$$

$x$ is the value of $X$ in $do(x)$, $x'$ refers to all values of $X$

---
exclude: true
## Instrumental Variables
Instrumental variables are those that only affect the treatment, with no relationship to anything else

---
exclude: false
# Identifying Covariate Specific Effects
Suppose a drug $D$ affects $C$, and that gender $G$ affects who takes the drug.

.pull-left[
What is the effect of $E$ on $C$ when we *observe* $D = d$?
```{r, echo=FALSE, out.width="70%"}
gaecdg <- dagitty('dag{
  A [pos="1,1"]
  E [pos="0,0"]
  C [pos="2,0"]
  D [pos="1, -1"]
  G [pos="2, -1"]
  A -> C A->E  E->C D->C G->D G->C}')
pgaecdg <- ggdag(gaecdg, text_size=11) + theme_dag_gray()
pgaecdg
```
]
.pull-right[
Estimate $P(C=c|do(E=e), D=d)$


Must use an adjustment set $S$, such that $S \cup D$ blocks all backdoor paths from $E$ to $C$.

No need to block backdoor path from $D$ to $C$.

For this graph, $S = \{A\}$ and and $S \cup D = \{A, D\}$

Regression: `lm(obsdata, C ~ E + A + D)`

Mathematically:

$$\begin{aligned} P(C=c&|do(E=e), D=d) = \\  &\sum_S P(C=c|E=e, D=d, S=s)P(S=s|D=d) \end{aligned}$$

]
???
Note that the summation in the adjustment formula is over $S$, not $S \cup D$. $S$ is the set of all the other variables needed to block the backdoor paths from $E$ to $C$.

As a notational note, when you see a statement such as $P(S=s)$ and $S$ contains more than one variable, it is just shorthand for $P(S_1=s_1, S_2=s_2...)$. In this example, because $S=\{A,D\}$ $P(S=s) = P(A=a, D=d)$.

We don't want to block the backdoor path from $D$ to $C$ here, because the query is asking for the effect of $E$ on $C$ when we intervene on $E$ and then in the resulting interventional distribution observe $D=d$.

In contrast, if we just want to identify $P(C=c|do(E=e))$ we have:

$$P(C=c|do(E=e)) = P(C=c|E=e, S=s)P(S=s)$$

where $S$ is a valid adjustment set that blocks the backdoor paths from $E$ to $C$.

The basic idea is that when we look at the interventional distribution (resulting from adjusting for $S$), we must stratify by $D$ to compute the covariate-specific effect, which means we are also adjusting for $D$. Thus we have to make sure that adjusting for $D$ does not open any backdoor paths for the effect from $E$ to $C$.

To see this, consider adding high blood pressure $B$ to the graph with the paths $E \rightarrow B \leftarrow D$, meaning that both $E$ and $D$ affect $B$. If we ask how $D$ affects $C$ among those with high blood pressure, we must adjust for $G$ to block the backdoor path $D \leftarrow G \rightarrow C$ and also $B$, a collider, which opens the backdoor path $D \rightarrow B \leftarrow E \rightarrow C$. As a result, we need to add $E$ to the adjustment set, giving

$$S = \{G, E\}$$

and

$$P(C=c|do(D=d), B=b) = P(C=c|D=d, B=b, G=g, E=e)P(B=b, G=g|B=b)$$

---
exclude: false
## Identifying the Effects of Multiple Interventions
What is the effect of simultaneously intervening on $E$ and $D$?
.pull-left[
```{r, echo=FALSE, out.width="70%"}
pgaecdg
```
]
.pull-right[
Estimate $P(C=c|do(E=e, D=d))$

Must use an adjustment set $S$ that blocks all backdoor paths from $E$ to $C$ and $D$ to $C$. For this graph:
```{r echo=FALSE}
s <- adjustmentSets(gaecdg, c("E", "D"), "C")
s <- s[[1]]
```

$$S = \{`r s`\}$$

Regression: `lm(C ~ E + A + D + G)`

Mathematically:

$$\begin{aligned} P(C=c&|do(E=e, D=d)) = \\  &\sum_S P(C=c|E=e, D=d, S=s)P(S=s) \end{aligned}$$
]

---
## Confounding is a Causal Concept

* The data alone are insufficient to identify confounders or control for confouding
* Without a causal model, observational data is limited to associational models that do not tell us how to intervene
  + This includes all common machine learning and statistical approaches
  + Example: If an associational predictive model finds a few highly predictive "risk factors" we cannot assume that intervening on those risk factors will decrease risk.
* Associational predictive models are still useful if we already have interventions that are known to work
  + If we know that a treatment works if given earlier in the course of a disease, better predictive models can tell us who to treat
  + Variables that increase the predictive power of an associational model may include effects of the outcome, such as symptoms of a disease you are trying to predict

---
class: center, middle, Large

Causal Inference Requires a Causal Model,

But...

How Do we Know if Our Model is Correct?

---
class: center, middle
# Part III: Model Evaluation

---
# Part III Learning Objectives: Model Evaluation
- Identify two different approaches for evaluating causal models
- Describe the advantages and disadvantages of each approach
- Define what constitutes an equivalence class of causal DAGs

---
# Two Approaches to Model Evaluation
Local fit tests

- Tests of model-implied conditional independence assumptions

Global fit tests

- Model chi-square: comparison of model vs. observed covariance matrices

- RMSEA: Root Mean Square Error of Approximation

- CFI: Comparative fit index

- SRMSR: Standaridzed Root Mean Square Residual
???
From: Kline RB. Principles and practice of structural equation modeling. Guilford publications; 2015 Nov 3.
---
## Model Evaluation: Local Fit Tests
Causal DAG's imply independence relations

.pull-left[
```{r, echo=FALSE, out.width="70%", fig.height=2}
ggdag(glinear, text_size=11) + theme_dag_gray()
```
```{r, echo=FALSE, out.width="70%", fig.height=3}
ggdag(gdiverging, text_size=11) + theme_dag_gray()
```
```{r, echo=FALSE, out.width="70%", fig.height=3}
ggdag(gconverging, text_size=11) + theme_dag_gray()
```
]
.pull-right[
**Linear Path (or Chain)**

  - $A \perp\!\!\!\perp C | E$



**Common Cause**

  - $E  \perp\!\!\!\perp  C | A$

$$\space$$
$$\space$$

**Common Effect**

  - $A  \perp\!\!\!\perp  E$
  - $A  \not\!\perp\!\!\!\perp  E | C$
]

---
## Example: Age, Gender, BMI and Cholesterol
.small[Causal DAG Implied by Gostynski M, Gutzwiller F, Kuulasmaa K, Döring A, Ferrario M, Grafnetter D, Pajak A. [**Analysis of the relationship between total cholesterol, age, body mass index among males and females in the WHO MONICA Project.**](https://www.nature.com/articles/0802714) *International Journal of Obesity.* 2004 Aug;28(8):1082.]

.pull-left[
```{r, echo=FALSE, out.width="70%"}
gmonica <- dagitty('dag{
    eage [pos="0,3"]
    sex [pos="2,3"]
    height [pos="0,2"]
    weight [pos="1, 2"]
    bmi [pos="0,0"]
    chol [pos="1,0"]
    chrx [pos="2,0"]
    chdt [pos="2,2"]
    high [pos="2,1"]
    eage -> {weight chol}
    sex -> {chol height weight}
    height -> {weight bmi}
    weight -> bmi
    bmi -> chol
    chrx -> chol
    chdt -> chol
    high -> {chrx chdt}
}')
pgmonica <- ggdag(gmonica) + theme_dag_gray()
pgmonica
```
]

.pull-right[
This Causal DAG implies the following conditional independencies

.scroll-box-8[
```{r, echo=FALSE, out.width="70%", comment=NA}
imp <- impliedConditionalIndependencies( gmonica, max.results = Inf )
imp
```
]

Given data on these variables, we can test for these statistical relationships.
]
???
The tests shown here form a basis set for this graph, even though I did not tell dagitty to return the basis set (using option type="basis.set"). Specifying this option gives a set that when passed to localTests returns a list with p-values that are not easily accessible for later tests. The 'piecewiseSEM' package appears to have better functions for calculating and then testing basis sets. It also has a function 'fisherC' to directly calculate the C statistic.

---
### Testing Implied Independencies with Observed Data

.pull-left[
```{r, echo=FALSE, out.width="70%"}
monica <- as_tibble(read.csv("monica-chol.csv"))
#monica <- select(monica, -one_of(c('eage', 'bmi', 'chol')))
pgmonica
```
]

.pull-right[
Independence tests for all implied independencies based on `r nrow(monica)` participants:


```{r, echo=FALSE, out.width="100%", comment=NA}
r <- localTests( gmonica, monica, type = "cis", tests=imp)
r$p.value <- p.adjust(r$p.value) # Perform Holm-Bonferroni correction
DT::datatable(r %>% select(estimate, p.value), fillContainer = FALSE, class = "nowrap hover", options = list(pageLength=3, dom = 'tp', order = list(list(2, "asc")))) %>% formatRound(c("estimate", "p.value"), 2)
```
]

???
The table is sorted in increasing order by `p.values`

`p.values` were corrected for multiple comparisons using the Holm-Bonferroni method

The null hypothesis is that the estimate is 0, which would indicate independence, so if $p < .05$ we can reject the null, meaning that the indpendence implication is unlikely to be true.

See Figure 1 in http://johannes-textor.name/papers/2017-ije.pdf for details (Textor J, van der Zander B, Gilthorpe MS, Liśkiewicz M, Ellison GT. Robust causal inference using directed acyclic graphs: the R package ‘dagitty’. *International journal of epidemiology.* 2016 Dec 1;45(6):1887-94.)

Tests for independence are linear in the dagitty package. Distance Correlation can detect non-linear independence For example:

```{r, eval=FALSE}
# Too slow to run each time slides are processed
library(energy) # provides test for marginal independence using distance correlation 

# Tests show probability of independence
# Low p values indicate dependence

#  Distance correlation is computationally intensive $O(n^2)$ and memory intensive, so use a random sample of data
monicasamp <- sample_n(monica, 1000)
dcor.ttest(monica$eage, monica$sex) # marginally independent

# Independence under partial correlation (whether distance or linear) is not identical to conditional independence. 
# The CDCSIS package (below) tests for conditional independence using distance correlation, whereas the energy package (above) only tests for partial distance correlation. 
# If all variables are multivariate gaussian, the partial linear correlation is 0 if and only if  X and Y are conditionally linearly independent given Z, but otherwise this is not true. 
# See the following paper for more details on the relationship between conditional independence of partial correlation: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-842X.2004.00360.x

# See section 4.3 of (http://real.mtak.hu/44934/1/pdcov.pdf) Székely GJ, Rizzo ML. Partial distance correlation with methods for dissimilarities. The Annals of Statistics. 2014;42(6):2382-412.

# Test independence under partial correlation independence hypothesis using distance correlation. 
pdcor.test(monicasamp$bmigrp, monicasamp$weight, monicasamp$bmi, R=99) # Dependent although model implies independence

library(cdcsis) # Provides cdcsis and cdcsis.test for conditional distance correlation

# compute conditional distance correlation 
cdcsis(as.matrix(monicasamp$bmigrp), monicasamp$weight, as.matrix(monicasamp$bmi), threshold = 1) # must give threshold here as 1 for this data
# Test conditional independence hypothesis using distance correlation
cdcov.test(as.matrix(monicasamp$bmigrp), monicasamp$weight, as.matrix(monicasamp$bmi)) # Dependent at the .05 level, but not .01
```

---
# Coefficients where $p < .05$

.pull-left[
```{r, echo=FALSE, fig.height=6}
r <- r[r$p.value<0.05,] 
pr <- ggplot(r, aes(x=estimate, y=rownames(r))) + geom_point() + geom_errorbarh(mapping=aes(y=rownames(r), xmin=!!sym("2.5%"), xmax = !!sym("97.5%")), height=0) + theme_bw(base_size = 15)
pr
```
]

.pull-right[
The data suggests that $eage$ and $height$, and $eage$ and $high$ are *not* marginally independent

This contradicts our model

Any ideas as to why $eage$ and $high$ are not marginally independent?
]

???
$high$ indicates that a patient was told that he or she had high cholesterol. This means that  $high$ is dependent on a past measurement of high cholesterol, which itself is likely tied to past age, weight, etc. A more accurate model would have these past values along with edges leading to the present values.
---
## Add $eage \rightarrow height$ and $eage \rightarrow high$ and Retest

.pull-left[
```{r, echo=FALSE, out.width="70%"}
gmonica <- dagitty('dag{
    eage [pos="0,3"]
    sex [pos="2,3"]
    height [pos="0,2"]
    weight [pos="1, 2"]
    bmi [pos="0,0"]
    chol [pos="1,0"]
    chrx [pos="2,0"]
    chdt [pos="2,2"]
    high [pos="2,1.5"]
    eage -> {weight chol height high}
    sex -> {chol height weight}
    height -> {weight bmi}
    weight -> bmi
    bmi -> chol
    chrx -> chol
    chdt -> chol
    high -> {chrx chdt}
}')
pgmonica <- ggdag(gmonica, layout="circle") + theme_dag_gray()
pgmonica
```
]

.pull-right[
Coefficients where $p < .05$

```{r, echo=FALSE, fig.height=6}
r <- localTests( gmonica, monica, "cis" )
r$p.value <- p.adjust(r$p.value)
rlowp <- r[r$p.value<0.05,] 
pr <- ggplot(rlowp, aes(x=estimate, y=rownames(rlowp))) + geom_point() + geom_errorbarh(mapping=aes(y=rownames(rlowp), xmin=!!sym("2.5%"), xmax = !!sym("97.5%")), height=0) +theme_bw(base_size = 15)
pr
```
]
???
There are still many discrepencies here

---
exclude: false
## Summarizing Independence Tests: Fisher's C
- Statistical test of the set of independence test results

$$C = -2\sum \ln p_i$$
 - $p_i$ are the $p$-values of all tests of conditional independence for the **basis set**
- **Basis Set**: the smallest set of independence relations implied by a DAG
- Fisher's C requires a basis set constructed by listing "each unique pair of nonadjacent
vertices $(v_i, v_j)$ and then forming the conditioning set of this pair by listing all other vertices in the DAG that are parents of either $v_i$ or of $v_j$.".fn[1]
```{R, echo=FALSE}
# dagitty does not return the basis set required for Fisher's C
# so we define functions based on Shipley, 2000 that return
# the correct basis set
library(dagitty)
basisSet <- function(dag, v) {
  vars <- unique(tidy_dagitty(gmonica)$data$name)
  vars2 <- vars[!vars %in% c(adjacentNodes(gmonica, v), v)]
  lapply(vars2, function(v2) list(X = v, Y=v2, Z=unique(c(parents(dag, v), parents(dag, v2)))))
}
shipleyBasisSet <- function(dag) {
  vars <- unique(tidy_dagitty(gmonica)$data$name)
  result <- unlist(lapply(vars, function(v) basisSet(gmonica, v)), recursive=FALSE)
  return(result)
}

gmonicaBasisSet <- shipleyBasisSet(gmonica)

r = localTests(gmonica, monica, tests=gmonicaBasisSet)
```

- $C$ has a $\chi^2$ distribution with $2k$ degrees of freedom, where $k$ is the number of elements in the basis set

 - $p \geq .05$ is good
 - If $p < .05$ reject model, consider alternatives: more links, different structure, etc.

- For the WHO MONICA DAG: $C$ = `r round(-2*sum(log(r$p.value)), 2)`, $p$ = `r round(1 - pchisq(-2*sum(log(r$p.value)), 2 * length(r$p.value)), 2)`

.footnote[[1] Shipley, B. A new inferential test for path models based on directed acyclic graphs. *Structural Equation Modeling*. 2000 Jun 1;7(2):206-18.]
???
See Shipley: A New Inferential Test for Path Models Based on Directed Acyclic Graphs: https://pdfs.semanticscholar.org/23e5/29d073a66c314f19900e5e1e7fc510453812.pdf

The basis set required for this test (according to Shipley is): "To obtain the basis set Bø from a DAG, simply list each unique pair of nonadjacent
vertices (vi, vj) and then form the conditioning set of this pair by listing all other vertices in the DAG that are parents of either vi or of vj."

To calculate $C$ and the $\chi^2$ $p$-value in R, given that the $p$-values for the WHO MONICA DAG are in `r$p.value`:

```{r}
C = -2*sum(log(r$p.value))
k = length(r$p.value) # Number of tests in basis set
p = 1 - pchisq(C, 2 * k)
cat("p = ", p)
```

The implied independencies shown in earlier slides happen to be a  basis set for this graph (which is not guaranteed for all graphs when the function is called without type="basis.set"). 'dagitty' can produce a basis set using the option type="basis.set" as shown here, but when supplied with the resulting implications, its localTests function did not return a value that allows easy extraction of the p values, so I used the basis set shown on the previous slide.

There can be more than one basis set for a given DAG, but not 
all basis sets are mutually independent in finite data samples.

Section 3.4 of 

Shipley B. Cause and correlation in biology: a user's guide to path analysis, structural equations and causal inference with R. Cambridge University Press; 2016 Apr 18. discusses these issues and includes a description of Shipley's method for producing a basis set that is independent in finite samples. This is also described in the paper:

Shipley B. A new inferential test for path models based on directed acyclic graphs. Structural Equation Modeling. 2000 Jun 1;7(2):206-18. https://pdfs.semanticscholar.org/23e5/29d073a66c314f19900e5e1e7fc510453812.pdf

---
# Limitations
- Most packages only detect (lack of) linear dependence
 - Independent variables are guaranteed to have $0$ correlation
 - But $0$ correlation does not imply independence when the relationship is nonlinear
 - A Partial correlation of $0$ (for conditional independence) implies conditional independence only when the joint distribution of the tested variables is multinormal
 
- Distance correlation (Szekely, et al. 2007) and conditional distance correlation (Wang, et al., 2015) can detect nonlinear dependence
 - But complexity is at worst $\mathcal{O}(n^2)$ and at best $\mathcal{O}(n\log n)$ (Chaudhuri and Hu, 2019)
 
???
 Székely GJ, Rizzo ML, Bakirov NK. Measuring and testing dependence by correlation of distances. The annals of statistics. 2007;35(6):2769-94.
 
  Wang X, Pan W, Hu W, Tian Y, Zhang H. Conditional distance correlation. Journal of the American Statistical Association. 2015 Oct 2;110(512):1726-34.
 
 Chaudhuri A, Hu W. A fast algorithm for computing distance correlation. Computational Statistics & Data Analysis. 2019 Jul 1;135:15-24.
 
---
## Global Fit Tests

```{r, echo=FALSE, message=FALSE, warning = FALSE}
model <- '
  # regressions
  height ~ eage + sex
  weight ~ height + eage + sex
  bmi ~ weight + height
  hychol ~ bmi + sex + eage + chdt + chrx
  chdt ~ high
  chrx ~ high
  '
fit <- sem(model, data=monica, std.ov=TRUE) # , ordered = c("hychol")
```

.pull-left-65[
.scroll-box-20[
All variables standardized before analysis

```{r, echo=FALSE, comment=NA}
summary(fit, fit.measures=TRUE)
```
]
]

.pull-right-35[
Recommended levels for global fit tests

- Model $\chi^2$ with df and $p$-value
 - prefer $p \geq .05$
- Comparative Fit Index
 - prefer CFI $> 0.90$
- Root mean square error of Approximation
 - prefer lower 90%CI to be $< .05$
- Standardized Root Mean Square Residual
 - prefer SRMR $< 0.10$
]
???
Pearl argues against global fit tests because they rarely tell you where your model might be wrong. In addition, many different models can often fit the same data. He also points out that for estimating specific effects, not all of the model needs to be correct. As long as the local independence tests are valid for the effect of interest, you can still use the effect estimates.
 
---
### Visualize Mismatch between Model and Data Covariance Matrices
This shows observed - model-implied correlations
```{r, echo=FALSE, out.width="70%", fig.height=4}
# switch this to corrr to remove low res plot: https://drsimonj.svbtle.com/exploring-correlations-in-r-with-corrr

# Function from: https://psu-psychology.github.io/r-bootcamp-2018/talks/lavaan_tutorial.html
plot_matrix <- function(matrix_toplot){
corrplot::corrplot(matrix_toplot, is.corr = FALSE,
               type = 'lower',
               order = "original",
               tl.col='black', tl.cex=.75)
}
plot_matrix(residuals(fit, type="cor")$cov) # 
#plot_matrix(residuals(fit))
```
---
exclude: false
## Examine Modification Indices
MI Indicates possible improvement in model test score
```{r, echo=FALSE}
modificationindices(fit, minimum.value = 20)
```

The highest MI involved adding an undirected edge between $chdt$ and $chrx$, indicating an unmeasured common cause.

---
## Age, Exercise, Cholesterol Model Implies No Independence Relations
.pull-left[
```{r, echo=FALSE, out.width="70%"}
ggdag(gobs, text_size=11)
```
]
.pull-right[
```{r}
print( impliedConditionalIndependencies( gobs ) )
```
* C is a collider between $A$ and $E$, but $A \rightarrow E$ **shields** the collider, rendering $A$ and $E$ dependent even when $C$ is not conditioned on

There are no implied independence relationships to test.
]

---
## Independence (d-separation) Equivalence
.pull-left[
Two Causal DAGs over the same set of variables are independence equivalent iff they have:

* the same adjacencies (edges between variables ignoring direction)
* the same unshielded colliders (any two nodes pointing into the collider are not adjacent)
  * In graphs 1 and 3, $C$ is a *shielded* collider, because $A$ and $E$ have an edge between them
  * If we deleted the edge between $A$ and $E$, then $C$ would be an unshielded collider

These graphs are all independence equivalent.  (Verma and Pearl, 1988)
]
.pull-right[
```{r, echo=FALSE, fig.height = 7}
ggdag_equivalent_dags(gobs) + theme_dag_gray()
```
]
???
This concept is explained in detail in the following lecture at around 28 minutes:

From https://www.youtube.com/watch?v=TISSNwWDfw4&list=PLO5mmwQolPRX858CyOOIHqYnmdzlHjIgS&index=8

(CCD Summer short course 2016:Day 2 Part 1)

---
exclude: true
### Completed Partially Directed Acyclic Graph (CPDAG)
.pull-left[
```{r, echo=FALSE, fig.height = 7}
# ggdag_equivalent_class(gobs) # Bug as of 5/24/19. See https://github.com/malcolmbarrett/ggdag/issues/8
# Use dagitty function instead
plot(equivalenceClass(gobs))
```
]
.pull-right[
A CPDAG represents an equivalance class of independence equivalent graphs

* A directed edge is shared by all graphs in the class
* An undirected edge is reversible as long as no cycles or additional unshielded colliders are introduced
]

---
exclude: true
### Global Tests of AEC Model are Uninformative
There are 5 free parameters and 0 degrees of freedom.

.scroll-box-18[
```{r, echo=FALSE}
model <- '
  # regressions
  E ~ A
  C ~ A + E
  '
fit <- sem(model, data=dobs)
summary(fit, standardized=TRUE)
```
]

---
exclude: true
### Age, Sex, Gender, BMI, and Cholesterol
.pull-left[
Hypothesized Graph

```{r, echo=FALSE, out.width="70%"}
pgmonica
```
]
.pull-right[
CPDAG

```{r, echo=FALSE, fig.height=4}
# ggdag_equivalent_class(gobs) # Bug as of 5/24/19. See https://github.com/malcolmbarrett/ggdag/issues/8
# Use dagitty function instead
plot(equivalenceClass(gmonica))
```
]

---
class: middle, center
# Part IV: Causal Discovery from Observational Data

---
# Part IV Learning Objectives: Causal Discovery from Observational Data

- Describe how conditional independence tests of a dataset constrains the possible causal DAGs that may have generated the data
- Describe the steps of the PC causal discovery algorithm

---
# Inferring Causal Structure from Data

.center[**Causal DAG `\\(\Rightarrow\\)` Independence Relations**]

.pull-left[
```{r, echo=FALSE, out.width="70%"}
pgmonica
```
]

.pull-right[
.scroll-box-16[
```{r, echo=FALSE, out.width="70%", comment=NA}
impliedConditionalIndependencies( gmonica, max.results = Inf )
```
]
]

.center[**Causal DAGS `\\(\Leftarrow\\)` Independence Relations**]

---
# PC Algorithm.fn[1]
1. Adjacency Phase
    - Determine which variables are adjacent
1. Orientation Phase
    - Orient direction of arrows if possible

.footnote[[1] Spirtes and Glymour, 1991]

---
### PC Algorithm: Adjacency Phase

$X$ and $Y$ are not adjacent if they are conditionally independent given any subset of variables that does not contain $X$ and $Y$

- Begin with a fully connected graph
- Remove edge between $X$ and $Y$ if $X \perp\!\!\!\perp Y | \text{ any set } S$
    - Begin with the empty set and test all connected pairs for marginal independence $X \perp\!\!\!\perp Y$
        - Remove edges for any pair that is independent 
    - For any remaining connected pair, test conditional independence vs. each single other variable, ie., $S$ of size 1
        - Remove edge between $X$ and $Y$ if $X \perp\!\!\!\perp Y | S \text{ of size 1}$
    - Repeat for $S$ of size 2, 3, etc.

---
background-image: url(PCAlignmentPhaseExample.png)
background-size: contain
.footnote[From https://www.ccd.pitt.edu/wp-content/uploads/2017/02/CCD_Summer_Course_2016_Day3_morning.pdf]

---
background-image: url(PCOrientationRules.png)
background-size: contain
.footnote[From https://www.ccd.pitt.edu/wp-content/uploads/2017/02/CCD_Summer_Course_2016_Day3_morning.pdf]
???
X -- Y -- Z, X and Z must not be adjacent. If Y is not in S, then X and Z, then Y is a collider for X and Z. If Y separates X and Z, then we cannot orient the edges (other than to say that Y is not a collider)

---
background-image: url(PCOrientationPhaseExample.png)
background-size: contain
.footnote[From https://www.ccd.pitt.edu/wp-content/uploads/2017/02/CCD_Summer_Course_2016_Day3_morning.pdf]

---
exclude: true
## PC Algorithm: Adjacency Phase Example
.pull-left[
Generative DAG
```{r, echo=FALSE}
knownDAG <- dagitty('dag{
  X1 [pos="0,2"]
  X2 [pos="0,0"]
  X3 [pos="1,1"]
  X4 [pos="2,1"]
  X1 -> X3
  X2 -> X3
  X3 -> X4}')
ggdag(knownDAG, text_size=4)  + theme_dag_blank()
```

Begin with: 
]
.pull-right[
Independencies

\begin{align}
X1 \perp\!\!\!\perp X2 & \\
X1 \perp\!\!\!\perp X4 &\mid X3 \\
X2 \perp\!\!\!\perp X4 &\mid X3 \\
\end{align}

```{r, echo=FALSE}
fcDAG <- dagitty('dag{
  X1 [pos="0,2"]
  X2 [pos="0,0"]
  X3 [pos="1,1"]
  X4 [pos="2,1"]
  X1 <-> X2
  X1 <-> X3
  X1 <-> X4
  X2 <-> {X3 X4}
  X3 <-> X4}')
ggdag(fcDAG, text_size=4) + theme_dag_blank()
```
]

---
exclude: true
### Results of the PC Algorithm (Spirtes and Glymour, 1991) On The WHO Monica Dataset 

.pull-left[
```{r, echo=FALSE, fig.height=6}
pcmonica <- monica %>% select(-bmigrp) %>% select(-hychol) # Remove bmigrp and hychol to match earlier model
suffStat <- list(C = cor(pcmonica), n = nrow(pcmonica))
pc.monica <- pc(suffStat, indepTest = gaussCItest, labels = names(pcmonica), alpha = 0.01, skel.method="stable", maj.rule = TRUE)
plot(pc.monica)
```
]

.pull-right[
Graph implied by [Nature paper]
```{r, echo=FALSE, fig.height=6}
pgmonica
```
]

---
## Results of PC Algorithm on the WHO Dataset
.pull-left[
```{r, echo=FALSE}
pcmonica <- monica %>% select(-bmigrp) %>% select(-hychol) # Remove bmigrp and hychol to match earlier model
res <- tetradrunner(algoId = "pc-all", df = pcmonica, alpha = .01, testId = "fisher-z", dataType = 'continuous', stableFAS = TRUE, colliderDiscoveryRule = 2)
graph_dot <- tetradrunner.tetradGraphToDot(res$graph)
grViz(graph_dot)
```
]

.pull-right[
Hypothesized DAG
```{r, echo=FALSE, fig.height=6}
pgmonica
```
]


???
The PC algorithm is one of the first causal discovery algorithms. As you can see from the inferred graph, many of the arrows are incorrect and some edges are not directional, meaning the algorithm could not orient the edge.

One limitation of PC (and most other algoirithms) is that the independence tests assume linear relationships. As a result, the algorithm will mistake nonlinearly dependent variables as independent.

See: Spirtes P, Glymour C. An algorithm for fast recovery of sparse causal graphs. Social science computer review. 1991 Apr;9(1):62-72.
---
## Adding Background Knowledge Improves Discovery

```{r, results="hide"}
# List of forbidden directed edges
forbid <- list(c('chol', 'bmi'), c('eage', 'bmi'), c('sex', 'bmi'))

# List of required directed edges
require <- list(c('height', 'bmi'), c('weight', 'bmi'))

# Forbid edges between eage and sex
forbiddenWithin <- c('sex', 'eage')
class(forbiddenWithin) <- 'forbiddenWithin'

# Add Temporal constraints
temporal <- list(forbiddenWithin, c('height', 'weight', 'chrx', 'chdt', 'high'), c('bmi', 'chol')) 

prior <- priorKnowledge(requiredirect = require, addtemporal = temporal, forbiddirect = forbid)
```
???
Since bmi is defined only by bmi, we forbid other variables from causing it. We could likewise forbid the two other defined variables in the dataset, eageg and hychol, but the algorithm works fine as is.

The temporal constraints indicate that lower tiers cannot have any directed edges to higher tiers. For instance, `bmi` and `chol` can have directed edges to `bmigrp` and `hychol`, as well as between themselves, but cannot have directed edges to the other variables that appear before them in the temporal list.
---

# PC Algorithm with Knowledge
.pull-left[
```{r, echo=FALSE}
pcmonica <- monica %>% select(-bmigrp) %>% select(-hychol) # Remove bmigrp and hychol to match earlier model
res <- tetradrunner(algoId = "pc-all", df = pcmonica, alpha = .01, testId = "fisher-z", dataType = 'continuous', stableFAS = TRUE, priorKnowledge = prior, verbose = TRUE, colliderDiscoveryRule = 2)
graph_dot <- tetradrunner.tetradGraphToDot(res$graph)
grViz(graph_dot)
```
]

.pull-right[
Causal discovery from observational data is difficult

Adding background knowledge helps, but there is no magic here.

The same algorithm, with different settings can produce a different graph.

There are many different algorithms to choose from. 
]

---
exclude: true
## Adding Knowledge: Required Edges

`TRUE` means the pair of variables must be connected by an edge

```{r, echo=FALSE, eval=FALSE}
# Required Edges
k <- matrix(FALSE, nrow = length(names(monica)), ncol = length(names(monica)))
rownames(k)<-names(monica)
colnames(k)<-names(monica)
k[["bmi", "bmigrp"]] <- TRUE # bmi must `cause` bmigrp
k[["bmigrp", "bmi"]] <- TRUE
k[["height", "bmi"]] <- TRUE
k[["bmi", "height"]] <- TRUE
k[["weight", "bmi"]] <- TRUE
k[["bmi", "weight"]] <- TRUE
k[["eage", "eageg"]] <- TRUE
k[["eageg", "eage"]] <- TRUE
k[["height", "weight"]] <- TRUE
k[["weight", "height"]] <- TRUE
k[["chol", "hychol"]] <- TRUE
k[["hychol", "chol"]] <- TRUE
as_tibble(k, rownames="id") %>% mutate_all(~cell_spec(.x, color = ifelse(.x == TRUE, "red", "gray"))) %>% knitr::kable(format="html", booktabs = FALSE, escape=FALSE, booktabs = TRUE) %>%
kable_styling(full_width = F)
```

---
exclude: true
## Adding Knowledge: Forbidden Edges

A $1$ means the row variable (e.g., `height`) must not have direct edge to the column variable (e.g., `bmi`)

```{r, echo=FALSE, eval=FALSE}
# Forbidden Edges
kf <- matrix(FALSE, nrow = length(names(monica)), ncol = length(names(monica)))
rownames(kf)<-names(monica)
colnames(kf)<-names(monica)
kf["sex",] <- c(TRUE,FALSE,TRUE,FALSE,FALSE,FALSE,TRUE,TRUE,TRUE)
kf["eage",] <- c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE,TRUE,TRUE)
kf["eageg",] <- c(TRUE,FALSE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE)
kf["chol",] <- c(FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE)
kf["weight",] <- c(FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,FALSE,TRUE,TRUE)
kf["height",] <- c(FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,FALSE,TRUE,TRUE)
kf["bmi",] <- c(TRUE,TRUE,TRUE,FALSE,FALSE,FALSE,FALSE,TRUE,FALSE)
kf["hychol",] <- c(TRUE,TRUE,TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE)
kf["bmigrp",] <- c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE,TRUE,TRUE)

as_tibble(kf, rownames="id") %>% mutate_all(~cell_spec(.x, color = ifelse(.x == TRUE, "red", "gray"))) %>% knitr::kable(format="html", booktabs = FALSE, escape=FALSE, booktabs = TRUE) %>%
kable_styling(full_width = F)
```
---
exclude: true
### Results of the PC Algorithm With Background Knowledge 

.pull-left[
```{r, echo=FALSE, fig.height=6, eval=FALSE}
suffStat <- list(C = cor(monica), n = nrow(monica))
pc.monica <- pc(suffStat, indepTest = gaussCItest, labels = names(monica), alpha = 0.05, skel.method="stable", fixedEdges = k, fixedGaps = kf)
plot(pc.monica)
```
]

.pull-right[
Graph implied by [Nature paper]
```{r, echo=FALSE, fig.height=6}
pgmonica
```
]


---
class: center, middle

## Part V: Counterfactual Reasoning

---
## Counterfactual Reasoning Example 
- Joe is $60$ years old, exercised $10$ units, and has total cholesterol of $140$.
 - What would Joe's cholesterol be if he had exercised $20$ units instead of $10$?

---
## Counterfactual Reasoning: Naïve Approach
- Joe is $60$ years old, exercised $10$ units, and has total cholesterol of $140$.
 - What would Joe's cholesterol be if he had exercised $20$ units instead of $10$?
 
- Use the estimated causally correct regression equation from observational data, with $A=60$ and $E=20$:

$$
\begin{align}
c &= `r round(Cintadjusted[1], 2)`+`r round(Cintadjusted[2], 2)`e + `r round(Cintadjusted[3], 2)`a \\
 &= `r round(Cintadjusted[1], 2)`+`r round(Cintadjusted[2], 2)`\cdot 20 + `r round(Cintadjusted[3], 2)`\cdot 60 \\
 &= `r round(Cintadjusted[1], 2) + round(Cintadjusted[2], 2) * 20 + round(Cintadjusted[3], 2)*60`
\end{align}
$$
- This is incorrect, because it ignores what we have learned about Joe by measuring what has already happened to him

 - This approach estimates $E(C|do(E=20, A=60))$
 - We need to estimate the expected value of Joe's cholesterol if he had exercised at $20$ units, given the previously observed evidence regarding Joe:
$$E(C_{E=20}|E=10, A=60, C=C_{E=10}=140)$$

---
## Causal DAGs as Structural Equations
A causal DAG is a set of structural equations with endogenous and exogenous variables.

.pull-left[
```{r tidy=FALSE, eval=FALSE}
n <- 10000
a <- rnorm(n, mean = 50, sd = 10)
e <- 0.3*a + rnorm(n)
c <- 0.5*a + -1*e + rnorm(n, mean = 100, sd = 5)
```

```{r, echo=FALSE, fig.height=4}
pgobs
```
]

.pull-right[
$$
\begin{align}
A &= U_a \\
E &= `r ea` A + U_e \\
C &= `r ca` A + `r ce` E + U_c
\end{align}
$$

$U = \{U_a, U_e, U_c\}$ is the set of exogenous variables
- Exogenous variables represent background factors with no explantory mechanism encoded in the model

$V = \{A, E, C\}$ is the set of endogenous variables

- The values of endogenous variables are defined through functions encoded in the model

Given a set of values for $U$, the model functions uniquely define the values of $V$.
]
???
ON the left is our known true model. On the right are the  structural equations that we estimated earlier from observational data, but with the intercepts replaced with exogenous variables.

---
## What would Joe's Cholesterol Have Been With E=20?
.pull-left[
**Step 1:** Use inferred equations and prior observations from Joe to determine Joe's background factors, $U$

Observed: $A = 60$, $E = 10$, $C = 140$

Substitute and solve:
$$
\begin{align}
A = 60 &= U_a \\
E = 10 &= `r ea` \cdot 60 + U_e \\
C = 140 &= `r ca` \cdot 60 + `r ce`\cdot 10 + U_c \\
\\
U_a &= 60 \\
U_e &= `r 10 - ea*60` \\
U_c &= `r 140-ca*60-ce*10`
\end{align}
$$
```{r, echo=FALSE}
# Save Uc for use below
Uc <- 140-ca*60-ce*10
```
]

.pull-right[
**Step 2:** Using inferred equations:

$$
\begin{align}
A &= U_a \\
E &= `r ea` A + U_e \\
C &= `r ca` A + `r ce` E + U_c
\end{align}
$$

Replace the function for $E$ with $20$, substitute Joe's values for $U$ and the remaining $V$ values to determine the counterfactual

$$E(C_{E=20}|E=10, A=60, C=C_{E=10}=140)$$
$$
\begin{align}
A &= 60 \\
E &= 20 \\
C_{E=20} &= `r ca` \cdot 60 + `r ce`\cdot 20 + `r Uc` = `r ca*60+ce*20+Uc`
\end{align}
$$
naïve version had $C = `r round(Cintadjusted[1], 2) + round(Cintadjusted[2], 2) * 20 + round(Cintadjusted[3], 2)*60`$
]

---
### Why the Difference?
The naïve estimate,

$$E(C|do(E=20, A=60)) = `r round(Cintadjusted[1], 2) + round(Cintadjusted[2], 2) * 20 + round(Cintadjusted[3], 2)*60`$$

does not reflect the fact the Joe's cholesterol, at $140$, is higher than expected for a person of Age 60 with 10 units of exercise:

$$E(C|do(E=10, A=60)) = `r round(Cintadjusted[1], 2) + round(Cintadjusted[2], 2) * 10 + round(Cintadjusted[3], 2)*60`$$

By using Joe's observed data to estimate $U$, then substituting $E=20$ for the original function $E = 0.3A+U_e$ the resulting structural equations estimate Joe's counterfactual result:

$$E(C_{E=20}|E=10, A=60, C=C_{E=10}=140) = `r ca*60+ce*20+Uc`$$

---
## Counterfactual Reasoning Requires a Causal Model
To compute the background factors $U$ for an individual we need the structural equations.

An incorrect model and set of equations will give an incorrect counterfactual result

A predictive associational model is not causal, so cannot be used for counterfactual reasoning:

- We should not use it to play what if games, such as what would your risk be if you exercised more?

---
# Use of Counterfactual Reasoning
- Average Causal Effect applies to populations, whereas counterfactual reasoning can answer questions about an individual
- Effect of Treatment on the Treated (ETT): 
$$ETT = E[Y_1 - Y_0|X = 1]$$
where $Y_1$ is outcome when treated, $Y_2$ is outcome when not treated, and $X = 1$ means all those treated

- Required for Mediation Analysis: what is the direct effect of a treatment, vs the indirect effect as mediated by other factors?
- Effect of additive interventions
 - What is the expected value of an outcome if we add 5mg to the population's existing dosage?

---
# Conclusion
Causal inference from observational data is difficult, but possible

Causal knowledge is required to reason about interventions--associative predictive or risk models are not enough

Analysis of observational health data requires causal inference

Causal Inference Remains Underutilized and Difficult to Master

- Despite considerable research, textbooks, software packages and tutorials, the area is difficult and confusing to learn--there is no single best introductory textbook

- Schisms in the field, such as between the Causal DAG, [Potential Outcomes](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/), and other groups (see [Andrew Gelman's blog](https://statmodeling.stat.columbia.edu/)) further add to the confusion.

- Awareness of causal inference among those who could benefit from it's application is still quite low, but that is slowly changing

---
# To Learn More
Non-technical intro
- Pearl J, Mackenzie D. [The book of why: the new science of cause and effect.](http://bayes.cs.ucla.edu/WHY/) Basic Books; 2018 May 15.

Best technical intro

- Pearl J, Glymour M, Jewell NP. [Causal inference in statistics: A primer.](http://bayes.cs.ucla.edu/PRIMER/) John Wiley & Sons; 2016 Mar 7.

Best Practical Intro (including potential outcomes and Causal DAGs)

- Hernán MA, Robins JM (2019). [Causal Inference.](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/) Boca Raton: Chapman & Hall/CRC, forthcoming. 

---
# Thanks
.pull-right[.pull-down[

<a href="mailto:todd.r.johnson@uth.tmc.edu">
`r icon::fa("paper-plane")` todd.r.johnson@uth.tmc.edu
</a>

<a href="https://github.com/tjohnson250/overview_causal_inference">
`r icon::fa("link")` https://github.com/tjohnson250/overview_causal_inference
</a>

<a href="http://twitter.com/johnsontoddr4">
`r icon::fa("twitter")` @johnsontoddr4
</a>

<a href="https://github.com/tjohnson250">
`r icon::fa("github")` @tjohnson250
</a>

]]




---
exclude: true
### Does Estrogen Replacement Therapy (ERT) Cause Uterine Cancer?
In the 1970's

- ERT widely prescribed to women after menopause
- Several studies found an association between ERT and uterine cancer
- ERT known to cause bleeding in some women
- Uterine cancer can cause bleeding

Two Hypotheses

- ERT causes uterine cancer
- ERT does not cause cancer, but causes bleeding, leading to a uterine exam, which increases the chance of finding undiagnosed uterine cancer

Yale researchers proposed limiting data analysis to women who have bled

Harvard researchers argued that such an analysis could find an association between ERT and Cancer even if ERT did not cause cancer

---
exclude: true
## Causal DAGs for each hypothesis
.pull-left[
```{r, echo=FALSE}
gERTspurious <- dagitty('dag{
  ERT [pos="0,1"]
  B [pos="1,1"]
  DC [pos="2,1"]
  UC [pos="0,0"]
  ERT -> B B -> DC UC -> B UC-> DC}')
ggdag(gERTspurious, text_size=4) #, use_labels="label")
```
]
---
exclude: true
## Adjust done at data collection or analysis (from Elwert)

---
exclude: true
## Causal DAGs and Joint Probability Distributions
.pull-left[
Any distribution generated by a Markovian causal DAG $M$, can be written as:

$$P(v_1,v_2,...,v_n) = \prod_{i} P(v_{i}|pa_i)$$

where $V_1, V_2,..., V_n$ are the variables in $M$, and $pa_i$ are the values of the parents of $V_i$ in M.
]

.pull-right[
```{r, echo=FALSE, fig.height=3}
ggdag(gobs, text_size=11)
```

$$P(a,e,c)=P(a)*P(e|a)*P(c|a, e)$$
]

For example, the probability that a patient has A=50, E = 5, and C = 120 is given by:

$$P(A=50,E=5,C=120)=P(a=50)*P(e=5|a=50)*P(c=120|a=50, e=5)$$


.whisper[Technical notes: A graph is Markovian if there are no cycles and all error terms jointly indpendent.]
  
---
exclude: true
# Overadjustment Bias


---
exclude: true
## Intervening Changes the Generative Model
.pull-left[
Dichotomous Observational Generative Model
```{r}
n <- 10000
g <- rbinom(n, 1, .5)
{{p_x_g <- c(1/3, 2/3)
x <- rbinom(n, 1, p_x_g[g])}} # Male: 1/3, Female: 2/3 
# gender  treatment p(heart attack)
#  male       0           .3
#  male       1           .4
# female      0           .05
# female      1           .07
p_h_gx = c(.3, .4, .05, .075)
h <- rbinom(n, 1, p_h_gx[g*2+x+1])
```
```{r, echo=FALSE, out.width="100%"}
dobs <- tibble(G = g, X = x, H=h)
library(dagitty)
library(ggdag)
coords <- list(
    x = c(x = 0, g = 1, h = 2),
    y = c(x = 0, g = 1, h = 0)
    )
gobs <- dagify(
    x ~ g,
    h ~ g + x,
    exposure = "x",
    outcome = "h",
    labels = c("h" = "Heart\n Attack", 
                  "g" = "Gender",
                  "x" = "Treatment"),
    coords = coords)

ggdag(gobs, text = FALSE, use_labels = "label")
```
]
.pull-right[
Dichotomous Interventional Generative Model
```{r}
n <- 10000
g <- rbinom(n, 1, .5) 
{{x <- rbinom(n, 1, .5)}} # Random chance
#          m¬t  mt  f¬t  ft
p_h_gx = c(.3, .4, .05, .075) # p(H|M,T)
h <- rbinom(n, 1, p_h_gx[g*2+x+1])
```
```{r, echo=FALSE, out.width="100%"}
dint <- tibble(G = g, X = x, H=h)
library(dagitty)
library(ggdag)
coords <- list(
    x = c(x = 0, g = 1, h = 2),
    y = c(x = 0, g = 1, h = 0)
    )
gint <- dagify(
    h ~ g + x,
    exposure = "x",
    outcome = "h",
    labels = c("h" = "Heart\n Attack", 
                  "g" = "Gender",
                  "x" = "Treatment"),
    coords = coords)

ggdag(gint, text = FALSE, use_labels = "label")
```
```{r}
# CPT of h given g and x
p_h = dint %>% group_by(G, X) %>% summarise(p_h = sum(H == 1)/n())
xtabs(p_h ~ G+X, p_h)
```
]

---
exclude: true
# CPT in R
P(h|g,x)
```{r}
# CPT of h given g and x, with marginals
p_h = dint %>% group_by(G, X) %>% summarise(p_h = sum(H == 1)/n())
addmargins(xtabs(p_h ~ G+X, p_h))
```

---
exclude: true
## GGPlotly graph
```{r, echo=FALSE}
set.seed(100)
d <- diamonds[sample(nrow(diamonds), 1000), ]
p <- ggplot( data = d, aes(x = carat, y = price, text = paste("Clarity:", clarity)) ) +
  geom_point() +
  facet_wrap(~ cut) +
 geom_smooth(aes(group=cut, color=cut, fill=cut), method='loess') 
p <- ggplotly(p)
frameWidget(p)
```