---
title: "課題11：動画サブスクの解約引き留めクーポンとuplift targeting"
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(tidyr)
library(fixest)
library(grf)

set.seed(20261)

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)
```

## この課題について

第11回の講義では、通信キャリアの解約引き留めクーポンRCTを題材に、「解約リスクが高い人ほど施策が効く」という直感が裏切られうること（Ascarza 2018 の教訓）、CATEの推定（T-learner、causal forest）、Qini曲線によるターゲティング評価、IPWによる政策価値評価を学んだ。

この課題では、**動画配信サブスクリプションの解約引き留めクーポンRCT**という別の設定で、同じ一連の技術を自分の手で実装する。

- **Part 1**：自分でデータ生成過程（DGP）を書き、20,000人分の顧客RCTデータをシミュレーションで作る。「解約リスク」と「クーポンの効果（uplift）」が非単調な関係を持つように設計する。
- **Part 2**：(i) ATEの推定、(ii) T-learnerによるCATE推定、(iii) causal forestによるCATE推定と2手法の比較、(iv) churnスコアtargeting vs upliftターゲティングのQini曲線と予算25%時点の比較、(v) 3つの候補ポリシーのIPW政策価値評価、(vi) 経営含意のとりまとめ、を行う。

配点は合計100点。手を動かせば2〜4時間程度で終わる分量である。各設問にヒントcalloutを用意しているので、詰まったら遠慮なく開いてほしい。

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

- 「解約リスクの高さ」と「施策への感応度（uplift）」が別の量であることを、自分で作ったDGPの中で確認できるようになる。
- T-learner（`feols`ベース）とcausal forest（`grf::causal_forest`）の両方でCATEを推定し、真値との対応を比較できるようになる。
- Qini曲線を自分の手で実装し、churnスコア順のターゲティングがランダム以下になりうることを確認する。
- IPW政策価値推定量を実装し、複数の候補ポリシーを信頼区間つきで比較できるようになる。
:::

## Part 1：RCTデータのシミュレーション設計（配点35点）

### モデルの設定

動画配信サービスの契約者20,000人を対象に、解約引き留めクーポン（月額の一部割引）を50%の確率でランダムに配布するRCTを実施したとする。共変量 $X$ は次の5つとする。

- `usage_hours`：週あたりの視聴時間
- `tenure_months`：契約継続月数
- `billing_trouble`：過去の請求トラブル（クレーム・問い合わせ）の有無（0/1）
- `satisfaction_score`：満足度アンケートに基づく潜在スコア（連続変数、平均0）
- `device_count`：登録デバイス数

クーポンなしでの解約確率は、次のロジスティックモデルに従うとする。

$$
\text{logit}\, P(\text{churn}=1 \mid X, D=0) = \beta_0 - \beta_1 \cdot \text{usage\_hours} - \beta_2 \cdot \text{tenure\_months} + \beta_3 \cdot \text{billing\_trouble} - \beta_4 \cdot \text{satisfaction\_score} + \beta_5 \cdot \text{device\_count}
$$

クーポンの真の効果（CATE）$\tau(x)$ は、解約確率を引き下げる量として定義し、次のような**非単調**な構造を持たせる。

- 解約リスクが中位の層で効果が最大になる（山型のボーナス項）。
- 過去に請求トラブルがあった層では平均効果が**負**になりやすい。これは層の全員がsleeping dogsであることを意味しない。
- 解約リスクが非常に高い層では効果がほぼ消える。

::: {.callout-tip}
## ヒント：DGPを書く順番

1. 5つの共変量を生成する（`runif`, `rgamma`, `rbinom`, `rnorm`, `sample`などを使い分ける）。
2. `billing_trouble`は、他の共変量（例えば`satisfaction_score`）とゆるく相関する形で生成すると、後の分析で「解約リスクとuplift」の非単調性が見えやすくなる。
3. クーポンなしの解約確率 `p_churn0` をロジスティックモデルで計算する。
4. まず効果index `tau_index` を、「中位リスクで山、`billing_trouble`で負、高リスクでペナルティ」という3項の和として作る。
5. `p_churn0_clipped` と `p_churn1 = clip(p_churn0_clipped - tau_index)` を作り、**クリップ後の実際のCATE**を `tau_true = p_churn0_clipped - p_churn1` と定義する。
6. 両潜在アウトカムを`rbinom()`で作り、`principal_stratum`を4類型に分ける。その後でランダム処置`D`に応じた観測値を作る。
7. 観測を60%の`train_df`と40%の`test_df`にランダム分割する。後者はモデル学習やルール選択に使わない。
:::

### 設問1-1（15点）：パラメータとDGPコードの実装

以下のパラメータを使って、20,000人分のデータを生成せよ。

- ベース解約確率のロジット切片：$\beta_0 = -0.9$
- `usage_hours`の係数：$-0.06$（`usage_hours ~ rgamma(shape=2.2, scale=2.5)`）
- `tenure_months`の係数：$-0.015$（`tenure_months ~ runif(1, 48)`）
- `billing_trouble`の係数：$+0.85$
- `satisfaction_score`の係数：$-0.30$（`satisfaction_score ~ rnorm(0, 1)`）
- `device_count`の係数：$+0.05$（`device_count`は1〜5の整数、`sample(1:5, ..., replace=TRUE)`）
- `billing_trouble`の生成：`satisfaction_score`が低いほど、`usage_hours`が短いほど発生しやすいロジスティックモデル（切片$-1.5$、`satisfaction_score`の係数$+0.3$、`usage_hours`の係数$-0.02$）で確率を作り、`pmin(pmax(p, 0.02), 0.7)`でクリップしてから`rbinom()`
- 真のCATEの山の中心：解約確率0.22、幅0.08、山の高さ0.14
- `billing_trouble`による効果indexのシフト：$-0.20$
- 高リスク層（`p_churn0 > 0.45`）でのペナルティ：$-0.04$

```{r}
#| eval: false
n <- 20000

