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

library(broom)
library(tidyverse)
library(patchwork)
library(stargazer)

なぜバイアスが問題なのか？

  • バイアスが存在すると、OLS などのパラメータ推定値が真の値からずれてしまう
  • ここでは、重回帰分析における変数選択に関連して、次の 3 つのケースを考える 　　
内容 問題点 バイアスの種類
✔ 必要な説明変数を含めない時 推定値にバイアスが生じる → 深刻 脱落変数バイアス (OVB)
✔ 入れるべきでない変数を含めた時 推定値にバイアスが生じる → 深刻 処置後変数バイアス
✔ 無関係な変数を含めた時 無益だが有害ではない
  • ここではこれら 3 つのケースを事例を交えながら解説する

1. 平均処置効果 \(ATE\) が推定できる条件

\(Y\) : 結果変数（健康状態）
\(D\) : 処置変数（0: 通院しない、1: 通院する）
\(e\) : 残差
  • 処置 \(D\) 以外に結果 \(Y\) に影響を与える変数が一切ない
    →　処置 \(D\) と残差 \(u\) の共分散がゼロ : \(Cov(D_i, e_i) = 0\)
    →　処置変数 \(D\) と残差 \(e\) が独立（＝両者は無関係）
    →　因果効果 \(\beta\) を推定できる

  • しかし、実際は「もともとの健康状態」という交絡変数（= 交絡因子） \(X\) が存在する
  • 交絡変数 (confounding variable)とは「脱落変数バイアスを生じさせるような変数」のこと
  • 「もともとの健康状態」という交絡変数 \(X\) をモデルに含めないで単回帰分析を実行すると
    → 交絡変数 \(X\) はモデルに含まれていないが、実際は処置 \(D\) と結果 \(Y\) の両方に影響を与えている
    →　交絡変数 \(X\) と 残差 \(e\) は二つまとめて、観察されていない残差 のような効果 \(u\)を発揮する
  • \(u\) の構成要素は交絡変数 \(X\) と 残差 \(e\)
  • 本来、残差は処置 \(D\) と相関してはいけない
  • しかし、\(u\) の構成要素の一つである交絡変数 \(X\) が \(D\) と相関してしまう
    →　処置 \(D\) と \(u\) の共分散がゼロではない : \(Cov(D_i, u_i) ≠ 0\)
    →　処置変数 \(D\) と \(u\) が独立（＝両者は無関係）は満たされない
    →　セレクションバイアスを引き起こす
    →　因果効果は推定できない

  • しかし、モデルの中に交絡変数 \(X\) を加えれば、セレクションバイアスを除去できる

2. 脱落変数バイアス

2.1 脱落変数バイアスとは

  • 「脱落変数バイアス」(Omitted Variable Bias: OVB) とは「必要な説明変数を含めない時のバイアス」

  • ここでは通院と健康状態の事例を使って解説する

  • 「通院するかしないか」というセレクションは「もともとの健康状態 \(X\)」によって生じると仮定

  • 「健康状態」\(Y\) を 「もともとの健康状態」\(X\) に回帰させる

  • ここで私たちが知りたいのは（平均処置効果: \(ATE\)） は \(\beta^l\) の値

  • 「もともとの健康状態」\(X\) は「結果」\(Y\) と「処置」\(D\) の両方に影響を与えている交絡変数
    →　モデルに含める必要あり

  • 「もともとの健康状態」\(X\) を含む左側のモデル＝正しいモデル

  • 「もともとの健康状態」\(X\) を含間内右側のモデル＝間違ったモデル

右側のモデル：含めるべき交絡変数 \(X\) を含めない「間違ったモデル」

  • モデル回帰式 (2) の中にある \(\alpha^s\) と \(\beta^s\) の右上にある \(s\) は「間違ったモデル」を表す
  • 便宜上、ここでは \(\gamma X_i + e_i\) をひとまとめにして \(u_i\) と表している

\[Y_i = \alpha^s + \beta^s D_i + u_i・・・(2)\]

左側のモデル：含めるべき交絡変数 \(X\) を含めた「正しいモデル」

  • モデル回帰式 (1) の中にある \(l\) は「正しいモデル」を表す

