課題7：メールクーポン実験のITT/TOT・CUPED・peeking診断

計量経済学II

作者

Kei Ikegami

この課題のゴール

第7回の講義では、メールクーポン実験（$N=2,000$）を例に、ランダム化推論、検出力・MDE、ITT/TOT、CUPED、peekingを一通り扱った。この課題では、舞台をECサイトのメールクーポン実験（$N=10,000$）に変え、自分の手でデータを生成し、講義で学んだ推定・推論の一式をこの新しいデータに適用する。

具体的には次を確認する。

単純比較（使用者vs非使用者）がTOTを大幅に過大評価すること。
ITT推定＋t検定＋ランダム化推論p値。
TOT：配信をIVにした2SLS。
CUPED/回帰調整によるSE削減率。
peekingによるType I errorの膨張曲線。
この実験のMDE（検出できる最小効果）。

配点は全体で100点。手を動かせば2〜4時間程度で終わる分量を想定している。ヒントのcalloutを活用してほしい。

設定：ECサイトのメールクーポン実験

あなたはECサイトのマーケティング分析チームにいるとしよう。$N=10,000$人の顧客を対象に、メールクーポンの配信実験を計画している。

配信$D_i \in \{0,1\}$：対象のちょうど50%を無作為に選ぶ完全無作為化で配信する。
実験前30日購買額$X_i$（共変量）：配信とは無関係に決まっている、実験開始前の顧客の購買傾向を表す変数。実験期間のアウトカムと相関$\rho \approx 0.6$を持つように設計する。
使用$U_i \in \{0,1\}$：配信された顧客のうち、実際にクーポンを使うかどうかは顧客の「能力」（購買意欲・ネットリテラシーなど、観測されない潜在変数）に依存する自己選択。配信されなかった顧客は使いようがない（one-sided noncompliance）。
実験期間購買額$Y_i$（アウトカム）：使用したことの効果（TOT）を受けて決まる。

真値は

\[ \text{true\_TOT} = 800\text{円} \]

（クーポンを実際に使った顧客の購買額を平均800円押し上げる）とする。ITT（配信そのものの効果）はTOTと使用率の積で決まるので、別途固定する必要はない（Part 1で自然に生成される）。

Part 1：シミュレーション設計（40点）

Q1（10点）：配信・能力・実験前購買額の生成

$N=10,000$の顧客について、次の変数を生成せよ。

配信 $D_i$：処置数を$N/2$に固定し、対象のちょうど50%をランダムに割り当てる。
能力（潜在変数）$\text{ability}_i$：標準正規分布 $N(0,1)$。
実験前30日購買額 $X_i$（pre_period）：能力と相関0.7を持つように、次の手順で生成する。
1. z_pre_indep <- rnorm(N, 0, 1)（能力とは独立なノイズ）
2. pre_z <- 0.7 * ability + sqrt(1 - 0.7^2) * z_pre_indep
3. pre_period <- 8000 + 2200 * pre_z

N <- 10000

D <- sample(c(rep(1, N / 2), rep(0, N / 2)))
ability <- rnorm(N, ___, ___)

z_pre_indep <- rnorm(N, 0, 1)
pre_z <- ___ * ability + sqrt(1 - ___^2) * z_pre_indep
pre_period <- 8000 + 2200 * ___

ヒント：なぜ能力から2段階でpre_periodを作るのか

pre_periodを単にabilityの線形関数として作ると、後で使用確率もpre_periodとほぼ同じ情報を持つ変数になってしまい、CUPEDの効果とTOTの自己選択バイアスが数値的に見分けにくくなる。z_pre_indepという独立なノイズを混ぜることで、「能力と実験前購買額は関連するが、完全には一致しない」という、より現実的な設定にしている。

Q2（10点）：使用の生成（one-sided noncompliance、能力による自己選択）

配信された顧客のみ、能力に応じた確率でクーポンを使用する。

\[ \text{use\_prob}_i = \text{clip}(0.35 + 0.28 \times \text{ability}_i,\ 0.02,\ 0.95) \]

