- このセクションで使っている R packages

library(DT)
library(gapminder)
library(gghighlight)
library(ggrepel)
library(tidyverse)
library(stargazer)

1. 検討するモデル

  • ここでは「14. セレクションバイアスと Rubin の因果モデル（理論）」で扱った、通院と健康状態の例 (Angrist and Pischke 2009) を使い、セルフセレクションによるセレクションバイアスをシミュレーションによって確かめる
  • ここでは下図のモデルを考える

  • 通常の調査・観察データではあり得ないことだが、ここでは「通院すること (D)」が「健康状態 (Y)」に 0.6 の因果効果があるとわかっているとする（通常は「真の因果関係」は神様だけが知っていることだが、ここでは、我々もこのことを知っていると想定して話を進める）

  • ここでは R を使って次の 3 つの変数を人工的に作り出す（観察数 = 10,000）

変数 変数名 詳細
目的変数 Y : 健康状態 (1 = 最悪　〜　5 = 最良)
説明変数 D : 1 = 通院する、0 = 通院しない
交絡変数 h_hidden : もともとの健康状態 (1 = 最悪、2 = 悪い、3 = 普通、4 = 良い、5 = 最良)

1.1 データの準備

  • シミュレーションの対象となる全体の人数 N を決める
  • ここでは 10,000人と設定
N <- 10000

  • もともとの健康状態 (1 = 最悪、2 = 悪い、3 = 普通、4 = 良い、5 = 最良)をランダムに決める

  • これは、通院する前の「隠れた」健康状態であり、通常の「調査・観察データ」では観測されないため我々には入手できない変数である

library(tidyverse)
set.seed(2021-03-23)
df1 <- tibble(h_hidden = sample(1:5, 
                               size = N, 
                               replace = TRUE,
                               prob = c(1, 2, 3, 2, 1)))
table(df1)
df1
   1    2    3    4    5 
1052 2190 3316 2290 1152
  • もともとの健康状態ごとにそれぞれ人数が振り分けられた
  • 1 = 最悪（1052人）、2 = 悪い（2190人）、3 = 普通（3316人）、4 = 良い（2290人）、5 = 最良（1152人）が人工的に作り出された
  • 想定したように 3 = 普通（3316人）が最も多い
names(df1)
[1] "h_hidden"
  • 「もともとの健康状態」を表す変数名は h_hidden である
  • 可視化してみる
df1 %>% 
  ggplot() +
  geom_bar(aes(x = h_hidden), fill = "deepskyblue") +
  labs(x = "健康状態", y = "人数") + 
  theme_bw(base_family = "HiraKakuProN-W3") # 文字化け対策

  • ここで「セルフセレクション」の基準を設定する
  • この隠れた健康状態に基づき、各個人が通院するかどうか決めるとする
  • 各自が「通院するかどうか」はこの「隠れた健康状態」次第で決まる
  • 他の条件が等しければ (ceteris paribus)、健康状態が悪いほど通院するはず
  • 健康状態ごとに通院する確率の平均値 \(μ\) が異なると考える
  • 各個人が通院する確率は、平均値 \(μ\) の正規分布からランダムに生成されると考える
  • 話を単純にするため、標準偏差は同じだと仮定
  • 健康状態が \(s\) の人が通院する確率 \(p_s\) は、\(p_s∼Normal(θ_s,σ)\) によって決まる
  • ただし、\(0≤p_s≤1\) になるように調整する
  • 例として、健康状態ごとの \(μ_s\) を \((μ_1,μ_2,μ_3,μ_4,μ_5) = (0.8, 0.6, 0.5, 0.3, 0.2)\) としてみる
    →　健康状態が最悪の人 (\(μ_1\)) が通院する確率は 80% に設定（つまり、\(μ_1 = 0.8\) ）
    →　健康状態が最良の人 (\(μ_5\)) が通院する確率は 20% に設定（つまり、\(μ_5 = 0.2\) ）
  • 標準偏差 \(σ\) は一律で 0.1 と設定
  • R で以下を実行して、各個人の通院確率を決める
mu <- c(0.8, 0.6, 0.5, 0.3, 0.2)
sigma <- 0.1
df1 <- df1 %>% 
  dplyr::mutate(prob = rnorm(n(), 
                             mean = mu[h_hidden], 
                             sd = sigma),
         prob = case_when(
           prob > 1 ~ 1,
           prob < 0 ~ 0,
           TRUE ~ prob
         ))
