「大人の教養・知識・気付き」を伸ばすブログ

一流の大人（ビジネスマン、政治家、リーダー…）として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。　→現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年1月21日改訂)。

MENU

効果検証入門（03/08）

経済学 機械学習

　計量経済学の知見をより深めるべく

効果検証入門〜正しい比較のための因果推論／計量経済学の基礎

効果検証入門〜正しい比較のための因果推論／計量経済学の基礎

Amazon

を基に因果推論と計量経済学を学んでいく。

前回

power-of-awareness.com

3.　介入効果を測るための回帰分析

　選択バイアスの影響を取り除きつつ介入効果を推測できる基本的な方法が回帰分析である。目的変数Y_iに対して

\begin{aligned} Y_i=\beta_0+\displaystyle{\sum_{i=1}^{k}\beta_i X_i}+\varepsilon_i \end{aligned}

で表す。回帰切片および回帰係数\beta_i,i=0,1,\cdots,kは母集団の特性値の推定量である。母集団においてY_iおよびその条件付期待値と回帰式の関係性は

\begin{aligned} Y_i=E[Y|X]+\varepsilon_i \end{aligned}

を満たす。

3.1　効果分析のための回帰分析

　前節で学んだように、介入効果は施策を行った場合と行わなかった場合の各結果の期待値の差分\tau=E[Y^{(1)}]-E[Y^{(0)}]で表される。効果分析のための回帰分析ではデータに存在する選択バイアスの影響を可能な限り減らし、その差分を推定する。
　効果分析のための回帰分析では

　　　 (1) 目的変数Y 介入による効果を測定する対象となる変数
　　　 (2) 介入変数Z 施策の有無を表す変数（ダミー変数）
　　　 (3) 共変量X 選択バイアスの発生源と分析者が想定する変数

と各変数を呼ぶ。*1。
　具体的には

\begin{aligned} Y=E[Y|\boldsymbol{X}]+\varepsilon&=(\beta_0,\beta_1,\cdots,\beta_k)\begin{bmatrix}1\\X_1\\\vdots\\X_k\end{bmatrix}+\beta_{k+1}Z+\varepsilon,\\ E[\varepsilon|\boldsymbol{X},Z]&=0,\\ \mathbb{Cov}[\varepsilon,\boldsymbol{X}]&=\boldsymbol{O},\\ \mathrm{Cov}[\varepsilon,Z]&=0 \end{aligned}

とモデル化する。

3.1.1　回帰分析による効果の推定

　この回帰モデルを前提とすれば介入効果は、

\begin{aligned} E[Y|\boldsymbol{X},Z=1]&=(\beta_0,\beta_1,\cdots,\beta_k)\begin{bmatrix}1\\X_1\\\vdots\\X_k\end{bmatrix}+\beta_{k+1}\cdot1,\\ E[Y|\boldsymbol{X},Z=0]&=(\beta_0,\beta_1,\cdots,\beta_k)\begin{bmatrix}1\\X_1\\\vdots\\X_k\end{bmatrix}+\beta_{k+1}\cdot0 \end{aligned}

を踏まえれば

\begin{aligned} E[Y|\boldsymbol{X},Z=1]-E[Y|\boldsymbol{X},Z=0]=&\{(\beta_0,\beta_1,\cdots,\beta_k)\begin{bmatrix}1\\X_1\\\vdots\\X_k\end{bmatrix}+\beta_{k+1}\cdot1\}\\ &-\{(\beta_0,\beta_1,\cdots,\beta_k)\begin{bmatrix}1\\X_1\\\vdots\\X_k\end{bmatrix}+\beta_{k+1}\cdot0\}\\ =&\beta_{k+1} \end{aligned}

が成り立つ。すなわち母数\beta_{k+1}の推定値が介入効果の値を示すことになり、それ以外は攪乱母数である。

3.1.2　回帰分析における有意差検定

　\beta_{k+1}の推定値が介入効果の推定値である以上、その有意性について検定を行う必要がある。古典的回帰モデルであれば、t検定を施せばよい*2。

##########################################
### 重回帰分析により介入効果を推定する ###
##########################################

### 推定実施
# バイアスのあるデータを利用
biased_reg <- lm(data = biased_data,formula = formula(spend ~ treatment + history))
summary(biased_reg) # RのVerが異なるので乱数を用いたサンプルではテキストと結果が違うのかも

summary(biased_reg)$coefficients

# β^ = 0.835711808, p-value = 7.820525e-08 なので恐らく効果あり

3.2　回帰分析におけるバイアス

　回帰分析で選択バイアスの影響を可能な限り避けて推定を行うには、共変量を正しく選択する必要がある。

　実データを用いて検討してみる。

### 回帰分析におけるバイアスを検討する

## RCTデータでの回帰とバイアスデータでの回帰を比較する
# (1) 単回帰の場合
rct_reg <- lm(data = male_df, formula = spend ~ treatment)
nonrct_reg <- lm(data = biased_data, formula = spend ~ treatment)