\[ U_i = D_i \times \text{Bernoulli}(\text{use\_prob}_i) \]

（配信されなかった顧客は自動的に$U_i=0$）

use_prob <- pmin(pmax(___ + ___ * ability, 0.02), 0.95)
U <- ifelse(D == 1, rbinom(N, 1, use_prob), 0)

cat("配信率:", mean(D), "\n")
cat("使用率(配信群内):", mean(U[D == 1]), "\n")

ヒント

使用率(配信群内)はだいたい35〜40%程度になっているはずである。極端に0%や100%に近い場合は、use_probの式を疑うこと。

Q3（10点）：実験期間購買額（アウトカム）の生成

実験期間購買額$Y_i$は、実験前購買額と相関$\rho=0.6$を持つ「処置前のベース水準」に、使用の効果（TOT=800円）を加えたものとする。

\[ \text{outcome\_base}_i = 8200 + 2200 \times (0.6 \times \text{pre\_z}_i + \sqrt{1-0.6^2} \times z_{\text{out},i}) \]

（$\text{pre\_z}_i$はQ1で作った変数を再利用する。$z_{\text{out},i}$は新しい独立なノイズ）

\[ Y_i = \max(\text{outcome\_base}_i,\ 0) + U_i \times 800 \]

true_tot <- 800

z_out_indep <- rnorm(N, 0, 1)
out_z <- ___ * pre_z + sqrt(1 - ___^2) * z_out_indep
outcome_base <- pmax(8200 + 2200 * out_z, 0)

Y <- outcome_base + U * ___

coupon_df <- data.frame(D = D, ability = ability, pre_period = pre_period, U = U, Y = Y)

ヒント：処置効果を除いた相関を確認する

$\text{corr}(\text{pre\_period}, Y)$をそのまま計算すると、使用$U_i$による処置効果が混ざるため、狙った0.6よりやや高めに出ることがある。処置効果を除いた相関を確認したい場合は、cor(pre_period, Y - U * true_tot)のように、処置効果を引いてから相関を取るとよい。

Q4（10点）：記述統計と可視化

配信率、使用率（配信群内）、$\text{corr}(\text{pre\_period}, Y)$を計算せよ。
pre_periodとYの散布図を描け（geom_smoothで回帰直線も重ねるとよい）。
配信群・非配信群それぞれの$Y$の平均・標準偏差を表にまとめよ。

cat("配信率:", mean(coupon_df$D), "\n")
cat("使用率(配信群内):", mean(coupon_df$U[coupon_df$D == 1]), "\n")
cat("corr(pre_period, Y):", cor(coupon_df$pre_period, coupon_df$Y), "\n")

# 散布図・記述統計表はここに追加

Part 2：推定・推論（60点）

Q5（10点）：単純比較（使用者vs非使用者）がTOTを過大評価することを確認する

配信群の中で、使用者と非使用者の平均$Y$の差を計算し、真値800円と比較せよ。過大評価の倍率も報告すること。

delivered <- coupon_df %>% filter(D == 1)
naive_diff <- mean(delivered$Y[delivered$U == 1]) - mean(delivered$Y[delivered$U == 0])

cat("単純比較(使用者-非使用者):", naive_diff, "\n")
cat("真値:", 800, "\n")
cat("過大評価の倍率:", naive_diff / 800, "\n")

ヒント：なぜ単純比較は過大評価するのか

使用するかどうかは能力に依存する自己選択である。能力の高い顧客ほどクーポンを使い、かつ能力の高い顧客はクーポンなしでも（outcome_baseが高いため）購買額が高い。この2つの効果が混ざって、単純比較はTOTを大幅に過大評価する。指定したseedではおおむね2〜3倍になるが、倍率そのものを正誤判定に使わず、コードと自己選択の向きを確認すること。

Q6（12点）：ITT推定＋t検定＋ランダム化推論p値