\[Y_i = \alpha^l + \beta^l D_i + \gamma^l X_i + e_i・・・(1)\]

  • 「もともとの健康状態」\(X\) を「処置」\(D\) に回帰する
  • 本来は「もともとの健康状態 \(X\)」が「通院するかどうか \(D\)」に影響すると仮定しているので、「\(D\) を \(X\) に回帰する」が正しい
  • しかし、ここでは \(X\) と \(D\) の相関を問題にしているので「\(X\) を \(D\) に回帰する」でも同じ結果が得られる
  • ここでは \(X\) を含まない式を得たいので「\(X\) を \(D\) に回帰する」を採用する

\[X_i = \alpha_0 + λD_i + ν_i・・・ (3)\]

  • モデル回帰式 \((1)\) から \(X\) を消去するため、\((1)\) に \((3)\) を代入する　　

\[Y_i = \alpha^l + \beta^l D_i + \gamma^l (\alpha_0 + λD_i + ν_i) + e_i ・・・(4)\]

\[Y_i = \alpha^l + \gamma^l\alpha_0 + (\beta^l + \gamma^lλ)D_i + e_i + \gamma^lν_i\]

  • 回帰式における \(D\) の係数 \(\beta^l + \gamma^lλ\) が「間違ったモデル」における「処置 \(D\)」の影響力
  • 以上の説明をまとめると次のようになる

\((\beta^l + \gamma^lλ)\) : 「間違ったモデル」(2) における傾き　= \(\beta ^s\)
\(\beta^l\) : 私たちが知りたい平均処置効果: ATE
\(\gamma^lλ\) : 脱落変数バイアス: OBV
\(\gamma^l\) : 交絡変数 \(X\) と結果 \(Y\) の共変関係：\(Y_i = \alpha^l + \beta^lD_i + \gamma^lX_i + e_i\)
\(λ\) : 交絡変数 \(X\) と処置 \(D\) の共変関係：\(X_i = \alpha_0 + \gamma^lD_i + ν\)
\(\alpha^l + \gamma^l\alpha_0\) :「間違ったモデル」における切片 = \(\alpha^s\)
\(e_i + \gamma^lν_i\) 「間違ったモデル」における残差 \(u_i\)
  • 「間違ったモデル」の \(D_i\) の傾き \((\beta^l + \gamma^lλ)\) には 私たちが知りたい値（平均処置効果: \(ATE\)） = \(\beta^l\)に余計なもの \(\gamma^lλ\)がついている
    →　これが脱落変数バイアス

まとめ ・ \(\gamma^l\) か \(λ\) のどちらかがゼロならバイアスは生じない
・ \(\gamma^l≠0\) かつ \(λ≠0\) のとき、\(X\) を交絡変数（共変量）と呼ぶ
・ 交絡変数とは「脱落変数バイアスを生じさせる変数」のこと
・ 交絡変数をコントロールしないと、脱落変数バイアスが生じる
→　セレクションバイアスが除去されずに残り、因果効果が正しく推定できない

2.2 脱落変数バイアスの事例

  • 「あつまれどうぶつの森」をプレイする時間と身長の関係は？
  • 「あつもり」をやると背が伸びない？

  • 理論的に考えると、ありえない！
  • しかし、回帰分析してみると・・・

データの準備

  • 学生20人の架空データ

  • 散布図を描いてみる

  • 「あつもりをやる時間」と「身長」には負の関係がある！

  • あつもりをやればやるほど背が伸びない？

  • 理論的には考えられない

  • 性別ごとに図を描き直してみる

  • 負の相関は消える！
  • 「性別」が「交絡変数」としてセレクションバイアスを生じさせていた

「何を説明したいのか？」に注意

  • 重回帰分析の目的が「予測」なのか「因果推論」なのかよく考えること
  • 分析の目的次第で、重回帰分析モデルにどのような変数を使うかが決まる
「予測」することが目的であれば

→　目的変数（説明したい変数）= 身長
→　説明変数（説明する変数）=「あつもり」をやる時間/1週間あたり
→　例えば「あつもりを週に 5 時間時間やる女性の身長は 165 cm」という予測は可能

「因果推論」することが目的であれば

→　目的変数（説明したい変数）= 「あつもり」をやる時間/1週間あたり
→　説明変数（説明する変数）= 身長
→　例えば「身長が 1 cm 高い男性は、週あたりあつもりをやる時間が 1 時間増える」という因果推論は可能
→　しかし「男性があつもりを週あたり 1 時間やると身長が 1 cm 増える」という因果推論は考えられない

  • どの変数がどの変数に「影響するか」どうかという問題は統計学や統計的因果推論では解決できない
    →　それぞれの専門分野における学問の蓄積をよく調べることが大切

  • 計量経済学や計量政治学などの社会科学では、説明したい変数（＝目的変数）を縦軸に、説明する変数（説明変数）を横軸にとるのが一般的

  • 上の例で、もし「身長」で「あつもり」を説明したいのであれば次のような散布図を描く

