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

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

MENU

効果検証入門(08/08)

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

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

6. 回帰不連続デザイン

6.1 厳格なルール下における選択バイアス

 介入が明確なルールによって決まっている状況を考える。このとき傾向スコアがルールを満たすか否かのみでフラグが立つため、傾向スコアも\mathrm{IPW}も利用できなくなる。そのための分析手法が回帰不連続デザインである。

6.1.1 回帰不連続デザインの仕組み

 ある変数X_iの値が閾値Aを超える場合に施策を行う場合を考える。このとき施策の実行を表すフラグZ_j


\begin{aligned}
Z_j=\begin{cases}
1,&X_i\geq A,\\
0,&X_i\lt A
\end{cases}
\end{aligned}

と与えられる。このようなルールで介入が割り振られる場合、介入発生の有無はX_iのみに依存することになる。変数X_iのように介入を決定する変数を割当変数 (running variable; forcing variable)という。介入の決定に用いる閾値をカットオフという。

6.1.2 選択バイアスの確認

 介入の有無で分けた2集団それぞれで算出した期待値の差分で効果を計測すると、


\begin{aligned}
\tau_{\mathrm{naive}}&=E[Y^{(1)}|Z=1]-E[Y^{(0)}|Z=0]\\
&=\tau+E[Y^{(0)}|Z=1]-E[Y^{(0)}|Z=0]\\
&=\tau+E[Y^{(0)}|X_i\geq A]-E[Y^{(0)}|X_i\lt A]
\end{aligned}

である。このようにE[Y^{(0)}|X_i\geq A]-E[Y^{(0)}|X_i\lt A]が選択バイアスとして乗ることになる。

6.2 回帰不連続デザイン

 介入がルールにより決定される場合、データの持つバイアスは割当変数のみにより発生する。したがって割当変数を制御することで選択バイアスを軽減することを考える。

6.2.1 線形回帰による分析

 まずは回帰分析による効果測定を考える。割当変数をXとして


\begin{aligned}
E[Y^{(0)}|X]&=\beta_0+\beta_1 X,\\
Y^{(1)}&=Y^{(0)}+\tau
\end{aligned}

とするとき、介入の効果を表す母数として\rhoを導入すれば、攪乱項をuとして


\begin{aligned}
Y=\beta_0+\beta_1 X+\rho Z+u
\end{aligned}

で与えられる。もし目的変数と割当変数が線形の関係にあればこれで選択バイアスを考慮することができる。

6.2.2 非線形回帰による分析

 もし割当変数と目的変数が非線形な関係を持つならば非線形性を考慮したモデルを考えることになる。たとえば


\begin{aligned}
Y&=\beta_0+f(X)+\rho Z,\\
f(X)&=\beta_1 X+\beta_2 X^2+\cdots+\beta_p X^p
\end{aligned}

を考える。他には交差項Z\times f(X)を加えることが検討できる。

6.2.3 ノンパラメトリックRDD

 利用データを閾値の前後のみに限定することで選択バイアスを小さくする方法としてノンパラメトリック\mathrm{RDD}がある。


\begin{aligned}
\hat{\tau}_{\mathrm{naive}}=\tau+E[Y^{(0)}|X\geq A]-E[Y^{(0)}|X\lt A]
\end{aligned}

ここでの選択バイアスE[Y^{(0)}|X\geq A]-E[Y^{(0)}|X\lt A]XAに近い値を取っている場合には0に近い値を取り得ると考えられる。
 これはカットオフ近辺でデータ分布が大きく変化している場合に生じている蓋然性が高くなる。

6.3 回帰不連続デザインの仮定

6.3.1 Continuity of Conditional Regression Functions

 \mathrm{RDD}では推定した効果が正しいと考えるためにContinuity of Conditional Regression Functionsと呼ばれる仮定を満たす必要がある。すなわちE[Y^{(0)}|X],E[Y^{(1)}|X]Xに対して連続であることを仮定する。別の介入があるときに不連続となる状況が発生しやすい。

6.3.2 non-manipulation

 分析の対象が自身の介入に関するステータスを調整できないという\mathrm{non}-\mathrm{manipulation}も仮定する。

6.3.3 LATEの妥当性

 \mathrm{RDD}による分析結果は、どのような手法においても実質的にはカットオフ近辺のデータから来ます。このとき得られた効果の推定値はカットオフ周辺における効果(局所平均処置効果(LATE))である。行われた介入の効果が割当変数の値によっては変化しないという仮定を検証する方法は存在しないため、分析者によるデータの解釈で当該仮定の現実性を評価することになる。
 ただし閾値設定の妥当性を評価する場合など、カットオフ近辺の効果を知りたい場合は\mathrm{LATE}が非常に重要な概念である。

6.4 Rによる実証分析

##########################
### 回帰不連続デザイン ###
##########################

library("tidyverse")
library("broom")

# データの読み込み
email_data <- read_csv("http://www.minethatdata.com/Kevin_Hillstrom_MineThatData_E-MailAnalytics_DataMiningChallenge_2008.03.20.csv")

# ルールによるメールの配信を行ったログを作成
## データの整形と割当変数の追加
male_data <- email_data %>%
  filter(segment %in% c("Mens E-Mail","No E-Mail")) %>%
  mutate(treatment = ifelse(segment == "Mens E-Mail", 1, 0),
         history_log = log(history))

## cut-off の値を指定
threshold_value <- 5.5

