# 第 49 章 有向無環圖 DAG & 可怕的因果 Causal Terror

The most newsworthy scientific studies are the least trustworthy. Maybe popular topics attract more and worse researchers, like flies drawn to the smell of honey?
Richard McElreath

set.seed(1914)
N <- 200 # num grant proposals
p <- 0.1 # proportion to select

# uncorrelated newsworthiness and trustworthiness
nw <- rnorm(N)
tw <- rnorm(N)

# select top 10 of combined scores
s <- nw + tw # total score
q <- quantile( s, 1-p ) # top 10% threshold

selected <- ifelse( s >= q, TRUE, FALSE)
cor(tw[selected], nw[selected])
## [1] -0.76800831
proposal <- data.frame(nw, tw, selected)
with(proposal, plot(nw[!selected], tw[!selected], col=c("black"),
xlim = c(-3,3.5), ylim = c(-3.5,3),
bty="n",
xlab = "newsworthiness",
ylab = "trustworthiness"))
points(proposal$nw[proposal$selected],
proposal$tw[proposal$selected],
col = c("blue"),
pch = 16)
abline(lm(proposal$nw[proposal$selected] ~ proposal$tw[proposal$selected]),
lty = 2, lwd = 2, col = c("blue"))
text(1, -2.8, "rejected")
text(2, 2.5, "selected", col = c("blue"))

## 49.1 多重共線性問題 multicollinearity

N <- 100                        # number of individuals
set.seed(909)
height <- rnorm(N, 10, 2)       # sim total height for each
leg_prop <- runif(N, 0.4, 0.5)  # leg as proportion of height
leg_left <- leg_prop * height + # sim left leg as proportion + error
rnorm( N, 0, 0.02 )
leg_right <- leg_prop * height + # sim right leg as proportion + error
rnorm( N, 0, 0.02 )
# combine into data frame
d <- data.frame(height, leg_left, leg_right)
##       height  leg_left leg_right
## 1  5.9314170 2.6794115 2.7092861
## 2  6.5129833 2.6764277 2.6800068
## 3  9.3466279 3.9271549 3.9849467
## 4  9.2330335 3.9641911 3.9933889
## 5 10.3571282 4.4275932 4.4187658
## 6 10.0889226 4.9566406 4.9718779

m6.1 <- quap(
alist(
height ~ dnorm( mu, sigma ),
mu <- a + bl * leg_left + br * leg_right,
a ~ dnorm( 10, 100 ),
bl ~ dnorm(2, 10),
br ~ dnorm(2, 10),
sigma ~ dexp( 1 )
), data = d
)
precis(m6.1)
##             mean          sd        5.5%      94.5%
## a     0.98127908 0.283955398  0.52746352 1.43509465
## bl    0.21185848 2.527037063 -3.82683482 4.25055177
## br    1.78367737 2.531250606 -2.26174999 5.82910472
## sigma 0.61710260 0.043434269  0.54768625 0.68651895

plot(precis(m6.1))

post <- extract.samples(m6.1)
plot(bl ~ br, post, col = col.alpha(rangi2, 0.1),
pch = 16)

sumblbr <- post$bl + post$br
dens(sumblbr, col = rangi2, lwd = 2, xlab = "sum of bl and br")

m6.2 <- quap(
alist(
height ~ dnorm( mu, sigma ),
mu <- a + bl * leg_left,
a ~ dnorm( 10, 100 ),
bl ~ dnorm(2, 10),
sigma ~ dexp( 1 )
), data = d
)

precis(m6.2)
##             mean          sd       5.5%      94.5%
## a     0.99793262 0.283646205 0.54461121 1.45125404
## bl    1.99206762 0.061157037 1.89432686 2.08980838
## sigma 0.61860375 0.043539981 0.54901845 0.68818905

### 49.1.1 哺乳動物奶質量數據中的共線性