2.3 脱落変数バイアスの事例（総選挙データ）

  • ここでは衆院選挙データ (1996-2021) を使って脱落変数バイアス OVB を確認する

2.3.1 データの準備

・衆議院議員総選挙の得票データ hr96-21.csv をダウンロード
→ RProject folder内の data フォルダに入れる
・ dataフォルダ内から read_csv で読み取り df というデータフレーム名をつける

df <- read_csv("data/hr96-21.csv")
names(df)
 [1] "year"          "pref"          "ku"            "kun"          
 [5] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "name"          "previous"     
[13] "age"           "exp"           "status"        "vote"         
[17] "voteshare"     "eligible"      "turnout"       "seshu_dummy"  
[21] "jiban_seshu"   "nojiban_seshu"

  • exppv （有権者一人あたり選挙費用）の作成

  • exp と eligible を使って、有権者一人あたりに使う選挙費用 (exppv) を作る

  • 分析で使う数値型変数 eligible, exp, turnout のデータ型をチェック

str(df$eligible)  
 chr [1:9660] "346774" "346774" "346774" "346774" "346774" "346774" ...
str(df$exp)
 chr [1:9660] "9828097" "9311555" "9231284" "2177203" "." "." "." ...
str(df$turnout)
 chr [1:9660] "49.2" "49.2" "49.2" "49.2" "49.2" "49.2" "49.2" "51.8" ...
  • データ型が文字型 (chr) なので、数値型 (numeric) に変換する　
df$exp <- as.numeric(df$exp) 
df$eligible <- as.numeric(df$eligible) 
df$turnout <- as.numeric(df$turnout)
  • mutate()関数を使って exppv を作る
df <- df %>% 
  mutate(exppv = exp / eligible)
  • 2012年の総選挙と必要な 4 つの変数を抜き出し df_2012 と名前を付ける
df_2012 <- df %>% 
  filter(year == 2012) %>% 
  select(voteshare, previous, exppv, turnout)
  • これら変数の記述統計は次のとおり
  • R-Markdown で html 表示する際にはチャンクオプションで {r, results = “asis”} と指定する
stargazer(as.data.frame(df_2012), 
          type ="html",
          digits = 2)
Statistic N Mean St. Dev. Min Max
voteshare 1,294 23.18 17.11 0.20 84.50
previous 1,294 1.12 2.17 0 14
exppv 1,280 17.18 13.50 0.001 94.59
turnout 1,294 59.37 3.48 50.20 69.80
変数名 詳細
voteshare 得票率 (%)
previous 当選回数
exppv 有権者一人あたりに費やす選挙費用（円）
turnout 小選挙区の投票率 (%)

2.3.2 必要な説明変数を含めない時のバイアス

  • ここでは次のような二つのモデルを考えてみる
  • 「正しいモデル」と「間違ったモデル」
  • 「正しいモデル」とは必要な説明変数（「選挙費用」と「当選回数」）が含まれているモデル
    →　ここでは「得票率」は「選挙費用」と「当選回数」の両方から影響を受けているものとする
  • 「間違ったモデル」とは必要な説明変数（「当選回数」） が含まれていないモデル
  • 「当選回数」は交絡変数
    →　「当選回数」は「選挙費用」と「得票率」の両方に影響を与えている
  • 実際は「選挙費用」と「当選回数」の両方が「得票率」に影響を与えている
  • ここでは本来必要な変数が全て含まれている「正しいモデル」と交絡変数の「当選回数」が欠落している「間違ったモデル」を作り、それぞれの推定結果を比較してみる

「当選回数」が欠落したケース

  • ここでは「当選回数」が交絡変数

  • 「当選回数」は「選挙費用」と「得票率」の両方に影響を与えていると考えられる

  • 当選回数の多いベテラン議員ほど資金に余裕があるので選挙にお金が使える（はず）←　理論から導かれる

  • 交絡変数である「当選回数」をモデルに含めないとセレクションバイアスが生じる

  • 「得票率」を \(Y\) 軸に指定した散布図を確認