DT::datatable(df1)
  • 健康状態別の通院確率の分布を図示する
hist_prob <- ggplot(df1, aes(x = prob)) + 
  geom_histogram(binwidth = 0.05, color = "black") + 
  facet_grid(rows = vars(h_hidden)) +
  labs(x = "通院する確率", y = "人数")  + 
  theme_bw(base_family = "HiraKakuProN-W3") # 文字化け対策
plot(hist_prob)

  • この確率に基づいて、ベルヌーイ試行で通院するかどうかを表す変数 Di∈{0,1} の値を決める
  • ランダムに選ばれた観察される処置の値（「通院する」集団と「通院しない」集団の人数）は次のとおり
df1 <- df1 %>% 
  mutate(D = rbinom(n(), 
                    size = 1, 
                    prob = prob))
  • 作成した値をチェックする
table(df1$D)

   0    1 
5334 4666
  • 「通院しない」集団が 5334人、「通院する」集団が 4666人割り振られた
  • 「通院する」4666人の割合を計算してみる
mean(df1$D)
[1] 0.4666

  • 10,000人中 46.66% （4666人）が「通院する」ことがわかる

  • ここで「通院が健康状態に与える平均処置効果 (ATE)」= beta を設定する

  • 単純化のため、ATE = ITE とする →　「処置効果はどの個人にとっても同じ」だと仮定する

  • 試しに、beta = 0.6 に設定してみる
    →　隠れた健康状態 1 から 5 の範囲内で「通院する」ことで 0.6 ポイント健康状態が改善されるという意味

  • 例）隠れた健康状態が 3 の人が通院すると健康状態が 3.6 になるということ　　

  • 実際には、健康状態は 1 から 5 の整数で観測される
    →　隠れた健康状態が 5 の 人が通院しても健康状態は 5 のままのはず

  • しかし、ここでは健康状態が 5.6 に改善されることを許容する

  • シミュレーションを単純化するための妥協

  • ここでは、便宜的に隠れた健康状態を「 1 から 5 の整数」と設定したが、現実的に考えると、人々の健康状態は「連続変数」のはず
    →　人々が健康状態を問われた時、実際の感覚として 5.2 だとすれば、その値に最も近い 5 を選ぶはず
    →　シミュレーションにとって、それほど大きな問題ではない

beta <- 0.6
  • この効果は、「通院した」人だけが加算されると設定して、隠れた健康状態 Y を計算する
df1 <- df1 %>% 
  mutate(Y = h_hidden + beta*D)
  • これで、必要なデータが揃った
DT::datatable(df1)

1.2 平均値の単純比較シミュレーション

  • 人工的に生成されたデータにはセレクションバイアスがないという前提で（つまり「通院する効果」(ATE) が計算できるという前提で）シミュレーションしてみる

  • 通院した人と行かなかった人の健康状態を単純に比較する

df1_D <- df1 %>% 
  mutate(D = factor(D,
                    label = c("通院しなかった", "通院した"))) %>% 
  group_by(D) %>% 
  summarize(health = mean(Y)) %>% 
  print()
# A tibble: 2 x 2
  D              health
  <fct>           <dbl>
1 通院しなかった   3.40
2 通院した         3.21
  • 「通院した」方が体調が悪いという結果が得られた
  • この結果を可視化してみる
go_not_go <- df1 %>% 
  mutate(D = factor(D,
                    label = c("通院しなかった", "通院した"))) %>% 
  group_by(D) %>% 
  summarize(health = mean(Y)) %>% 
  
  ggplot() +
  geom_bar(aes(x = D, y = health, fill = D), 
           stat = "identity", position = "dodge") + 
  labs(x = "通院の有無", y = "体調") +
  theme_minimal(base_family = "HiraKakuProN-W3")

go_not_go

  • 「通院しなかった」集団の健康状況は 3.40 で、「通院した」集団の健康状況は 3.21

得られた「通院する効果」(ATE) = 3.21 - 3.40 = − 0.19