data(milk)
d <- milk
d$K <- standardize( d$kcal.per.g )
d$F <- standardize( d$perc.fat )
d$L <- standardize( d$perc.lactose )

# kcal.per.g regressed on perc.fat
m6.3 <- quap(
alist(
K ~ dnorm( mu, sigma ),
mu <- a + bF * F,
a ~ dnorm( 0, 0.2 ),
bF ~ dnorm( 0, 0.5 ),
sigma ~ dexp(1)
), data = d
)

# kcal.per.g regressed on perc.lactose
m6.4 <- quap(
alist(
K ~ dnorm( mu, sigma ),
mu <- a + bL * L,
a ~ dnorm( 0, 0.2 ),
bL ~ dnorm( 0, 0.5 ),
sigma ~ dexp(1)
), data =  d
)

precis( m6.3 )
##                mean          sd        5.5%      94.5%
## a     1.5355260e-07 0.077251946 -0.12346338 0.12346368
## bF    8.6189697e-01 0.084260876  0.72723181 0.99656212
## sigma 4.5101794e-01 0.058707561  0.35719192 0.54484396
precis( m6.4 )
##                 mean          sd        5.5%       94.5%
## a      7.4388949e-07 0.066616333 -0.10646502  0.10646651
## bL    -9.0245498e-01 0.071328484 -1.01645167 -0.78845828
## sigma  3.8046526e-01 0.049582594  0.30122270  0.45970783

m6.5 <- quap(
alist(
K ~ dnorm( mu, sigma ),
mu <- a + bF * F + bL * L,
a ~ dnorm( 0, 0.2 ),
bF ~ dnorm( 0, 0.5 ),
bL ~ dnorm( 0, 0.5 ),
sigma ~ dexp( 1 )
), data = d
)
precis( m6.5 )
##                 mean          sd        5.5%       94.5%
## a     -3.1721359e-07 0.066035771 -0.10553823  0.10553760
## bF     2.4349835e-01 0.183578648 -0.04989579  0.53689248
## bL    -6.7808254e-01 0.183776698 -0.97179320 -0.38437188
## sigma  3.7674180e-01 0.049183935  0.29813637  0.45534722

pairs( ~ kcal.per.g + perc.fat + perc.lactose, data = d,
col = rangi2)

## 49.2 治療後偏倚 post-treatment bias

set.seed(71)
# number of plants
N <- 100

# simulate initial heights
h0 <- rnorm(N, 10, 2)

# assign treatments and simulate fungus and growth
treatment <- rep(0:1, each = N/2)
fungus <- rbinom( N, size = 1, prob = 0.5 - treatment * 0.4)
h1 <- h0 + rnorm( N, 5 - 3*fungus )

# compose a clean data frame
d <- data.frame( h0 = h0, h1 = h1, treatment = treatment, fungus = fungus)
##           h0        h1 treatment fungus
## 1  9.1363156 14.345788         0      0
## 2  9.1056256 15.623924         0      0
## 3  9.0428548 14.386665         0      0
## 4 10.8342908 15.837422         0      0
## 5  9.1641987 11.469124         0      1
## 6  7.6256722 11.107757         0      0
precis(d)
##                mean         sd       5.5%     94.5%    histogram
## h0         9.959780 2.10116231  6.5703278 13.078737 ▁▂▂▂▇▃▂▃▁▁▁▁
## h1        14.399204 2.68808705 10.6180024 17.933694     ▁▁▃▇▇▇▁▁
## treatment  0.500000 0.50251891  0.0000000  1.000000   ▇▁▁▁▁▁▁▁▁▇
## fungus     0.230000 0.42295258  0.0000000  1.000000   ▇▁▁▁▁▁▁▁▁▂

### 49.2.1 設定模型

\begin{aligned} h_{1,i} & \sim \text{Normal}(\mu_i, \sigma) \\ \mu_i & = h_{0,i} \times p \end{aligned}