exppv.vs <- ggplot(df_2012, aes(exppv, voteshare)) +
  geom_point() +
  labs(x = "選挙費用 (exppv)", 
       y = "得票率 (voteshare)",
       title = "得票率と選挙費用の散布図") +
  stat_smooth(method = lm) + # (method = lm, se = FALSE) → 95% 信頼区間が消える
  theme_minimal(base_family = "HiraKakuProN-W3") 
previous.vs <- ggplot(df_2012, aes(previous, voteshare)) +
  geom_point() +
  labs(x = "当選回数 (previous)", 
       y = "得票率 (voteshare)",
       title = "得票率と当選回数の散布図") +
  stat_smooth(method = lm) + # (method = lm, se = FALSE) → 95% 信頼区間が消える
  theme_minimal(base_family = "HiraKakuProN-W3")
  • ２つのグラフを並べて表示させる
library(patchwork)
exppv.vs + previous.vs

  • 予想どおり、どちらの説明変数も「得票率」とは正の相関がある

  • 交絡変数である「当選回数」と「選挙費用」の関係を調べてみる

exppv.previous <- ggplot(df_2012, aes(previous, exppv)) +
  geom_point() +
  labs(x = "当選回数 (previous)", 
       y = "選挙費用 (exppv)",
       title = "当選回数と選挙費用の散布図") +
  stat_smooth(method = lm) + # (method = lm, se = FALSE) → 95% 信頼区間が消える
  theme_minimal(base_family = "HiraKakuProN-W3") 
exppv.previous

  • これも予想どおり、正の相関がある

  • 次の二つの回帰分析を実行してみる

  • 「正しいモデル」：選挙費用 (exppv) と当選回数 (previous) をモデルに含める Model_1

  • 「間違ったモデル」：選挙費用 (exppv) だけをモデルに含める Model_2

Model_1 <- lm(voteshare ~ exppv + previous, data = df_2012)
Model_2 <- lm(voteshare ~ exppv, data = df_2012)
  • 分析結果を stargazer()を使って表示させる
    （チャンクオプションに {r, results = "asis"} と指定すること）
stargazer(Model_1, Model_2,
                     type = "html",
                     digits = 2,
          style = "apsr",
          dep.var.caption = "得票率(%)", # 応答変数の名前
          dep.var.labels = "正しいモデル・・・・・間違ったモデル", # 応答変数の名前
          covariate.labels = c("選挙費用 (exppv)", # 変数名を変更  
                   "当選回数 (previous)"),
          title = "(1): 正しいモデル & (2): 間違ったモデル") # タイトル
(1): 正しいモデル & (2): 間違ったモデル
正しいモデル・・・・・間違ったモデル
(1) (2)
選挙費用 (exppv) 0.62*** 0.83***
(0.03) (0.03)
当選回数 (previous) 3.06***
(0.16)
Constant 9.13*** 8.99***
(0.52) (0.58)
N 1,280 1,280
R2 0.55 0.43
Adjusted R2 0.55 0.43
Residual Std. Error 11.42 (df = 1277) 12.91 (df = 1278)
F Statistic 792.11*** (df = 2; 1277) 961.34*** (df = 1; 1278)
p < .1; p < .05; p < .01
  • 正しいモデル (1) と間違ったモデル (2) の重回帰式を書いてみる

正しいモデル (1)

\[得票率= 8.81 + 0.55選挙費用 + 3.17当選回数 + e_i\]

間違ったモデル (2)

\[得票率= 8.99 + 0.83選挙費用 + u_i\]

  • 正しいモデル (1) では「選挙費用」の係数は 0.55
  • 間違ったモデル (2) では「選挙費用」の係数は 0.83

結論 ・ モデルに本来入れるべき説明変数（ここでは「当選回数」）を含めないと、脱落変数バイアスのため「選挙費用」の係数が過大に評価される
・ ここでは 0.28パーセンテージ・ポイントだけ「選挙費用」が「得票率」に与える影響を過小評価されている

解説

  • 重回帰分析とは「他の説明変数を全てコントロールした状態で」個別の変数の変化が応答変数に与える影響を推定すること
  • ここでは「当選回数」という交絡変数をモデルに加えることで、セレクションバイアスを除去した推定ができた
  • 二つの説明変数（「選挙費用」と「当選回数」）の間には正の相関関係があるず
    → 「当選回数」を重ねたベテラン議員ほど選挙で使えるお金が多い（はず？）
  • 「当選回数」と「選挙費用」の両方を含んだ正しいモデルにおける「選挙費用」の係数 0.55 は、「当選回数」が同じであると仮定した場合に、「選挙費用」だけを変えることによって「得票率」がどれだけ変化するかを測定した値
  • しかし、間違ったモデルでは、「選挙費用」が「得票率」に与える影響 (0.83) には
    「選挙費用」→ 「得票率」(0.55) だけでなく、
    「当選回数」→ 「選挙費用」という経路の影響力 (0.28)も含まれている

