# CIS Primer Question 2.4.1

Here are my solutions to question 2.4.1 of Causal Inference in Statistics: a Primer (CISP).

We’ll use the following simulated dataset to verify our answers. The coefficients are chosen so that each variable has approximately unit variance.

```
set.seed(29490)
N <- 10000 # sample size
sigma <- 1 # variance of nodes
e <- 0.1 # variance of errors
df <- tibble(
id = 1:N,
z1 = rnorm(N, 0, sigma),
z2 = rnorm(N, 0, sigma),
z3 = sqrt(1/2) * (z1 + z2) + rnorm(N, 0, e),
x = sqrt(1/2) * (z1 + z3) + rnorm(N, 0, e),
w = x + rnorm(N, 0, e),
y = sqrt(1/3) * (w + z3 + z2) + rnorm(N, 0, e)
)
```

## Part a

The nodes \(\{W, Z_1\}\), \(\{W, Z_2\}\), and \(\{W, Z_3\}\) are each d-separated by \(X\).

The nodes \(\{X, Y\}\) are d-separated by \(\{W, Z_1, Z_3\}\).

The nodes \(\{X, Z_2\}\) are d-separated by \(\{Z_1, Z_3\}\).

The nodes \(\{Y, Z_1\}\) are d-separated by \(\{W, Z_2, Z_3\}\).

The nodes \(\{Z_1, Z_2\}\) are d-separated by \(\emptyset\).

In the above statements, the former nodes are all conditionally independent given the d-separating set. In particular, the only two variables that are unconditionally independent are \(Z_1\), \(Z_2\).

```
part_a <- list(
w_z1 = formula(w ~ 1 + z1 + x),
w_z2 = formula(w ~ 1 + z2 + x),
w_z3 = formula(w ~ 1 + z3 + x),
x_y = formula(x ~ 1 + y + w + z1 + z3),
x_z2 = formula(x ~ 1 + z2 + z1 + z3),
y_z1 = formula(y ~ 1 + z1 + w + z2 + z3),
z1_z2 = formula(z1 ~ 1 + z2)
) %>%
map(lm, df) %>%
map_dfr(broom::tidy, .id = 'model') %>%
separate(model, c('source', 'target'), '_', remove = FALSE) %>%
filter(source == term | target == term) %>%
transmute(
model,
source,
target,
term,
lower = estimate - 2 * std.error,
estimate,
upper = estimate + 2 * std.error
)
```

## Part b

The only conditioning sets used in part a that involve \(Z_2\) were for separating \(\{Y, Z_1\}\). It’s doesn’t appear to me possible to find a d-separating set that doesn’t contain \(Z_2\). For example, the only way to block the chain \(Z_1 \rightarrow Z_3 \rightarrow Y\) is to condition on \(Z_3\), which in turn unblocks the path \(Z_1 \rightarrow Z_3 \leftarrow Z_2 \rightarrow Y\), which can only be blocked by conditioning on \(Z_2\).

## Part c

Two nodes are independent conditional on all other nodes if and only if there is a path between them consisting purely of colliders.

Conditional on all other nodes:

- \(\{W, Z_1\}\) are independent
- \(\{W, Z_2\}\) are not independent
- \(\{W, Z_3\}\) are not independent
- \(\{X, Y\}\) are independent
- \(\{X, Z_2\}\) are independent
- \(\{Y, Z_1\}\) are independent
- \(\{Z_1, Z_2\}\) are not independent

```
part_b <- list(
w = formula(w ~ 0 + x + y + z1 + z2 + z3),
x = formula(x ~ 0 + w + y + z1 + z2 + z3),
y = formula(y ~ 0 + w + x + z1 + z2 + z3),
z1 = formula(z1 ~ 0 + w + x + y + z2 + z3),
z2 = formula(z2 ~ 0 + w + x + y + z1 + z3),
z3 = formula(z3 ~ 0 + w + x + y + z1 + z2)
) %>%
map(lm, df) %>%
map_dfr(broom::tidy, .id = 'response') %>%
# filter(term != '(Intercept)') %>%
transmute(
response,
term,
lower = estimate - 2 * std.error,
estimate,
upper = estimate + 2 * std.error
)
```

## Part d

No node can be independent of its parents or children, so the minimal set for each node must include its parents and children. The minimal sets that render each variable independent of all other variables are:

- \(W\): \(\{X, Y\}\)
- \(X\): \(\{W, Z_1, Z_3 \}\)
- \(Y\): \(\{W, Z_2, Z_3 \}\)
- \(Z_1\): \(\{ X, Z_2, Z_3 \}\)
- \(Z_2\): \(\{ W, Y, Z_1, Z_3\}\)
- \(Z_3\): \(\{ W, X, Y, Z_1, Z_3 \}\)

## Part e

I originally thought that measuring just the two root nodes, \(Z_1\) and \(Z_2\) would be sufficient since every other node is a function of those two. However, comparing models via ANOVA suggests otherwise.

```
e0 <- lm(y ~ z1 + z2, data = df)
e_full <- lm(y ~ w + x + z1 + z2 + z3, data = df)
anova(e0, e_full)
```

```
Analysis of Variance Table
Model 1: y ~ z1 + z2
Model 2: y ~ w + x + z1 + z2 + z3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 9997 263.24
2 9994 101.40 3 161.84 5317.2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Since \(Y\) is a direct function of \(W, Z_2, Z_3\), restricting ourselves to those three measurements would certainly be just as good. ANOVA agrees with us on that one.

```
e1 <- lm(y ~ w + z2 + z3, data = df)
anova(e1, e_full)
```

```
Analysis of Variance Table
Model 1: y ~ w + z2 + z3
Model 2: y ~ w + x + z1 + z2 + z3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 9996 101.41
2 9994 101.40 2 0.012199 0.6012 0.5482
```

Indeed, the coefficients for the remaining variables are not statistically significantly different from 0.

`summary(e1)`

```
Call:
lm(formula = y ~ w + z2 + z3, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.40148 -0.06922 -0.00039 0.06870 0.35607
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0007178 0.0010075 -0.712 0.476
w 0.5793390 0.0058715 98.670 <2e-16 ***
z2 0.5796524 0.0043199 134.182 <2e-16 ***
z3 0.5728814 0.0100160 57.196 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1007 on 9996 degrees of freedom
Multiple R-squared: 0.9965, Adjusted R-squared: 0.9965
F-statistic: 9.364e+05 on 3 and 9996 DF, p-value: < 2.2e-16
```

I’m not sure how to prove that this set of three measurments is minimal though.

## Part f

Although \(Z_2\) has no direct causes, it does have direct effects on \(Z_3\) and \(Y\).

```
f_full <- lm(z2 ~ w + x + y + z1 + z3, data = df)
f0 <- lm(z2 ~ y + z3, data = df)
anova(f0, f_full)
```

```
Analysis of Variance Table
Model 1: z2 ~ y + z3
Model 2: z2 ~ w + x + y + z1 + z3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 9997 2784.93
2 9994 118.92 3 2666 74681 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

However, the above comparison suggests that conditioning on just those direct effects is not as good as conditioning on everything. This is due to the fact that \(Z_3\) and \(Y\) are colliders, which open up information flow between \(Z_2\) and \(W\), \(Z_1\). Adding those two to the conditioning set opens up no new paths to \(Z_2\), so we are finished.

```
f1 <- lm(z2 ~ y + z1 + z3, data = df)
f2 <- lm(z2 ~ w + y + z1 + z3, data = df)
anova(f0, f1, f2, f_full)
```

```
Analysis of Variance Table
Model 1: z2 ~ y + z3
Model 2: z2 ~ y + z1 + z3
Model 3: z2 ~ w + y + z1 + z3
Model 4: z2 ~ w + x + y + z1 + z3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 9997 2784.93
2 9996 140.02 1 2644.92 222271.83 <2e-16 ***
3 9995 118.92 1 21.09 1772.51 <2e-16 ***
4 9994 118.92 1 0.00 0.05 0.823
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

The ANOVA agrees.

## Part g

First of all, the model comparison suggests that the prediction quality does improve if we add \(W\).

```
g0 <- lm(z2 ~ z3, data = df)
g1 <- lm(z2 ~ z3 + w, data = df)
anova(g0, g1)
```

```
Analysis of Variance Table
Model 1: z2 ~ z3
Model 2: z2 ~ z3 + w
Res.Df RSS Df Sum of Sq F Pr(>F)
1 9998 4956.5
2 9997 543.6 1 4412.9 81149 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

This is because conditioning on the collider \(Z_3\) opens up the path

\[ Z_2 \rightarrow Z_3 \rightarrow Z_1 \rightarrow X \rightarrow W \]

which associates \(Z_2\) and \(W\).