• $$h_{0,i}$$ 是在時間 $$t = 0$$ 時的植物高度；
• $$h_{1,i}$$ 是在時間 $$t = 1$$ 時植物的高度；
• $$p$$ 是比例係數，也就是 $$h_{1,i}$$$$h_{0,i}$$ 之間的比值，$$p = \frac{h_{1,i}}{h_{0,i}}$$。如果 $$p = 1$$ 說明在時間 $$t = 1$$ 時植物並沒有比在時間 $$t = 0$$ 時有長高。

$\beta \sim \text{Log-Normal}(0, 0.25)$

sim_p <- rlnorm(10000, 0, 0.25)
precis(sim_p)
##          mean         sd       5.5%     94.5%    histogram
## sim_p 1.03699 0.26298937 0.67068301 1.4963969 ▁▁▃▇▇▃▁▁▁▁▁▁
dens( sim_p, xlim = c(0,3), adj = 0.1)

m6.6 <- quap(
alist(
h1 ~ dnorm(mu, sigma),
mu <- h0 * p ,
p ~ dlnorm( 0, 0.25 ),
sigma ~ dexp(1)
), data  = d
)
precis(m6.6)
##            mean          sd      5.5%     94.5%
## p     1.4266259 0.017609916 1.3984818 1.4547699
## sigma 1.7932857 0.125172622 1.5932357 1.9933357

$$p$$ 的事後概率分佈均值是 1.43，也就是預估平均每單位時間植物會長高大約 40%。接下來如果加入另外兩個變量，治療組，和是否出現菌落。我們會把這兩個變量對植物施加的影響使用線性回歸模型的方式加到 $$p$$ 上去：

\begin{aligned} h_{1, i} & \sim \text{Normal}(\mu_i, \sigma) \\ \mu_i & = h_{0,i} \times p \\ p & = \alpha + \beta_T T_i + \beta_F F_i \\ \alpha & \sim \text{Log-Normal}(0, 0.25) \\ \beta_T & \sim \text{Normal}(0, 0.5) \\ \beta_F & \sim \text{Normal}(0, 0.5) \\ \sigma & \sim \text{Exponential}(1) \end{aligned}

m6.7 <- quap(
alist(
h1 ~ dnorm(mu, sigma),
mu <- h0 * p ,
p <- a + bt * treatment + bf * fungus,
a ~ dlnorm( 0, 0.2 ),
bt  ~ dnorm(0, 0.5),
bf  ~ dnorm(0, 0.5),
sigma ~ dexp(1)
), data = d
)
precis(m6.7)
##                mean          sd         5.5%        94.5%
## a      1.4813914677 0.024510693  1.442218647  1.520564289
## bt     0.0024122219 0.029869649 -0.045325247  0.050149691
## bf    -0.2667189147 0.036547721 -0.325129231 -0.208308598
## sigma  1.4087974415 0.098620704  1.251182509  1.566412374

m6.8 <- quap(
alist(
h1 ~ dnorm(mu, sigma),
mu <- h0 * p,
p <- a + bt * treatment,
a ~ dlnorm(0, 0.25),
bt ~ dnorm(0, 0.5),
sigma ~ dexp(1)
), data = d
)
precis(m6.8)
##              mean          sd        5.5%     94.5%
## a     1.381693513 0.025197600 1.341422882 1.4219641
## bt    0.083654229 0.034313231 0.028815059 0.1384934
## sigma 1.746295968 0.121908002 1.551463435 1.9411285

plant_dag <- dagitty("dag{
H_0 -> H_1
F -> H_1
T -> F
}")
coordinates( plant_dag ) <- list(x = c(H_0 = 0 , T = 2, F = 1.5, H_1 = 1),
y = c(H_0 = 0 , T = 0, F = 0, H_1 = 0))
drawdag(plant_dag)

# define our coordinates
dag_coords <-
tibble(name = c("H0", "H1", "M", "F", "T"),
x    = c(1, 2, 2.5, 3, 4),
y    = c(2, 2, 1, 2, 2))