→　間違ったモデルでは「選挙費用」の係数を正しく推定できていない

まとめ ・必要な説明変数をもモデルに入れない　→　有害
・しかし、私たちは「世界の真の状態」を知ることはできない

解決策

→「関係していそうな説明変数は、できるだけ全てモデルに入れる」
→「コントロール変数」としてモデルに多くの変数を含める根拠
・しかし観測数 (N) が少ないときにコントロール変数を入れすぎる
→ 自由度を消費してしまうので注意が必要

3. 不必要な変数を含めた時のバイアス

  • ここでも投票率を説明するモデルを考えてみる
  • 左が「正しいモデル」で右が「間違ったモデル」だとしよう

  • ここでは次のような二つのモデルを考えてみる
  • 「正しいモデル」と「間違ったモデル」
  • 「正しいモデル」とは必要な説明変数（「選挙費用」と「当選回数」）が含まれているモデル
    →　ここでは「得票率」は「選挙費用」と「当選回数」の両方から影響を受けているものとする
  • 「間違ったモデル」とは必要な説明変数以外の変数で、応答変数に影響を与えないと思われる変数（ここでは「投票率」）が含まれるモデルと仮定
    →　「投票率」は交絡変数とは考えられない
  • ここでは本来必要な変数が全て含まれている「正しいモデル」と不必要な変数を含んだ「間違ったモデル」を作り、それぞれの推定結果を比較してみる

投票率が得票率に影響を与えないと思われる理由

・投票率が「特定の政党の候補者の得票率」に影響を与えることはあり得る
・例えば、経済的不況後の選挙では、政権与党全体の得票率は下がるはず
→　政権与党に所属する候補者の得票率は下がるず
・しかし、ここでは「候補者全員の得票率」が分析対象
・特別の理由がない限り「投票率」が候補者全員の「得票率」に影響を与えるとは考えられない

  • 分析に使うデータフレーム df_2012 の型を確認
str(df_2012)
tibble [1,294 × 4] (S3: tbl_df/tbl/data.frame)
 $ voteshare: num [1:1294] 40.7 31.8 19.3 8.2 44.8 31.9 15.2 8.1 36.7 34.9 ...
 $ previous : num [1:1294] 0 0 0 0 5 0 0 0 0 5 ...
 $ exppv    : num [1:1294] 40.76 10.35 31.75 3.27 19.81 ...
 $ turnout  : num [1:1294] 52.7 52.7 52.7 52.7 56.8 56.8 56.8 56.8 56.8 56.8 ...

  • 全て数値型 (numeric) なので問題なし

  • 「得票率」を \(Y\) 軸、「投票率」を \(X\) に指定した散布図を確認

ggplot(df_2012, aes(turnout, voteshare)) +
  geom_point() +
  labs(x = "投票率 (turnout)", 
       y = "得票率 (voteshare)",
       title = "投票率と得票率の散布図") +
  stat_smooth(method = lm) + # (method = lm, se = FALSE) → 95% 信頼区間が消える
  theme_minimal(base_family = "HiraKakuProN-W3")

  • 予想どおり、ほとんど無相関

  • 次の二つの回帰分析を実行する

  1. 正しいモデル：投票率 (turnout) をモデルに含めない Model_1
  2. 間違ったモデル：投票率 (turnout) をモデルに含めた Model_3
Model_1 <- lm(voteshare ~ exppv + previous,   
              data = df_2012)  
Model_3 <- lm(voteshare ~ exppv + previous + turnout,   
              data = df_2012)
  • 分析結果を stargazer()を使って表示させる
    （チャンクオプションに {r, results = "asis"} と指定すること）
stargazer(Model_1, Model_3,
                     type = "html",
                     digits = 2,
          style = "apsr",
          dep.var.caption = "得票率(%)", # 応答変数の名前
          dep.var.labels = "正しいモデル・・・・・間違ったモデル", # 応答変数の名前
          covariate.labels = c("選挙費用 (exppv)", # 変数名を変更  
                   "当選回数 (previous)",
                   "投票率 (turnout)"),
          title = "(1): 正しいモデル & (2): 間違ったモデル") # タイトル