rct_reg_coef <- summary(rct_reg) %>% tidy()
nonrct_reg_coef <- summary(nonrct_reg) %>% tidy()

nonrct_reg_coef[2,2] - rct_reg_coef[2,2] # この差分0.105が選択バイアス

# (2) 重回帰：共変量recency, channel, historyを加える
rct_mreg <- lm(data = male_df,
               formula = formula(spend ~ treatment + recency + channel + history))
nonrct_mreg <- lm(data = biased_data,
                  formula = formula(spend ~ treatment + recency + channel + history))

rct_mreg_coef <- rct_mreg %>% tidy()
nonrct_mreg_coef <- nonrct_mreg %>% tidy()


print(paste0("RCTデータによる介入効果推定値は",round(rct_reg_coef[2,2],4),"であった。"))
print(paste0("バイアスデータによる介入効果推定値は",round(nonrct_reg_coef[2,2],4),"であった。"))
print(paste0("バイアスデータによる共変量を加えた介入効果推定値は",round(nonrct_mreg_coef[2,2],4),"であった。"))

# RCTデータによる介入効果の推定値は0.7698であった。
# バイアスデータによる介入効果の推定値は0.875であった。
# バイアスデータによる共変量を加えた介入効果の推定値は0.8058であった。
3.2.1　脱落変数バイアス

　選択バイアスの影響をより小さくするには「目的変数Yと介入変数Zに対して相関の有る変数を加えるべき」である。
　共変量の有無で相違する2つのモデル

\begin{aligned} A:&Y_i=\alpha_0+\alpha_1 Z_i+\varepsilon_i,\\ B:&Y_i=\beta_0+\beta_1 Z_i+\beta_2 X_{\mathrm{omit},i}+\eta_i \end{aligned}

を考える。なお後者は選択バイアスが除外されたものとする。
　このとき

\begin{aligned} \varepsilon_i=\beta_2 X_{\mathrm{omit},i}+\eta_i \end{aligned}

であり、\beta_2 X_{\mathrm{omit},i}が選択バイアスに相当することが分かる。このように本来は必要であるもののモデルに欠落している変数X_{\mathrm{omit},i}を脱落変数という。
　このときモデルAにおいて介入効果を示す\alpha_1は

\begin{aligned} \alpha_1=\beta_1+\gamma_1\beta_2,\gamma_1\in\mathbb{R} \end{aligned}

で表されることが知られている。\beta_1はモデルBにて推定される効果を示しており、他方でこの\gamma_1\beta_2は脱落変数バイアスと呼ばれる。
　\beta_2はモデルBにおいて推定されたX_{\mathrm{omit},i}とY_iの相関に相当するもの(感応度)であり、\gamma_1はX_{\mathrm{omit},i}にZ_iを回帰させた回帰係数であり、これらの相関に相当するもの(感応度)と言える。すなわち

\begin{aligned} X_{\mathrm{omit},i}=\gamma_1 Z_i+e_i \end{aligned}

と書ける。

3.2.2　脱落変数バイアスのシミュレーション
########################
### 脱落変数バイアス ###
########################
library(broom)

# モデル式をベクトルとして生成
formula_vc <- c(spend ~ treatment + recency + channel,
                spend ~ treatment + recency + channel + history,
                history ~ treatment + recency + channel)
names(formula_vc) <- paste("reg",LETTERS[1:3],sep = "_") # LETTERはアルファベットとして内部で定義されている

# モデル式をデータフレーム化
models <- formula_vc %>% enframe(name = "model_index", value = "formula")

# 回帰分析を一括で実行
df_models <- models %>%
  mutate(model = map(.x = formula, .f = lm, data = biased_data)) %>%
  mutate(lm_result = map(.x = model, .f = tidy))

df_results <- df_models %>%
  mutate(formula = as.character(formula)) %>%
  select(formula, model_index, lm_result) %>%
  unnest(cols = c(lm_result))

# モデルA,B,cにおけるtreatment項の回帰係数の推定値を抽出
treatment_coef <- df_results %>%
  filter(term == "treatment") %>%
  pull(estimate)
history_coef <- df_results %>%
  filter(model_index == "reg_B",
         term == "history") %>%
  pull(estimate)

# 脱落変数バイアスを計算
OVB <- history_coef * treatment_coef[3]
coef_gap <- treatment_coef[1] - treatment_coef[2]
OVB # beta_2 + gamma_1
coef_gap # alpha_1 - beta_1
# 0.01842479で両者は一致
3.2.3　脱落変数バイアスから読み取れること

　Z,Yの両方に関係のあるような変数を交絡因子といい、バイアスの少ない効果推定を行うには交絡因子を喰わることが重要である。
　脱落変数バイアスが小さいという状態には、

  • (1)Y,Zが本来関係が薄い場合、
  • (2)既にモデルに含まれている変数の影響によって小さくなっている場合

