---
title: "課題7：メールクーポン実験のITT/TOT・CUPED・peeking診断"
subtitle: "計量経済学II"
author: "Kei Ikegami"
lang: ja
format:
  html:
    toc: true
    toc-depth: 3
    toc-location: left
    number-sections: false
    embed-resources: true
    code-copy: true
    theme: cosmo
    fig-width: 8
    fig-height: 4.8
execute:
  warning: false
  message: false
---

```{r}
#| include: false
library(ggplot2)
library(dplyr)
library(fixest)

set.seed(20267)

theme_lecture <- function() {
  fam <- if (Sys.info()[["sysname"]] == "Darwin") "Hiragino Sans" else ""
  theme_minimal(base_size = 13, base_family = fam)
}
theme_set(theme_lecture())
knitr::opts_chunk$set(fig.align = "center", dpi = 144)
options(scipen = 8)
```

::: {.callout-important}
## この課題のゴール

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

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

1. 単純比較（使用者vs非使用者）がTOTを大幅に過大評価すること。
2. ITT推定＋t検定＋ランダム化推論p値。
3. TOT：配信をIVにした2SLS。
4. CUPED/回帰調整によるSE削減率。
5. peekingによるType I errorの膨張曲線。
6. この実験の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`

```{r}
#| label: q1-skeleton
#| eval: false
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 * ___
```

::: {.callout-tip}
## ヒント：なぜ能力から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$）

```{r}
#| label: q2-skeleton
#| eval: false
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")
```

::: {.callout-tip}
## ヒント

使用率(配信群内)はだいたい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
$$

```{r}
#| label: q3-skeleton
#| eval: false
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)
```

::: {.callout-tip}
## ヒント：処置効果を除いた相関を確認する

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

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

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

```{r}
#| label: q4-skeleton
#| eval: false
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円と比較せよ。過大評価の倍率も報告すること。

```{r}
#| label: q5-skeleton
#| eval: false
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")
```

::: {.callout-tip}
## ヒント：なぜ単純比較は過大評価するのか

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

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

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

```{r}
#| label: q6-skeleton
#| eval: false
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)
```

::: {.callout-tip}
## ヒント：計算時間について

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

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

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

```{r}
#| label: q7-skeleton
#| eval: false
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)
```

::: {.callout-tip}
## ヒント：係数の取り出し方

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

```{r}
#| label: coef-helper
#| eval: false
get_coef_for <- function(model, varname) {
  cf <- coef(model)
  hit <- grepl(varname, names(cf), fixed = TRUE)
  cf[hit][1]
}
```
:::

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

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

```{r}
#| label: q8-skeleton
#| eval: false
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を計算し、分散削減率を出す
```

::: {.callout-tip}
## ヒント：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回）を比較し、図または表にまとめよ。

```{r}
#| label: q9-skeleton
#| eval: false
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程度に抑えると計算時間が現実的になる
```

::: {.callout-tip}
## ヒント：計算時間の目安

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

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

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

```{r}
#| label: q10-skeleton
#| eval: false
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行以内で報告文を書け。

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

::: {.callout-warning}
## 提出前のチェック

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