(1): 正しいモデル & (2): 間違ったモデル
正しいモデル・・・・・間違ったモデル
(1) (2)
選挙費用 (exppv) 0.62*** 0.62***
(0.03) (0.03)
当選回数 (previous) 3.06*** 3.06***
(0.16) (0.16)
投票率 (turnout) -0.02
(0.09)
Constant 9.13*** 10.44*
(0.52) (5.46)
N 1,280 1,280
R2 0.55 0.55
Adjusted R2 0.55 0.55
Residual Std. Error 11.42 (df = 1277) 11.43 (df = 1276)
F Statistic 792.11*** (df = 2; 1277) 527.70*** (df = 3; 1276)
p < .1; p < .05; p < .01
  • 正しいモデル (1) と間違ったモデル (3) の重回帰式を書いてみる

正しいモデル (1)

\[得票率= 8.81 + 0.55選挙費用 + 3.17当選回数 + e_i\]

間違ったモデル (2)

\[得票率= 10.1 + 0.55選挙費用 + 3.17当選回数 - 0.02 投票率 + u_i\]

  • 正しいモデルでも間違ったモデルでも「選挙費用」の係数 (0.55) は同じ
  • 間違ったモデルにおける「投票率」の係数は −0.02 で統計的に有意ではない

結論 モデルに本来入れるべきではない説明変数（投票率）を含めても、無益であっても有害ではない

解説

  • 「投票率」と「当選回数」が関係があるとは想定できない
  • しかし、より多くの「選挙費用」を使うと有権者が選挙動員されるため「投票率」を上げることはあり得る
  • もし「投票率」が応答変数なら、このことは考慮すべき
  • しかし、ここでは「得票率」が応答変数だから「選挙費用」→「投票率」の関係は論じない

4. 処置後変数バイアス

  • 処置後変数バイアス (post treatment variable bias)
  • 処置の影響を受けた変数を含むバイアス
  • 「どの因果効果を知りたいか」が重要！
  • 「あつもり」を「性別」と「身長」に回帰する

　　

  • 「あつもり」が「身長」に与える処置効果（因果効果）を知りたいのなら（上の左側のグラフ）
    →　このモデルで問題なし

  • 交絡変数である「性別」をモデルに含めればセレクションバイアスを除去した「身長」が「あつもり」に与える処置効果 \(\gamma\) を推定できる

  • しかし「性別」が「あつもり」に与える因果効果を知りたいのなら（上の右側のグラフ）
    →　問題がある

4.1 知りたいことが「身長」→「あつもり」

  • データの準備

  • 学生20人の架空データ atumori.csv をダウンロード
    → RProject folder内の dataフォルダに入れる
    ・ data フォルダ内から read_csv で読み取り df というデータフレーム名をつけ

  • 回帰式は次のようになる

　　

  • 実際に架空データを使って重回帰分析してみる
  • gender を男性 (male) なら 1、女性 (female) なら 0 に変更
df <- mutate(df, gender = ifelse(gender == "male", 1, 0))

データの準備

  • 学生20人の架空データ atumori.csv をダウンロード
    → RProject folder内の data フォルダに入れる
    ・ dataフォルダ内から read_csv で読み取り df というデータフレーム名をつける
df <- read_csv("data/atsumori.csv")
  • 散布図を描いてみる
df %>% 
  ggplot(aes(height, atsumori)) + 
  geom_point() +
  stat_smooth(method = lm, 
              se = FALSE) + 
  labs(x = "身長 (cm)", y = "「あつもり」をやる時間/ 1 週間あたり (hour)",
       title = "「あつもり」と「身長」の関係")  +
  theme_bw(base_family = "HiraKakuProN-W3")

  • 「あつもりをやる時間」と「身長」には負の関係がある！
  • 身長が高いほどあつもりをやらない？
  • 理論的には考えられない
  • 性別ごとに図を描き直してみる
df %>% 
  ggplot(aes(height, atsumori)) + 
  geom_point(aes(color = gender)) + 
  geom_smooth(method = lm, 
              se = FALSE, 
              aes(color = gender)) +
 labs(x = "身長 (cm)", y = "「あつもり」をやる時間/ 1 週間あたり (hour)",
       title = "「あつもり」と「身長」の関係：男女別") +
  scale_color_discrete(name = "性別",
                       labels = c("女性", "男性")) +
  theme_bw(base_family = "HiraKakuProN-W3")

  • 負の相関は消える！

  • 重回帰分析をしてみる