# save our DAG
dag <-
dagify(F ~ M + T,
H1 ~ H0 + M,
coords = dag_coords)

# plot
dag %>%
ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
geom_dag_point(aes(color = name == "M"),
alpha = 1/2, size = 6.5, show.legend = F) +
geom_point(x = 2.5, y = 1,
size = 6.5, shape = 1, stroke = 1, color = "orange") +
geom_dag_text(color = "black") +
geom_dag_edges() +
scale_color_manual(values = c("steelblue", "orange")) +
theme_dag()
set.seed(71)
N <- 1000
h0 <- rnorm(N, 10, 2)
treatment <- rep(0:1, each = N/2)
M <- rbern(N)
fungus <- rbinom( N, size = 1, prob = 0.5 - treatment * 0.5 + 0.4 * M)
h1 <- h0 + rnorm( N, 5 + 3 * M)
d2 <- data.frame( h0 = h0, h1 = h1, treatment = treatment, fungus = fungus)
precis(d2)
##                mean         sd       5.5%     94.5%      histogram
## h0        10.123143 1.99233421  6.9383408 13.428682 ▁▁▂▂▅▇▇▅▃▂▁▁▁▁
## h1        16.627464 2.64352784 12.3803465 20.692560      ▁▁▃▇▇▇▂▁▁
## treatment  0.500000 0.50025019  0.0000000  1.000000     ▇▁▁▁▁▁▁▁▁▇
## fungus     0.461000 0.49872610  0.0000000  1.000000     ▇▁▁▁▁▁▁▁▁▇
# incorrectly included fugus
m6.7 <- quap(
alist(
h1 ~ dnorm(mu, sigma),
mu <- h0 * p ,
p <- a + bt * treatment + bf * fungus,
a ~ dlnorm( 0, 0.2 ),
bt  ~ dnorm(0, 0.5),
bf  ~ dnorm(0, 0.5),
sigma ~ dexp(1)
), data = d2
)
precis(m6.7)
##              mean          sd        5.5%       94.5%
## a     1.510467777 0.014175426 1.487812708 1.533122846
## bt    0.066779619 0.015124914 0.042607085 0.090952152
## bf    0.160976680 0.015183931 0.136709826 0.185243533
## sigma 2.109488889 0.047098698 2.034216072 2.184761705
# the correct model to see the treatment effect
m6.8 <- quap(
alist(
h1 ~ dnorm(mu, sigma),
mu <- h0 * p,
p <- a + bt * treatment,
a ~ dlnorm(0, 0.25),
bt ~ dnorm(0, 0.5),
sigma ~ dexp(1)
), data = d2
)
precis(m6.8)
##               mean           sd         5.5%        94.5%
## a      1.625584764 0.0096235572  1.610204461 1.6409650674
## bt    -0.016667992 0.0136202768 -0.038435825 0.0050998411
## sigma  2.222812199 0.0496210313  2.143508207 2.3021161910

## 49.3 對撞因子偏倚 collider bias

grant_dag <- dagitty("dag{
T -> S
N -> S
}")
coordinates( grant_dag ) <- list(x = c(T = 0.5, S = 1, N = 1.5),
y = c(T = 0, S = 0, N = 0))
drawdag(grant_dag)

### 49.3.1 虛假的傷心對撞因子 collider of false sorrow

marriage_dag <- dagitty("dag{
H -> M
A -> M
}")
coordinates( marriage_dag ) <- list(x = c(H = 0.5, M = 1, A = 1.5),
y = c(H = 0, M = 0, A = 0))
drawdag(marriage_dag)

