はじめに
過去に動的離散選択モデルについてまとめています。
今回はサブスクのプラン選択のシミュレーションデータに対して、モデルの当てを行ってみます。
動的離散選択モデルについて
動的離散選択モデルとは
動的離散選択モデル（Dynamic Discrete Choice Model）は、個人や企業が「将来の利益」を考慮しながら、複数の選択肢から行動を選ぶ意思決定を表現するモデルである。 労働経済学、マーケティング、教育分野などで広く利用されている。
例えば、ある個人が時点 において状態
のもとで行動
を選択するとする。
行動には効用
が存在し、個人は将来の効用も含めた期待値を最大化するように行動すると考える。
ここで は将来効用を現在価値へ割り引く割引率である。
ベルマン方程式
将来の状態は不確実であるため、状態遷移確率を用いて期待値を計算する。
このとき、価値関数 はベルマン方程式として再帰的に表現される。
つまり、現在の効用と将来価値の期待値を合わせて最適な行動を選択する。
サブスクプラン選択を題材に
今回は動的離散選択モデルをサブスクリプションのプラン変更の分析に適用するパターンを試してみます。
はじめに、プラン選択行動における効用関数は以下のように考える。
ここでは特定のプラン、
は利用度数、
は現在のプラン、
は構造パラメータ、
は各プラン特徴リスト。
各プランには利用度数があることを仮定し、各プランでの利用度数には上限があると設定します。 そのため、現状の利用度数または最大度数が現状利用している度数となる。 上位プランに上げたことによる最大度数が増加して利用度数が自然と増えるという動きは表現していない。
第4項はプランスイッチングコストを表現しており、現状と異なるプランを選択した場合にはだけのコストがかかることを意味する。
移行前と移行後のプランのペアでコストが変化する可能性があるが、今回は均一なコストがかかるとしている。
この効用値を用い、効用が高い選択肢を選択していくことになる。
# 効用を計算する関数定義 flow_util <- function(p_d, u, p_c, theta, plan = PLAN) { eff_usage <- pmin(u, plan$capacity[p_d]) # 上限内の有効利用量 theta["value_usage"] * eff_usage + theta["value_support"] * plan$support[p_d] - theta["price"] * plan$price[p_d] - theta["switch_cost"] * (p_d != p_c) }
今回のシミュレーションでは、4プラン存在するとして、各パラメータの設定は以下のようにしています。
library(ggplot2) library(dplyr) library(tidyr) library(scales) set.seed(42) # 定数・プラン定義・真のパラメータ定義 ------------------------ PLAN_NAMES <- c("Free", "Simple", "Pro", "Max") N_PLANS <- 4L N_USAGE <- 60L # 利用度グリッド (1〜60) BETA <- 0.95 # 割引因子 EULER_GAMMA <- 0.5772156649 # プラン特徴の定義 ------------------------ PLAN <- list( price = c( 0, 9, 19, 39), # 月額料金（ドル） capacity = c( 12, 25, 45, 60), # 利用上限グリッド値 support = c(0.0, 0.3, 0.7, 1.0) # サポート品質スコア ) # 真の構造パラメータ定義 ------------------------ THETA_TRUE <- c( value_usage = 1.5, # 有効利用量 1 単位あたりの価値 value_support = 5.0, # サポート品質スコア price = 1.1, # 価格影響 switch_cost = 12.0 # プラン変更の摩擦コスト ) # ---- 利用度遷移 (Δu ∈ {-2,-1,0,1,2,3}) 定義 ------------------------ # 高プランほど利用度が増加しやすい設定 DELTA_VALS <- c(-2L, -1L, 0L, 1L, 2L, 3L) TRANS_PROBS <- list( Free = c(0.07, 0.18, 0.35, 0.30, 0.07, 0.03), Simple = c(0.04, 0.12, 0.38, 0.28, 0.13, 0.05), Pro = c(0.02, 0.08, 0.30, 0.35, 0.17, 0.08), Max = c(0.02, 0.05, 0.23, 0.35, 0.23, 0.12) ) # 利用度遷移行列の構築 ---------------------- build_usage_trans <- function(n, probs, deltas) { P <- matrix(0.0, n, n) for (i in seq_len(n)) { for (k in seq_along(deltas)) { j <- min(max(i + deltas[k], 1L), n) P[i,j] <- P[i,j] + probs[k] } } P } P_usage <- lapply(TRANS_PROBS, build_usage_trans, n = N_USAGE, deltas = DELTA_VALS)
これらの定義に従い、サンプルデータを生成してみます。
# データ生成シミュレーション ------------------- simulate_data <- function(theta, n_users = 300, n_periods = 72) { res <- value_fn_iter(theta) rows <- vector("list", n_users * n_periods) idx <- 1L for (usr in seq_len(n_users)) { u <- sample(3:20, 1) p_c <- sample(1:N_PLANS, 1, prob = c(0.5, 0.25, 0.15, 0.10)) for (t in seq_len(n_periods)) { d <- sample(1:N_PLANS, 1, prob = res$ccp[u, , p_c]) rows[[idx]] <- list(user = usr, period = t, usage = u, current_plan = p_c, choice = d) idx <- idx + 1L delta <- sample(DELTA_VALS, 1, prob = TRANS_PROBS[[d]]) u <- min(max(u + delta, 1L), N_USAGE) p_c <- d } } df <- as.data.frame(do.call(rbind, lapply(rows, as.data.frame))) df[] <- lapply(df, unlist) df$choice_name <- factor(PLAN_NAMES[df$choice], levels = PLAN_NAMES) df$current_name <- factor(PLAN_NAMES[df$current_plan], levels = PLAN_NAMES) df } cat("データ生成中...\n") df <- simulate_data(THETA_TRUE) cat(sprintf("観測数: %d (ユーザー数: %d)\n", nrow(df), length(unique(df$user))))
生成されたプラン選択行動における各プランの選択割合は以下のようになった。
> print(round(prop.table(table(df$choice_name)), 3)) Free Simple Pro Max 0.249 0.328 0.347 0.076
プラン遷移行列は以下のようになった。
PLAN_COLORS <- c(Free = "#9CA3AF", Simple = "#2563EB", Pro = "#0F6E56", Max = "#DC2626") trans_df <- df |> count(current_name, choice_name) |> group_by(current_name) |> mutate(prob = n / sum(n)) |> ungroup() p1 <- trans_df |> ggplot() + aes(x = choice_name, y = current_name, fill = prob) + geom_tile(color = "white", linewidth = 1) + geom_text(aes(label = percent(prob, accuracy = 0.1)), size = 4.5) + scale_fill_gradient(low = "#EFF6FF", high = "#1D4ED8", labels = percent) + labs(title = "プラン遷移確率（観測データ）", #subtitle = "行 = 現プラン、列 = 次期プラン", x = "次期プラン（選択）", y = "現プラン", fill = "遷移確率") + theme_bw(base_size = 13, base_family = "HiraKakuProN-W3") p1
プランごとの利用度数の分布は以下のようになった。
p2 <- df |> ggplot() + aes(x = usage, fill = choice_name) + geom_histogram(bins = 30, alpha = 0.85, position = "identity") + scale_fill_manual(values = PLAN_COLORS) + facet_wrap(~choice_name, ncol = 1) + labs(title = "選択プラン別の利用度分布", #subtitle = "高プランは利用度の高いユーザーに選ばれやすい", x = "利用度グリッド", y = "観測数", fill = "プラン") + theme_bw(base_size = 12, base_family = "HiraKakuProN-W3") + theme(legend.position = "none") p2
推定コード
NFXP推定に必要な関数の定義。
# 価値関数反復の定義 value_fn_iter <- function(theta, plan = PLAN, beta = BETA, tol = 1e-12, max_iter = 6000) { u_idx <- seq_len(N_USAGE) EV <- matrix(0.0, N_USAGE, N_PLANS) # 将来価値の事前計算キャッシュ（EV固定時） calc_future <- function(EV) { fv <- matrix(0.0, N_USAGE, N_PLANS) for (d in 1:N_PLANS) fv[, d] <- beta * as.vector(P_usage[[d]] %*% EV[, d]) fv } for (iter in seq_len(max_iter)) { fv <- calc_future(EV) EV_new <- matrix(0.0, N_USAGE, N_PLANS) for (p_c in 1:N_PLANS) { v_mat <- matrix(0.0, N_USAGE, N_PLANS) for (d in 1:N_PLANS) v_mat[, d] <- flow_util(d, u_idx, p_c, theta, plan) + fv[, d] mx <- apply(v_mat, 1, max) EV_new[, p_c] <- log(rowSums(exp(v_mat - mx))) + mx + EULER_GAMMA } if (max(abs(EV_new - EV)) < tol) { EV <- EV_new; break } EV <- EV_new } # CCP[u, d, p_c] の計算 fv <- calc_future(EV) ccp <- array(NA_real_, dim = c(N_USAGE, N_PLANS, N_PLANS), dimnames = list(NULL, PLAN_NAMES, PLAN_NAMES)) for (p_c in 1:N_PLANS) { v_mat <- matrix(0.0, N_USAGE, N_PLANS) for (d in 1:N_PLANS) v_mat[, d] <- flow_util(d, u_idx, p_c, theta, plan) + fv[, d] mx <- apply(v_mat, 1, max) log_denom <- log(rowSums(exp(v_mat - mx))) + mx for (d in 1:N_PLANS) ccp[, d, p_c] <- exp(v_mat[, d] - log_denom) } list(EV = EV, ccp = ccp) } # 外側ループ用の対数尤度 log_lik <- function(params, data, plan = PLAN) { theta <- c( value_usage = exp(params[1]), value_support = exp(params[2]), price = exp(params[3]), switch_cost = exp(params[4]) ) res <- tryCatch(value_fn_iter(theta, plan = plan), error = function(e) NULL) if (is.null(res)) return(1e10) u <- data$usage d <- data$choice p <- data$price p_c <- data$current_plan ll <- sum(vapply(seq_len(nrow(data)), function(i) log(max(res$ccp[u[i], d[i], p_c[i]], 1e-300)), numeric(1))) if (!is.finite(ll)) return(1e10) -ll }
上記の推定関数を用いて実際に推定を行う。
opt <- optim( par = log(c(1.2, 4.0, 1.0, 4.0)), fn = log_lik, data = df, method = "Nelder-Mead", control = list(maxit = 8000, reltol = 1e-11) ) theta_est <- setNames(exp(opt$par), c("value_usage", "value_support", "price","switch_cost"))
推定結果は以下のようになった。 真値に近い値を取得することができている。
> for (nm in names(THETA_TRUE)) + message(sprintf(" %-15s: 真値 = %5.3f, 推定値 = %5.3f", + nm, THETA_TRUE[nm], theta_est[nm])) value_usage : 真値 = 1.500, 推定値 = 1.434 value_support : 真値 = 5.000, 推定値 = 4.628 price : 真値 = 1.100, 推定値 = 1.047 switch_cost : 真値 = 12.000, 推定値 = 11.490
反実仮想分析
サブスクリプションのプランの構成要素を変更した際に、ユーザーの反応傾向がどのように変化するのかという反実仮想を求めて分析することができる。
いくつかのプラン要素の変更を行った時の各プランの契約割合を算出する。
はじめに、反実仮想分析用関数を作成する。
推定された係数パラメーターを用いて、各プランの条件(特徴)を変更した場合の選択割合を取得して、その平均選択確率を算出する。
# 反実仮想分析用関数の作成 run_counterfactural <- function(label, plan_cf, res_base) { res_cf <- value_fn_iter(theta_est, plan = plan_cf) cat(sprintf("\n■ %s\n", label)) # 全利用度・全現プランの平均選択確率 base_avg <- apply(res_base$ccp, 2, mean) cf_avg <- apply(res_cf$ccp, 2, mean) tab <- data.frame( プラン = PLAN_NAMES, ベース選択率 = round(base_avg, 4), 変更後選択率 = round(cf_avg, 4), 変化幅 = round(cf_avg - base_avg, 4) ) print(tab, row.names = FALSE) invisible(list(base = res_base, cf = res_cf)) } res_base <- value_fn_iter(theta_est)
一つ目はプランの値上げ。
# CF1: Simple を 9 → 15 ドルに値上げ plan_cf1 <- PLAN plan_cf1$price[2] <- 15 run_counterfactural("Simple を $9 → $15 に値上げ", plan_cf1, res_base)
結果は以下となり、Simpleプランの契約割合が大幅に低下し、Proプランの契約が増加する結果となった。
■ Simple を $9 → $15 に値上げ プラン ベース選択率 変更後選択率 変化幅 Free 0.1729 0.2133 0.0404 Simple 0.2156 0.0011 -0.2144 Pro 0.5373 0.7113 0.1740 Max 0.0743 0.0743 0.0000
次にProプランの利用上限を変更してみる。
# CF2: Pro の利用上限を 45 → 60 に拡張（容量増強） plan_cf2 <- PLAN plan_cf2$capacity[3] <- 60 run_counterfactural("Pro の利用上限を 45 → 60 に拡張", plan_cf2, res_base)
Proプランの契約割合が増え、Simple, Maxプランの契約数が減っている。 Proプランの価値が上がり、Maxプランの契約の意味が薄れてしまった影響であると考えられる。
■ Pro の利用上限を 45 → 60 に拡張 プラン ベース選択率 変更後選択率 変化幅 Free 0.1729 0.1404 -0.0324 Simple 0.2156 0.1647 -0.0509 Pro 0.5373 0.6949 0.1576 Max 0.0743 0.0000 -0.0743
3つ目として、Freeプランの機能制限の上限を引き下げた場合を確認した。
# CF3: Free の機能制限強化（上限 12 → 6） plan_cf3 <- PLAN plan_cf3$capacity[1] <- 6 run_counterfactural("Free の利用上限を 12 → 6 に縮小", plan_cf3, res_base)
Freeプランの契約割合が低下し、SimpleとProの割合が増加した。
■ Free の利用上限を 12 → 6 に縮小 プラン ベース選択率 変更後選択率 変化幅 Free 0.1729 0.0828 -0.0901 Simple 0.2156 0.2890 0.0734 Pro 0.5373 0.5539 0.0166 Max 0.0743 0.0743 0.0000
最後にSimpleプランを値上げした場合の利用度ごとのプラン選択確率がどのように変化するかを確認する。
# ---- 図3: 反実仮想 — Simple 値上げの効果 ---- res_cf1 <- value_fn_iter(theta_est, plan = plan_cf1) cf_compare <- expand.grid(usage = 1:N_USAGE, choice = PLAN_NAMES, current = PLAN_NAMES, stringsAsFactors = FALSE) |> rowwise() |> mutate( base = res_base$ccp[usage, match(choice, PLAN_NAMES), match(current, PLAN_NAMES)], cf = res_cf1$ccp[ usage, match(choice, PLAN_NAMES), match(current, PLAN_NAMES)], diff = cf - base ) |> ungroup() |> mutate(choice = factor(choice, levels = PLAN_NAMES), current = factor(current, levels = PLAN_NAMES)) # Simpleを現プランとするユーザーへの影響に絞る p3 <- cf_compare |> filter(current == "Simple") |> pivot_longer(c(base, cf), names_to = "scenario", values_to = "prob") |> mutate(scenario = recode(scenario, base = "ベースライン ($9)", cf = "値上げ後 ($15)")) |> ggplot() + aes(x = usage, y = prob, color = choice, linetype = scenario)+ geom_line(linewidth = 0.7) + scale_color_manual(values = PLAN_COLORS) + scale_linetype_manual(values = c("solid", "dashed")) + labs(title = "反実仮想: Simple を $9 → $15 に値上げ（現プラン = Simple ユーザー）", x = "利用度グリッド", y = "選択確率", color = "選択プラン", linetype = "シナリオ") + theme_bw(base_size = 12, base_family = "HiraKakuProN-W3") + theme(legend.position = "bottom") p3
値上げによりSimpleプランの契約は Free へのダウングレードと Pro へのアップグレードになり、Simpleプランの契約が皆無になるという結果だった。 そのため、このようなプラン変更を行う場合は、プラン自体を4プランから3プランへ変更するのが望ましいと考えられる。 そして、Freeプランの利用度数の上限は17程度が望ましいと考えられます。