## ルールによる介入を再現したデータを作成
### 割当変数がcut-offより大きければが配信されたデータのみ残す
### 逆の場合には配信されなかったデータのみ残す
### 割当変数を0.1単位で区切ったグループわけの変数も追加しておく
rdd_data <- male_data %>%
  mutate(history_log_grp = round(history_log/0.1,0)*0.1) %>%
  filter(((history_log > threshold_value) &
            (segment == "Mens E-Mail")) |
           (history_log <= threshold_value) &
           (segment == "No E-Mail"))

# RCTデータとRDDデータの傾向の比較
## 割当変数とサイト来訪率のプロット(RCTデータ)
male_data %>%
  mutate(history_log_grp = round(history_log/0.1,0)*0.1) %>%
  group_by(history_log_grp, segment) %>%
  summarise(visit = mean(visit),
            N = n()) %>%
  filter(N > 10) %>%
  ggplot(aes(y = visit,
             x = history_log_grp,
             shape = segment,
             size = N)) +
  geom_point() +
  ylim(0,NA) +
  theme_bw() +
  xlab("log(history)") +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "bottom",
        plot.margin = margin(1,2,1,1, "cm"))

## 割当変数とサイト来訪率のプロット(RDDデータ)
rdd_data %>%
  group_by(history_log_grp, segment) %>%
  summarise(visit = mean(visit),
            N = n()) %>%
  filter(N > 10) %>%
  ggplot(aes(y = visit,
             x = history_log_grp,
             shape = segment,
             size = N)) +
  geom_point() +
  geom_vline(xintercept = 5.5, linetype = 2) +
  ylim(0,NA) +
  theme_bw() +
  xlab("log(history)") +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "bottom",
        plot.margin = margin(1,2,1,1, "cm"))

# 集計による分析
## RCTデータでの比較
rct_data_table <- male_data %>%
  filter(history_log > 5, history_log < 6) %>%
  group_by(treatment) %>%
  summarise(count = n(),
            visit = mean(visit))
## RDDデータでの比較
rdd_data_table <- rdd_data %>%
  group_by(treatment) %>%
  summarise(count = n(),
            visit_rate = mean(visit))

# 回帰分析による分析
## 線形回帰による分析
rdd_lm_reg <- rdd_data %>%
  mutate(treatment = ifelse(segment == "Mens E-Mail", 1, 0)) %>%
  lm(data = ., formula = visit ~ treatment + history_log) %>%
  summary() %>%
  tidy() %>%
  filter(term == "treatment")

## 非線形回帰による分析
library("rddtools")
nonlinear_rdd_data <- rdd_data(y = rdd_data$visit,
                               x = rdd_data$history_log,
                               cutpoint = 5.5)

nonlinear_rdd_ord4 <- rdd_reg_lm(rdd_object=nonlinear_rdd_data, order=4)
nonlinear_rdd_ord4
plot(nonlinear_rdd_ord4)

# 分析に使うデータの幅と分析結果のプロット
bound_list <- 2:100/100
result_data <- data.frame()
for(bound in bound_list){
  out_data <- rdd_data %>%
    filter(between(history_log, 5.5 - bound, 5.5 + bound)) %>%
    group_by(treatment) %>%
    summarise(count = n(),
              visit_rate = mean(visit),
              sd = sd(visit))
  
  late <- out_data$visit_rate[2] - out_data$visit_rate[1]
  N <- sum(out_data$count)
  se <- sqrt(sum(out_data$visit_rate^2))/sqrt(N)
  result_data <- rbind(result_data, data.frame(late, bound, N, se))
}
result_data %>%
  ggplot(aes(y = late,
             x = bound)) +
  geom_ribbon(aes(ymax = late + 1.96*se,
                  ymin = late - 1.96*se), fill = "grey70") +
  geom_line() +
  theme_bw()

# nonparametric RDD
## ライブラリの読み込み
library("rdd")

## non-parametric RDDの実行
rdd_result <- RDestimate(data = rdd_data,
                         formula = visit ~ history_log,
                         cutpoint = 5.5)

## 結果のレポート
summary(rdd_result)

## 結果のプロット
plot(rdd_result)

## manipulat
DCdensity(runvar = rdd_data %>% pull(history_log),
          cutpoint = 5.5,
          plot = FALSE)

7. 因果推論のビジネスへの利用

 因果推論を「どのような環境」で「どのように使うか」。

7.1 因果推論を活用できる環境

 大前提として「より正しい情報がより多くのビジネス上の価値をもたらす」という条件が必要になる。

7.2 より正しい意思決定をするために

 因果推論を活用できる環境があったとして、計測方法という技術のみならずどのプロセス内でその技術を使っていくべきかを考える必要がある。

  1. 施策・介入の目的を明確にすること
  2. 手法の仮定が満たされる状況を作ること

 目的が不明瞭だと、手当たり次第の計測が必要となり労力がかかりすぎる。またチェリー・ピッキングが発生し得、さらにその効果が本質的か否かではなく効果の計測可否に注目し意味のない分析を行うことになり得る。
 因果推論の手法は実験コストが払えない場合や実験ができない状況でも利用できる点が利点である一方で、利用に当たっての仮定を満たす必要がある。そのため、その仮定を満たすことにコストが発生し得るものの、如何にそれを満たすような環境を構築できるかがポイントである。

今回のまとめ

  • ルールで介入が割り振られる場合、介入発生の有無はX_iのみに依存することになる。変数X_iのように介入を決定する変数を割当変数 (running variable; forcing variable)という。介入の決定に用いる閾値をカットオフという。
  • 介入が明確なルールによって決まっている状況での分析手法が回帰不連続デザインである。
プライバシーポリシー お問い合わせ