1. 每年有20名實驗對象出生，且他們擁有符合均一分佈特徵的幸福感。
2. 每年，實驗對象年齡自然會增加一歲。然而幸福感並不會因年齡的增加而增加或減少。
3. 當實驗對象18歲時，有些人會結婚。結婚本身的比值 (odds) 則於該實驗對象的幸福感成一定的比例關係 (proportional)。
4. 當一名實驗對象結婚了以後，她/他保持結婚的狀態，不會離婚。
5. 年齡到65歲之後，該實驗對象離開本次研究。
d <- sim_happiness( seed = 1977, N_years = 1000)
precis(d)
##                     mean         sd       5.5%      94.5%     histogram
## age        3.3000000e+01 18.7688832  4.0000000 62.0000000 ▇▇▇▇▇▇▇▇▇▇▇▇▇
## married    3.0076923e-01  0.4587690  0.0000000  1.0000000    ▇▁▁▁▁▁▁▁▁▃
## happiness -1.0000698e-16  1.2144211 -1.7894737  1.7894737      ▇▅▇▅▅▇▅▇

d %>%
mutate(married = factor(married,
labels = c("unmarried", "married"))) %>%

ggplot(aes(x = age, y = happiness)) +
geom_point(aes(color = married), size = 1.75) +
scale_color_manual(NULL, values = c("grey85", "forestgreen")) +
scale_x_continuous(expand = c(.015, .015)) +
theme(panel.grid = element_blank())

$\mu_i = \alpha_{\text{MID}[i]} + \beta_A A_i$

• $$\text{MID}[i]$$ 是實驗對象 $$i$$ 是否已經結婚的索引變量 (index variable)，當它等於1時表示單身，等於2時表示已婚。
• 這其實是我們人為地給已婚者和單身者的幸福感和年齡之間關係的直線設定了各自的截距。

