計量経済学応用
ランダム化比較試験（RCT）
矢内 勇生
2021-04-13（改訂：2022-04-09, 2023-10-06）
準備
必要なパッケージを読み込み、作図用の日本語フォントを設定する。
pacman::p_load(tidyverse,
patchwork)
if (.Platform$OS.type == "windows") {# Windows
library(fontregisterer)
my_font <- "Yu Gothic"
} else if (capabilities("aqua")) { # macOS
my_font <- "HiraginoSans-W3"
} else { # Unix/Linux
my_font <- "IPAexGothic"
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))このトピックで使うRコードの説明
このトピックでは新しい関数は使わない。 理解できないコードがあれば、遠慮なくSlack で質問を。
ランダム化のシミュレーション
母集団の設定
ある処置 \(D\) の平均処置効果が正であると仮定して、母集団のデータを作る。 共変量として二値の性別（社会科学におけるジェンダー [gender] は2種類より多いが、ここでは説明の単純化のために二値とする）と、年齢（20歳から65歳）を考える。
set.seed(2023-10-6)
N <- 1e6
y_1_m <- rnorm(N, mean = 10, sd = 4)
y_1_f <- rnorm(N, mean = 12, sd = 4)
y_0_m <- rnorm(N, mean = 0, sd = 2)
y_0_f <- rnorm(N, mean = 4, sd = 2)
age <- sample(20:65, size = N, replace = TRUE)
df_pop <- tibble(id = 1:N,
age = age) |>
mutate(male = rbinom(n(), size = 1, prob = 0.5),
y_1 = ifelse(male, y_1_m, y_1_f) + 0.05 * (age - mean(age)),
y_0 = ifelse(male, y_0_m, y_0_f) + 0.05 * (age - mean(age)),
ITE = y_1 - y_0)作成した母集団の特徴を確認する。
summary(df_pop)## id age male y_1
## Min. : 1 Min. :20.0 Min. :0.0000 Min. :-9.544
## 1st Qu.: 250001 1st Qu.:31.0 1st Qu.:0.0000 1st Qu.: 8.187
## Median : 500000 Median :43.0 Median :1.0000 Median :11.000
## Mean : 500000 Mean :42.5 Mean :0.5003 Mean :11.000
## 3rd Qu.: 750000 3rd Qu.:54.0 3rd Qu.:1.0000 3rd Qu.:13.815
## Max. :1000000 Max. :65.0 Max. :1.0000 Max. :30.553
## y_0 ITE
## Min. :-9.9064 Min. :-13.832
## 1st Qu.:-0.1375 1st Qu.: 5.918
## Median : 1.9958 Median : 9.007
## Mean : 1.9958 Mean : 9.005
## 3rd Qu.: 4.1277 3rd Qu.: 12.094
## Max. :14.6178 Max. : 30.756ITE（個体処置効果）の分布を確認しておこう。
pop_ite <- ggplot(df_pop, aes(x = ITE,
after_stat(density))) +
geom_histogram(color = "black") +
labs(y = "確率密度",
title = "母集団における処置効果")
plot(pop_ite)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
母集団における平均処置効果を計算しておこう。
## 平均処置効果
pop_ATE <- mean(df_pop$ITE) |>
print() ## [1] 9.004586## 男性における平均処置効果
pop_ATE_m <- df_pop |>
filter(male == 1) |>
with(mean(ITE)) |>
print()## [1] 10.0053## 女性における平均処置効果
pop_ATE_f <- df_pop |>
filter(male == 0) |>
with(mean(ITE)) |>
print()## [1] 8.002555以下では、この母集団に対して異なる方法でRCTを実施した場合にどのような結果が得られるか、シミュレーションで確かめよう。
ベルヌーイ実験
ベルヌーイ実験 (Bernoulli experiments)
をRで実行しよう。まず、シミュレーションを繰り返し行うために、実験を実行する関数を作る。標本サイズ（1回の実験における被験者数）
を sample_size で、各ベルヌーイ試行で処置に割り付ける確率を
prob で指定できるようにする。
bern_experiment <- function(population,
sample_size = 100,
prob = 0.5,
info = FALSE) {
df <- slice_sample(.data = population, n = sample_size)
df <- df |>
mutate(D = rbinom(sample_size, size = 1, prob = prob))
ATE <- with(df, mean(y_1[D == 1]) - mean(y_0[D == 0]))
if (info) cat('処置群の個体数:',
sum(df$D),
'; 統制群の個体数:',
sample_size - sum(df$D),
'\n')
return(ATE)
}この関数を使ってベルヌーイ実験を1度実行してみる。sample_size = 100、prob = 0.5
とする（既定値なので指定しなくてよい）。
bern_experiment(population = df_pop, info = TRUE)## 処置群の個体数: 49 ; 統制群の個体数: 51## [1] 10.0504もう1度やってみる。
bern_experiment(population = df_pop, info = TRUE)## 処置群の個体数: 48 ; 統制群の個体数: 52## [1] 9.557958もう1回やってみる。
bern_experiment(population = df_pop, info = TRUE)## 処置群の個体数: 46 ; 統制群の個体数: 54## [1] 9.822413当然ながら、実験結果は毎回変わる。
実験を1,000回繰り返し、結果を図示してみよう。母集団における本当の平均処置効果を赤い縦線で示す。 標本サイズは100である。
bern_100 <- replicate(1000, bern_experiment(df_pop))
hist_bern_100 <- tibble(ATE = bern_100) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数", title = "ベルヌーイ実験")
plot(hist_bern_100)
1つひとつの実験は必ずしも ATE を正確に推定しているわけではないが、平均すれば正しい効果を推定できていることがわかる。「正しい」ATEが9 であるのに対し、RCT による推定の平均値は 9 である。
処置に割り付ける確率は、0.5 でなくても良い。
bern_experiment(population = df_pop, prob = 0.8, info = TRUE)## 処置群の個体数: 82 ; 統制群の個体数: 18## [1] 8.670659bern_experiment(population = df_pop, prob = 0.2, info = TRUE)## 処置群の個体数: 21 ; 統制群の個体数: 79## [1] 9.06003bern_experiment(population = df_pop, prob = 0.1, info = TRUE)## 処置群の個体数: 14 ; 統制群の個体数: 86## [1] 11.11127処置に割り付ける確率を0.8にして、実験を1,000回繰り返し、その結果を可視化しよう。
bern_p80 <- replicate(1000, bern_experiment(df_pop, prob = 0.8))
hist_bern_p80 <- tibble(ATE = bern_p80) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数")
plot(hist_bern_p80)
標本サイズ（被験者の数）を変えるとどうなるだろうか。 標本サイズが30の場合：
bern_30 <- replicate(1000, bern_experiment(df_pop, sample_size = 30))
hist_bern_30 <- tibble(ATE = bern_30) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数")
plot(hist_bern_30)
標本サイズが50の場合：
bern_50 <- replicate(1000, bern_experiment(df_pop, sample_size = 50))
hist_bern_50 <- tibble(ATE = bern_50) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数")
plot(hist_bern_50)
標本サイズが500の場合：
bern_500 <- replicate(1000, bern_experiment(df_pop, sample_size = 500))
hist_bern_500 <- tibble(ATE = bern_500) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数")
plot(hist_bern_500)
標本サイズが小さくても、平均的には正しい因果効果を推定できている。 しかし、標本サイズが小さいと、推定のばらつき（標準誤差）が大きくなる。
sd(bern_30)## [1] 1.349286sd(bern_50)## [1] 1.066696sd(bern_100)## [1] 0.7201027sd(bern_500)## [1] 0.3334029標準誤差が大きいと、過小推定や過大推定が多くなることがわかる。実験を1回（2回、3回）しか行わないなら、標本サイズを十分な大きさにしないと、過小推定や過大推定をしてしまうかもしれない。
完全ランダム化実験
次に、完全ランダム化実験 (completely random experiments) を行う。
sample_size と、処置に割付ける個体の数
n_treated
を選べるようにする。ただし、sample_size > n_treated
である。
cr_experiment <- function(population,
sample_size = 100,
n_treated = ceiling(sample_size / 2),
info = FALSE) {
if (sample_size <= n_treated | n_treated <= 0)
stop('n_treated must be a positive integer smaller than sample_size.')
df <- slice_sample(.data = population, n = sample_size) |>
mutate(D = sample(rep(c(1, 0), c(n_treated, sample_size - n_treated)),
size = sample_size,
replace = FALSE))
ATE <- with(df, mean(y_1[D == 1]) - mean(y_0[D == 0]))
if (info) cat('処置群の個体数:',
n_treated,
'; 統制群の個体数:',
sample_size - n_treated,
'\n')
return(ATE)
}sample_size = 100、n_treatd = 50
（どちらも既定値）にして、1回実験を行ってみよう。
cr_experiment(df_pop, info = TRUE)## 処置群の個体数: 50 ; 統制群の個体数: 50## [1] 8.901491もう1度やってみる。
cr_experiment(df_pop, info = TRUE)## 処置群の個体数: 50 ; 統制群の個体数: 50## [1] 8.520713この実験を1,000回繰り返し、結果を可視化しよう。
cr_100 <- replicate(1000, cr_experiment(df_pop))
hist_cr_100 <- tibble(ATE = cr_100) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数", title = "完全ランダム化実験")
plot(hist_cr_100)
やはり、平均的にはうまく推定できているようだ。ATEの推定値の平均値は、9.01 である。
完全ランダム化実験で処置に割付ける数は、標本サイズの半分でなくても良い。
cr_experiment(df_pop, sample_size = 100, n_treated = 70, info = TRUE)## 処置群の個体数: 70 ; 統制群の個体数: 30## [1] 9.121827cr_experiment(df_pop, sample_size = 100, n_treated = 25, info = TRUE)## 処置群の個体数: 25 ; 統制群の個体数: 75## [1] 8.165093サンプルサイズ100のうち、処置に20個体を割り付ける実験を1,000回繰り返し、結果を可視化しよう。
cr_prop20 <- replicate(1000, cr_experiment(df_pop, sample_size = 100, n_treated = 20))
hist_cr_prop20 <- tibble(ATE = cr_prop20) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数")
plot(hist_cr_prop20)
標本サイズ（被験者数）の20％のみを処置に割付けても、平均的には正しいATEを推定できる。 つまり、無作為化比較試験でサンプルを2群に分けるとき、サンプルを半々に分ける必要はないことがわかる（必要がないだけで、半々に分けたほうが効率はよい）。
ブロック別ランダム化実験
続いて、ブロック別ランダム化実験 (stratified randomized experiments)
を実施してみよう。
ここでは、性別（男性ダミー変数）によってブロックを作る。
処置に割り付ける数は、各群における処置の割合で指定できるようにする。
例えば、男性ブロックでは7割を処置、女性ブロックでは4割を処置するなら、prop = c(0.7, 0.4)
と指定できるようにする。
strat_experiment <- function(population,
sample_size = 100,
prop = c(0.5, 0.5),
info = FALSE) {
if (min(prop) <= 0 | max(prop) >= 1) stop("prop must be in (0, 1).")
df <- slice_sample(.data = population, n = sample_size)
Male <- df |> filter(male == 1)
Female <- df |> filter(male == 0)
n_treated <- ceiling(c(nrow(Male), nrow(Female)) * prop)
ATE_male <- cr_experiment(Male,
sample_size = nrow(Male),
n_treated = n_treated[1])
ATE_female <- cr_experiment(Female,
sample_size = nrow(Female),
n_treated = n_treated[2])
male_prop <- mean(df$male)
ATE <- male_prop * ATE_male + (1 - male_prop) * ATE_female
if (info) {
cat('男性ブロック 処置群の個体数',
n_treated[1],
'; 統制群の個体数:',
nrow(Male) - n_treated[1],
'\n')
cat('女性ブロック 処置群の個体数',
n_treated[2],
'; 統制群の個体数:',
nrow(Female) - n_treated[2],
'\n')
}
return(ATE)
}sample_size = 100、prop = c(0.5, 0.5)
（どちらも既定値）として、1回実験してみる。
strat_experiment(df_pop, info = TRUE)## 男性ブロック 処置群の個体数 29 ; 統制群の個体数: 29
## 女性ブロック 処置群の個体数 21 ; 統制群の個体数: 21## [1] 8.904447もう1度やってみる
strat_experiment(df_pop, info = TRUE)## 男性ブロック 処置群の個体数 30 ; 統制群の個体数: 29
## 女性ブロック 処置群の個体数 21 ; 統制群の個体数: 20## [1] 9.343559この実験を1,000回繰り返し、結果を可視化しよう。
strat_100 <- replicate(1000, strat_experiment(df_pop))
hist_strat_100 <- tibble(ATE = strat_100) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数", title = "ブロック別ランダム化実験")
plot(hist_strat_100)
やはり、平均的にはうまく推定できているようだ。ATEの推定値の平均値は、8.98 である。
処置する個体の数（割合）は、ブロックごとに変えても良い。
strat_experiment(df_pop, prop = c(0.3, 0.15), info = TRUE)## 男性ブロック 処置群の個体数 15 ; 統制群の個体数: 33
## 女性ブロック 処置群の個体数 8 ; 統制群の個体数: 44## [1] 6.993871strat_experiment(df_pop, prop = c(0.8, 0.1), info = TRUE)## 男性ブロック 処置群の個体数 39 ; 統制群の個体数: 9
## 女性ブロック 処置群の個体数 6 ; 統制群の個体数: 46## [1] 7.788258この設定の実験を1,000回繰り返してみよう。
strat_b <- replicate(1000, strat_experiment(df_pop, prop = c(0.8, 0.2)))
hist_strat_b <- tibble(ATE = strat_b) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数")
plot(hist_strat_b)
ブロックごとに処置に割り付ける割合が違っても、平均的には正しい推定ができる。
ペア毎のランダム化実験
4番目の方法として、ペア毎のランダム化試験 （対ランダム化実験, paired randomized experiments）を実施しよう。 サンプルの中で、性別が同じで、年齢が近い者同士をペアにして、ペア毎に一方を処置に、他方を統制に割り付ける。 ペアを組む都合上、標本サイズは偶数に制限し、サンプル内の男女の数がそれぞれ偶数になるように調整する。
paired_experiment <- function(population, sample_size = 100) {
## sample_size が奇数の場合はエラー
if (sample_size %% 2 != 0) stop("sample_size must be an even number.")
n_male <- sample_size / 2
if (n_male %% 2 != 0) n_male <- n_male + sample(c(-1, 1), size = 1)
n_female <- sample_size - n_male
Male <- population |>
filter(male == 1) |>
slice_sample(n = n_male) |>
arrange(age) |>
mutate(D = as.vector(replicate(n() / 2, sample(0:1, size = 2))))
Female <- population |>
filter(male == 0) |>
slice_sample(n = n_female) |>
arrange(age) |>
mutate(D = as.vector(replicate(n() / 2, sample(0:1, size = 2))))
df <- bind_rows(Male, Female)
ATE <- with(df, mean(y_1[D == 1]) - mean(y_0[D == 0]))
return(ATE)
}標本サイズ sample_size = 100
（既定値）で1回実験してみる。
paired_experiment(df_pop)## [1] 9.189225もう1回やってみよう。
paired_experiment(df_pop)## [1] 8.195619この実験を1,000回繰り返し、結果を可視化しよう（処理にそれなりに時間がかかるので注意）。
pair_100 <- replicate(1000, paired_experiment(df_pop))
hist_pair_100 <- tibble(ATE = pair_100) |>
ggplot(aes(x = ATE)) +
geom_histogram(binwidth = 0.2, color = "black") +
geom_vline(xintercept = pop_ATE, color = "tomato") +
labs(y = "度数", title = "ペア毎のランダム化実験")
plot(hist_pair_100)
やはり、平均的にはうまく推定できているようだ。ATEの推定値の平均値は、9.01 である。
各ペアで処置に割り付ける確率は、0.5 でなくても良い。自分で確かめてほしい。
4つのRCT を比較する
最後に、それぞれのRCTの推定精度を比較しよう。どれも標本サイズは100で、処置の確率・割合・数が全体の半分（程度）になるようにしたものを比べる。
1,000回の実験をシミュレートしたヒストグラムを再掲する。
MIN <- min(c(bern_100, cr_100, strat_100, pair_100)) - 0.2
MAX <- max(c(bern_100, cr_100, strat_100, pair_100)) + 0.2
plot(
((hist_bern_100 + xlim(MIN, MAX) + ylim(0, 150)) |
(hist_cr_100 + xlim(MIN, MAX) + ylim(0, 150))) /
((hist_strat_100 + xlim(MIN, MAX) + ylim(0, 150)) |
(hist_pair_100 + xlim(MIN, MAX) + ylim(0, 150))))
標準誤差を比較する。
sd(bern_100)## [1] 0.7201027sd(cr_100)## [1] 0.7234095sd(strat_100)## [1] 0.6316186sd(pair_100)## [1] 0.6247247このように、処置のパタンの数が最も少ない（かつブロック[ペア]がよ り同質的な）ペア毎のランダム化実験の標準誤差が最も小さく、推定精度が高いことがわかる。
経済実験におけるランダム化について、さらに詳細な説明は List, Sadoff, and Wagner (2011) を参照。また、RCT全般については、実験デザイン（小谷浩示先生）と実験経済学（林良平先生）の授業で詳しく扱うはずなので、それらの授業まだ受講していない人はぜひ受講してほしい。
トピック3の課題
- 上で行ったシミュレーションを、条件を変えて試し、どのような場合に実験がうまくいき、どのような条件だと失敗しやすいか研究しなさい。また、面白い結果や納得できない結果が出たら、Slack でシェアし、受講生同士で議論しなさい。
今回の課題は提出不要だが、必ず自分で実行してほしい。