usage_hours <- ____________________
tenure_months <- ____________________
device_count <- ____________________
satisfaction_score <- ____________________

# billing_troubleの生成（satisfaction_score, usage_hoursに依存）
logit_billing <- ____________________
p_billing <- 1 / (1 + exp(-logit_billing))
billing_trouble <- rbinom(n, 1, pmin(pmax(p_billing, 0.02), 0.7))

# クーポンなしの解約確率
logit_churn0 <- ____________________
p_churn0 <- 1 / (1 + exp(-logit_churn0))
risk_score <- pmin(pmax(p_churn0 + rnorm(n, 0, 0.02), 0.001), 0.999)

# クーポンのランダム割り当て
D <- rbinom(n, 1, 0.5)

# CATEを作る効果index
mid_risk_bump <- ____________________
tau_index <- ____________________

# 潜在アウトカムと観測値の生成
p_churn0_clipped <- pmin(pmax(p_churn0, 0.001), 0.999)
p_churn1 <- pmin(pmax(p_churn0_clipped - tau_index, 0.001), 0.999)
tau_true <- ____________________

churn_potential_0 <- rbinom(n, 1, p_churn0_clipped)
churn_potential_1 <- rbinom(n, 1, p_churn1)
stay_potential_0 <- 1 - churn_potential_0
stay_potential_1 <- 1 - churn_potential_1
principal_stratum <- case_when(
  stay_potential_0 == 0 & stay_potential_1 == 1 ~ "persuadable",
  stay_potential_0 == 1 & stay_potential_1 == 1 ~ "sure thing",
  stay_potential_0 == 0 & stay_potential_1 == 0 ~ "lost cause",
  stay_potential_0 == 1 & stay_potential_1 == 0 ~ "sleeping dog"
)
churn_obs <- ifelse(D == 1, churn_potential_1, churn_potential_0)
stay_obs <- 1 - churn_obs

set.seed(202611)
train_idx <- sample(seq_len(n), floor(0.60 * n))
sample_role <- ifelse(seq_len(n) %in% train_idx, "train", "test")