Model_1 <- lm(atsumori ~ height + gender, data = df)
  • 分析結果を stargazer を使って示す
  • チャンクオプションに ```{r, results = "asis"} と設定すること
stargazer(Model_1,
                     type = "html",
                     digits = 2,
          style = "apsr",
          dep.var.caption = "あつもりをする時間/1週間 (hour)", # 応答変数の名前
          covariate.labels = c("身長", # 変数名を変更  
                   "性別"),
          title = "「あつもりをする時間/1週間 (hour)」と「身長」の関係") # タイトル
「あつもりをする時間/1週間 (hour)」と「身長」の関係
atsumori
身長 0.14
(1.11)
性別 -27.82***
(8.31)
Constant 11.46
(188.36)
N 20
R2 0.74
Adjusted R2 0.70
Residual Std. Error 8.75 (df = 17)
F Statistic 23.62*** (df = 2; 17)
p < .1; p < .05; p < .01
  • 回帰式は次のようになる

\[あつもり_i = 11.46 - 27.82性別 + 0.14身長_i + e_i\]

解釈 ・ 「あつもり」と「身長」の間には関係はない

4.2 知りたいことが「性別」→「あつもり」（理論）

  • 「性別」が「あつもり」に与える因果効果を知りたいのなら、右側のモデルは問題あり
  • その理由：処置後変数「身長」も「あつもり」に影響しているから

  • 「性別」だけを含めた「正しいモデル」では、「性別」は「あつもり」に \(\beta\) という 影響を与えている

  • しかし「性別」と「身長」両方の変数を入れた「間違ったモデル」では「性別」が「あつもり」に与える影響は 2 つある：\(\beta^w\) と \(\gamma^w λ\) の二つ

  • 「性別」→「身長」→「あつもり」という経路で、「性別」が「あつもり」に \(\gamma ^wλ\) という影響を与えていることに注意

  • \(β^w\) だけで判断すると「性別」が「あつもり」に与える影響を \(\gamma ^wλ\) 分だけ過小評価（もしくは過大評価）してしまう

  • 「性別」の処置効果の一部 \(\gamma^w λ\) は「身長」を通じて「あつもり」に伝わる

　　

\(\beta\) : 私たちが知りたい「性別」の因果効果（平均処置効果）: ATE
\(\beta^w\) : 間違ったモデルで得られる「性別」の回帰係数
\(\gamma^w λ\) :処置後変数バイアス
\(\gamma^w\) と \(λ\) の符号が同じ場合 : バイアスにより過小推定
\(\gamma^w\) と \(λ\) の符号が異なる場合 : バイアスにより過大推定
\(\gamma^w\) または \(λ\) がゼロの場合 : バイアスは生じない

回帰分析（予測 vs. 因果効果） ・ 回帰分析を行う際には「予測したいのか」それとも「因果効果を確認したい」のかを明らかにする
・ 機械学習のように、例えば「あつもり」を「予測」したいのであれば「性別」を含め、できるだけ多くの変数をモデルに入れてコントロールし、大きな寄与率 \(R^2\) を求める
・ 「性別」が「あつもり」に与える処置効果を確認したいのなら、処置後変数である「身長」はモデルには入れない
・ 処置効果を確認したいのなら、寄与率 \(R^2\) の値は気にしなくてもいい

4.3 知りたいことが「性別」→「あつもり」 （実践）

  • 実際のデータを使って回帰分析し、処置後変数バイアスを確認してみる

　　

　　

  • 正しいモデル（上図の左側のモデル）を確認してみる
right_model <- lm(atsumori ~ gender, data = df)

broom::tidy(right_model)
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)     35.2      2.69     13.1  1.25e-10
2 gender男性     -26.9      3.81     -7.07 1.36e- 6

結果 ・「正しいモデル」と比べると「間違ったモデル」は「性別 \(X\)」が「あつもり \(Y\)」に与える因果効果を 0.924 ポイントだけ過大評価していることがわかる

Rを使って上の作業を実行してみる　　

  • 散布図を描いてみる
fig_1 <- df %>% 
  ggplot(aes(gender, atsumori)) + 
  geom_point() +
  stat_smooth(method = lm, 
              se = FALSE) + 
  labs(x = "性別（女性=0、男性=1)", y = "あつもりする時間 /週あたり (hour)",
       title = "「あつもり」と「性別」の関係")  +
  theme_bw(base_family = "HiraKakuProN-W3")
fig_2 <- df %>% 
  ggplot(aes(height, atsumori)) + 
  geom_point() +
  stat_smooth(method = lm, 
              se = FALSE) + 
  labs(x = "身長 (cm)", y = "あつもりする時間 /週あたり (hour)",
       title = "「あつもり」と「身長」の関係")  +
  theme_bw(base_family = "HiraKakuProN-W3")
fig_1 + fig_2

  • 「性別」と「あつもり」には負の関係がある
  • 「あつもり」と「身長」には負の関係がある
  • 男性ほどあつもりをしない
  • 重回帰分析をしてみる
Model_2 <- lm(atsumori ~ height + gender, data = df)
  • 分析結果を stargazer を使って示す
  • チャンクオプションに ```{r, results = "asis"} と設定すること