しかし、実際の「通院する効果」(ATE) は 0.6 のはず

  • 本来の「通院する効果」(0.6) を過小評価している
  • これはシミュレーションなので、本来は計算できないはずのセレクションバイアスも計算できる
  • セレクションバイアスはフォーマルには次のように表せる
    （詳細は「15.1. セレクションバイアスと Rubin の因果モデル（理論）」を参照）

  • 求めたいセレクションバイアスは次の値である

\[E[Y_i(0)|D = 1] - E[Y_i(0)|D = 0]\]

  • \(E[Y_i(0)|D = 1]\) の値は
e1 <- df1 %>% 
  filter(D == 1) %>% 
  pull(h_hidden) %>% 
  mean() %>% 
  print()
[1] 2.606515
  • \(E[Y_i(0)|D = 0]\) の値は
e0 <- df1 %>% 
  filter(D == 0) %>% 
  pull(h_hidden) %>% 
  mean() %>% 
  print()
[1] 3.40045
  • 従って、セレクションバイアス \(E[Y_i(0)|D = 1] - E[Y_i(0)|D = 0]\) の値は
s_bias <- e1 - e0  
s_bias
[1] -0.7939347
  • セルフセレクションの影響で、負のセレクションバイアス (−0.7939347) が生じている
  • ATE とセレクションバイアスを足した値は
beta + s_bias
[1] -0.1939347
  • この値は単純比較による効果の推定値に一致する
diff(df1_D$health)
[1] -0.1939347

シミュレーションの結果✔　セレクションバイアスのあるデータを使って「通院の効果」を推定した結果、本来 0.6 のはずなのに − 0.19 と過小評価された

1.3 単回帰分析によるシミュレーション

  • 上のモデルに関して x 軸を D (「通院しない = 0」「通院する = 1」)、y 軸を Y（体調）に設定した散布図を描いてみる
df1 %>% 
  ggplot(aes(D, Y)) + 
  geom_point() +
  stat_smooth(method = lm)  + 
  geom_jitter(width = 0.04)

  • 「通院 (D)」は「健康状態 (Y)」に 0.6 の因果効果を与えているはず（なぜなら、我々がそのようなデータを意図的に作成したのだから）

  • しかし、散布図の回帰直線の傾きはマイナス
    →　なぜなら、セレクションバイアスがあるから（我々はセレクションバイアスの存在も知っている）

  • 単回帰分析を行ってみる

fit_1 <- lm(Y ~ D, data = df1)

summary(fit_1)

Call:
lm(formula = Y ~ D, data = df1)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4005 -0.6065  0.3935  0.5996  2.3935 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.40045    0.01482 229.398   <2e-16 ***
D           -0.19393    0.02170  -8.937   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.083 on 9998 degrees of freedom
Multiple R-squared:  0.007925,  Adjusted R-squared:  0.007826 
F-statistic: 79.87 on 1 and 9998 DF,  p-value: < 2.2e-16
  • 「1.2. 平均値の単純比較シミュレーション」と同様の推定結果が得られた
  • D の傾きが -.019 => 「通院する」と「健康状態」が 0.19 だけ悪くなる

得られた「通院する効果」(ATE) = − 0.19

しかし、実際の「通院する効果」(ATE) は 0.6 のはず

  • 「通院する」ことが健康状況に与える効果は有意水準 0.001 で統計的に有意 (p値 = 2e-16)
    →　統計的有意性は因果推論の役に立たない

1.4 単回帰を繰り返すシミュレーション

  • 前のセクションの単回帰分析を繰り返し実行できるように関数にまとめる
  • 返り値は、単純比較によって得られる「バイアスを含んだ因果効果」　　
sim_hospital <- function(beta = 0.6, # 「通院」の因果効果を 0.6 に設定
                         N = 1e4, # 繰り返し回数を 10000 と指定
                         mu = c(0.8, 0.6, 0.5, 0.3, 0.2),  # もともとの健康状態ごとに各個人が通院する確率を指定
                         sigma = rep(0.1, 5)) {   # sigma は、健康状態ごとに変えても良いことにする
  if (length(sigma == 1)) sigma <- rep(sigma, 5)
  
  h_hidden <- sample(1:5, size = N, replace = TRUE,
                    prob = c(1, 2, 3, 2, 1))
  prob <-  rnorm(N, mean = mu[h_hidden], sd = sigma[h_hidden])
  prob <- case_when(
    prob > 1 ~ 1,
    prob < 0 ~ 0,
    TRUE ~ prob
  )
  D <- rbinom(N, size = 1, prob = prob)
  Y <- h_hidden + beta * D
  fit <- lm(Y ~ D)
  return(coef(fit)[2])
}
  • beta = 0.6 と設定し、この関数を 1 度だけ実行してみる