d2 <- d[ d$age > 17, ] # adults only d2$A <- (d2$age - 18) / (65 - 18 ) 經過上述的數據處理，我們使得年齡變量 A 的範圍控制在 0-1 之間，其中 0 代表 18 歲，1 代表 65 歲。幸福感則是一個範圍在 -2, 2 之間的數值。這樣的話，假定年齡和幸福感之間呈現的是極其強烈的正關係，那麼這最極端的斜率也就是 $$(2 - (-2)) / 1 = 4$$。所以，一個較為合理的斜率的先驗概率分佈，可以是95%的斜率取值分佈在小於極端斜率之內的範圍。其次是為截距 $$\alpha$$ 設定合理的先驗概率分佈。因為 $$\alpha$$ 本身是年齡等於零，也就是18歲時的幸福感，我們需要這個數據能夠覆蓋所有可能的幸福感取值，-2，2 之間。那麼，標準正（常）態分佈是一個不錯的選擇。 d2$mid <- d2$married + 1 # construct the marriage status index variable m6.9 <- quap( alist( happiness ~ dnorm(mu, sigma), mu <- a[mid] + bA * A, a[mid] ~ dnorm(0, 1), bA ~ dnorm(0, 2), sigma ~ dexp(1) ), data = d2 ) precis(m6.9, depth = 2) ## mean sd 5.5% 94.5% ## a[1] -0.23508769 0.063489864 -0.33655675 -0.13361862 ## a[2] 1.25855173 0.084959885 1.12276943 1.39433404 ## bA -0.74902744 0.113201124 -0.92994470 -0.56811018 ## sigma 0.98970796 0.022557999 0.95365592 1.02576000 看，這個模型似乎很確定，年齡和幸福感是呈現負關係的。對比一下沒有假如婚姻狀態的變量的模型給出的估計結果： m6.10 <- quap( alist( happiness ~ dnorm(mu, sigma), mu <- a + bA * A, a ~ dnorm( 0, 1 ), bA ~ dnorm(0, 2), sigma ~ dexp(1) ), data = d2 ) precis(m6.10) ## mean sd 5.5% 94.5% ## a 1.6492480e-07 0.076750151 -0.12266140 0.12266173 ## bA -2.7286203e-07 0.132259761 -0.21137692 0.21137637 ## sigma 1.2131876e+00 0.027660803 1.16898030 1.25739491 m6.10 才是正確的模型。它正確的給出了年齡和幸福感之間並無關係的結果。模型 m6.9 錯誤地把對撞因子 – 婚姻狀況作為預測變量之一加入了模型中，而婚姻狀況在這個數據背景下，同時是年齡和幸福感的結果（common consequence of age and happiness）。結果就會像 m6.9 那樣，出現年齡和幸福感之間虛假的負關係 (false negative association between the two causes)。m6.9 告訴我們看起來似乎年齡的增長和幸福感呈現負關係，這種僅僅在模型中給出的變量之間的關係應該嚴格來說只能算是一種“統計學的關係 (statistical association)”，不能算是真實的因果關係 (causal association)。正如圖 49.11 所顯示的。當已知實驗對象是已婚或者未婚，實驗對象的年齡似乎能告訴我們他/她的幸福度。看綠色點的部分，這些人都是已婚者，年齡越大，越多人結婚，那麼這個已婚人群的幸福度數值就會平均被拉低。相似的，看空白點的部分，這些人都是未婚者，年齡越大，其中幸福感較強的人都結婚而加入了已婚人群陣營，那麼剩下的人就會感覺幸福度越來越低。所以，當把人群分成未婚和已婚兩個部分的話，這兩個人群中的幸福度都隨著年齡增加而呈現下降趨勢。但是，我們知道，這並不是真實的因果關係。 碰撞因子偏倚本身會出現在當模型中的預測變量加入了某個同時是結果和某一個預測變量的結果的變量（common consequence）。 ### 49.3.2 對撞因子偏倚另一實例（未測量變量造成的碰撞偏倚） 當我們知道並了解了這樣的對撞因子偏倚現象之後，DAG圖非常有助於幫助我們避免陷入這樣的困境。但是，最可怕的其實是未知（未測量）變量可能造成的對撞偏倚。 假設我們想通過數據了解父母親的受教育程度 (P)，祖父母的受教育程度 (G)，和子女的學習成績 (C)，之間的關係，特別是 P, G 對 C 可能的貢獻或者影響。我們很容易就能理解圖 49.12 中所表示的這三個變量之間應該存在的相互因果關係。祖父母的受教育程度很顯然也會對子女的學習成績造成影響 ($$G \rightarrow C$$)。 dag6.3.2 <- dagitty( "dag{ G -> P; P -> C; G -> C }" ) coordinates(dag6.3.2) <- list( x=c(G=0,P=1,C=1) , y=c(G=0,P=0,C=1) ) drawdag( dag6.3.2 ) 如果此時我們發現，可能有一個未被測量的變量可能同時影響父母 (P)，和子女 (C)，但是對祖父母的教育程度並無直接影響，例如父母和子女生活的社區的環境 (U)，通常祖父母不一定和子女和孫輩生活在一個社區。那麼這個未觀察的社區變量就很可能會造成一個碰撞偏倚： dag6.3.2_u <- dagitty( "dag{ G -> P; P -> C; G -> C ; U -> P ; U -> C U [unobserved]}" ) coordinates(dag6.3.2_u) <- list( x=c(G=0,P=1,C=1,U=1.3) , y=c(G=0,P=0,C=1,U=0.5) ) drawdag( dag6.3.2_u ) 如果因果關係 49.13 成立的話，那麼 P 就是 G 和 U 之間的對撞因子。如果此時我們建立 $$G \rightarrow C$$ 模型時把 $$P$$ 加入預測變量中，即便我們並沒有測量這個 U 變量，對撞因子偏倚也變得不可避免。 我們實際使用計算機模擬來理解這一過程。數據模擬符合下列條件 1. P 是 G, U 共同影響的結果 (common consequence)。 2. C 則是 G, U, P 三者共同預測的結果。 3. G 和 U 不受任一變量的影響（完全獨立）。 N <- 200 # number of grandparent-parent-child traids b_GP <- 1 # direct effect of G on P b_GC <- 0 # direct effect of G on C b_PC <- 1 # direct effect of P on C b_U <- 2 # direct effect of U on P and C 上面的模擬參數其實類似於回歸模型中的回歸係數。注意到我們假設祖父母對孫輩學習成績的影響是 0。接下來我們用這些回歸係數來採集一些符合上述設定的隨機樣本數據： set.seed(1) U <- 2 * rbern( N, 0.5 ) - 1 G <- rnorm( N ) P <- rnorm( N, b_GP * G + b_U * U) C <- rnorm( N, b_PC * P + b_GC * G + b_U * U ) d <- data.frame(C = C, P = P, G = G, U = U) 那麼，該怎樣設定一個模型來分析祖父母對孫輩學習成績的影響呢？我們認為祖父母的教育程度會通過父母親這條通路，影響到子女的學習成績，所以模型中應該要控制 P。那麼下面的代碼建立的是一個簡單的線性回歸模型，結果變量是 C，預測變量是 P 和 G，注意這裡我們假裝沒有測量到，也不知道 U 變量的存在： m6.11 <- quap( alist( C ~ dnorm( mu, sigma ), mu <- a + b_PC * P + b_GC * G, a ~ dnorm(0, 1), c(b_PC, b_GC) ~ dnorm(0, 1), sigma ~ dexp(1) ), data = d ) precis(m6.11) ## mean sd 5.5% 94.5% ## a -0.11747518 0.099195743 -0.27600914 0.041058772 ## b_PC 1.78689146 0.044553553 1.71568628 1.858096644 ## b_GC -0.83895371 0.106140454 -1.00858666 -0.669320767 ## sigma 1.40948908 0.070111393 1.29743753 1.521540626 看模型 m6.11 給出的父母子女的學習成績的影響是多麼的顯著。甚至比我們設定的關係還要大兩倍。這並不奇怪，因為這裡給出的一部分 $$P \rightarrow C$$ 的關係應該歸功於 $$U$$，但是該模型本身不知道 $$U$$ 的存在。但是更令人感到驚訝的是，本來設定的祖父母不影響子女學習成績的關係在這個模型中變得顯著呈現負相關，這是有悖常識的。再怎麼說祖父母的受教育程度越高也不能給子女的學習成績帶來負面的影響才是。這個數學模型本身並沒有錯，但是它不能被賦予一個因果關係的解釋。這裏的對撞因子偏倚是如何形成的呢？ d$G_st <- standardize(d$G) d$C_st <- standardize(d$C) d$P_st <- standardize(d$P) P_45 <- quantile(d$P_st, prob = .45)
P_60 <- quantile(d\$P_st, prob = .60)