ITT（配信群と非配信群の$Y$の平均差）を計算せよ。
t.test(Y ~ D, data = coupon_df)でt検定を実行し、p値を報告せよ。
「全顧客で配信効果がゼロ」というFisherのsharp nullに対して、実際の完全無作為化と同じ配信人数を保ったまま処置ラベル$D$を2,000回シャッフルし、Monte Carlo p値を計算せよ。平均効果ゼロというweak nullを扱うWelchのt検定とは帰無仮説が違うことを明記した上で、数値を比較せよ。

itt_hat <- mean(coupon_df$Y[coupon_df$D == 1]) - mean(coupon_df$Y[coupon_df$D == 0])

t_result <- t.test(Y ~ D, data = coupon_df)

n_perm <- 2000
perm_diffs <- numeric(n_perm)
for (b in 1:n_perm) {
  D_shuffled <- sample(coupon_df$D)
  perm_diffs[b] <- mean(coupon_df$Y[D_shuffled == 1]) - mean(coupon_df$Y[D_shuffled == 0])
}
p_perm <- (1 + sum(abs(perm_diffs) >= abs(itt_hat))) / (n_perm + 1)

ヒント：計算時間について

$N=10,000$で2,000回のシャッフルを行うと数秒〜十数秒かかることがある。sample(coupon_df$D)でラベルだけをシャッフルすればよい。回数を減らすとMonte Carlo誤差は大きくなるため、p値の末尾が一致することではなく、検定対象の違いと結論が大きく矛盾しないかを確認する。

Q7（12点）：TOT推定：配信をIVにした2SLS

ランダム割り当て・first stage・排除制約・monotonicity・one-sided noncomplianceの各仮定を確認し、Wald推定量（この設定ではTOT = ITT ÷ 使用率）を計算せよ。
fixest::feols()で、配信$D$を操作変数として使用$U$を内生変数にした2SLSを実行し、TOTを推定せよ。Wald推定量と一致することを確認せよ。
Q5の単純比較、Q7のTOT推定値、真値800円を並べた比較表を作れ。

usage_rate <- mean(coupon_df$U[coupon_df$D == 1])
wald_tot <- itt_hat / usage_rate

iv_fit <- feols(Y ~ 1 | U ~ D, data = coupon_df, vcov = "hetero")
summary(iv_fit)

ヒント：係数の取り出し方

fixestのバージョンによって、内生変数の係数名が接頭辞付きになることがある。第3回・第7回講義ノートで使った、名前の部分一致で探すヘルパー関数を再利用するとよい。

get_coef_for <- function(model, varname) {
  cf <- coef(model)
  hit <- grepl(varname, names(cf), fixed = TRUE)
  cf[hit][1]
}

Q8（10点）：CUPED/回帰調整によるSE削減率

pre_periodを共変量にしたCUPED調整を行い、調整前後のATE（ITT）推定値とSEを比較せよ。
分散削減率とSE削減率を計算し、母集団の最適係数を使う理論式$\rho^2$の標本対応物$\text{corr}(\text{pre\_period}, Y)^2$と比較せよ。有限標本で必ず一致する値ではないことにも触れよ。
同じ検出力を達成するのに必要なサンプルサイズが何%減るかを報告せよ。

theta_hat <- cov(coupon_df$Y, coupon_df$pre_period) / var(coupon_df$pre_period)
pre_mean <- mean(coupon_df$pre_period)
coupon_df$Y_cuped <- coupon_df$Y - theta_hat * (coupon_df$pre_period - pre_mean)

# 調整前・調整後のSEを計算し、分散削減率を出す

ヒント：ANCOVAとの関係で検算する

CUPEDによる調整は、feols(Y ~ D + pre_period, data = coupon_df)という回帰（ANCOVA）の$D$の係数・SEともほぼ同じ結果になるはずである。余裕があれば、この回帰を実行してCUPEDの結果と付き合わせてみるとよい（完全に一致はしないが、非常に近い値になる）。

Q9（10点）：peekingシミュレーション（Type I error膨張曲線）

