課題11：動画サブスクの解約引き留めクーポンとuplift targeting

計量経済学II

作者

Kei Ikegami

この課題について

第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を用意しているので、詰まったら遠慮なく開いてほしい。

重要この課題のゴール
  • 「解約リスクの高さ」と「施策への感応度（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であることを意味しない。
  • 解約リスクが非常に高い層では効果がほぼ消える。
ヒントヒント：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\)
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)を計算し、値の符号と大きさについてコメントせよ。

ヒントヒント

講義ノートの該当図とほぼ同じコードで描ける。負の相関になっていれば、「解約リスクが高いほどクーポンが効く」という直感が、この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はターゲティング方針の最終評価まで触らない。

ヒントヒント

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)を計算せよ。

ヒントヒント：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 に渡すこと。

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

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

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

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

  1. 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枚に重ねて描け。
  2. 配布割合25%の時点における3方針の累積増分継続者数を表にまとめよ。
  3. churnスコア順のターゲティングが、25%時点でランダムと比べてどうなっているか（上回るか、下回るか、ほぼ同じか）を1〜2行でコメントせよ。
ヒントヒント：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の人にクーポンを配布する。
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)
}
ヒントヒント：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行程度で返信メモを書け。何が問題で、代わりにどうすべきかを簡潔に述べること。

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

提出物

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