with(d, plot(G_st[(P_st < P_45 | P_st > P_60) & U == -1],
C_st[(P_st < P_45 | P_st > P_60) & U == -1],
col = c("black"),
xlim = c(-3,3), ylim = c(-2.5, 2.3),
bty="n",
main = "Parents in 45th to 60th centiles",
xlab = "grandparent education (G)",
ylab = "grandchild education (C)"))
with(d, points(G_st[(P_st >= P_45 & P_st <= P_60) & U == -1],
C_st[(P_st >= P_45 & P_st <= P_60) & U == -1],
col = c("black"),
pch = 16))
with(d, points(G_st[(P_st < P_45 | P_st > P_60) & U == 1],
C_st[(P_st < P_45 | P_st > P_60) & U == 1],
col = rangi2))
with(d, points(G_st[(P_st >= P_45 & P_st <= P_60) & U == 1],
C_st[(P_st >= P_45 & P_st <= P_60) & U == 1],
col = rangi2,
pch = 16))
with(d, abline(lm(C_st[(P_st >= P_45 & P_st <= P_60)] ~
G_st[(P_st >= P_45 & P_st <= P_60)]),
lty = 2, lwd = 1, col = c("blue")))
text(-1.5, 2.1, "Good neighborhoods", col = rangi2)