sub_df <- data.frame(
  usage_hours = usage_hours, tenure_months = tenure_months,
  billing_trouble = billing_trouble, satisfaction_score = satisfaction_score,
  device_count = device_count, risk_score = risk_score,
  D = D, churn = churn_obs, stay = stay_obs,
  stay_potential_0 = stay_potential_0, stay_potential_1 = stay_potential_1,
  principal_stratum = principal_stratum, tau_true = tau_true,
  p_churn0 = p_churn0_clipped, sample_role = sample_role
)
train_df <- sub_df %>% filter(sample_role == "train")
test_df <- sub_df %>% filter(sample_role == "test")
head(sub_df)
```

### 設問1-2（10点）：非単調性の可視化

`risk_score`（横軸）と`tau_true`（縦軸）の散布図を、`geom_smooth()`（平滑化線）付きで描け。また、`cor(risk_score, tau_true)`を計算し、値の符号と大きさについてコメントせよ。

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

講義ノートの該当図とほぼ同じコードで描ける。負の相関になっていれば、「解約リスクが高いほどクーポンが効く」という直感が、このDGPでは成り立っていないことの数値的な確認になる。
:::

### 設問1-3（10点）：記述統計

以下を表にまとめよ。

- 5つの共変量それぞれの平均・標準偏差。
- 処置群・対照群それぞれの解約率（`churn`の平均）。
- `billing_trouble`が1の人と0の人それぞれの、平均`tau_true`。
- `tau_true < 0`と`tau_true >= 0`の各層で`principal_stratum`の構成比を計算し、「負のCATE」と「sleeping dogという個人類型」がなぜ別物かを1〜2行で説明せよ。

## Part 2：推定・推論（配点65点）

### 設問2-1（10点）：ATEの推定

RCTデータなので、単純な群間差でATEが不偏に推定できる。ここでは後のcausal forestと対象標本をそろえるため、`train_df`だけを使う。`stay`（継続したかどうか）について処置群と対照群の平均の差を計算し、その標準誤差・95%信頼区間も求めよ（`t.test()`を使ってよい）。`mean(train_df$tau_true)`と比較せよ。`test_df`はターゲティング方針の最終評価まで触らない。

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

`t.test(stay ~ D, data = train_df)`で群間差の検定と信頼区間が一度に得られる。ただし符号（処置群−対照群）に注意すること。
:::

### 設問2-2（10点）：T-learnerによるCATE推定

`train_df`だけを使い、処置群・対照群それぞれで`stay ~ usage_hours + tenure_months + billing_trouble + satisfaction_score + device_count`を`feols()`で推定せよ。その2つのモデルを`test_df`に適用し、CATE推定値 `tau_hat_tlearner` を作れ。`test_df`で`cor(tau_hat_tlearner, tau_true)`を計算せよ。

::: {.callout-tip}
## ヒント：T-learnerの実装手順

1. `mu1_fit <- feols(stay ~ ..., data = train_df %>% filter(D == 1))`
2. `mu0_fit <- feols(stay ~ ..., data = train_df %>% filter(D == 0))`
3. `predict(mu1_fit, newdata = test_df)`と`predict(mu0_fit, newdata = test_df)`の差を取る。最終評価標本は学習に使わない。
:::

### 設問2-3（15点）：causal forestによるCATE推定と比較

`train_df`で`grf::causal_forest(X, Y, W)`を推定し、`predict(cf, X_test)$predictions`で`test_df`のCATE推定値を作れ。RCTの既知の割当確率0.5を`W.hat` に渡すこと。

(a) `cor(tau_hat_cf, tau_true)`を計算し、T-learnerの相関係数と比較する表を作れ。
(b) `tau_true`（横軸）と`tau_hat_tlearner`・`tau_hat_cf`（縦軸、facetで2パネル）の散布図を、45度線付きで描け。
(c) `average_treatment_effect(cf)`でATEを推定し、設問2-1で求めた`train_df`内の単純な群間差、`mean(train_df$tau_true)`と並べた比較表を作れ。3行がすべて同じtrain標本を対象にしていることを表のラベルにも明記せよ。

```{r}
#| eval: false
X_train <- as.matrix(train_df %>% select(____________________))
X_test <- as.matrix(test_df %>% select(____________________))
Y_train <- ____________________
W_train <- ____________________