stargazer(Model_2,
                     type = "html",
                     digits = 2,
          style = "apsr",
          dep.var.caption = "あつもりをする時間/1週間 (hour)", # 応答変数の名前
          covariate.labels = c("身長 (cm)", # 変数名を変更  
                   "性別"),
          title = "「あつもりをする時間/1週間 (hour)」と「身長」の関係") # タイトル
「あつもりをする時間/1週間 (hour)」と「身長」の関係
atsumori
身長 (cm) 0.14
(1.11)
性別 -27.82***
(8.31)
Constant 11.46
(188.36)
N 20
R2 0.74
Adjusted R2 0.70
Residual Std. Error 8.75 (df = 17)
F Statistic 23.62*** (df = 2; 17)
p < .1; p < .05; p < .01
  • 参考：broom() 関数を使って、簡単に結果を表示することもできる
broom::tidy(Model_2)
# A tibble: 3 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   11.5      188.      0.0609 0.952  
2 height         0.140      1.11    0.126  0.901  
3 gender男性   -27.8        8.31   -3.35   0.00380
  • 回帰式は次のとおり

\[あつもり_i = 11.5 - 27.8性別 + 0.14身長_i + e_i・・・(1)\]

  • 「身長」を「性別」に回帰してみる
Model_3 <- lm(height ~ gender, data = df)
stargazer(Model_3, 
                     type = "html",
                     digits = 2,
          style = "apsr",
          dep.var.caption = "身長",
          title = "「身長」と「性別」の関係") # タイトル
「身長」と「性別」の関係
height
gender男性 6.60***
(0.83)
Constant 169.70***
(0.59)
N 20
R2 0.78
Adjusted R2 0.77
Residual Std. Error 1.86 (df = 18)
F Statistic 63.03*** (df = 1; 18)
p < .1; p < .05; p < .01
  • 回帰式は次のとおり

\[身長_i = 169.7 + 6.6性別 + ν_i・・・(2)\]

  • (2)を(1)に代入して回帰式から「身長」を消すと

\[あつもり_i = 11.5 - 27.8性別 + 0.14身長_i + e_i\\ = 11.5 - 27.8性別 + 0.14(169.7 + 6.6性別 + ν_i) + e_i\\ = 35.22 + (-27.8 + 0.14・6.6)性別 + 0.14ν_i + e_i\\ = 35.22 -26.9 + 0.14ν_i + u_i\\ = 35.22 -26.9 + e_i\\\]

  • 「性別」が「あつもり」に与える因果効果は -26.9 だとわかる
\(\beta = -26.9\) : 私たちが知りたい「性別」の因果効果（平均処置効果）: ATE
\(\beta^w = -27.8\) : 間違ったモデルで得られる「性別」の回帰係数
\(\gamma^w λ = 0.924\) :処置後変数バイアス
\(\gamma^w\) と \(λ\) の符号が同じ場合 : バイアスにより過小推定
\(\gamma^w\) と \(λ\) の符号が異なる場合 : バイアスにより過大推定
\(\gamma^w\) または \(λ\) がゼロの場合 : バイアスは生じない

解釈 ・「性別」が「あつもり」に与える因果効果 \(\beta = -26.9\)
・処置後変数である「身長」をモデルに加えたために「性別」の回帰係数は \(\beta^w = -27.8\) と表示される
・\(\gamma^w\) と \(λ\) の符号が異なる
→　処置後変数バイアスのため「性別」が「あつもり」に与える因果効果を 0.924 ポイントだけ過大推定してしまう