m6.12 <- quap(
alist(
C ~ dnorm( mu, sigma ),
mu <- a + b_PC * P + b_GC * G + b_U * U,
a ~ dnorm(0, 1),
c(b_PC, b_GC, b_U) ~ dnorm(0, 1),
sigma ~ dexp(1)
), data = d
)

precis(m6.12)
##               mean          sd        5.5%         94.5%
## a     -0.121975097 0.071925877 -0.23692654 -0.0070236546
## b_PC   1.011611028 0.065972577  0.90617411  1.1170479480
## b_GC  -0.040813729 0.097287159 -0.19629740  0.1146699415
## b_U    1.996489923 0.147704623  1.76042941  2.2325504390
## sigma  1.019599112 0.050801756  0.93840809  1.1007901302

## 49.4 直面混雜效應

### 49.4.1 兩條通路

dag_6.1 <- dagitty("dag{
U [unobserved]
X -> Y
X <- U <- A -> C -> Y
U -> B <- C
}")
coordinates(dag_6.1) <- list( x=c(X = 0, Y = 1, U = 0,
A = 0.5, B = 0.5, C = 1),
y=c(X = 0, Y = 0, U = -1,
A = -1.3, B = -0.7, C = -1) )
drawdag(dag_6.1)

1. $$X \leftarrow U \leftarrow A \rightarrow C \rightarrow Y$$
2. $$X \leftarrow U \rightarrow B \leftarrow C \rightarrow Y$$
3. $$X \rightarrow Y$$

adjustmentSets( dag_6.1, exposure = "X", outcome = "Y")
## { C }
## { A }

### 49.4.2 華夫餅的後門

dag_6.2 <- dagitty("dag{
A -> D
A -> M -> D
A <- S -> M
S -> W -> D
}")
coordinates(dag_6.2) <- list( x=c(A = 0, D = 1, M = 0.5,
S = 0, W = 1),
y=c(A = 0, D = 0, M = -0.5,
S = -1, W = -1))
drawdag(dag_6.2)

49.16 中，

• S 是所在州是否屬於南方州；
• A 是所在州的結婚年齡中位數；
• M 是所在州的結婚率；
• W 是所在州的瓦夫餅餐廳數量；
• D 是所在州的離婚率。

• 南方州的話，結婚年齡中位數較低 $$S \rightarrow A$$
• 南方州的話，結婚率較高 $$S \rightarrow M$$
• 上述 $$S \rightarrow M$$ 的關係同時還經過 $$A$$，也就是 $$S \rightarrow A \rightarrow M$$
• 南方州的華夫餅餐廳數量較多 $$S \rightarrow W$$
• 年齡中位數 $$A$$，結婚率 $$M$$，還有華夫餅餐廳數量 $$W$$ 都直接影響離婚率，$$A\rightarrow D; M\rightarrow D; W \rightarrow D$$

1. $$W \leftarrow S \rightarrow A \rightarrow D$$
2. $$W \leftarrow S \rightarrow M \rightarrow D$$
3. $$W \leftarrow S \rightarrow A \rightarrow M \rightarrow D$$

adjustmentSets( dag_6.2, exposure = "W", outcome = "D")
## { A, M }
## { S }

impliedConditionalIndependencies( dag_6.2 )
## A _||_ W | S
## D _||_ S | A, M, W
## M _||_ W | S

1. 當控制了是否是南方州($$S$$)之後，華夫餅餐廳$$W$$ 和結婚年齡中位數 $$A$$ 之間相互獨立。
2. 當同時控制了$$A$$$$M$$$$W$$之後，$$S$$ 和離婚率之間相互獨立。
3. 當控制了 $$S$$ 之後，$$W$$ 和 結婚率 $$M$$ 之間相互獨立。