sim_hospital(beta = 0.6)
         D 
-0.2280425
  • シミュレーションを1000回繰り返してみよう。
sim_1000 <- replicate(1000, sim_hospital())
  • シミュレーションの結果の初めの部分だけを表示させる
head(sim_1000)
         D          D          D          D          D          D 
-0.2079258 -0.2356270 -0.2424896 -0.1953436 -0.1960446 -0.1928219 
p_s2 <- tibble(beta = sim_1000) %>% 
  ggplot(aes(x = beta, y = after_stat(density))) +
  geom_histogram(color = "black", binwidth = 0.01) +
  labs(x= "単純比較から得られる「効果」", y = "確率密度") + 
  theme_bw(base_family = "HiraKakuProN-W3") +
  geom_vline(xintercept = 0.6, color = "tomato") # 本当の因果効果 (0.6) に縦線を引く
plot(p_s2)

  • 結果の記述統計を表示させてみる 　
summary(sim_1000)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.2702 -0.2139 -0.2001 -0.1999 -0.1865 -0.1383
  • 当初設定したように、本当の「通院の因果効果」は 0.6
  • しかし、1,000回のシミュレーションで推定された因果効果の平均値は-0.199
  • 中央値は -0.200 であり、因果効果を大幅に過小推定
  • セルフセレクションによって、セレクションバイアスが生じているのがその原因

→　セレクションバイアスは因果推論の敵

  • ここでは「通院の因果効果」をプラスと設定
  • 1000回のシミュレーションを実行
    →　「通院の因果効果」はマイナスという結果
    その理由　→　セレクションバイアスの存在
  • 「目的変数 (Y)」と「説明変数 (D)」の両方に影響を与えている「交絡変数」(h_hidden) をモデルに含めずに回帰分析した結果、実際とはかけ離れた因果効果を得た
  • セレクションバイアスが生じるような状況では、単純比較（単回帰）によって因果効果を推定することはできない

セレクションバイアスへの対処法 →　重回帰分析を行う

1.5 重回帰分析によるシミュレーション

  • このセクションで取り扱った通院と健康状況の関係を図示してみる

  • 単回帰分析 (fit_1) と重回帰分析 (fit_2) の結果を比べてみる
fit_1 <- lm(Y ~ D, data = df1)

fit_2 <- lm(Y ~ D + h_hidden, data = df1)
stargazer(fit_1, fit_2, type = "html")
Dependent variable:
Y
(1) (2)
D -0.194*** 0.600***
(0.022) (0.000)
h_hidden 1.000***
(0.000)
Constant 3.400*** 0.000***
(0.015) (0.000)
Observations 10,000 10,000
R2 0.008 1.000
Adjusted R2 0.008 1.000
Residual Std. Error 1.083 (df = 9998) 0.000 (df = 9997)
F Statistic 79.866*** (df = 1; 9998) 317,563,416,726,861,527,786,985,619,456.000*** (df = 2; 9997)
Note: p<0.1; p<0.05; p<0.01
  • それぞれの回帰式は次のとおり

単回帰分析の結果

\[Y = 3.4 - 0.194D\]

  • 単回帰分析では「通院 (D)」の因果効果を（0.6 と設定した因果効果を）誤って −0.19 と推定している

重回帰分析の結果
\[Y = 0 + 0.6D + 1h.hidden\]

  • しかし、交絡変数である「もともとの健康状態 (h_hidden」をモデルに加えると
    　→「通院 (D)」の因果効果を（0.6 と設定した因果効果を）正しく 0.6 と推定できる

2. Exercise

  • ここでは同一のサンプル結果を得られるよう特定のノイズを set.seed(2021-03-23)と設定したが、この設定をせず、「1.4 単回帰を繰り返すシミュレーション」に習ってシミュレーションを実行して可視化しなさい
  • その際、実行するたびに当初設定した「本当の通院の因果効果 (0.6)」からどの程度、ずれるか確かめなさい　　
参考文献