講義ノートのsimulate_peeking()関数（コピーして使ってよい）を使い、この実験のサンプルサイズ設計（$N_{\max}=10{,}000$程度を想定）のもとで、「何人ごとに検定して止めるか」を変えながらType I errorがどう変わるかを確認せよ。少なくとも4パターン（例：1回だけ、10回、50回、100回）を比較し、図または表にまとめよ。

simulate_peeking <- function(n_max, check_every, n_sims, p0 = 0.10) {
  false_positive_count <- 0
  checkpoints <- seq(check_every, n_max, by = check_every)
  for (s in 1:n_sims) {
    d_full <- rbinom(n_max, 1, 0.5)
    y_full <- rbinom(n_max, 1, p0)
    stopped_significant <- FALSE
    for (n_cur in checkpoints) {
      d_cur <- d_full[1:n_cur]; y_cur <- y_full[1:n_cur]
      n1 <- sum(d_cur == 1); n0 <- sum(d_cur == 0)
      if (n1 < 5 || n0 < 5) next
      p1_hat <- mean(y_cur[d_cur == 1]); p0_hat <- mean(y_cur[d_cur == 0])
      p_pool <- mean(y_cur)
      se <- sqrt(p_pool * (1 - p_pool) * (1 / n1 + 1 / n0))
      if (se == 0) next
      z <- (p1_hat - p0_hat) / se
      if (abs(z) > 1.96) { stopped_significant <- TRUE; break }
    }
    if (stopped_significant) false_positive_count <- false_positive_count + 1
  }
  false_positive_count / n_sims
}

# n_sims は1000程度に抑えると計算時間が現実的になる

ヒント：計算時間の目安

n_max=10000、n_sims=1000程度であれば、通常のノートPCで数十秒程度で完了する。4パターンすべてを合わせても2〜3分以内に収まるはずである。もし遅い場合はn_simsを減らしてよい（結論の方向性は変わらない）。

Q10（3点）：この実験のMDE

Q6で得たITTのSE（あるいは調整前のY全体の分散）を使い、この実験規模でのMDE（検出力80%、$\alpha=5\%$）を円単位で計算せよ。
MDEを平均購買額に対する%で表現せよ。
真のTOT（800円）×使用率と比較して、この実験がITTを検出するのに十分な規模だったかコメントせよ。

n1 <- sum(coupon_df$D == 1); n0 <- sum(coupon_df$D == 0)
se_itt <- sqrt(var(coupon_df$Y[coupon_df$D == 1]) / n1 + var(coupon_df$Y[coupon_df$D == 0]) / n0)

z_alpha <- qnorm(0.975)
z_beta <- qnorm(0.80)
mde <- (z_alpha + z_beta) * se_itt

cat("MDE:", mde, "円\n")
cat("平均購買額に対する%:", mde / mean(coupon_df$Y) * 100, "%\n")

Q11（3点）：経営含意（3行で）

あなたはECサイトのマーケティング責任者に、この実験結果を報告しなければならない。次の3点に触れながら、3行以内で報告文を書け。

このクーポン施策を「配信し続けるべきか」の判断材料としてITTとTOTのどちらを使うべきか。
単純比較（使用者vs非使用者）だけで判断するとどんな誤りを犯すか。
CUPEDやサンプルサイズ設計の観点から、次回以降の実験でどんな改善が考えられるか。

提出前のチェック

チャンクを上から順に実行して、未定義変数によるエラーが出ないか確認したか。
itt_hat、naive_diff（Q5）、wald_tot・2SLS推定値（Q7）を明確に区別して報告しているか。
Q9のpeekingシミュレーションで、覗く機会を増やすほどType I errorが全体として上昇することを確認したか（Monte Carlo誤差で隣り合う推定値が逆転することはある）。
Q10のMDEが、真のITT（true_tot * usage_rateとだいたい一致するはず）と比べてどちらが大きいか（この実験でITTを検出できる規模だったか）にコメントしているか。