がある。バイアスを発生し得る変数X_{\mathrm{omit}}がデータとして得られないとしても常識的に考えてそれとY,Zとそれぞれ正・負のどちらの相関があり得るかを考えるだけでもバイアスの方向性を検討できる。

3.2.4　CIA

　効果検証のための回帰分析における理想的な共変量の選択はモデルに含まれていない変数による脱落変数バイアスがすべて0になるような状態である。このような状態では共変量がY^{(0)},Y^{(1)}と条件付き独立になる。これを\mathrm{CIA（Conditional\ Independent}\mathrm{ Assumption）}という。

3.2.5　変数の選び方とモデル評価

　モデルの構築方法は以下の手順で行われる：

　　　 (1) 介入の割り当てがどのように決定されるかを想定する。
　　　 (2) 想定を表現できるような共変量を選択する。
　　　 (3) 選択した共変量とYとの関係性を考慮して関数形を決定する。

　推定された効果の妥当性を主張するためには、作ったモデルの共変量が\mathrm{CIA}を満たしていると考える必要がある。しかし問題点が2点ある：

　 (1) バイアス評価が出来ない 得られた効果の推定値がどの程度バイアスを有しているかを評価する方法が無い。
　 (2) 必要な共変量がデータ内に存在しない場合がある 手持ちのデータセット内に存在する変数の実ではバイアスを十分に減らせない可能性がある。

これらについてはより高度な分析手法を用いない限り、分析者の初見で判断するしかない。
　手持データに含まれない変数が選択バイアスを起こしているかを判断する方法の1つに\mathrm{Sensitivity\ analysis}がある。分析者が重要だと認識している共変量以外の共変量をモデルから除いて推定したときに効果の推定値が大きく変動しないかを確認するというものである。変動が小さければ回帰の結果が他の変数による影響を受けにくく、仮にデータセットに含まれていないような偏巣を含めても大きく変化しないことを示唆する。

3.2.6　Post Treatment Bias

　選択バイアスを減らす可能性があるからといって共変量を増やし過ぎて介入の影響を受け得る変数をモデルに含めると、回帰分析の結果がゆがむ可能性がある。これは介入の影響を受けた変数(=介入変数と高い相関を持つ共変量)を加えることで介入の効果がその変数の推定値に吸着されてしまうことで発生する。このような介入によって影響を受けた変数を分析に加えることで発生するバイアスを\mathrm{Post\ Treatment\ Bias}という。

3.3　回帰分析における諸論点

3.3.1　予測と効果推定

　予測能力が低いために回帰分析を利用すべきでないとの議論がある。しかし「モデルのデータに対する説明能力や未知の標本に対する予測能力を高めることと“効果検証について有用である”という保証はない」ため、予測能力と効果検証は別問題である。

3.3.2　制限被説明変数

　Y=0,1やY\geq0など目的変数(被説明変数)の取る値に制約がある場合、制限被説明変数という。このときには線形回帰以外を利用することも一考の余地がある。

3.3.3　対数を用いた線形回帰

　説明変数および目的変数としてそれらの(原指標の)対数階差を用いて回帰することも考えられる。これは対数の差分が両者の変化率に近似し得ることが背景にある。介入効果が比率で表されるときや、共変量と目的変数の関係が比率で扱われる場合に対数変換は使い得る。

3.3.4　多重共線性

　回帰モデルに含まれる説明変数間に強い相関が持つとき多重共線が発生しているという。この場合、推定量の標準誤差が異常に大きくなり、推定結果が大きく歪む。実際、変数kが他のある変数と多重共線を起こしており、それらの相関をr_kとすれば

\begin{aligned} \bar{X}_k&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}x_{k,i}},\\ V[\hat{\beta}_k]&=\displaystyle{\frac{\sigma^2}{(1-r_k\displaystyle{\sum_{i=1}^{n}(X_{k,i}-\bar{X}_{k})^2})}} \end{aligned}

と書ける。r_k\rightarrow 1ならばV[\hat{\beta}_k]\rightarrow\infty出が成り立つ。このように多重共線性の問題点は母数の推定量の標準誤差が信頼できない点である。

今回のまとめ

  • 選択バイアスの影響を取り除きつつ介入効果を推測できる基本的な方法が回帰分析である。説明変数に介入効果を表すダミー変数（介入変数）を加える。
  • 介入効果の推定量は介入変数の回帰係数の推定量に等しい。
  • 選択バイアスの影響をより小さくするには「目的変数Yと介入変数Zに対して相関の有る変数を加えるべき」である。

次回

power-of-awareness.com

*1:回帰分析の文脈では介入変数及び共変量を総称して説明変数と呼ぶ。

*2:t検定についてはたとえばpower-of-awareness.comを参照。

プライバシーポリシー お問い合わせ