はじめに
こんにちは、事業会社で働いているデータサイエンティストです。
さて、今年もアドベントカレンダーの季節がやってきましたね！
今回は弊社アドベントカレンダーの記事として、因果推論的な考え方がどこまで有効なのか、そしてどこに限界があるのかについて、私なりに整理してみたいと思います。
近年、因果推論の考え方がアカデミアだけでなく、ビジネスの現場にも広く浸透してきたのは、とても良い流れだと思います。政府や企業は、単に状況を受動的に観察するだけでなく、何らかの施策や介入を行い、その「効果」を知ろうとするようになりました。
たとえば、政府（正確には中央銀行ですが）が為替市場に介入するケースや、WEB 系企業の企画職がアプリのタブ配置を変更してユーザー行動の変化を見る、といった例は、因果推論的な発想が自然に活用されているわかりやすい場面でしょう。
一方で、私たちは常に「介入できる立場」にいるわけではありません。ここで少し立ち止まって考えてみたいと思います。
たとえば、人事担当者にとって「職歴」と「入社後の業績」の因果関係は、本当にそれほど重要でしょうか。もちろん学術的には興味深いテーマですが、現実問題として、人事担当者が中途採用候補者の職歴そのものを変更できるわけではありません。
別の例を挙げてみます。仮に、X 国の Y 国に対する外交姿勢が、「X 国で与党が A 党になると厳しくなる」という相関関係を、データから確認できたとしましょう。因果推論の観点から見れば、それが厳密な因果関係なのかを検証したくなるのは、ごく自然な発想だと思います。
しかし、ここで少し視点を変えて考えてみます。例えば、実際には A 党そのものが原因なのではなく、失業率が上昇すると、Y 国に対する世論が厳しくなる一方で、失業者に手厚い政策を掲げる A 党に票が集まりやすくなる、という構造が背後にあったとしましょう。
このような分析結果が得られたとして、それは「A 党が当選したら、X 国との外交にはより注意すべきだ」という Y 国の大臣の態度を、どの程度変えるべきでしょうか。さらに言えば、Y 国の大臣は、何らかの施策によって X 国の失業率を改善することが、本当に可能なのでしょうか。もしそうでなければ、A 党と X 国の外交政策が本当の意味での因果関係かどうかを探求することに、果たして意味あるのでしょうか？
このように考えてみると、因果関係を厳密に特定することと、実際の意思決定に役立つことは、必ずしも同じではないのではないか、という疑問が自然と浮かんできます。
改めて振り返ると、人事担当者の例も、Y 国の大臣の例も、究極的には「将来をどう予測するか」という点に関心があるだけだと言えます。しかし、その予測に用いたい変数に対して、人は常に介入できるわけではありません。時間を遡って候補者の職務経歴を変えることは物理的に不可能ですし、他国の内政に踏み込んで失業率を下げることは政治的にも許されません。こうした意味で、予測に用いる変数への介入が不可能である場面は、実は少なくないのです。
候補者自身と X 国の政策担当者・中央銀行だったら、実際の介入を実施する可能性がより高いため、おそらくこの因果関係に興味あると思います。したがって、ただの予測がしたいのか、それとも介入の効果が知りたいのかは分析の消費者の立場によって全く異なります。予測タスクに無理やり因果推論の考え方を持ち込むと、壊滅的な結果になることを示すのが、本記事の目的です。
本記事では、こうした問題意識を出発点に、「因果推論は非常に強力な道具である一方で、すべての問いに対する万能な答えではないのではないか」という点について、操作変数を予測に用いるシナリオのシミュレーションを通じて考えていきたいと思います。
参考になったのはこちらの論文です：
The Sense and Non-Sense of Holdout Sample Validation in the Presence of Endogeneity
予測タスクにおける操作変数法について
特に社会科学をバックグラウンドに持つ方であれば、計量経済学の授業で次のような説明を聞いたことがあるかもしれません。すなわち、回帰モデル
$$
Y = X \beta + \epsilon
$$
において、
$$
E(\epsilon X) \ne 0
$$
すなわち $X$ と $\epsilon$ が独立でない場合、いわゆる内生性が生じ、パラメータ $\beta$ の推定にはバイアスが含まれます。その結果、一致性という統計学的に望ましい性質も失われ、標本サイズをいくら大きくしても、$\beta$ の値を正しく推定できない、ということを学んだのではないでしょうか。
本記事で取り上げてきた例に当てはめると、例えば
$$
入社後の業績_{i} = \alpha + \beta \cdot 職務経歴_{i} + \epsilon_{i}
$$
や
$$
外交政策の厳しさ = \alpha + \beta \cdot A党が与党_{i} + \epsilon_{i}
$$
といったモデルは、それぞれ
$$
E(\epsilon \cdot 職務経歴) \ne 0
$$
および
$$
E(\epsilon \cdot A党が与党) \ne 0
$$
となっている可能性が高そうです。職務経歴は測定が不可能な本人の優秀さに相関するし、A 党が与党かどうかは、例えば外国では把握しづらい X 国の失業率と相関します。こうした状況で「ではどうすればよいのか」と考えると、多くの場合、計量経済学の教科書や因果推論の入門書で紹介されるように、誤差項 $\epsilon$ と無相関な操作変数 $Z$ を探し、
$$
E(\epsilon \cdot Z) = 0
$$
という条件のもとで、
$$
E\left(\left(入社後の業績 - (\alpha + \beta \cdot 職務経歴) \right) \cdot Z\right) = 0
$$
といったモーメント条件を用いて $\beta$ を推定する、という発想に至るのではないでしょうか。
さて、ここで一つ問いを投げかけてみます。
確かに、操作変数法によって推定された係数は、一致性という意味で統計学的に望ましい性質を備えています。では一方で、介入がない中で予測がしたいので、その操作変数推定量 $\hat{\beta_{iv}}$ を用いて行う予測と、内生性の影響を受けている可能性がある OLS 推定量 $\hat{\beta_{ols}}$ を用いた予測とでは、介入がない中で実際にはどちらの予測精度が高くなるのでしょうか？
シミュレーション
ここでは
$$
y = -0.5 \cdot x + \epsilon_{2}
$$
$$
x = 0.9 \cdot z + \epsilon_{1}
$$
$$
\begin{bmatrix}
\epsilon_1 \
\epsilon_2
\end{bmatrix}'
\sim Normal \left(
\begin{bmatrix}
0 \
0
\end{bmatrix}',
\begin{bmatrix}
1 & 0.8 \\
0.8 & 1
\end{bmatrix}
\right)
$$
という、操作変数を利用すべき典型的な状況でデータを10000件生成します。R言語で操作変数法を利用する際は ivreg パッケージがおすすめです。Python の場合はお馴染みの statsmodels がおすすめです。本記事では R 言語を用いてシミュレーションします。
ここで $y$ に対する $x$ の効果を推定しようとすると、$\epsilon_{1}$ と $\epsilon_{2}$ の間に相関があるため、$x$ の係数である本来の値 -0.5 を OLS では適切に推定できません。計量経済学の教科書に忠実に従うのであれば、このような場合には $z$ を操作変数として用い、操作変数法による推定を行うべきだ、ということになります。
set.seed(123)
# 学習データ
df <- 50000 |>
mvtnorm::rmvnorm(c(0, 0), sigma = matrix(c(1, 0.8, 0.8, 1), nrow = 2)) |>
tibble::tibble(
e = _
) |>
dplyr::mutate(
e_1 = e[,1],
e_2 = e[,2],
z = rnorm(dplyr::n()),
x = 0.9 * z + e_1,
y = - 0.5 * x + e_2
)
# 検証データ
pred_observe_df <- 50000 |>
mvtnorm::rmvnorm(c(0, 0), sigma = matrix(c(1, 0.8, 0.8, 1), nrow = 2)) |>
tibble::tibble(
e = _
) |>
dplyr::mutate(
e_1 = e[,1],
e_2 = e[,2],
z = rnorm(dplyr::n()),
x = 0.9 * z + e_1,
y = - 0.5 * x + e_2
)
では、OLS と操作変数でそれぞれモデルを推定します：
m_lm <- df |>
lm(y ~ x, data = _)
m_iv <- df |>
ivreg::ivreg(y ~ x | z, data = _)
まず OLS の推定結果を見ると
> summary(m_lm)
Call:
lm(formula = y ~ x, data = df)
Residuals:
Min 1Q Median 3Q Max
-3.0872 -0.5466 0.0000 0.5443 3.4431
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0005051 0.0036139 0.14 0.889
x -0.0598993 0.0026802 -22.35 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8081 on 49998 degrees of freedom
Multiple R-squared: 0.009891, Adjusted R-squared: 0.009871
F-statistic: 499.5 on 1 and 49998 DF, p-value: < 2.2e-16
うーん、係数の値が -0.0598993 で、正解の -0.5 からかなり離れていますね、、、
一方で操作変数推定量を確認すると：
> summary(m_iv)
Call:
ivreg::ivreg(formula = y ~ x | z, data = df)
Residuals:
Min 1Q Median 3Q Max
-4.231205 -0.677212 0.006117 0.673143 4.011416
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.002101 0.004485 0.468 0.639
x -0.500236 0.004943 -101.201 <2e-16 ***
Diagnostic tests:
df1 df2 statistic p-value
Weak instruments 1 49998 41359 <2e-16 ***
Wu-Hausman 1 49997 40346 <2e-16 ***
Sargan 0 NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.003 on 49998 degrees of freedom
Multiple R-Squared: -0.5246, Adjusted R-squared: -0.5247
Wald test: 1.024e+04 on 1 and 49998 DF, p-value: < 2.2e-16
おっ、係数の推定値が -0.500236！ちゃんと真の値にかなり近いですね。これは素直に嬉しい結果です。
となると、「よし、この推定値をそのまま予測に使おう！やっぱり因果推論は最高だ😁」と言いたくなるかもしれません。ですが、ここで少し立ち止まって考えてみましょう。
今、同じデータ生成過程で作成された pred_observe_df に対して、予測を実施して誤差率を確認しましょう：
> pred_observe_lm <- predict(m_lm, newdata = pred_observe_df)
>
> pred_observe_iv <- predict(m_iv, newdata = pred_observe_df)
>
> pred_observe_df |>
dplyr::bind_cols(
pred_observe_lm = pred_observe_lm,
pred_observe_iv = pred_observe_iv
) |>
dplyr::mutate(
ape_lm = abs((y - pred_observe_lm)/y),
ape_iv = abs((y - pred_observe_iv)/y),
uno = 1
) |>
dplyr::summarise(
mape_observe_lm = mean(ape_lm),
mape_observe_iv = mean(ape_iv),
.by = uno
)
# A tibble: 1 × 3
uno mape_observe_lm mape_observe_iv
<dbl> <dbl> <dbl>
1 1 1.50 6.33
ええええええええええええ！！！！！！操作変数で推定した係数の方が真の値に近いのに、なんで誤差率が OLS の 4 倍くらいなの！？！？！？？！？！絶対たまたまでしょう、あり得ないわ💦
はい、そんなあなたに、大規模なシミュレーションも用意したので、頻度論の分布で見ていきましょう。
大規模シミュレーション
ここでは、同じシミュレーションを 10 万回回しましょう。OLS は R 言語の関数名にちなんで lm、操作変数は iv とします。
future::plan(future::multisession, workers = 16)
simul_df <- 100000 |>
seq_len() |>
furrr::future_map(
\(i){
df <- 50000 |>
mvtnorm::rmvnorm(c(0, 0), sigma = matrix(c(1, 0.8, 0.8, 1), nrow = 2)) |>
tibble::tibble(
e = _
) |>
dplyr::mutate(
e_1 = e[,1],
e_2 = e[,2],
z = rnorm(dplyr::n()),
x = 0.9 * z + e_1,
y = - 0.5 * x + e_2
)
pred_observe_df <- 50000 |>
mvtnorm::rmvnorm(c(0, 0), sigma = matrix(c(1, 0.8, 0.8, 1), nrow = 2)) |>
tibble::tibble(
e = _
) |>
dplyr::mutate(
e_1 = e[,1],
e_2 = e[,2],
z = rnorm(dplyr::n()),
x = 0.9 * z + e_1,
y = - 0.5 * x + e_2
)
m_lm <- df |>
lm(y ~ x, data = _)
m_iv <- df |>
ivreg::ivreg(y ~ x | z, data = _)
pred_observe_lm <- predict(m_lm, newdata = pred_observe_df)
pred_observe_iv <- predict(m_iv, newdata = pred_observe_df)
pred_observe_df |>
dplyr::bind_cols(
pred_observe_lm = pred_observe_lm,
pred_observe_iv = pred_observe_iv
) |>
dplyr::mutate(
ape_lm = abs((y - pred_observe_lm)/y),
ape_iv = abs((y - pred_observe_iv)/y),
uno = i
) |>
dplyr::summarise(
mape_observe_lm = mean(ape_lm),
mape_observe_iv = mean(ape_iv),
.by = uno
) |>
dplyr::bind_cols(
beta_lm = coef(m_lm)[2],
beta_iv = coef(m_iv)[2]
)
},
.options = furrr::furrr_options(seed = 1),
.progress = TRUE
) |>
dplyr::bind_rows()
では、早速 OLS と IV の誤差率の分布を確認しましょう：
simul_df |>
tidyr::pivot_longer(!uno, names_to = "method", values_to = "value") |>
dplyr::filter(stringr::str_detect(method, "mape_observe")) |>
dplyr::filter(dplyr::between(value, 0, 15)) |>
ggplot2::ggplot() +
ggplot2::geom_density(ggplot2::aes(x = value, fill = method), alpha = 0.5) +
ggplot2::labs(
x = "MAPE"
)
IV が完全に負けました、、、、、念の為 IV が OLS に勝った回数を確認しましょう；
> simul_df |>
dplyr::mutate(
iv_win = dplyr::case_when(
mape_observe_iv < mape_observe_lm ~ 1,
TRUE ~ 0
)
) |>
dplyr::pull(iv_win) |>
sum()
[1] 0
10 万回中 0 回、負け方としてすごすぎません？
一方で、係数の分布を見ると：
simul_df |>
tidyr::pivot_longer(!uno, names_to = "method", values_to = "value") |>
dplyr::filter(stringr::str_detect(method, "beta")) |>
ggplot2::ggplot() +
ggplot2::geom_density(ggplot2::aes(x = value, fill = method), alpha = 0.5) +
ggplot2::geom_vline(xintercept = -0.5, linetype = "dashed", color = "red", alpha = 0.6, linewidth = 1.5)
操作変数の方が係数の正解値である -0.5 に近いのに、なぜ予測タスクで負けたでしょうか？
操作変数が捨てた情報
操作変数法が予測で不利になった理由は、実はそれほど複雑ではありません。もう一度、データ生成過程を見てみましょう。
$$
y = -0.5 \cdot x + \epsilon_{2}
$$
$$
x = 0.9 \cdot z + \epsilon_{1}
$$
$$
\begin{bmatrix}
\epsilon_1 \
\epsilon_2
\end{bmatrix}'
\sim Normal \left(
\begin{bmatrix}
0 \
0
\end{bmatrix}',
\begin{bmatrix}
1 & 0.8 \\
0.8 & 1
\end{bmatrix}
\right)
$$
ご覧の通り、$\epsilon_{1}$ と $\epsilon_{2}$ の間にはかなり強い相関があります。これは、$x$ が単なる「ノイズまみれの説明変数」ではなく、$\epsilon_{1}$ を通じて $y$ の誤差項 $\epsilon_{2}$ と相関した情報を含んでおり、結果として $y$ を予測する上で有用な情報を持っていることを意味しています。
ところが、操作変数法は「$\epsilon$ と相関している情報をすべて排除する」ことを目的とした推定手法です。その結果、$x$ のうち $z$ によって説明できる成分だけが残され、$\epsilon_{1}$ に由来する成分——すなわち $\epsilon_{2}$ と相関している、予測にとっては有益な情報——は意図的に捨てられてしまいます。
この点で、操作変数法は因果効果をきれいに取り出す代わりに、予測に役立つかもしれない情報も一緒に切り落としていると言えます。一方、OLS は内生性という問題を抱えつつも、結果として $\epsilon_{2}$ と相関した情報を $x$ を通じて間接的に利用しており、それが予測精度の高さにつながっています。
つまりここで観察されているのは、因果推論としては正しい推定が、予測としては必ずしも最適ではないという、ごく自然ではあるものの、直感に反しやすい現象なのです。
結論
本記事では、「操作変数法は因果効果を正しく推定できるにもかかわらず、予測タスクでは壊滅的に弱くなりうる」という、一見すると直感に反する現象を、シミュレーションを通じて確認しました。
ここで重要なのは、操作変数法が「間違っている」わけでは決してない、という点です。むしろその逆で、操作変数法は「介入したときに何が起きるか」という問いに対して、非常に誠実で、理論的にも美しい答えを与えてくれます。係数の一致性という性質が保証されていることは、因果効果を議論する上で極めて重要です。
しかし同時に、本記事が示したのは次の事実です。
「介入しない世界での予測」という目的においては、
因果推論的に「正しい」推定が、必ずしも「役に立つ」とは限らない。
操作変数法は、内生性を生む情報を「ノイズ」とみなし、それを徹底的に排除します。しかし、予測の観点から見ると、その「ノイズ」はしばしば未来を当てるための重要なシグナルでもあります。OLS が予測で勝った理由は、まさにそのシグナルを（理論的に正しくない形であっても）利用していたからでした。
ここから得られる教訓は、決して「因果推論は役に立たない」というものではありません。むしろ重要なのは、
-
自分はいま何を知りたいのか
-
その分析結果を、誰が、どのような意思決定に使うのか
-
その意思決定主体は、実際に介入できるのか
という問いを、分析の最初にきちんと立てることです。
介入が可能で、その効果を知りたいのであれば、因果推論は極めて強力な武器になります。一方で、介入が不可能で、純粋に将来を予測したいのであれば、因果推論の作法をそのまま持ち込むことは、むしろ有害ですらありえます。
因果推論と予測は、どちらが「上」でも「下」でもありません。それぞれが答えている問いが、そもそも違うのです。
本記事が、その分析は、誰のための、どんな問いに答えているのか？ を改めて考えるきっかけになれば、とても嬉しく思います。
メリークリスマス、そして良いアドベントカレンダーを🎄
最後に、私たちと一緒に働きたい方はぜひ下記のリンクもご確認ください：