cf <- causal_forest(X_train, Y_train, W_train,
                    W.hat = rep(0.5, nrow(train_df)))
test_df$tau_hat_cf <- predict(cf, X_test)$predictions
```

::: {.callout-tip}
## ヒント：grfのAPIは最小限に

`causal_forest()`と`predict()`、`average_treatment_effect()`以外のgrfの機能（tuning、variable importance、best linear projectionなど）は今回は使わなくてよい。この3つの関数だけで十分に今回の課題は完結する。
:::

### 設問2-4（15点）：Qini曲線と予算25%時点の比較

講義ノートと同じ考え方で、以下を実装せよ。

(a) `test_df`に限定し、`risk_score`順、`tau_hat_cf`順、ランダムの3方針について、配布割合を横軸、累積増分継続者数（$(\overline{Y}_{D=1,\text{top-}k} - \overline{Y}_{D=0,\text{top-}k}) \times k$）を縦軸としたQini型の累積uplift曲線を1枚に重ねて描け。
(b) 配布割合25%の時点における3方針の累積増分継続者数を表にまとめよ。
(c) churnスコア順のターゲティングが、25%時点でランダムと比べてどうなっているか（上回るか、下回るか、ほぼ同じか）を1〜2行でコメントせよ。

::: {.callout-tip}
## ヒント：Qini曲線の関数化

講義ノートの`qini_curve()`関数を流用してよい。ランダムのベンチマークにも`test_df`の群間差を使う。文献によってQini/uplift curveの正規化は異なるが、ここでは上記の定義に固定する。
:::

### 設問2-5（10点）：3つの候補ポリシーのIPW評価

次の3つの候補ポリシーを`test_df`に適用し、IPW政策価値と95%信頼区間を計算し、1つの図で比較せよ。

1. 全員にクーポンを配布する。
2. `risk_score`上位25%にクーポンを配布する。
3. `tau_hat_cf > 0`の人にクーポンを配布する。

```{r}
#| eval: false
ipw_policy_value <- function(policy_mask, D, Y, e = 0.5) {
  p_realized <- ifelse(D == 1, e, 1 - e)
  matched <- ____________________
  contrib <- ifelse(matched, Y / p_realized, 0)
  n <- length(Y)
  v_hat <- mean(contrib)
  se_hat <- sd(contrib) / sqrt(n)
  data.frame(v_hat = v_hat, se = se_hat,
             ci_low = v_hat - 1.96 * se_hat, ci_high = v_hat + 1.96 * se_hat)
}
```

::: {.callout-tip}
## ヒント：matchedの中身

`matched`は「実際の割り当て`D_i`が、評価したいポリシー`policy_mask[i]`と一致しているかどうか」を表す論理値ベクトルである。`policy_mask`は`TRUE`/`FALSE`（あるいは`1`/`0`）で「そのポリシーがこの人にクーポンを配るかどうか」を表しているので、`(D == 1) == policy_mask`のように書けばよい。

ここの通常のSEは、ルールが`test_df`から独立に固定されていることを条件とする。閾値やモデルをtestで選び直すと選択バイアスが入る。その場合は学習・validation・最終testの3分割、またはnested cross-fittingが必要である。OOB予測だけで通常のSEが保証されるわけではない。
:::

### 設問2-6（5点）：経営含意3行

あなたはこの動画配信サービスのCRM部門のアナリストである。上司から「解約リスクが高い会員に絞ってクーポンを配りたいので、リスクスコア上位のリストを月次で送ってほしい」と言われた。この課題の分析結果を踏まえて、**3行程度**で返信メモを書け。何が問題で、代わりにどうすべきかを簡潔に述べること。

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

(1) リスクスコアと施策効果が別物であること、(2) 価値・コスト・予算から候補ルールを作り、独立testで一律政策も含めて比較すること、(3) 運用時にランダムholdoutを残すこと、の3点を意識するとよい。
:::

## 提出物

- 本ファイル（`assignment11.qmd`）をレンダリングしたHTMLファイル。
- 全てのコードチャンクが上から順に実行可能であること（`eval: false`のチャンクは自分で`eval: true`に直すか、コードをコピーして実行すること）。