FacebookのMMMのOSS「Robyn」のチュートリアルさわってみた

目次
・Robynとは
・とりあえずチュートリアルやってみる
・向き合い方
・参考情報

Robynとは

Robyn（ロビン）はFacebook（META）が開発しているMarketing-Mix-Modeling（以降、MMM）のオープンソース（https://facebookexperimental.github.io/Robyn/）です。主にR言語で開発されています。（Python版は目下開発中らしいです。）

MMMは、マーケティングの広告投資の予算を、効果を最大化するためにどこに配分するかを数理的なモデルで決めようとするものです。
そもそも、MMMにそれほど明るくない方もおられると思いますが、その際は、『データ活用のための数理モデリング入門』の100ページ目に目的やら簡単な概要が載っていますので、参考にされると良いと思います。

また、RobynにはRobyn Open Source MMM Usersというユーザーグループがあるようです。そこそこ活発に運営されているようです。

私は以前、このコミュニティーのイベント（Facebook Project Robyn Open Source MMM 2021 Community Summit）があったので聴講しました。英語が苦手なので何度も聞きなおしましたが。。

この会では、世界中のマーケターがRobynを使ってプロモーションの分析をしているのを知れました。彼らはRobynのアーリーアダプターなんだろうなと思いました。

とりあえずチュートリアルやってみる

とりあえず、Robynがどういうツールなのか知るためにチュートリアル（https://facebookexperimental.github.io/Robyn/docs/quick-start/）をやってみることにしました。

まずは、Rを最新版にします。今回は「R version 4.2.0 (2022-04-22) — “Vigorous Calisthenics”」にしました。

# package install
install.packages("Robyn")

1 2	# package install install.packages("Robyn")

結構な数の依存パッケージがインストールされているようです。スクロールバーがめちゃ小さくなりました。時間もそれなりにかかるようです。

## RobynではPythonのライブラリを使う
install.packages("reticulate")
library(reticulate)
virtualenv_create("r-reticulate")
py_install("nevergrad", pip = TRUE)
use_virtualenv("r-reticulate", required = TRUE)

## RobynではPythonのライブラリを使う

install.packages("reticulate")

library(reticulate)

virtualenv_create("r-reticulate")

py_install("nevergrad", pip = TRUE)

use_virtualenv("r-reticulate", required = TRUE)

Robynパッケージを読み込み、サンプルデータを呼び出します。

library(Robyn)

data("dt_simulated_weekly")
head(dt_simulated_weekly)

data("dt_prophet_holidays")
head(dt_prophet_holidays)

library(Robyn)

data("dt_simulated_weekly")

head(dt_simulated_weekly)

data("dt_prophet_holidays")

head(dt_prophet_holidays)

モデルを作る際の入力変数の設定をします。途中でProphetが使われているようです。確かに、Prophetは時系列の分析にちょうど適したライブラリではあります。

さて、指定する入力変数は結構多くて、以下の通りです。
・データセット
・従属変数
・従属変数のタイプ
・Prophetでの周期性などのオプション
・国
・競合の情報やイベントなどのコンテキスト変数
・ペイドメディアの支出
・ペイドメディアのインプレッションやクリックなど
・オーガニックな変数
・コンテキスト変数のなかでオーガニックなもの
・期初、期末
・広告の残存効果

実務でMMMを使っているものからすると馴染み深いですので、「ああ、このデータですね」とスッと頭に入ってきます。コードでいろいろ個別にやるととっ散らかるので、こういうカタマリで処理を実行できるのはいいですね。MMMが初めての方は、扱うデータセットの概要をよく調べてから入力変数にするようにしたほうがいいと思います。

## 入力変数の設定

InputCollect <- robyn_inputs(
  dt_input = dt_simulated_weekly
  , dt_holidays = dt_prophet_holidays
  , date_var = "DATE"  # "2020-01-01"みたいにする
  , dep_var = "revenue" # 従属変数
  , dep_var_type = "revenue" # "revenue" (ROI) or "conversion" (CPA)
  , prophet_vars = c("trend", "season", "holiday") # "trend","season", "weekday" & "holiday"
  , prophet_country = "DE" # input one country. dt_prophet_holidays includes 59 countries by default
  , context_vars = c("competitor_sales_B", "events") # e.g. competitors, discount, unemployment etc
  , paid_media_spends = c("tv_S","ooh_S",	"print_S"	,"facebook_S", "search_S") # mandatory input
  , paid_media_vars = c("tv_S", "ooh_S"	,	"print_S"	,"facebook_I" ,"search_clicks_P") # mandatory.
  # paid_media_vars must have same order as paid_media_spends. Use media exposure metrics like
  # impressions, GRP etc. If not applicable, use spend instead.
  , organic_vars = c("newsletter") # marketing activity without media spend
  , factor_vars = c("events") # specify which variables in context_vars or organic_vars are factorial
  , window_start = "2016-11-23"
  , window_end = "2018-08-22"
  , adstock = "geometric" # geometric, weibull_cdf or weibull_pdf.
)

print(InputCollect)

## 入力変数の設定

InputCollect <- robyn_inputs(

dt_input = dt_simulated_weekly

, dt_holidays = dt_prophet_holidays

, date_var = "DATE" # "2020-01-01"みたいにする

, dep_var = "revenue" # 従属変数

, dep_var_type = "revenue" # "revenue" (ROI) or "conversion" (CPA)

, prophet_vars = c("trend", "season", "holiday") # "trend","season", "weekday" & "holiday"

, prophet_country = "DE" # input one country. dt_prophet_holidays includes 59 countries by default

, context_vars = c("competitor_sales_B", "events") # e.g. competitors, discount, unemployment etc

, paid_media_spends = c("tv_S","ooh_S", "print_S" ,"facebook_S", "search_S") # mandatory input

, paid_media_vars = c("tv_S", "ooh_S" , "print_S" ,"facebook_I" ,"search_clicks_P") # mandatory.

# paid_media_vars must have same order as paid_media_spends. Use media exposure metrics like

# impressions, GRP etc. If not applicable, use spend instead.

, organic_vars = c("newsletter") # marketing activity without media spend

, factor_vars = c("events") # specify which variables in context_vars or organic_vars are factorial

, window_start = "2016-11-23"

, window_end = "2018-08-22"

, adstock = "geometric" # geometric, weibull_cdf or weibull_pdf.

)

print(InputCollect)

次に、ハイパーパラメータの設定を行います。『StanとRでベイズ統計モデリング (Wonderful R)』の作法に従うならば、あまりハイパーパラメータを恣意的に決めて収束させるのはやりたくないですが、明らかに符号がおかしいとかの制約は付けてもいいのかなと思います。

## ハイパーパラメータの設定
hyper_names(adstock = InputCollect$adstock, all_media = InputCollect$all_media)

hyperparameters <- list(
  facebook_S_alphas = c(0.5, 3)
  ,facebook_S_gammas = c(0.3, 1)
  ,facebook_S_thetas = c(0, 0.3)
  
  ,print_S_alphas = c(0.5, 3)
  ,print_S_gammas = c(0.3, 1)
  ,print_S_thetas = c(0.1, 0.4)
  
  ,tv_S_alphas = c(0.5, 3)
  ,tv_S_gammas = c(0.3, 1)
  ,tv_S_thetas = c(0.3, 0.8)
  
  ,search_S_alphas = c(0.5, 3)
  ,search_S_gammas = c(0.3, 1)
  ,search_S_thetas = c(0, 0.3)
  
  ,ooh_S_alphas = c(0.5, 3)
  ,ooh_S_gammas = c(0.3, 1)
  ,ooh_S_thetas = c(0.1, 0.4)
  
  ,newsletter_alphas = c(0.5, 3)
  ,newsletter_gammas = c(0.3, 1)
  ,newsletter_thetas = c(0.1, 0.4)
)

## ハイパーパラメータの設定

hyper_names(adstock = InputCollect$adstock, all_media = InputCollect$all_media)

hyperparameters <- list(

facebook_S_alphas = c(0.5, 3)

,facebook_S_gammas = c(0.3, 1)

,facebook_S_thetas = c(0, 0.3)

,print_S_alphas = c(0.5, 3)

,print_S_gammas = c(0.3, 1)

,print_S_thetas = c(0.1, 0.4)

,tv_S_alphas = c(0.5, 3)

,tv_S_gammas = c(0.3, 1)

,tv_S_thetas = c(0.3, 0.8)

,search_S_alphas = c(0.5, 3)

,search_S_gammas = c(0.3, 1)

,search_S_thetas = c(0, 0.3)

,ooh_S_alphas = c(0.5, 3)

,ooh_S_gammas = c(0.3, 1)

,ooh_S_thetas = c(0.1, 0.4)

,newsletter_alphas = c(0.5, 3)

,newsletter_gammas = c(0.3, 1)

,newsletter_thetas = c(0.1, 0.4)

)

ハイパーパラメータを設定したら、アルゴリズムを実行します。裏側でベイズ推定をしていることから、結構時間がかかります。Prophetを動かすということはStanを動かしていることと同義ですから。

InputCollect <- robyn_inputs(InputCollect = InputCollect, hyperparameters = hyperparameters)
print(InputCollect)

## Run all trials and iterations. Use ?robyn_run to check parameter definition
OutputModels <- robyn_run(
  InputCollect = InputCollect # feed in all model specification
  #, cores = NULL # default
  #, add_penalty_factor = FALSE # Untested feature. Use with caution.
  , iterations = 2000 # recommended for the dummy dataset
  , trials = 5 # recommended for the dummy dataset
  , outputs = FALSE # outputs = FALSE disables direct model output
)
print(OutputModels)

InputCollect <- robyn_inputs(InputCollect = InputCollect, hyperparameters = hyperparameters)

print(InputCollect)

## Run all trials and iterations. Use ?robyn_run to check parameter definition

OutputModels <- robyn_run(

InputCollect = InputCollect # feed in all model specification

#, cores = NULL # default

#, add_penalty_factor = FALSE # Untested feature. Use with caution.

, iterations = 2000 # recommended for the dummy dataset

, trials = 5 # recommended for the dummy dataset

, outputs = FALSE # outputs = FALSE disables direct model output

)

print(OutputModels)

イタレーションごとのモデルの目的関数の事後分布を可視化します。徐々に収束してそうに見えます。

## Check MOO (multi-objective optimization) convergence plots
OutputModels$convergence$moo_distrb_plot

1 2	## Check MOO (multi-objective optimization) convergence plots OutputModels$convergence$moo_distrb_plot

OutputModels$convergence$moo_cloud_plot

1	OutputModels$convergence$moo_cloud_plot

モデルがいろいろと求まったので、パレート最適な組み合わせの計算をします。

## Calculate Pareto optimality, cluster and export results and plots. See ?robyn_outputs
OutputCollect <- robyn_outputs(
  InputCollect, OutputModels
  , pareto_fronts = 3
  # , calibration_constraint = 0.1 # range c(0.01, 0.1) & default at 0.1
  , csv_out = "pareto" # "pareto" or "all"
  , clusters = TRUE # Set to TRUE to cluster similar models by ROAS. See ?robyn_clusters
  , plot_pareto = TRUE # Set to FALSE to deactivate plotting and saving model one-pagers
  , plot_folder = robyn_object # path for plots export
)
print(OutputCollect)

## Calculate Pareto optimality, cluster and export results and plots. See ?robyn_outputs

OutputCollect <- robyn_outputs(

InputCollect, OutputModels

, pareto_fronts = 3

# , calibration_constraint = 0.1 # range c(0.01, 0.1) & default at 0.1

, csv_out = "pareto" # "pareto" or "all"

, clusters = TRUE # Set to TRUE to cluster similar models by ROAS. See ?robyn_clusters

, plot_pareto = TRUE # Set to FALSE to deactivate plotting and saving model one-pagers

, plot_folder = robyn_object # path for plots export

)

print(OutputCollect)

パレート最適な組み合わせで返された複数のモデルから一つを選びます。

print(OutputCollect)
select_model <- "1_143_4" # select one from above
ExportedModel <- robyn_save(
  robyn_object = robyn_object # model object location and name
  , select_model = select_model # selected model ID
  , InputCollect = InputCollect
  , OutputCollect = OutputCollect
)
print(ExportedModel)

print(OutputCollect)

select_model <- "1_143_4" # select one from above

ExportedModel <- robyn_save(

robyn_object = robyn_object # model object location and name

, select_model = select_model # selected model ID

, InputCollect = InputCollect

, OutputCollect = OutputCollect

)

print(ExportedModel)

選んだモデルの係数などを確認します。

> print(ExportedModel)
Exported file: MyRobyn.RDS
Exported model: 1_143_4

Media Summary for Selected Model:
rn      coef mean_spend mean_response   roi_mean total_spend total_response
1: facebook_S  23432.93  136111.15      11556.96 0.08490823     5988890       594478.3
2:      ooh_S 237278.85  262577.85      90572.07 0.34493417    10240536      5956088.7
3:    print_S 450597.23   77589.33      79774.15 1.02815874     2793216      3123007.2
4:   search_S 185857.45   47618.18      28802.09 0.60485482     3666600      2440274.7
5:       tv_S 551227.13  256198.38     110879.89 0.43278921    10247935      5859543.8
roi_total
1: 0.09926351
2: 0.58161884
3: 1.11806864
4: 0.66554157
5: 0.57177799

> print(ExportedModel)

Exported file: MyRobyn.RDS

Exported model: 1_143_4

Media Summary for Selected Model:

rn coef mean_spend mean_response roi_mean total_spend total_response

1: facebook_S 23432.93 136111.15 11556.96 0.08490823 5988890 594478.3

2: ooh_S 237278.85 262577.85 90572.07 0.34493417 10240536 5956088.7

3: print_S 450597.23 77589.33 79774.15 1.02815874 2793216 3123007.2

4: search_S 185857.45 47618.18 28802.09 0.60485482 3666600 2440274.7

5: tv_S 551227.13 256198.38 110879.89 0.43278921 10247935 5859543.8

roi_total

1: 0.09926351

2: 0.58161884

3: 1.11806864

4: 0.66554157

5: 0.57177799

この選んだモデルをもとに、最適なアロケーションを計算します。

# Run the "max_historical_response" scenario: "What's the revenue lift potential with the
# same historical spend level and what is the spend mix?"
AllocatorCollect1 <- robyn_allocator(
  InputCollect = InputCollect
  , OutputCollect = OutputCollect
  , select_model = select_model
  , scenario = "max_historical_response"
  , channel_constr_low = 0.7
  , channel_constr_up = c(1.2, 1.5, 1.5, 1.5, 1.5)
  , export = TRUE
  , date_min = "2016-11-21"
  , date_max = "2018-08-20"
)
print(AllocatorCollect1)

# Run the "max_historical_response" scenario: "What's the revenue lift potential with the

# same historical spend level and what is the spend mix?"

AllocatorCollect1 <- robyn_allocator(

InputCollect = InputCollect

, OutputCollect = OutputCollect

, select_model = select_model

, scenario = "max_historical_response"

, channel_constr_low = 0.7

, channel_constr_up = c(1.2, 1.5, 1.5, 1.5, 1.5)

, export = TRUE

, date_min = "2016-11-21"

, date_max = "2018-08-20"

)

print(AllocatorCollect1)

推定した、選んだモデルでの最適な広告のアロケーション結果を出力します。予算を削った方がいい広告経路、増やした方がいい広告経路などが示されます。

> print(AllocatorCollect1)
Model ID: 1_143_4
Scenario: Maximum Historical Response
Media Skipped (coef = 0): None
Relative Spend Increase: 0% (+0)
Total Response Increase (Optimized): 31.2%
Window: 2016-11-21:2018-08-20 (92 weeks)

Allocation Summary:

- facebook_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 17.4% -> Optimized = 12.2%
Mean Response: 11,557 -> Optimized = 8,016
Mean Spend (per time unit): 136.1K -> Optimized = 95.28K [Delta = -30%]

- ooh_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 33.7% -> Optimized = 24.3%
Mean Response: 90,572 -> Optimized = 80,605
Mean Spend (per time unit): 262.6K -> Optimized = 189.6K [Delta = -28%]

- print_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 9.95% -> Optimized = 14.9%
Mean Response: 79,774 -> Optimized = 140,513
Mean Spend (per time unit): 77.59K -> Optimized = 116.4K [Delta = 50%]

- search_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 6.1% -> Optimized = 9.16%
Mean Response: 28,802 -> Optimized = 57,873
Mean Spend (per time unit): 47.62K -> Optimized = 71.43K [Delta = 50%]

- tv_S:
Optimizable Range (bounds): [-30%, 20%]
Mean Spend Share (avg): 32.8% -> Optimized = 39.4%
Mean Response: 110,880 -> Optimized = 134,921
Mean Spend (per time unit): 256.2K -> Optimized = 307.4K [Delta = 20%]

> print(AllocatorCollect1)

Model ID: 1_143_4

Scenario: Maximum Historical Response

Media Skipped (coef = 0): None

Relative Spend Increase: 0% (+0)

Total Response Increase (Optimized): 31.2%

Window: 2016-11-21:2018-08-20 (92 weeks)

Allocation Summary:

- facebook_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 17.4% -> Optimized = 12.2%

Mean Response: 11,557 -> Optimized = 8,016

Mean Spend (per time unit): 136.1K -> Optimized = 95.28K [Delta = -30%]

- ooh_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 33.7% -> Optimized = 24.3%

Mean Response: 90,572 -> Optimized = 80,605

Mean Spend (per time unit): 262.6K -> Optimized = 189.6K [Delta = -28%]

- print_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 9.95% -> Optimized = 14.9%

Mean Response: 79,774 -> Optimized = 140,513

Mean Spend (per time unit): 77.59K -> Optimized = 116.4K [Delta = 50%]

- search_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 6.1% -> Optimized = 9.16%

Mean Response: 28,802 -> Optimized = 57,873

Mean Spend (per time unit): 47.62K -> Optimized = 71.43K [Delta = 50%]

- tv_S:

Optimizable Range (bounds): [-30%, 20%]

Mean Spend Share (avg): 32.8% -> Optimized = 39.4%

Mean Response: 110,880 -> Optimized = 134,921

Mean Spend (per time unit): 256.2K -> Optimized = 307.4K [Delta = 20%]

続いて、支出の上限を決めた上での、7日間でのアロケーションを行います。

AllocatorCollect2 <- robyn_allocator(
  InputCollect = InputCollect
  , OutputCollect = OutputCollect
  , select_model = select_model
  , scenario = "max_response_expected_spend"
  , channel_constr_low = c(0.7, 0.7, 0.7, 0.7, 0.7)
  , channel_constr_up = c(1.2, 1.5, 1.5, 1.5, 1.5)
  , expected_spend = 1000000 # Total spend to be simulated
  , expected_spend_days = 7 # Duration of expected_spend in days
  , export = TRUE
)

AllocatorCollect2 <- robyn_allocator(

InputCollect = InputCollect

, OutputCollect = OutputCollect

, select_model = select_model

, scenario = "max_response_expected_spend"

, channel_constr_low = c(0.7, 0.7, 0.7, 0.7, 0.7)

, channel_constr_up = c(1.2, 1.5, 1.5, 1.5, 1.5)

, expected_spend = 1000000 # Total spend to be simulated

, expected_spend_days = 7 # Duration of expected_spend in days

, export = TRUE

)

こちらが、出力した結果です。

> print(AllocatorCollect2)
Model ID: 1_143_4
Scenario: Maximum Response with Expected Spend
Media Skipped (coef = 0): None
Relative Spend Increase: 28.2% (+1.1M in 7 days)
Total Response Increase (Optimized): 38.8%
Window: 2016-11-21:2018-08-20 (92 weeks)

Allocation Summary:

- facebook_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 17.4% -> Optimized = 12%
Mean Response: 11,557 -> Optimized = 10,238
Mean Spend (per time unit): 136.1K -> Optimized = 119.7K [Delta = -12%]

- ooh_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 33.7% -> Optimized = 38.5%
Mean Response: 90,572 -> Optimized = 102,837
Mean Spend (per time unit): 262.6K -> Optimized = 385.1K [Delta = 47%]

- print_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 9.95% -> Optimized = 11.6%
Mean Response: 79,774 -> Optimized = 140,513
Mean Spend (per time unit): 77.59K -> Optimized = 116.4K [Delta = 50%]

- search_S:
Optimizable Range (bounds): [-30%, 50%]
Mean Spend Share (avg): 6.1% -> Optimized = 7.14%
Mean Response: 28,802 -> Optimized = 57,873
Mean Spend (per time unit): 47.62K -> Optimized = 71.43K [Delta = 50%]

- tv_S:
Optimizable Range (bounds): [-30%, 20%]
Mean Spend Share (avg): 32.8% -> Optimized = 30.7%
Mean Response: 110,880 -> Optimized = 134,921
Mean Spend (per time unit): 256.2K -> Optimized = 307.4K [Delta = 20%]

> print(AllocatorCollect2)

Model ID: 1_143_4

Scenario: Maximum Response with Expected Spend

Media Skipped (coef = 0): None

Relative Spend Increase: 28.2% (+1.1M in 7 days)

Total Response Increase (Optimized): 38.8%

Window: 2016-11-21:2018-08-20 (92 weeks)

Allocation Summary:

- facebook_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 17.4% -> Optimized = 12%

Mean Response: 11,557 -> Optimized = 10,238

Mean Spend (per time unit): 136.1K -> Optimized = 119.7K [Delta = -12%]

- ooh_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 33.7% -> Optimized = 38.5%

Mean Response: 90,572 -> Optimized = 102,837

Mean Spend (per time unit): 262.6K -> Optimized = 385.1K [Delta = 47%]

- print_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 9.95% -> Optimized = 11.6%

Mean Response: 79,774 -> Optimized = 140,513

Mean Spend (per time unit): 77.59K -> Optimized = 116.4K [Delta = 50%]

- search_S:

Optimizable Range (bounds): [-30%, 50%]

Mean Spend Share (avg): 6.1% -> Optimized = 7.14%

Mean Response: 28,802 -> Optimized = 57,873

Mean Spend (per time unit): 47.62K -> Optimized = 71.43K [Delta = 50%]

- tv_S:

Optimizable Range (bounds): [-30%, 20%]

Mean Spend Share (avg): 32.8% -> Optimized = 30.7%

Mean Response: 110,880 -> Optimized = 134,921

Mean Spend (per time unit): 256.2K -> Optimized = 307.4K [Delta = 20%]

続いて、特定の広告経路の目的関数に対しての影響度が支出に応じてどう変わっていくか、つまりサチっているかどうかを見てみます。

## QA optimal response
# Pick any media variable: InputCollect$all_media
select_media <- "search_S"
# For paid_media_spends set metric_value as your optimal spend
metric_value <- AllocatorCollect1$dt_optimOut[channels == select_media, optmSpendUnit]
# # For paid_media_vars and organic_vars, manually pick a value
# metric_value <- 10000

if (TRUE) {
  optimal_response_allocator <- AllocatorCollect1$dt_optimOut[
    channels == select_media, optmResponseUnit]
  optimal_response <- robyn_response(
    robyn_object = robyn_object,
    select_build = 0,
    media_metric = select_media,
    metric_value = metric_value)
  plot(optimal_response$plot)
  if (length(optimal_response_allocator) > 0) {
cat("QA if results from robyn_allocator and robyn_response agree: ")
cat(round(optimal_response_allocator) == round(optimal_response$response), "( ")
cat(optimal_response$response, "==", optimal_response_allocator, ")\n")
}
}

## QA optimal response

# Pick any media variable: InputCollect$all_media

select_media <- "search_S"

# For paid_media_spends set metric_value as your optimal spend

metric_value <- AllocatorCollect1$dt_optimOut[channels == select_media, optmSpendUnit]

# # For paid_media_vars and organic_vars, manually pick a value

# metric_value <- 10000

if (TRUE) {

optimal_response_allocator <- AllocatorCollect1$dt_optimOut[

channels == select_media, optmResponseUnit]

optimal_response <- robyn_response(

robyn_object = robyn_object,

select_build = 0,

media_metric = select_media,

metric_value = metric_value)

plot(optimal_response$plot)

if (length(optimal_response_allocator) > 0) {

cat("QA if results from robyn_allocator and robyn_response agree: ")

cat(round(optimal_response_allocator) == round(optimal_response$response), "( ")

cat(optimal_response$response, "==", optimal_response_allocator, ")\n")

}

支出に関して、目的関数がサチっているかどうかを見てみます。

新しいデータで、現在のモデルをアップデートします。

# Run ?robyn_refresh to check parameter definition
Robyn <- robyn_refresh(
  robyn_object = robyn_object
  , dt_input = dt_simulated_weekly
  , dt_holidays = dt_prophet_holidays
  , refresh_steps = 1
  , refresh_mode = "manual"
  , refresh_iters = 1000 # 1k is estimation. Use refresh_mode = "manual" to try out.
  , refresh_trials = 3
  , clusters = FALSE
)

# Export this refreshed model you wish to export
last_refresh_num <- sum(grepl('listRefresh', names(Robyn))) + 1 # Pick any refresh.
# Here's the final refresh using the model recommended by least combined normalized nrmse and decomp.rssd
ExportedRefreshModel <- robyn_save(
  robyn_object = robyn_object
  , select_model = Robyn[[last_refresh_num]]$OutputCollect$selectID
  , InputCollect = Robyn[[last_refresh_num]]$InputCollect
  , OutputCollect = Robyn[[last_refresh_num]]$OutputCollect
)

# Run ?robyn_refresh to check parameter definition

Robyn <- robyn_refresh(

robyn_object = robyn_object

, dt_input = dt_simulated_weekly

, dt_holidays = dt_prophet_holidays

, refresh_steps = 1

, refresh_mode = "manual"

, refresh_iters = 1000 # 1k is estimation. Use refresh_mode = "manual" to try out.

, refresh_trials = 3

, clusters = FALSE

)

# Export this refreshed model you wish to export

last_refresh_num <- sum(grepl('listRefresh', names(Robyn))) + 1 # Pick any refresh.

# Here's the final refresh using the model recommended by least combined normalized nrmse and decomp.rssd

ExportedRefreshModel <- robyn_save(

robyn_object = robyn_object

, select_model = Robyn[[last_refresh_num]]$OutputCollect$selectID

, InputCollect = Robyn[[last_refresh_num]]$InputCollect

, OutputCollect = Robyn[[last_refresh_num]]$OutputCollect

)

アップデートした場合、先ほどと同様に、推定結果や予算に応じたアロケーションを出力します。

# Run ?robyn_allocator to check parameter definition
AllocatorCollect <- robyn_allocator(
  robyn_object = robyn_object
  #, select_build = 1 # Use third refresh model
  , scenario = "max_response_expected_spend"
  , channel_constr_low = c(0.7, 0.7, 0.7, 0.7, 0.7)
  , channel_constr_up = c(1.2, 1.5, 1.5, 1.5, 1.5)
  , expected_spend = 2000000 # Total spend to be simulated
  , expected_spend_days = 14 # Duration of expected_spend in days
)
print(AllocatorCollect)

# Get response for 80k from result saved in robyn_object
Spend1 <- 60000
Response1 <- robyn_response(
  robyn_object = robyn_object
  #, select_build = 1 # 2 means the second refresh model. 0 means the initial model
  , media_metric = "search_S"
  , metric_value = Spend1)
Response1$response/Spend1 # ROI for search 80k
Response1$plot

# Run ?robyn_allocator to check parameter definition

AllocatorCollect <- robyn_allocator(

robyn_object = robyn_object

#, select_build = 1 # Use third refresh model

, scenario = "max_response_expected_spend"

, channel_constr_low = c(0.7, 0.7, 0.7, 0.7, 0.7)

, channel_constr_up = c(1.2, 1.5, 1.5, 1.5, 1.5)

, expected_spend = 2000000 # Total spend to be simulated

, expected_spend_days = 14 # Duration of expected_spend in days

)

print(AllocatorCollect)

# Get response for 80k from result saved in robyn_object

Spend1 <- 60000

Response1 <- robyn_response(

robyn_object = robyn_object

#, select_build = 1 # 2 means the second refresh model. 0 means the initial model

, media_metric = "search_S"

, metric_value = Spend1)

Response1$response/Spend1 # ROI for search 80k

Response1$plot

以上、チュートリアルを行いましたが、過去にMMMを実務で使ったことがあるものとしては、Robynはかなりオートマチックなツールだなぁと思いました。時系列のベイズモデリングに対してProhpetに感じた感情と似ているかもしれません。
パレート最適なものを見つけたり、アロケーションをどうするかを決めたりする関数までもが用意されており、適切にモデルを作成することさえできれば、データサイエンティストの業務時間をかなり削減することができると思います。
ただ、残存効果をカスタマイズしたり、独自のモデルをやる自由度はある程度犠牲にしていると思うので、当てはまりにこだわる場合、これまで通りStanなどで独自にアルゴリズムを書くこともあってしかるべきかなと思います。

参考情報

・https://facebookexperimental.github.io/Robyn/
・データ活用のための数理モデリング入門
・Robyn Open Source MMM Users
・Facebook Project Robyn Open Source MMM 2021 Community Summit
・https://facebookexperimental.github.io/Robyn/docs/quick-start/

Bayesian Methods for Media Mix Modeling with Carryover and Shape Effects[A4一枚まで備忘録]

A4用紙1枚にまとめるイメージでメモを残そうという取り組みです。

今回はメディアミックスモデリングに関して仕事があった際に参照できるメモの一つとして残します。
先日、Python/STAN Implementation of Multiplicative Marketing Mix Modelに参考文献としてBayesian Methods for Media Mix Modeling with Carryover and Shape Effectsの存在を知りました。時系列の重回帰の一種じゃないかと思っていたんですが、天下のGoogle様がどのようにマーケティングの成果を測ろうとしているのか、この論文から知ってみようと思います。

ちなみに、Python/STAN Implementation of Multiplicative Marketing Mix Modelの写経したものはMarketingMixModelingに残しています。

このコードでマーケティングにおけるチャネルごとのROAS(Return On Advertising Spend)を推定することができます。

以下の図はコードを回すことで得られるものですが、チャネルごとのROASの事後分布となります。

この図もコードを回すことで得られるものですが、広告などの支出に関する売上の関数をフィットしたものになります。精度高くこれらの関数を求めれば、過剰に広告に投資している可能性のあるチャネルを見つけることが可能になります。

Abstract

ROASやmROASに関するアトリビューション分析をベイズモデリングで行い、事後分布から効果を推定している。
広告関連のデータが、売上などに対して、数ヶ月後に遅れて効果が現れるキャリーオーバー効果の想定、支出を増やしても効果がなくなる収獲低減の想定などをモデルに置いている。
モデルの識別に関してはBICを評価基準として用いている。

Introduction

Media mix models(MMM)はメディア支出が売上にどのように影響を与えるのかを理解するために、最適なメディア投資を行うための支出の配分を決めるために使われる。
扱うデータとしては売上、物価、プロダクトの分布、様々なメディアの支出、マクロ経済や天候、季節性、競合の情報などの外部要因が含まれる。
RCTが難しいため、あるいは傾向スコアなどのモデリングもデータが少なく難しいため、回帰分析が行われることが多い。
線形回帰だと、広告がサチっている効果（広告が飽和している状態）や収獲低減の効果も表現できない。
広告にはキャリーオーバー効果と言って、売上への効果が遅れて現れるという可能性が広く信じられている。
事前分布に広告に関する様々なノウハウを反映させることで、データ数に比して推定するパラメータの数が多い分析に向きあう。

Model Specification

モデル構築には週次データが使われるケースが多い。
売上のデータ、メディア支出のデータ、コントロール変数として売上と関係してそうなデータが扱われる。コントロール変数は業界によって異なる。
キャリーオーバー効果
$$ adstock(x_{t-L+1, m}, \dots, x_{t, m} ; w_{m} ,L) = \frac{\sum_{t=0}^{L-1} w_m(l) x_{t-l, m} }{\sum_{l=0}^{L-1} w_m(l)} $$

$ w_m(l) $は非負のウェイトでラグ期間の関数になっている。
$ L $はキャリーオーバー効果の最大期間。
$ x_{t, m} $はある期間のあるメディアへの支出。
論文ではL=13とされている。
ウェイトの関数について、メディアの売上に対する効果は徐々に減っていくが、効果が遅れて生じることもあるので、効果のピークをずらした表現も扱われている。
ウェイトの関数で最初にピークがくる場合、

$$ w_{m}^{g} (l; \alpha_{m}) = \alpha_{m}^{l}, \\ l = 0, \dots, L-1 , \\ 0 < \alpha_m < 1 $$

ウェイトの関数で遅れてピークがくる場合、

$$ w_{m}^{d} (l; \alpha_{m}, \theta_m) = \alpha_{m}^{(l – \theta)^2}, \\ l = 0, \dots, L-1 , \\ 0 < \alpha_m < 1 ,\\ 0 \leq \theta_m \leq L -1 $$

となる。
$ \alpha_m $はメディアの効果の、ある期から次の期への維持率。
$ \theta_m $は遅れてピークが現れる程度に関するパラメータ。
形状効果
$$Hill(x_{t,m};\mathcal{K}_m, \mathcal{S}_m) = \frac{1}{1+ \left (
\frac{x_{t, m}}{\mathcal{K}_m} \right )^{- \mathcal{S}_m} }, \\ x_{t, m} \geq 0 $$

$ x_{t, m} $はある期間のあるメディアへの支出。
$ \mathcal{S}_m $は形状に関するパラメータで傾きに関わるもの。
$ \mathcal{K}_m $は形状に関するパラメータでサチるポイントに関わるもの。
分子に$ \mathcal{K}_m ^{\mathcal{S}_m} $を足して引くことで、

$$ \beta_{m} Hill_{m} (x_{t, m}) = \beta_{m} – \frac{\mathcal{K}_m^{\mathcal{S}_m}\beta_{m}}{x_{t,m}^{\mathcal{S}_m}+\mathcal{K}_m^{\mathcal{S}_m}} $$

に変換できる。
なお、$ \beta_m $は各メディアの回帰係数となる。これでもって各メディア支出の売上に対する効果を見れる。パラメータの大きさによって様々な形状をとりうる。
キャリーオーバー効果と形状効果をモデルで一緒に表現する。
シンプルにするためにメディア間のシナジー効果を無視している。

$$ y_t = \tau + \sum_{m=1}^{M} \beta_m Hill(x_{t, m }^{*}; \mathcal{K}_m,\mathcal{S}_m) + \sum_{c=1}^{C} \gamma_{c}z_{t, c} + \epsilon_t $$

$ y_t $は売上。
$ x_{t, m }^{*} = adstock(x_{t-L+1, m}, \dots, x_{t, m} ; w_{m}, L) $は前述の広告支出のキャリーオーバー効果を表現したもの。
$ \tau $は広告支出によらないベースラインの売上。
$ z_{t,c} $はコントロール変数。
$ \gamma_c $はコントロール変数の係数。
$ \epsilon_t $はホワイトノイズ。

Estimating the Bayesian Model

階層ベイズに拡張することもできる。
ハイパーパラメータはメディアの効果が非負であるとか、割合は0から1だとか色々と事前情報として設定する。

Attribution Metrics and Optimal Media Mix

経路ごとのROSAやmROASを計算してメディアの最適化をしたい。
メディアに投資することでの予測売上と、投資なかった場合の売上の差分でもって成果を測る。
ROASの定義式
$$ ROAS_{m} = \frac{\sum_{t_{0} \leq t \leq t_{1} + L – 1 } \hat Y_{t}^{m} (x_{t-L+1, m}, \dots, x_{t, m}; \Phi) – \hat Y_{t}^{m} (\tilde{x}_{t-L+1, m}, \dots, \tilde{x}_{t, m}; \Phi) }{\sum_{t_{0} \leq t \leq t_{1}} x_{t, m} } $$

$ \hat{Y} $は売上の予測。
$ \tilde{x} $はメディアの支出の変化。
mROASの定義式
$$ mROAS_{m} = \frac{\sum_{t_{0} \leq t \leq t_{1} + L – 1 } \hat Y_{t}^{m} (\tilde{\tilde{x}}_{t-L+1, m}, \dots, \tilde{\tilde{x}}_{t, m}; \Phi) – \hat Y_{t}^{m} (x_{t-L+1, m}, \dots, x_{t, m}; \Phi) }{ 0.01 \times \sum_{t_{0} \leq t \leq t_{1}} x_{t, m} } $$

$ \hat{Y} $は売上の予測。
$ \tilde{\tilde{x}} $はメディアの支出の1%の変化。
ROASの分布は事後分布と上述の定義式から得られる。
メディアミックスの最適化
予算という制約付きの極値問題を解く。ラグランジュ未定乗数法を用いる。

Application to a Simulated Data Set

シミュレーションしたデータを使っている。2年間の週次のデータで、売上、3つのメディア、一つのコントロール変数からなる。

Impact of Sample Size on Estimation Accuracy

単一のデータセットのみで評価することが難しいので、それぞれが独立した500個のデータセットを生成し、それぞれ推定をしている。

Choice of Priors

サンプルサイズが小さい場合、事後分布が事前分布の影響を大きく受けるため、事前分布に関して色々と考察している。
$ \beta $に関する事前分布。設定次第で結構変わる。メディアの支出にかかってくるので分析に与える影響は大きい。
$ \mathcal{K}_m $に関する事前分布。要はサチるポイントに関するもの。設定次第でグラフの挙動が結構変わる。

Application to Real Data and Model Selection

シャンプーの広告主の実際のデータを使って分析。
2.5年の週次データで、反応変数として売上、メディアデータとしてTV・雑誌・ディスプレイ・YouTube・検索、コントロール変数として小売のデータ（平均オンスあたり価格）・重み付けされた全商品の流通量や広告量が扱われている。
BICを頼りにモデルを選択した。BICは事後分布のデータから計算できる。

Conclusion

メディアミックスモデルはモデルの識別によって結果が大きく変わる。
ベイズモデルでキャリーオーバー効果や形状効果を考慮したモデルを求め、それらの事後分布からROASやmROASを計算した。
データ数が少ないとバイアスの伴った推定となる。事前情報の選定にも慎重になるべき。
BICでモデルを識別した結果、キャリーオーバー効果や形状効果を想定した設定が選ばれた。
残差に自己相関があったため、今回のモデルはまだ不完全で改善の余地がある。

2019年に読んだデータ分析系の本の振り返り（21+1冊）

はじめに

2020年、あけましておめでとうございます。年末に自分自身を振り返ろうと思ったのですが、結局データ分析と勉強しかしていないわけで、書籍を振り返ろうと思うに至りました。私の知り合いのデータサイエンティストはだいたい全冊持っているであろうと思われますが、良い本だと思うので思い出していただければ幸いです。

1.『ベイズモデリングの世界』（岩波書店）

基本的に階層ベイズモデルを使って、個体ごとの異質性を考慮した分析手法が提案されています。前半はオムニバス形式で様々な先生がモデルの適用について執筆されており、後半では伊庭先生による階層ベイズモデルの講義になっています。途中でスタイン統計量による縮小推定の話があげられ、柔軟なモデリングのためには「階層化した方が少なくとも望ましい推定量が得られる」という数学的証明を捨てることもやむを得ないと書かれています。

2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』（コロナ社）

この本はトピックモデルの教科書というよりも、ベイズ推定の教科書という側面が強い印象があります。途中で出てくる数式は流し読みするのは難しく、最低2冊以上のノートが別途必要になると思います。一度でもLDAのパラメータを導出してみたいという方には良い教科書だと思います。疑似コードが提供されているので、それをもとにRやPythonでコーディングしていけば、一番シンプルなLDAが非常に短い行で実行できてしまうことに驚かれるかもしれません。人間が手を動かして推定アルゴリズムを導出しているからこそ、短いコードで済むということを実感できるはずです。

3.『構造的因果モデルの基礎』（共立出版）

グラフィカルなアプローチで因果推論を扱っている書籍です。Judea Pearl流の因果推論アプローチについて記すことを目的に書かれています。基礎と書かれていますが決して簡単ではありません。ただ、扱われる数学のレベルとしては確率と線形代数がわかれば大丈夫だと思われます。余談ではありますが、1章の相関関係と因果関係の事例紹介で「おむつとビールの話」が都市伝説ではなくきちんと記事としてWall Street Journalという雑誌に掲載されていたことが明らかにされています。

4.『現場で使える！PyTorch開発入門深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』（翔泳社）

PyTorchを触ったことがないが、深層学習の手法について知っている層を対象とした本です。6章まではGoogleのColabで動かせるのでGoogleに課金することなく深層学習による回帰、CNN、GAN、RNN、Encoder-Decoderモデル、ニューラル行列因子分解をPyTorchで試すことができます。写経したものはこちら。転移学習や高解像度化や画像生成、文章のクラス分類、文書生成、機械翻訳などもできるので、PyTorchでこれくらいの量をコーディングしたらこれくらいのことができるのかという学びや、他の人の書いたPyTorchコードを読みやすくなるなどの便益は十分にあると思いました。

5.『作ってわかる! アンサンブル学習アルゴリズム入門』（シーアンドアール研究所）

会社で行っているPythonもくもく会用に買った本で、scikit-learnを使わずに機械学習のアルゴリズム（アンサンブル系）をコーディングするための本です。pythonのコードについて逐次、細かい解説が行われているわけではないので、1行1行自分でコメントを加えながら写経をしていけば力が付くという本かなと思われます。sklearnはそれはそれで素晴らしいですが、こういう本でフルスクラッチで修行できるのはいいですね。

6.『数理統計学―基礎から学ぶデータ解析』（内田老鶴圃）

統計検定1級を合格された方のブログで紹介されていた教科書です。理系の大学生レベルの数学知識があれば、数理統計学の基礎を学べると思います。中心極限定理の証明や、様々な分布の期待値や分散、様々な分布の性質について数式を用いてしっかり理解することができます。数式もほどよく端折られているので、無論ですがノートが数冊必要になります。各章毎にある練習問題も解くことで力が付くと思います。日本の大学の授業の教科書がこれだったらジェノサイド（再履修者の大量発生）が起きるんだろうなと思ってしまった。

7.『44の例題で学ぶ統計的検定と推定の解き方』（オーム社）

統計の検定に関してだけ扱った珍しい本です。第3部までは統計学の普通の教科書ですが、それ以降であらゆる検定の例題が44件も載せられています。パラメトリックな検定から、ノンパラメトリックな検定まで幅広く扱われています。一番気にいっているのは仮説検定法の分類の表です。これさえあれば、どのデータに対してどの検定を行えばいいかが一目瞭然です。

8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』（技術評論社）

機械学習の原理を手早く数式を交えて学べる本です。かゆいところに手が届いていると言うか、既出の教科書では捨象されがちな、条件付き確率における2変数以上の条件づけでの表現に紙面を割いていたりしてくれるのが嬉しいです。ある程度数学の話はわかるが、だいぶ忘れているビジネスパーソンには大変にありがたいコンテンツと言えると思います。ベイズ線形回帰に関しても行列を用いた、わかりやすい導出方法が紹介されています。またコラムで紹介されている、測度論にどう向き合えばいいかの著者の見解は参考になります。

9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)』

R言語とstanを用いてベイズ統計学を入門レベルから学べる本です。各トピックごとにそれなりの紙面が割かれています。例題も豊富にあるので、線形回帰・MCMC・情報量基準・階層ベイズモデルまで、ベイズ統計学を基礎から応用までしっかりと学べると思います。youtubeで著者の講義も配信されているので、留学気分を味わえます。

10.『scikit-learnとTensorFlowによる実践機械学習』（オライリージャパン）

2019年に日本で開かれたML SummitでTFの開発者がおすすめしていた教科書です。前半部分で機械学習の入門から応用までをわかりやすい説明で学ぶことができます。数式は少ないですが、図とソースコード（Python）がちりばめられており、手を動かして理解を進めることができます。後半部分はTensorFlowを用いた深層学習の基礎を同様に手を動かして学ぶことができます。ただ、TFのバージョンも変わってきているので前半の説明をアテにして読むのも良いと思います。

11.『AIアルゴリズムマーケティング自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)』

マーケティングへのデータサイエンスの適用に関する珍しい書籍です。ソースコードはついていないですが、業務で使う際のアイデアが手に入ることもあります。一般的な回帰、生存時間分析、オークション、アトリビューション分析、アップリフトモデリング以外にも、情報検索やレコメンデーションやトピックモデルなどマーケティングながら学際的なトピックも扱われています。レコメンドなどで使われる、ランク学習に関して詳しく書かれた書籍をあまり知らないので、この本はその点においてもありがたい本でもあります。

12.『入門統計的因果推論』（朝倉書店）

ほぼ全ての章でグラフィカルなアプローチで因果推論を扱っています。例題も豊富なので、一つ一つ丁寧にやれば理解が捗ります。おそらく、例題の多さを含め一番丁寧にd分離性、do演算子、バックドア基準、フロントドア基準に関する説明をしてくれている本なのかなと思いました。グラフでの因果推論に関して初めての人でも、確率さえ知っていれば読み進めることができるはずです。また、途中で操作変数法の紹介もされ、経済学出身者としては読みやすい。ただ、傾向スコアのくだりや、DIDなどのくだりはあまり出てきません。あと、やってないですが章末の練習問題に対するSolution Manualが提供されているようです。

13.『実践ベイズモデリング -解析技法と認知モデル-』（朝倉書店）

ベイズモデリングを様々な事例に適用する方法がオムニバス形式で記された本です。ワイブル分布や異質性を考慮した二項分布、無制限複数選択形式のアンケートデータに対する手法、トピックモデル、項目反応理論などが扱われています。マーケティングの実務で使える事例が多いように感じました。こちらはサポートサイトでRコードとstanコードが提供されています。あと、appendixにあるプレート表現の見方も参考になります。

14.『機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)』

機械学習などで用いるベイズ推論を扱った教科書です。入門とありますが、入門者は書かれた数式をそのまま見ていても頭に入らないのではないでしょうか。手を動かしてなんぼの本だと思います。ノート2冊は絶対に必要です。たぶん、数式の展開を丁寧に記すと倍以上の厚みの本になると思います。各々のモデルに関してグラフィカルモデルが記されているのや、サンプルコードとしてGitHubにJuliaで書かれたソースコードが提供されているのも良いです。

15.『その問題、数理モデルが解決します』（ベレ出版）

物語形式で、様々な問題に対して数理モデリングのアプローチが紹介されています。途中でマッチング理論やゲーム理論やオークションなども登場することから、経済学出身者も喜ぶ内容かもしれません。社会人になってからナッシュ均衡という言葉が書かれた本は中々出会って来なかった。

16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』（CCCメディアハウス）

2013年と結構古い本ですが、データ分析を様々な事象に対して適用した事例紹介本です。アップリフトモデリングへの言及もあり、こういったものに関して日本は何年も遅れてブームが来るんだなという実感を与えてくれた本でもありました。appendixに分析事例が147個ほどあげられているのも参考になります。

17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』（北大路書房）

主にstanを用いたベイズモデリングによる分析事例が1と2で38本もオムニバス形式で載っています。ほとんどの事例で階層ベイズモデルが扱われています。2では若干マーケティングに近い内容の題材も扱われ、データサイエンティストの人にも嬉しい内容かもしれません。もちろんデータとstanとRのコードがサポートサイトで提供されています。

18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』（共立出版）

状態空間モデルで時系列予測を行うための手法が記されている本です。RのKFASパッケージが全面に渡って扱われています。トレンドを考慮したり、カレンダー効果を追加したり、共変量を追加したりなど様々なアプローチが紹介されコードも伴っているわけですから、業務でも抜群に役に立ちました。

19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』（SBクリエイティブ）

自分のいる会社で最低限の数学がわかると思われる若いメンバーに買ってもらうように言っている本です。微積分・線形代数だけでなく、カルシュ・キューン・タッカー条件（最適化数学）に関しても扱ってくれているので、ここで出てくる数学がわかれば大体の論文に立ち向かえると思います。さらに、Pythonの基礎もこれで学ぶことができるので一石二鳥な素敵な本ですね。また、最後の方でスクラッチでアルゴリズムを書くパートがあり、こちらも勉強になります。

20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

機械学習における前処理の指針を与えてくれる本です。Pythonのコードが提供されています。例えばですが、「テストデータにだけある、新しい単語は取り除いてしまえばいい」などの細かいアドバイスが何気に嬉しいです。「Effectコーディング」「特徴量ハッシング」「ビンカウンティング」「バックオフ」「leakage-proof統計量」などは読むまで知らないところだったので勉強になりました。

21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』（オライリージャパン）

データ分析の仕事をする上で最低限必要な知識を幅広く抑えることができる本です。数式は少なく、ところどころ出てくるコードはR言語です。参考文献などがブログだったりするため厳密さがめちゃあるわけではないですが、業務で使う分には問題ないと思います。分類問題において、AUCなどの評価指標だけでなく、予測値自体の探索的分析のすすめなどが書かれており、参考になりました。また、特徴量エンジンとしてのk-NN法の話も面白いと思いました。

[+α]『プログラマのためのGoogle Cloud Platform入門サービスの全体像からクラウドネイティブアプリケーション構築まで』（翔泳社）

Google Cloud Platformを初めて触るデータ分析者にはちょうど良い本です。説明もわかりやすいので、いきなりアカウントを作ってドキュメントを解読するよりかは戸惑いは減るはずです。この本を土台に、GCS・GCEを駆使してML系のAPIを呼び出して使うなどの最低限の操作は私でもできるようになりました。GCPの画面や機能もどんどん変わっていくので書籍を買ってもアレなんですが、歴史的な背景も若干記述されているので、それはそれで勉強になります。ただ、エンジニアにこの本を買うべきか聞いた際にネガティブな意見があったのですが、たぶん現役プログラマからすると簡単過ぎるからなんだろうなと思います。

終わりに

2019年もぼちぼち勉強できましたが、2020年もこれまで同様にノートとペンを大事にする勉強を続けていき、コーディングも分析ももっともっと数をこなして会社や社会に求められるようなデータ分析官を目指していこうと思います。あぁ、英会話などの勉強をする時間を作るのが難しい。

[Stan]ロジスティック回帰の階層ベイズモデルとk-foldsクロスバリデーション

はじめに

stanは意思決定のための分析などでのパラメータ推定に使うことが多く、機械学習のために扱うことはありませんでした。ただ、過去にリク面などでお話したデータサイエンティストの方はstanで機械学習していて、クロスバリデーションもしているとの発言をされていました。
先日、記事を漁っていたらstanでクロスバリデーションを行うためのコードを書いている方を見つけたので、その方のコードをもとに大人のirisであるwineデータを用いて、質の高いワインかどうかを分類するために階層ベイズでのロジスティック回帰モデルを回してみたいと思います。

データについて

UCI Machine Learning Repositoryにある、赤ワインの評価と成分のデータです。データに関する説明はワインの味（美味しさのグレード）は予測できるか？（１）で丁寧になされていますので、ご確認ください。今回は6点以上であれば1を、そうでなければ0を取るものを教師データとしています。

分析方針

今回は階層ベイズモデルを扱うことから、グループごとにロジスティック回帰のパラメータが異なるという仮定を置きます。そのために、citric.acidというデータを3つのカテゴリデータに変換して、それをグループとして扱います。モデルでは、今回のデータセットの変数を全て回帰項として使います。もちろん、回帰用の式からはcitric.acidは除外します。
全体の80%を訓練データに、20%をテストデータとして、10foldsクロスバリデーションでのstanによる予測結果の平均AUCを評価指標とします。最後に、テストデータを用いた予測のAUCを確かめます。また、階層ベイズモデルではないモデルでの10foldsクロスバリデーションでのAUCとも比較します

分析概要

・データ整形
・訓練データとテストデータの分割
・クロスバリデーション用のデータの作成
・stanの実行
・クロスバリデーション結果の出力
・テストデータでの予測
・非階層モデルとの比較

全体のコード以下のリンクにあります。
kick_logistic_regression_allowing_k_hold_cross_validation_hierachical.R
logistic_regression_allowing_k_fold_cross_validation_hierachical.stan

データ整形

階層ベイズで扱うグループをcitric.acidから作っています。

library(tidyverse)
library(rstan)
library(GGally)
library(shinystan)
library(pbmcapply)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("dataset/winequality-red.csv")

# Visualization -----------------------------------------------------------
ggpairs(wine_dataset)


# Making group -----------------------------------------------------------
wine_dataset <- wine_dataset %>% mutate(citric_acid_group = if_else(citric.acid < 0.2, 1,
                                                                    if_else(citric.acid < 0.4, 2, 3)))

wine_dataset <- wine_dataset %>% select(-citric.acid)

library(tidyverse)

library(rstan)

library(GGally)

library(shinystan)

library(pbmcapply)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("dataset/winequality-red.csv")

# Visualization -----------------------------------------------------------

ggpairs(wine_dataset)

# Making group -----------------------------------------------------------

wine_dataset <- wine_dataset %>% mutate(citric_acid_group = if_else(citric.acid < 0.2, 1,

if_else(citric.acid < 0.4, 2, 3)))

wine_dataset <- wine_dataset %>% select(-citric.acid)

訓練データとテストデータの分割

ワインの質に関するバイナリーデータをこちらで作成し、80%を訓練データに、20%をテストデータに分割しています。

# Split Train and Test ----------------------------------------------------
smp_size <- floor(0.8 * nrow(wine_dataset))

## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(wine_dataset)), size = smp_size)

train <- wine_dataset[train_ind, ]
test <- wine_dataset[-train_ind, ]

y <- if_else(wine_dataset$quality > 5, 1, 0)
x <- as.matrix(wine_dataset %>% select(-quality))

y_train <- y[train_ind]
y_test <- y[-train_ind]

x_train <- x[train_ind,]
x_test <- x[-train_ind,]

x_train_group <- x_train[,ncol(x_train)]
x_test_group <- x_test[,ncol(x_test)]

x_train <- x_train[,1:(ncol(x_train)-1)]
x_test <- x_test[,1:(ncol(x_test)-1)]

x_train <- scale(x_train)
x_test <- scale(x_test)

# Split Train and Test ----------------------------------------------------

smp_size <- floor(0.8 * nrow(wine_dataset))

## set the seed to make your partition reproducible

set.seed(123)

train_ind <- sample(seq_len(nrow(wine_dataset)), size = smp_size)

train <- wine_dataset[train_ind, ]

test <- wine_dataset[-train_ind, ]

y <- if_else(wine_dataset$quality > 5, 1, 0)

x <- as.matrix(wine_dataset %>% select(-quality))

y_train <- y[train_ind]

y_test <- y[-train_ind]

x_train <- x[train_ind,]

x_test <- x[-train_ind,]

x_train_group <- x_train[,ncol(x_train)]

x_test_group <- x_test[,ncol(x_test)]

x_train <- x_train[,1:(ncol(x_train)-1)]

x_test <- x_test[,1:(ncol(x_test)-1)]

x_train <- scale(x_train)

x_test <- scale(x_test)

クロスバリデーション用のデータの作成

こちらのコードでは任意の数でクロスバリデーション用のデータを作成し、stanで扱う訓練用データのlistに追加しています。
また、参考にしているブログより転用したstan_kfoldという関数を定義しています。k分割した際のstanの推定結果をリストに格納するための関数です。

N <- length(y_train) # sample size
n_fold <- 10 # number of folds
K <- ncol(x_train) #n umber of predictors

# Create cross validation data --------------------------------------------
# create 10 folds of data
hh <- kfold_split_random(n_fold, N) #hh index the fold ID of each data point
holdout_k <- matrix(0, nrow = N, ncol = n_fold)
for(i in 1:N) holdout_k[i, hh[i]] <- 1

# turn into a list
holdout_k <- split(holdout_k,rep(1:ncol(holdout_k),each=nrow(holdout_k)))


# the basic data object
data_m <- list(N=N,
               K=K,
               M = length(unique(x_train_group)),
               citric_acid_group = x_train_group,
               X=x_train,
               y=y_train
               )

# create a list of data list
data_l <- rep(list(data_m),n_fold)
# add the holdout index to it
for(i in 1:n_fold) data_l[[i]]$holdout <- holdout_k[[i]]


# Define function ---------------------------------------------------------
# function to parrallelize all computations
# need at least two chains !!!
stan_kfold <- function(file, list_of_datas, chains, cores,...){
  library(pbmcapply)
  badRhat <- 1.1 # don't know why we need this?
  n_fold <- length(list_of_datas)
  model <- stan_model(file=file)
  # First parallelize all chains:
  sflist <- 
    pbmclapply(1:(n_fold*chains), mc.cores = cores, 
               function(i){
                 # Fold number:
                 k <- ceiling(i / chains)
                 s <- sampling(model, data = list_of_datas[[k]], 
                               chains = 1, chain_id = i)
                 return(s)
               })
  
  # Then merge the K * chains to create K stanfits:
  stanfit <- list()
  for(k in 1:n_fold){
    inchains <- (chains*k - (chains - 1)):(chains*k)
    #  Merge `chains` of each fold
    stanfit[[k]] <- sflist2stanfit(sflist[inchains])
  }  
  return(stanfit) 
}

N <- length(y_train) # sample size

n_fold <- 10 # number of folds

K <- ncol(x_train) #n umber of predictors

# Create cross validation data --------------------------------------------

# create 10 folds of data

hh <- kfold_split_random(n_fold, N) #hh index the fold ID of each data point

holdout_k <- matrix(0, nrow = N, ncol = n_fold)

for(i in 1:N) holdout_k[i, hh[i]] <- 1

# turn into a list

holdout_k <- split(holdout_k,rep(1:ncol(holdout_k),each=nrow(holdout_k)))

# the basic data object

data_m <- list(N=N,

K=K,

M = length(unique(x_train_group)),

citric_acid_group = x_train_group,

X=x_train,

y=y_train

)

# create a list of data list

data_l <- rep(list(data_m),n_fold)

# add the holdout index to it

for(i in 1:n_fold) data_l[[i]]$holdout <- holdout_k[[i]]

# Define function ---------------------------------------------------------

# function to parrallelize all computations

# need at least two chains !!!

stan_kfold <- function(file, list_of_datas, chains, cores,...){

library(pbmcapply)

badRhat <- 1.1 # don't know why we need this?

n_fold <- length(list_of_datas)

model <- stan_model(file=file)

# First parallelize all chains:

sflist <-

pbmclapply(1:(n_fold*chains), mc.cores = cores,

function(i){

# Fold number:

k <- ceiling(i / chains)

s <- sampling(model, data = list_of_datas[[k]],

chains = 1, chain_id = i)

return(s)

})

# Then merge the K * chains to create K stanfits:

stanfit <- list()

for(k in 1:n_fold){

inchains <- (chains*k - (chains - 1)):(chains*k)

# Merge `chains` of each fold

stanfit[[k]] <- sflist2stanfit(sflist[inchains])

}

return(stanfit)

}

stanの実行

こちらのstanのコードでは、M個のグループごとにパラメータが異なるというモデルを書いています。modelブロックの途中でholdoutを入れることで一部のデータを推定に使わないようにしています。

data {
  int<lower=0> N;   // number of data items
  int<lower=0> K;   // number of predictors
  int<lower=0> M;   // number of group
  row_vector[K] X[N];
  int<lower=0,upper=1> y[N];
  int<lower=0,upper=M> citric_acid_group[N];
  int<lower=0, upper=1> holdout[N];
  // index whether the observation should be held out (1) or used (0)
}

parameters {
  real mu[K];
  real<lower=0> sigma[K];
  vector[K] beta[M];
}

model {
  for (k in 1:K) {
    mu[k] ~ normal(0, 100);
    sigma[k] ~ inv_gamma(1, 1);

    for (m in 1:M)
      beta[m,k] ~ normal(mu[k], sigma[k]);
  }
  for (n in 1:N){
    if(holdout[n] == 0){
      target += bernoulli_lpmf( y[n] | inv_logit(X[n] * beta[citric_acid_group[n]]));
  }
 }
}

data {

int<lower=0> N; // number of data items

int<lower=0> K; // number of predictors

int<lower=0> M; // number of group

row_vector[K] X[N];

int<lower=0,upper=1> y[N];

int<lower=0,upper=M> citric_acid_group[N];

int<lower=0, upper=1> holdout[N];

// index whether the observation should be held out (1) or used (0)

}

parameters {

real mu[K];

real<lower=0> sigma[K];

vector[K] beta[M];

}

model {

for (k in 1:K) {

mu[k] ~ normal(0, 100);

sigma[k] ~ inv_gamma(1, 1);

for (m in 1:M)

beta[m,k] ~ normal(mu[k], sigma[k]);

}

for (n in 1:N){

if(holdout[n] == 0){

target += bernoulli_lpmf( y[n] | inv_logit(X[n] * beta[citric_acid_group[n]]));

}

こちらはstanをキックするためのコードです。いつもと違い、先程定義したstan_kfoldを用いています。

# Kick the stan code ------------------------------------------------------
# run the functions
ss <- stan_kfold(file="model/logistic_regression_allowing_k_fold_cross_validation_hierachical.stan",
                 data_l,
                 chains=4,
                 cores=2)

# Kick the stan code ------------------------------------------------------

# run the functions

ss <- stan_kfold(file="model/logistic_regression_allowing_k_fold_cross_validation_hierachical.stan",

data_l,

chains=4,

cores=2)

クロスバリデーション結果の出力

以下は、k個ずつ手に入ったクロスバリデーションでの推定結果から、各パラメータの平均値を計算し、ロジスティック回帰モデルで２値の予測を行い、平均AUCを計算するコードです。

# Calculate Mean AUC ------------------------------------------------------
## hierachical model
set.seed(123)
cv_mean_auc <- NULL
for (i in 1:n_fold){
  ext_fit <- extract(ss[[i]]) # choose 1 chunk
  
  # グループによって推定したパラメータが違う
  coef_list <- NULL
  group_list <- data_l[[i]]$citric_acid_group[data_l[[i]]$holdout > 0]
  for(j in 1:length(group_list)){
    coef_list <- rbind(coef_list, colMeans(ext_fit$beta[,group_list[j],]))
  }
  
  lin_comb <- rowSums(data_l[[i]]$X[data_l[[i]]$holdout > 0, ] * coef_list)
  prob <- 1/(1 + exp(-lin_comb))
  pred_value <- rbinom(sum(data_l[[i]]$holdout), 1, prob)
  
  # Syntax (response, predictor):
  auc = pROC::auc(data_l[[i]]$y[data_l[[i]]$holdout > 0], pred_value)[1]
  cv_mean_auc <- append(cv_mean_auc, auc)
}

cv_mean_auc
mean(cv_mean_auc)

# Calculate Mean AUC ------------------------------------------------------

## hierachical model

set.seed(123)

cv_mean_auc <- NULL

for (i in 1:n_fold){

ext_fit <- extract(ss[[i]]) # choose 1 chunk

# グループによって推定したパラメータが違う

coef_list <- NULL

group_list <- data_l[[i]]$citric_acid_group[data_l[[i]]$holdout > 0]

for(j in 1:length(group_list)){

coef_list <- rbind(coef_list, colMeans(ext_fit$beta[,group_list[j],]))

}

lin_comb <- rowSums(data_l[[i]]$X[data_l[[i]]$holdout > 0, ] * coef_list)

prob <- 1/(1 + exp(-lin_comb))

pred_value <- rbinom(sum(data_l[[i]]$holdout), 1, prob)

# Syntax (response, predictor):

auc = pROC::auc(data_l[[i]]$y[data_l[[i]]$holdout > 0], pred_value)[1]

cv_mean_auc <- append(cv_mean_auc, auc)

}

cv_mean_auc

mean(cv_mean_auc)

平均AUCは0.675となりました。すごくいいわけではないですが、手抜きモデルとしてはまずまずと言ったところでしょうか。

テストデータでの予測

以下のコードで最初に分けていたテストデータでの予測結果を返します。

# Prediction --------------------------------------------------------------
# Choose best model in cross-validation
ext_fit <- extract(ss[[1]]) # choose 1 chunk

# グループによって推定したパラメータが違う
beta_post <- NULL
group_list <- x_test_group

for(j in 1:length(group_list)){
  beta_post <- rbind(beta_post, colMeans(ext_fit$beta[,group_list[j],]))
}

lin_comb <- rowSums(x_test * beta_post)
prob <- 1/(1 + exp(-lin_comb))
pred_value <- rbinom(nrow(x_test), 1, prob)


# Syntax (response, predictor):
auc = pROC::auc(y_test, pred_value)[1]
auc

# Prediction --------------------------------------------------------------

# Choose best model in cross-validation

ext_fit <- extract(ss[[1]]) # choose 1 chunk

# グループによって推定したパラメータが違う

beta_post <- NULL

group_list <- x_test_group

for(j in 1:length(group_list)){

beta_post <- rbind(beta_post, colMeans(ext_fit$beta[,group_list[j],]))

}

lin_comb <- rowSums(x_test * beta_post)

prob <- 1/(1 + exp(-lin_comb))

pred_value <- rbinom(nrow(x_test), 1, prob)

# Syntax (response, predictor):

auc = pROC::auc(y_test, pred_value)[1]

auc

実行の結果、AUCは0.665と、クロスバリデーションでの平均AUCと比べてあまり下がりませんでした。

非階層モデルとの比較

非階層モデルでも同様に10foldsクロスバリデーションの平均AUCを計算しました。非階層モデルよりもAUCが1%ポイントくらいは高いようです。

> mean(cv_mean_auc)
[1] 0.6745282

> mean(cv_mean_auc_normal)
[1] 0.6640103

> mean(cv_mean_auc)

[1] 0.6745282

> mean(cv_mean_auc_normal)

[1] 0.6640103

おわりに

現時点において、stanでの柔軟なモデリングを機械学習に活かす作法について紹介されている文献はあまりなく、選手人口はどれくらいいるのか気になるところですが、発見したブログのやり方でクロスバリデーションをカジュアルに行えるので、より多くの方がstanでの機械学習にチャレンジしうるものだなと思いました。ただ、このレベルの階層ベイズだとrstanarmで簡単にできてしまうので、より深く分析してモデルをカスタムしていきたいですね。

参考情報

[1]Lionel Hertzog (2018), “K-fold cross-validation in Stan,datascienceplus.com”
[2]Alex Pavlakis (2018), “Making Predictions from Stan models in R”, Medium
[3]Richard McElreath (2016), “Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)”, Chapman and Hall/CRC
[4]松浦健太郎 (2016), 『StanとRでベイズ統計モデリング (Wonderful R)』, 共立出版
[5]馬場真哉 (2019), 『実践Data Scienceシリーズ RとStanではじめるベイズ統計モデリングによるデータ分析入門』, 講談社

[Stan]生存時間分析のコードと便利なデータセットについて

はじめに

仕事で生存時間分析を使うことは結構あるのですが、マーケティングの良いデータセットがない印象でブログにしにくいと感じていました。また、Stanでの生存時間分析の事例もあまり把握していません。そこで使えそうなデータセットやStanのコードを探して、そのデータに対して生存時間分析を適用してみたいと思います。

目次
・生存時間分析とは
・生存時間分析で使えるデータ
・生存時間分析をマーケティングで使う際の用途
・先行研究
・生存時間分析で使えるデータセット
・Stanでの実行例
・おわりに
・参考文献

生存時間分析とは

生存時間分析は、ある時点から興味のあるイベント（マーケティングだと解約など）が発生するまでの期間を分析対象としています。データを手に入れた時点で、すでに解約して、真の累積の契約期間が判明している人と、解約しておらず今後いつ解約するかわからない中での累積の契約期間が残されている人のようなデータを扱うことが多いです。ここでの後者をcensoring（打ち切り）されたデータと呼びます。

生存時間分析をマーケティングで使う際の用途

ブログなどを読み漁る限り、以下の用途で生存時間分析を活用できるようです。

顧客のサービス離脱率の予測、離脱原因の特定
顧客がマーケティングキャンペーンに反応するまでの期間の長さ
故障率の予測、故障原因の特定

先行研究

Stanを用いた分析事例は、調べた限りですが以下のモデルがあるようです。

指数分布のモデル
Weibull（ワイブル）分布による比例ハザードモデル
ハザードの対数値についてのランダムウォークモデル
2階差分のマルコフ場モデル（生存時間の確率分布は正規分布）
1階差分のランダムウォークモデル（生存時間の確率分布は正規分布）

生存時間分析で使えるデータセット

事例を調べる過程で見つけた、生存時間分析に適したデータセットは以下の通りです。

RのMASSパッケージに含まれているgehan
Rのsurvivalパッケージに含まれているleukemia（白血病）
Rのsurvivalパッケージに含まれているveteran
Pythonのlifelineパッケージに含まれているrossi（逮捕された人が再逮捕されるまでの期間と共変量）
kaggleで提供されているEmployee Attrition Can you forecast employee attrition?のデータセット
Princeton Divorce Studyのデータ
IBMが提供している、Using Customer Behavior Data to Improve Customer Retentionという電話会社の解約に関すCRMデータ

どうやら、マーケティング、HR、離婚、再犯と幅広いデータがオープンソースで手に入るようです。

Stanでの実行例

今回は、「Using Customer Behavior Data to Improve Customer Retention」のデータセットを用いて、先行研究のソースコードにより生存時間分析をしてみようと思います。データは電話会社の顧客の解約に関するもので、様々な顧客の履歴データなどが用意されています。
先行研究のソースコードはWeibull分布を想定した比例ハザードモデルです。今回は決済の電子化の有無と離脱の関係を確かめてみます。なお、今回の打ち切りデータは契約期間となります。

まずはStanのコードはこちらです。Xobs_bgに説明変数が来るようにデータを用意しておきます。

functions { 
   vector sqrt_vec(vector x) { 
     vector[dims(x)[1]] res; 
  
     for (m in 1:dims(x)[1]){ 
       res[m] = sqrt(x[m]); 
     } 
  
     return res; 
   } 
  
   vector bg_prior_lp(real r_global, vector r_local) { 
     r_global ~ normal(0.0, 10.0); 
     r_local ~ inv_chi_square(1.0); 
  
     return r_global * sqrt_vec(r_local); 
   } 
 } 
  
data { 
   int<lower=0> Nobs; 
   int<lower=0> Ncen; 
   int<lower=0> M_bg; 
   vector[Nobs] yobs; 
   vector[Ncen] ycen; 
   matrix[Nobs, M_bg] Xobs_bg; 
   matrix[Ncen, M_bg] Xcen_bg; 
 } 
  
transformed data { 
   real<lower=0> tau_mu; 
   real<lower=0> tau_al; 
  
   tau_mu = 30.0; 
   tau_al = 30.0; 
 } 
  
parameters { 
   real<lower=0> tau_s_bg_raw; 
   vector<lower=0>[M_bg] tau_bg_raw; 
  
   real alpha_raw; 
   vector[M_bg] beta_bg_raw; 
  
   real mu; 
 } 
  
transformed parameters { 
   vector[M_bg] beta_bg; 
   real alpha; 
  
   beta_bg = bg_prior_lp(tau_s_bg_raw, tau_bg_raw) .* beta_bg_raw; 
   alpha = exp(tau_al * alpha_raw); 
 } 
  
model { 
   yobs ~ weibull(alpha, exp(-(mu + Xobs_bg * beta_bg)/alpha)); 
   target += weibull_lccdf(ycen | alpha,
                            exp(-(mu + Xcen_bg * beta_bg)/alpha)); 
  
   beta_bg_raw ~ normal(0.0, 1.0); 
   alpha_raw ~ normal(0.0, 1.0); 
  
   mu ~ normal(0.0, tau_mu); 
 } 
  
generated quantities { 
     real yhat_uncens[Nobs + Ncen]; 
     real log_lik[Nobs + Ncen]; 
     real lp[Nobs + Ncen]; 
  
     for (i in 1:Nobs) { 
         lp[i] = mu + Xobs_bg[i,] * beta_bg; 
         yhat_uncens[i] = weibull_rng(alpha,
                                      exp(-(mu + Xobs_bg[i,] * beta_bg)/alpha)); 
         log_lik[i] = weibull_lpdf(yobs[i] | alpha,
                                    exp(-(mu + Xobs_bg[i,] * beta_bg)/alpha)); 
     } 
     for (i in 1:Ncen) { 
         lp[Nobs + i] = mu + Xcen_bg[i,] * beta_bg; 
         yhat_uncens[Nobs + i] = weibull_rng(alpha,
                                              exp(-(mu + Xcen_bg[i,] * beta_bg)/alpha)); 
         log_lik[Nobs + i] = weibull_lccdf(ycen[i] | alpha,
                                              exp(-(mu + Xcen_bg[i,] * beta_bg)/alpha)); 
     } 
 }

functions {

vector sqrt_vec(vector x) {

vector[dims(x)[1]] res;

for (m in 1:dims(x)[1]){

res[m] = sqrt(x[m]);

}

return res;

}

vector bg_prior_lp(real r_global, vector r_local) {

r_global ~ normal(0.0, 10.0);

r_local ~ inv_chi_square(1.0);

return r_global * sqrt_vec(r_local);

}

data {

int<lower=0> Nobs;

int<lower=0> Ncen;

int<lower=0> M_bg;

vector[Nobs] yobs;

vector[Ncen] ycen;

matrix[Nobs, M_bg] Xobs_bg;

matrix[Ncen, M_bg] Xcen_bg;

}

transformed data {

real<lower=0> tau_mu;

real<lower=0> tau_al;

tau_mu = 30.0;

tau_al = 30.0;

}

parameters {

real<lower=0> tau_s_bg_raw;

vector<lower=0>[M_bg] tau_bg_raw;

real alpha_raw;

vector[M_bg] beta_bg_raw;

real mu;

}

transformed parameters {

vector[M_bg] beta_bg;

real alpha;

beta_bg = bg_prior_lp(tau_s_bg_raw, tau_bg_raw) .* beta_bg_raw;

alpha = exp(tau_al * alpha_raw);

}

model {

yobs ~ weibull(alpha, exp(-(mu + Xobs_bg * beta_bg)/alpha));

target += weibull_lccdf(ycen | alpha,

exp(-(mu + Xcen_bg * beta_bg)/alpha));

beta_bg_raw ~ normal(0.0, 1.0);

alpha_raw ~ normal(0.0, 1.0);

mu ~ normal(0.0, tau_mu);

}

generated quantities {

real yhat_uncens[Nobs + Ncen];

real log_lik[Nobs + Ncen];

real lp[Nobs + Ncen];

for (i in 1:Nobs) {

lp[i] = mu + Xobs_bg[i,] * beta_bg;

yhat_uncens[i] = weibull_rng(alpha,

exp(-(mu + Xobs_bg[i,] * beta_bg)/alpha));

log_lik[i] = weibull_lpdf(yobs[i] | alpha,

exp(-(mu + Xobs_bg[i,] * beta_bg)/alpha));

}

for (i in 1:Ncen) {

lp[Nobs + i] = mu + Xcen_bg[i,] * beta_bg;

yhat_uncens[Nobs + i] = weibull_rng(alpha,

exp(-(mu + Xcen_bg[i,] * beta_bg)/alpha));

log_lik[Nobs + i] = weibull_lccdf(ycen[i] | alpha,

exp(-(mu + Xcen_bg[i,] * beta_bg)/alpha));

}

続いて、このStanコードをキックするためのRのソースコードです。元のデータが7043件と多すぎるのでランダムサンプリングしています。サンプリング数を8000、チェイン数を4にして実行します。（なお、可視化のソースコードもあるので結構長くなっていますので。頑張ってスクロールしてください。）

library(tidyverse)
library(rstan)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())


# Data Import -------------------------------------------------------------
dataset <- read_csv("dataset/WA_Fn-UseC_-Telco-Customer-Churn.csv")
dataset$Churn <- as.factor(dataset$Churn)
dataset$censored <- if_else(dataset$Churn == "Yes", 0, 1)

dataset_obs <- sample_n(dataset %>% filter(censored == 0), 300)
dataset_cens <- sample_n(dataset %>% filter(censored == 1), 75)

dataset <- rbind(dataset_cens, dataset_obs)

# Kick Stan Code ----------------------------------------------------------
stan_data <- list(
                ## 離脱のイベントが計測された顧客
                Nobs = sum(dataset$censored == 0),
                ## 途中で打ち切られた顧客
                Ncen = sum(dataset$censored == 1),
                ## 共変量の数
                M_bg = 1,
                ## 離脱イベントが計測された顧客の契約期間
                yobs = dataset$tenure[dataset$censored == 0],
                ## 途中で打ち切られた顧客の契約期間
                ycen = dataset$tenure[dataset$censored == 1],
                ## 離脱のイベントが計測された顧客の共変量
                Xobs_bg = matrix(as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 0]),
                ## 途中で打ち切られた顧客の共変量
                Xcen_bg = matrix(as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 1])
              )


fit <- rstan::stan(file = "model/weibull_fit.stan",
                   data = stan_data,
                   iter = 8000,
                   chains = 4,
                   seed = 1234,
                   control = list(max_treedepth = 15,adapt_delta=0.99)
                   )

# diagnose ----------------------------------------------------------------

fit

summary_table <- data.frame(summary(fit)$summary)
ggplot(data = data.frame(Rhat = summary_table$Rhat), aes(Rhat)) + geom_histogram()

rstan::traceplot(fit, par = c("alpha","mu","beta_bg"))

bayesplot::mcmc_acf(as.matrix(fit), pars = c("alpha","mu","beta_bg[1]"))

bayesplot::mcmc_areas(as.matrix(fit), pars = c("alpha","mu","beta_bg[1]"), prob = 0.95)


# visualization ----------------------------------------------------------------

draws <- tidybayes::tidy_draws(fit)
draws


treatment_assignment <- c(as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 0],
                          as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 1])
treatment_assignment_df <- data_frame(obs = 1:nrow(dataset),treatment = treatment_assignment)
treatment_assignment_df


draws_yhat_uncens <- draws %>%
                      select(.chain, .iteration, .draw, starts_with("yhat_uncens")) %>%
                      gather(key = key, value = yhat_uncens, starts_with("yhat_uncens")) %>%
                      separate(col = key, sep = "uncens", into = c("key","obs")) %>%
                      select(-key) %>%
                      ## Avoid using regular expressions with square brackets (syntax highlighter broke).
                      ## https://stringr.tidyverse.org/articles/stringr.html
                      mutate(obs = as.integer(str_sub(obs, 2, -2))) %>%
                      left_join(y = treatment_assignment_df)
draws_yhat_uncens


ggplot(data = draws_yhat_uncens,
       mapping = aes(x = yhat_uncens, color = factor(treatment))) +
  geom_density(n = 512*10) +
  coord_cartesian(xlim = c(0,160)) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.key = element_blank(),
        plot.title = element_text(hjust = 0.5),
        strip.background = element_blank())



## Constructor for treatment-specific survival function
construct_survival_function <- function(alpha, mu, beta, x) {
  function(t) {
    sigma_i <- exp(-1 * (mu + beta * x) / alpha)
    exp(- (t / sigma_i)^alpha)
  }
}

## Random functions
survival_functins <- draws %>%
                        select(.chain, .iteration, .draw, alpha, mu, `beta_bg[1]`) %>%
                        ## Simplify name
                        rename(beta = `beta_bg[1]`) %>%
                        ## Construct realization of random functions
                        mutate(`S(t|1)` = pmap(list(alpha, mu, beta), function(a,m,b) {construct_survival_function(a,m,b,1)}),
                               `S(t|0)` = pmap(list(alpha, mu, beta), function(a,m,b) {construct_survival_function(a,m,b,0)}))
survival_functins



times <- seq(from = 0, to = 160, by = 0.1)
times_df <- data_frame(t = times)

## Try first realizations
survival_functins$`S(t|1)`[[1]](times[1:10])


survival_functins$`S(t|0)`[[1]](times[1:10])


## Apply all realizations
survival <- survival_functins %>%
              mutate(times_df = list(times_df)) %>%
              mutate(times_df = pmap(list(times_df, `S(t|1)`, `S(t|0)`),
                                     function(df, s1, s0) {df %>% mutate(s1 = s1(t),
                                                                         s0 = s0(t))})) %>%
              select(-`S(t|1)`, -`S(t|0)`) %>%
              unnest() %>%
              gather(key = treatment, value = survival, s1, s0) %>%
              mutate(treatment = factor(treatment,
                                        levels = c("s1","s0"),
                                        labels = c("Yes","No")))

## Average on survival scale
survival_mean <- survival %>%
                    group_by(treatment, t) %>%
                    summarize(survival_mean = mean(survival),
                              survival_95upper = quantile(survival, probs = 0.975),
                              survival_95lower = quantile(survival, probs = 0.025))

ggplot(data = survival,
       mapping = aes(x = t, y = survival, color = treatment, group = interaction(.chain,.draw,treatment))) +
  geom_line(size = 0.1, alpha = 0.02) +
  geom_line(data = survival_mean,
            mapping = aes(y = survival_mean, group = treatment)) +
  geom_line(data = survival_mean,
            mapping = aes(y = survival_95upper, group = treatment),
            linetype = "dotted") +
  geom_line(data = survival_mean,
            mapping = aes(y = survival_95lower, group = treatment),
            linetype = "dotted") +
  facet_grid(. ~ treatment) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.key = element_blank(),
        plot.title = element_text(hjust = 0.5),
        strip.background = element_blank())



## Average on parameter space
average_parameters <- draws %>%
                        summarize(alpha = mean(alpha),
                                  mu = mean(mu),
                                  beta = mean(`beta_bg[1]`))
average_parameters


average_params_survival1 <- with(average_parameters,
                                              construct_survival_function(alpha, mu, beta, 1))
average_params_survival0 <- with(average_parameters,
                                              construct_survival_function(alpha, mu, beta, 0))
average_params_survival <-
                            data_frame(t = seq(from = 0, to = 160, by = 0.1),
                                       s1 = average_params_survival1(t),
                                       s0 = average_params_survival0(t)) %>%
                            gather(key = treatment, value = survival, -t) %>%
                            mutate(treatment = factor(treatment,
                                                      levels = c("s1","s0"),
                                                      labels = c("Yes","No")))

average_params_survival %>%
  ggplot(mapping = aes(x = t, y = survival, color = treatment, group = treatment)) +
  geom_line() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.key = element_blank(),
        plot.title = element_text(hjust = 0.5),
        strip.background = element_blank())



ggplot(data = survival,
       mapping = aes(x = t, y = survival, color = treatment, group = interaction(.chain,.draw,treatment))) +
  geom_line(size = 0.1, alpha = 0.02) +
  geom_line(data = survival_mean,
            mapping = aes(y = survival_mean, group = treatment)) +
  geom_line(data = average_params_survival,
            mapping = aes(group = treatment),
            linetype = "dotted") +
  facet_grid(. ~ treatment) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.key = element_blank(),
        plot.title = element_text(hjust = 0.5),
        strip.background = element_blank())

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

library(tidyverse)

library(rstan)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

# Data Import -------------------------------------------------------------

dataset <- read_csv("dataset/WA_Fn-UseC_-Telco-Customer-Churn.csv")

dataset$Churn <- as.factor(dataset$Churn)

dataset$censored <- if_else(dataset$Churn == "Yes", 0, 1)

dataset_obs <- sample_n(dataset %>% filter(censored == 0), 300)

dataset_cens <- sample_n(dataset %>% filter(censored == 1), 75)

dataset <- rbind(dataset_cens, dataset_obs)

# Kick Stan Code ----------------------------------------------------------

stan_data <- list(

## 離脱のイベントが計測された顧客

Nobs = sum(dataset$censored == 0),

## 途中で打ち切られた顧客

Ncen = sum(dataset$censored == 1),

## 共変量の数

M_bg = 1,

## 離脱イベントが計測された顧客の契約期間

yobs = dataset$tenure[dataset$censored == 0],

## 途中で打ち切られた顧客の契約期間

ycen = dataset$tenure[dataset$censored == 1],

## 離脱のイベントが計測された顧客の共変量

Xobs_bg = matrix(as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 0]),

## 途中で打ち切られた顧客の共変量

Xcen_bg = matrix(as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 1])

)

fit <- rstan::stan(file = "model/weibull_fit.stan",

data = stan_data,

iter = 8000,

chains = 4,

seed = 1234,

control = list(max_treedepth = 15,adapt_delta=0.99)

)

# diagnose ----------------------------------------------------------------

fit

summary_table <- data.frame(summary(fit)$summary)

ggplot(data = data.frame(Rhat = summary_table$Rhat), aes(Rhat)) + geom_histogram()

rstan::traceplot(fit, par = c("alpha","mu","beta_bg"))

bayesplot::mcmc_acf(as.matrix(fit), pars = c("alpha","mu","beta_bg[1]"))

bayesplot::mcmc_areas(as.matrix(fit), pars = c("alpha","mu","beta_bg[1]"), prob = 0.95)

# visualization ----------------------------------------------------------------

draws <- tidybayes::tidy_draws(fit)

draws

treatment_assignment <- c(as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 0],

as.numeric(dataset$PaperlessBilling == "Yes")[dataset$censored == 1])

treatment_assignment_df <- data_frame(obs = 1:nrow(dataset),treatment = treatment_assignment)

treatment_assignment_df

draws_yhat_uncens <- draws %>%

select(.chain, .iteration, .draw, starts_with("yhat_uncens")) %>%

gather(key = key, value = yhat_uncens, starts_with("yhat_uncens")) %>%

separate(col = key, sep = "uncens", into = c("key","obs")) %>%

select(-key) %>%

## Avoid using regular expressions with square brackets (syntax highlighter broke).

## https://stringr.tidyverse.org/articles/stringr.html

mutate(obs = as.integer(str_sub(obs, 2, -2))) %>%

left_join(y = treatment_assignment_df)

draws_yhat_uncens

ggplot(data = draws_yhat_uncens,

mapping = aes(x = yhat_uncens, color = factor(treatment))) +

geom_density(n = 512*10) +

coord_cartesian(xlim = c(0,160)) +

theme_bw() +

theme(axis.text.x = element_text(angle = 90, vjust = 0.5),

legend.key = element_blank(),

plot.title = element_text(hjust = 0.5),

strip.background = element_blank())

## Constructor for treatment-specific survival function

construct_survival_function <- function(alpha, mu, beta, x) {

function(t) {

sigma_i <- exp(-1 * (mu + beta * x) / alpha)

exp(- (t / sigma_i)^alpha)

}

## Random functions

survival_functins <- draws %>%

select(.chain, .iteration, .draw, alpha, mu, `beta_bg[1]`) %>%

## Simplify name

rename(beta = `beta_bg[1]`) %>%

## Construct realization of random functions

mutate(`S(t|1)` = pmap(list(alpha, mu, beta), function(a,m,b) {construct_survival_function(a,m,b,1)}),

`S(t|0)` = pmap(list(alpha, mu, beta), function(a,m,b) {construct_survival_function(a,m,b,0)}))

survival_functins

times <- seq(from = 0, to = 160, by = 0.1)

times_df <- data_frame(t = times)

## Try first realizations

survival_functins$`S(t|1)`[[1]](times[1:10])

survival_functins$`S(t|0)`[[1]](times[1:10])

## Apply all realizations

survival <- survival_functins %>%

mutate(times_df = list(times_df)) %>%

mutate(times_df = pmap(list(times_df, `S(t|1)`, `S(t|0)`),

function(df, s1, s0) {df %>% mutate(s1 = s1(t),

s0 = s0(t))})) %>%

select(-`S(t|1)`, -`S(t|0)`) %>%

unnest() %>%

gather(key = treatment, value = survival, s1, s0) %>%

mutate(treatment = factor(treatment,

levels = c("s1","s0"),

labels = c("Yes","No")))

## Average on survival scale

survival_mean <- survival %>%

group_by(treatment, t) %>%

summarize(survival_mean = mean(survival),

survival_95upper = quantile(survival, probs = 0.975),

survival_95lower = quantile(survival, probs = 0.025))

ggplot(data = survival,

mapping = aes(x = t, y = survival, color = treatment, group = interaction(.chain,.draw,treatment))) +

geom_line(size = 0.1, alpha = 0.02) +

geom_line(data = survival_mean,

mapping = aes(y = survival_mean, group = treatment)) +

geom_line(data = survival_mean,

mapping = aes(y = survival_95upper, group = treatment),

linetype = "dotted") +

geom_line(data = survival_mean,

mapping = aes(y = survival_95lower, group = treatment),

linetype = "dotted") +

facet_grid(. ~ treatment) +

theme_bw() +

theme(axis.text.x = element_text(angle = 90, vjust = 0.5),

legend.key = element_blank(),

plot.title = element_text(hjust = 0.5),

strip.background = element_blank())

## Average on parameter space

average_parameters <- draws %>%

summarize(alpha = mean(alpha),

mu = mean(mu),

beta = mean(`beta_bg[1]`))

average_parameters

average_params_survival1 <- with(average_parameters,

construct_survival_function(alpha, mu, beta, 1))

average_params_survival0 <- with(average_parameters,

construct_survival_function(alpha, mu, beta, 0))

average_params_survival <-

data_frame(t = seq(from = 0, to = 160, by = 0.1),

s1 = average_params_survival1(t),

s0 = average_params_survival0(t)) %>%

gather(key = treatment, value = survival, -t) %>%

mutate(treatment = factor(treatment,

levels = c("s1","s0"),

labels = c("Yes","No")))

average_params_survival %>%

ggplot(mapping = aes(x = t, y = survival, color = treatment, group = treatment)) +

geom_line() +

theme_bw() +

theme(axis.text.x = element_text(angle = 90, vjust = 0.5),

legend.key = element_blank(),

plot.title = element_text(hjust = 0.5),

strip.background = element_blank())

ggplot(data = survival,

mapping = aes(x = t, y = survival, color = treatment, group = interaction(.chain,.draw,treatment))) +

geom_line(size = 0.1, alpha = 0.02) +

geom_line(data = survival_mean,

mapping = aes(y = survival_mean, group = treatment)) +

geom_line(data = average_params_survival,

mapping = aes(group = treatment),

linetype = "dotted") +

facet_grid(. ~ treatment) +

theme_bw() +

theme(axis.text.x = element_text(angle = 90, vjust = 0.5),

legend.key = element_blank(),

plot.title = element_text(hjust = 0.5),

strip.background = element_blank())

Rhatは全て1.05以下になっています。

traceplotを見る限り、重なり合っているので問題なさそうです。

各パラメータごとの自己相関係数です。こちらも問題なさそうです。

推定したパラメータの分布です。

横軸は推定した継続期間です。決済の電子化をしていない消費者は、契約期間の短い際の確率密度が低い傾向があるようです。

どうやら離脱率に関して決済の電子化をしていない消費者は、そうでない消費者よりも低いようです。

こちらは95%で取りうる範囲をそれぞれプロットしたものです。

おわりに

Stanで生存時間分析を行うという事例はそんなに多くはないものの、業界の長たちが良いコードを作成してくれていました。また、面白そうなデータセットも見つけることができました。このようなデータがもっと広まっていけば、マーケティングにおける生存時間分析がより活発に行われるのかもしれません。

参考文献

[1] 豊田秀樹 (2017) 『実践ベイズモデリング -解析技法と認知モデル-』朝倉書店
[2]生存時間解析入門
[3]比例ハザードモデルはとってもtricky！
[4]Stanで生存時間解析(Weibull 回帰)
[5]生存時間分析をStanで実行してみた
[6]階層ベイズ生存解析を用いたwebサイトの訪問者分析に関するStanでの実装
[7]生存時間分析 – ハザード関数に時間相関の制約を入れる
[8]Bayesian Survival Analysis 1: Weibull Model with Stan
[9]Bayesian Inference With Stan ~062~
[10]生存時間解析について – 概要編
[11]Survival Analysis for Employee Attrition ※kaggleで提供されているHR系のデータをサバイバル分析に用いている。
[12]Survival Analysis with R※Random Forests Modelによる生存時間の推定がなされている。
[13]Survival Analysis with R and Aster ※服役後の犯罪に関する分析や、離婚の分析などをしている。
[14]Survival Analysis of Mobile Prepaid Customers Using the Weibull Distribution（ダウンロード注意）

[Stan]項目反応理論(IRT)の段階反応モデルでbaysemのアンケートデータの分析をしてみる

はじめに

stanのユーザーガイドを見ていて、項目反応理論（IRT）についての章があり気になりました。勉強会のLTなどで手法の名前をちらっと聞いたことはあったのですが、使い道について調べていませんでした。ビジネスにおける実活用もしやすそうだと思ったので、カジュアルに分析して備忘録として残したいと思います。

目次
・項目反応理論（Item Response Theory：IRT）とは
・ビジネスでの適用可能性について
・データ
・モデルの推定
・結果の解釈
・おわりに

項目反応理論（Item Response Theory：IRT）とは

関西学院大学の教授のブログによると、

項目反応理論とは、テストについての計量モデルで、問題に対する正解・不正解のデータから、問題の特性や、回答者の学力を推定するためのモデルです。

とあります。また、Wikipediaによると、TOEFLの問題の評価のために使われているそうです。

主に、バイナリーと順序変数のモデルがあるようで、以下の母数がモデルに想定されています。どちらもほぼ同じです。

回答が2値変数のモデル

2母数のロジスティックモデル
- 特性値（例えば、広告配信の満足度とか）
- 識別度母数（項目特性曲線の傾き）
- 困難度母数（項目特性曲線の切片）
- 定数

回答が順序変数のモデル（まずい < まぁまぁ < おいしい）

段階反応モデル
- 特性値（例えば、広告配信の満足度とか）
- 識別度母数（項目特性曲線の傾き）
- 困難度母数（項目特性曲線の切片）

※項目特性曲線は横軸に特性値、縦軸に質問の正答率を取ったものです。

ビジネスでの適用可能性について

顧客のアンケート結果の解釈
- 異質な集団間の得点を比較可能
- 異なる尺度間の得点を比較可能（昔のアンケートだと5段階、今のアンケートは7段階などの状況はビジネスデータでありうる。）
人事評価のバイアスの統制
- 採用面接時の個人特性の正当な評価
アンケート項目の項目削減によるアンケートコストの低減
- 各アンケート項目が理解されたかどうかを分析し、一つ一つのアンケート項目の精度を高める

データ

今回扱うデータはbaysemパッケージに入っているデータセットです。Yellow Pagesの広告プロダクトにおける満足度サーベイの回答データで、全ての回答は1から10のスケールで点数が付けられています（1がPoorで10がExcellent）。質問数は10個で、回答数は1811件です。

各質問の内容（baysemパッケージのドキュメントに載っていました。）
q1：全体の満足度

価格について
q2：競争的な価格設定
q3：昨年と同じ広告の最小値に対しての価格の引き上げ
q4：消費者の数に対しての適切な価格設定

効果について
q5：広告の購入の潜在的な影響
q6：広告を通じて自身のビジネスへの集客ができたか
q7：多くの消費者にリーチしたかどうか
q8：年間を通じて消費者に対する長期での露出があったか
q9：多くの家計やビジネスを必要としている人に届いたかどうか
q10：ビジネスを必要としている地理上のエリアに届いたかどうか

今回のIRT適用における特性値は、「広告プロダクトに関する満足度の傾向」としてみたいと思います。

モデルの推定

今回は教科書にならって以下の段階反応モデルを用います。

ここでaは識別力（広告の満足度が高まりやすいかどうか）、bは境界パラメータ（回答カテゴリ間の境界値）、θは特性（回答者がどれだけ広告に満足しているか）を表しています。Dは定数項で、以下では1とおいています。cはアンケートの回答のカテゴリ番号です。今回の例では10段階の評価が入ることになります。最後に、uは反応を、jは質問の番号を表しています。

実践ベイズモデリング -解析技法と認知モデル-

こちらの本のサポートサイトからダウンロードできるzipファイルにstanのコードやRコードがありますので、そちらを利用しています。

モデルですが、以下のような設定となっています。

data{
    int ni; // the number of record
    int nj; // the number of item
    int nc; // the number of grade
    real D; // Constant
    int<lower=1,upper=10> y[ni,nj]; // data
}

parameters{
  vector<lower=0,upper=5>[nj] a;
    ordered[nc-1] ba[nj];
    vector<lower=-4,upper=4>[ni] theta;
}

transformed parameters{
    real b[nj,nc];
    vector<lower=0,upper=1>[nc-1] pa[ni,nj];
    simplex[nc] p[ni,nj];
    for (j in 1:nj){
        for (c in 1:nc){
            if (c ==1){
                b[j,c] = ba[j,c];
            }else if (c ==nc){
                b[j,c] = ba[j,c-1];
            }else{
                b[j,c] = (ba[j,c-1]+ba[j,c])/2;
            }
        }
    }
    for (i in 1:ni){
        for (j in 1:nj){
            for (c in 1:nc-1){
                pa[i,j,c] = 1/(1+exp(-D*a[j]*(theta[i] - ba[j,c])));
            }       
        }
    }
    for (i in 1:ni){
        for (j in 1:nj){
            for(c in 1:nc){
                if (c==1){
                    p[i,j,c] = 1-pa[i,j,c];
                }else if(c==nc){
                    p[i,j,c] = pa[i,j,c-1];
                }else{
                    p[i,j,c] = pa[i,j,c-1] - pa[i,j,c];
                }
            }
        }
    }
}

model{
    for (i in 1:ni){
        theta[i] ~ normal(0,1);
        for (j in 1:nj){
            y[i,j] ~ categorical(p[i,j]);
        }
    }
    for (j in 1:nj){
        a[j] ~ lognormal(0,sqrt(0.5));
        for (c in 1:nc-1){
            ba[j,c] ~ normal(0,2);
        }
    }
}

generated quantities{
    real bg[nj,nc];
    bg = b;
}

data{

int ni; // the number of record

int nj; // the number of item

int nc; // the number of grade

real D; // Constant

int<lower=1,upper=10> y[ni,nj]; // data

}

parameters{

vector<lower=0,upper=5>[nj] a;

ordered[nc-1] ba[nj];

vector<lower=-4,upper=4>[ni] theta;

}

transformed parameters{

real b[nj,nc];

vector<lower=0,upper=1>[nc-1] pa[ni,nj];

simplex[nc] p[ni,nj];

for (j in 1:nj){

for (c in 1:nc){

if (c ==1){

b[j,c] = ba[j,c];

}else if (c ==nc){

b[j,c] = ba[j,c-1];

}else{

b[j,c] = (ba[j,c-1]+ba[j,c])/2;

}

for (i in 1:ni){

for (j in 1:nj){

for (c in 1:nc-1){

pa[i,j,c] = 1/(1+exp(-D*a[j]*(theta[i] - ba[j,c])));

}

for (i in 1:ni){

for (j in 1:nj){

for(c in 1:nc){

if (c==1){

p[i,j,c] = 1-pa[i,j,c];

}else if(c==nc){

p[i,j,c] = pa[i,j,c-1];

}else{

p[i,j,c] = pa[i,j,c-1] - pa[i,j,c];

}

model{

for (i in 1:ni){

theta[i] ~ normal(0,1);

for (j in 1:nj){

y[i,j] ~ categorical(p[i,j]);

}

for (j in 1:nj){

a[j] ~ lognormal(0,sqrt(0.5));

for (c in 1:nc-1){

ba[j,c] ~ normal(0,2);

}

generated quantities{

real bg[nj,nc];

bg = b;

}

こちらをキックするためのRコードです。

library(bayesm)
library(rstan)
library(shinystan)
library(bayesplot)
library(tidyverse)
library(gridExtra)

# Data Import -------------------------------------------------------------
data("customerSat")
dataset <- customerSat

# Kick Stan model ---------------------------------------------------------
ni <- nrow(dataset) # 分析対象者の数
nj <- ncol(dataset) # 項目数
nc <- length(table(as.factor(dataset$q1))) # 回答の種類
D <- 1 # 定数項

stan_data <- list(y = dataset,
                  nj = nj,
                  ni = ni,
                  nc= nc,
                  D = D)

par <- c("theta","ba","a","b")
war <- 2500
ite <- 5000
see <- 1234
dig <- 2
cha <- 4

fit <- stan(file = "model/graded_response_model.stan",
                data = stan_data,
                pars = par,
                verbose = F,
                seed = see,
                chains = cha,
                warmup = war,
                iter = ite)

# Diagnose ----------------------------------------------------------------

traceplot(fit)

print(fit, pars = par, digits_summary = dig)

summary_table <- data.frame(summary(fit)$summary)
ggplot(data = data.frame(Rhat = summary_table$Rhat), aes(Rhat)) + geom_histogram()

library(bayesm)

library(rstan)

library(shinystan)

library(bayesplot)

library(tidyverse)

library(gridExtra)

# Data Import -------------------------------------------------------------

data("customerSat")

dataset <- customerSat

# Kick Stan model ---------------------------------------------------------

ni <- nrow(dataset) # 分析対象者の数

nj <- ncol(dataset) # 項目数

nc <- length(table(as.factor(dataset$q1))) # 回答の種類

D <- 1 # 定数項

stan_data <- list(y = dataset,

nj = nj,

ni = ni,

nc= nc,

D = D)

par <- c("theta","ba","a","b")

war <- 2500

ite <- 5000

see <- 1234

dig <- 2

cha <- 4

fit <- stan(file = "model/graded_response_model.stan",

data = stan_data,

pars = par,

verbose = F,

seed = see,

chains = cha,

warmup = war,

iter = ite)

# Diagnose ----------------------------------------------------------------

traceplot(fit)

print(fit, pars = par, digits_summary = dig)

summary_table <- data.frame(summary(fit)$summary)

ggplot(data = data.frame(Rhat = summary_table$Rhat), aes(Rhat)) + geom_histogram()

処理時間としては、2014年末モデルのMacbook Proのcorei5、メモリ8GBで数時間程度でした。（正確な時間はわかりませんが、寝て起きたら計算が終わっていました。）

どうやら収束してそうです。

Rhatも1.1未満におさまっています。

結果の解釈

# histogram of theta
theta <- rstan::extract(fit)$theta %>% apply(2,mean)
ggplot(data = data.frame(theta_mean = theta), aes(theta_mean)) + geom_histogram()

# histogram of theta

theta <- rstan::extract(fit)$theta %>% apply(2,mean)

ggplot(data = data.frame(theta_mean = theta), aes(theta_mean)) + geom_histogram()

まず、推定した特性値の値のユーザーごとの平均値を求めて、ヒストグラムを描いてみます。どうやら、上限周辺にやたらと高い評価をしてそうなユーザーがいるようです。

最後に、項目特性曲線を質問ごとに、そして回答ごとに描いてみようと思います。

# Visualization -----------------------------------------------------------
alpha <- rstan::extract(fit)$a %>% apply(2,mean)

for (i in 1:nj){
  eval(parse(text=paste0("beta_",i," <- rstan::extract(fit)$b[1:10000,",i,",1:nc] %>% apply(2,mean)")))
}

for (i in 1:nj){
  eval(parse(text=paste0("ggdf_",i," <- data.frame(matrix(ncol = nc,nrow = length(theta))) ; colnames(ggdf_",i,") <- 1:nc")))
}

## probability
for(i in 1:nj){
  for (j in 1:nc){
    eval(parse(text=paste0("ggdf_",j,"[,",i,"] <- 1/(1+exp(-alpha[",i,"]*(theta-beta_",j,"[",i,"])))")))
    eval(parse(text=paste0("ggdf_",j,"$theta <- theta")))
  }
}

## gather
for (i in 1:nj){
  eval(parse(text=paste0("ggdf_gt_",i," <- ggdf_",i," %>% tidyr::gather(key=var,value,-theta,factor_key=TRUE)")))
}

## ggplot
for (i in 1:nj){
  eval(parse(text=paste0("p",i," <- ggplot(data = ggdf_gt_",i,", aes(x = theta, y = value, colour = var)) + geom_line() + ggtitle(\"Q",i,"\")")))
}

# 2×3でグラフを描画
grid.arrange(p1, p2, p3, p4, p5, nrow = 3)
grid.arrange(p6, p7, p8, p9, p10,  nrow = 3)

# Visualization -----------------------------------------------------------

alpha <- rstan::extract(fit)$a %>% apply(2,mean)

for (i in 1:nj){

eval(parse(text=paste0("beta_",i," <- rstan::extract(fit)$b[1:10000,",i,",1:nc] %>% apply(2,mean)")))

}

for (i in 1:nj){

eval(parse(text=paste0("ggdf_",i," <- data.frame(matrix(ncol = nc,nrow = length(theta))) ; colnames(ggdf_",i,") <- 1:nc")))

}

## probability

for(i in 1:nj){

for (j in 1:nc){

eval(parse(text=paste0("ggdf_",j,"[,",i,"] <- 1/(1+exp(-alpha[",i,"]*(theta-beta_",j,"[",i,"])))")))

eval(parse(text=paste0("ggdf_",j,"$theta <- theta")))

}

## gather

for (i in 1:nj){

eval(parse(text=paste0("ggdf_gt_",i," <- ggdf_",i," %>% tidyr::gather(key=var,value,-theta,factor_key=TRUE)")))

}

## ggplot

for (i in 1:nj){

eval(parse(text=paste0("p",i," <- ggplot(data = ggdf_gt_",i,", aes(x = theta, y = value, colour = var)) + geom_line() + ggtitle(\"Q",i,"\")")))

}

# 2×3でグラフを描画

grid.arrange(p1, p2, p3, p4, p5, nrow = 3)

grid.arrange(p6, p7, p8, p9, p10, nrow = 3)

質問1~10に関して、10段階の回答ごとの項目反応曲線を以下に描いています。上まで戻るのが面倒なので、質問内容を再掲します。

q1：全体の満足度
q2：競争的な価格設定
q3：昨年と同じ広告の最小値に対しての価格の引き上げ
q4：消費者の数に対しての適切な価格設定
q5：広告の購入の潜在的な影響
q6：広告を通じて自身のビジネスへの集客ができたか
q7：多くの消費者にリーチしたかどうか
q8：年間を通じて消費者に対する長期での露出があったか
q9：多くの家計やビジネスを必要としている人に届いたかどうか
q10：ビジネスを必要としている地理上のエリアに届いたかどうか

これらの傾向から、9〜10点を獲得するにはある程度は特性値が高まる必要がある質問としては、q1〜q6のように見えます。価格や購買など自身のビジネスに直結しそうな質問が多い印象です。逆にふわっとした質問であるq7~q10は特性値が低くても9〜10点を取れる可能性が高い傾向があります。

おわりに

Stanのユーザーガイドを読むことで、普段自分が業務で扱っているアプローチなどが如何に限定的であることが実感できました。今回はIRTのアンケートデータへの適用事例を知れ、そこから様々な文献や便利なコードに至ることができました。社内のアンケートデータへの適用は面白そうだと思いますので業務で使ってみようと思います。

参考情報

[1] 豊田秀樹 (2017) 『実践ベイズモデリング -解析技法と認知モデル-』朝倉書店
[2] Yoshitake Takebayashi (2015) 「項目反応理論による尺度運用」 SlideShare
[3] 持主弓子・今城志保 (2011) 「IRTの組織サーベイへの応用」
[4] 清水裕士 (2017) 「項目反応理論をStanで実行するときのあれこれ」 Sunny side up!
[5] 清水裕士 (2016) 「Stanで多次元項目反応理論」 Sunny side up!
[6] 小杉考司 (2013) 「項目反応理論について」
[7] Daniel C. Furr et al. (2016) “Two-Parameter Logistic Item Response Model”
[8] daiki hojo (2018) “Bayesian Sushistical Modeling” Tokyo.R#70
[9] abrahamcow (2017) 「［RStan］項目反応理論の応用でフリースタイルダンジョン登場ラッパーの強さをランキングしてみた」

R Advent Calendar 2018 一発屋芸人の検索トレンドの分析

はじめに

昨年のR Advent Calendarはポケモンのデータをrvestでスクレイピングして、レアポケモンがどのような種族値における特徴があるのかを探ったり、経験値が必要な割に種族値が低い「コスパの悪いポケモン」などを見つけました。
今年のR Advent Calendarでは、年末年始といえば一発屋芸人のテレビなどでの露出が多くなることから、一発屋芸人の検索トレンドのデータを手に入れて分析してみたいと思います。

分析工程

・データの収集
・データの整形
・可視化
・分析

データの収集

こちらのサイト（流行した一発屋芸人一覧／年代流行）に一発屋の芸人さんが列挙されていました。私は普段テレビを見ないので大体の芸人さんがわからないです。

Googleトレンドから、芸人名に関するGoogle検索の時系列データを収集します。

非常に残酷なデータだなと思いました。

ただ、一つ弁護すると、Googleトレンドはレベルではなくピークを1として標準化した数値をデータとして提供してくれていますので、
ピークが著しく高ければ、今の水準が低くてもそこそこ検索されている可能性はあるとだけ言っておきます。

本当の検索回数が必要な場合は、Google Adwords（検索連動型広告）のアカウントの開設とともに検索ボリューム取得APIなどの申請が必要なので、正確なデータが必要な場合は会社として取り組んだほうが良いと思います。個人では厳しいです。

データの整形

各芸人さん（総勢21名）の検索トレンドデータのピークの6ヶ月前までのデータと6ヶ月後のデータまでの合計1年間の検索トレンドを各々抽出してみようと思います。
GoogleトレンドのデータはCSVでダウンロードできますので、そのCSVを読み込み、トレンドのデータを文字列から数値にし、ピークの前後12ヶ月ずつのデータを抽出します。
そうすることで、一発屋芸人のピークの前後に関するデータを作ります。（ただし、今朝、Google Trendのデータを取得できるgtrendsRというパッケージがR bloggerで紹介されていました。APIないはずなんですが、URLの工夫か裏でSelenium動かしていたりするんですかね。）

library(tidyverse)
library(directlabels)
library(TSclust)
library(gghighlight)

# データの準備 ------------------------------------------------------------------
trend_dataset <- data.frame()

for (i in 1:21) {
  #データの読込
  trend_data <- read_csv(file = paste0("multiTimeline_",i,".csv"),skip = 2)
  trend_data <- trend_data %>% mutate( gsub(": (日本)","",colnames(trend_data)[2]) )
  colnames(trend_data) <- c("month", "trend", "keyword")
  #1未満のデータをゼロにする
  trend_data <- trend_data %>% mutate(trend = as.numeric(replace(trend, trend=="1 未満", 0)))
  #ピークの月の前後12ヶ月を抽出
  trend_data <- trend_data[(which.max(trend_data$trend)-12):(which.max(trend_data$trend)+12),]
  trend_dataset <- trend_dataset %>% rbind(trend_data)
}

#キーワードごとにインデックスをふる

trend_dataset <- trend_dataset %>%
                        group_by(keyword) %>%
                        mutate(period = 1:n())

library(tidyverse)

library(directlabels)

library(TSclust)

library(gghighlight)

# データの準備 ------------------------------------------------------------------

trend_dataset <- data.frame()

for (i in 1:21) {

#データの読込

trend_data <- read_csv(file = paste0("multiTimeline_",i,".csv"),skip = 2)

trend_data <- trend_data %>% mutate( gsub(": (日本)","",colnames(trend_data)[2]) )

colnames(trend_data) <- c("month", "trend", "keyword")

#1未満のデータをゼロにする

trend_data <- trend_data %>% mutate(trend = as.numeric(replace(trend, trend=="1 未満", 0)))

#ピークの月の前後12ヶ月を抽出

trend_data <- trend_data[(which.max(trend_data$trend)-12):(which.max(trend_data$trend)+12),]

trend_dataset <- trend_dataset %>% rbind(trend_data)

}

#キーワードごとにインデックスをふる

trend_dataset <- trend_dataset %>%

group_by(keyword) %>%

mutate(period = 1:n())

可視化

作成したデータを実際にプロットしてみます。

old = theme_set(theme_gray(base_family="HiraKakuProN-W3"))

ggplot(data = trend_dataset,
            aes(x = period, y = trend, color=keyword))  + geom_line()

old = theme_set(theme_gray(base_family="HiraKakuProN-W3"))

ggplot(data = trend_dataset,

aes(x = period, y = trend, color=keyword)) + geom_line()

一発屋にも盛り上がり方に違いがあるようですね。

時系列クラスタリングの適用

多様な盛り上がり方があることから、TSclustというライブラリを使って時系列クラスタリングを行い、トレンドに関しての分類的なものを得たいと思います。
今回初めて使うのですが、参考文献によると様々な類似性指標を指定して、時系列ごとの類似性を計算するようです。ピアソン相関係数のようなシンプルなものもあれば、ユークリッド距離のものやFrechet距離とかいう聞いたことないものまで幅広く用意されています。今回はシンプルにピアソン相関係数にしてみます。そして、類似性指標を出してから、そのまま階層クラスタリングを行います。

trend_dataset_spread <- trend_dataset %>% select(-month) %>% tidyr::spread(key = keyword, value = trend)

# COR距離で距離行列を作成
d <- diss(trend_dataset_spread %>% select(-period), "COR")

#デフォルトの設定で階層クラスタリング
h <- hclust(d)

#階層クラスタリングの結果の可視化
par(cex=0.6)
par(family = "HiraKakuProN-W3")
plot(h, hang = -1)

trend_dataset_spread <- trend_dataset %>% select(-month) %>% tidyr::spread(key = keyword, value = trend)

# COR距離で距離行列を作成

d <- diss(trend_dataset_spread %>% select(-period), "COR")

#デフォルトの設定で階層クラスタリング

h <- hclust(d)

#階層クラスタリングの結果の可視化

par(cex=0.6)

par(family = "HiraKakuProN-W3")

plot(h, hang = -1)

こちらが、TSclustのdiss関数を用いて計算した時系列データごとの距離を、階層クラスタリングにより描いたデンドログラムです。

この分類だけ見ても、芸人さんを知らない私からすると何も共感がありませんので、先程のクラスタリング結果をもとに可視化をしてみます。
そこで、Tokyo.Rで知らない人はいないであろう、yutaniさんの作られたgghighlightを使ってみようと思います。

ただ、日本語のラベルの表示がうまくいかなかったので、芸人さんの名前をGoogleSpreadSheetのGoogle翻訳関数（GOOGLETRANSLATE）で英訳しておきます。

（Anyway bright YasumuraやThick slice Jasonは結構キャッチーなのでは？）

# クラスタ数は3とする
data.frame(cutree(h, 3))
clusters <- data.frame(cluster_number=cutree(h, 3))
clusters$keyword <- rownames(clusters)
rownames(clusters ) <- NULL
trend_dataset_withcluster <- trend_dataset %>% left_join(clusters, by = "keyword")

#英訳したデータの読み込みと結合
rename_keywordlist <- read_csv("rename_keywordlist.csv")
trend_dataset_withcluster <- trend_dataset_withcluster %>% left_join(rename_keywordlist, by = "keyword")


gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),
                 predicate = max(cluster_number) == 1)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),
                 predicate = max(cluster_number) == 2)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),
                 predicate = max(cluster_number) == 3)

# クラスタ数は3とする

data.frame(cutree(h, 3))

clusters <- data.frame(cluster_number=cutree(h, 3))

clusters$keyword <- rownames(clusters)

rownames(clusters ) <- NULL

trend_dataset_withcluster <- trend_dataset %>% left_join(clusters, by = "keyword")

#英訳したデータの読み込みと結合

rename_keywordlist <- read_csv("rename_keywordlist.csv")

trend_dataset_withcluster <- trend_dataset_withcluster %>% left_join(rename_keywordlist, by = "keyword")

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),

predicate = max(cluster_number) == 1)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),

predicate = max(cluster_number) == 2)

gghighlight_line(trend_dataset_withcluster, aes(period, trend, colour = keyword_en),

predicate = max(cluster_number) == 3)

まずはクラスター1

比較的短期でピークに達し、すぐに検索されなくなる、一発屋の名に相違ない傾向を持ったクラスターのように思われます。「日本エレキテル連合」とか「楽しんご」とか「8.6秒バズーカ」とかです。

続いてクラスター2

急激にピークに達するものの、ややしぶとく残り続けるような一発屋のクラスターなのかなと思います。「レイザーラモンHG」とか「厚切りジェイソン」とか「ピコ太郎」とか「世界のナベアツ」です。

そしてクラスター3

3人の芸人さんしか属していないですね。クラスターの数は2個でもよかったかもしれない。段階的にピークに達し、一気に落とされるという一発屋のクラスターのようです。「とにかく明るい安村」とか「藤崎マーケット」とか「すぎちゃん」とかです。

様々な傾向の一発屋さんがいるのがわかりました。

トレンドの推定

今回扱っているデータは芸人さんの数×時点のデータの多変量時系列となります。都合の良いものはないかと考えていましたが、古典的なVARではサンプルサイズ的にかなり苦しいと思い、Stanによるダイナミックパネルデータ分析などの事例はないか漁っていましたが、なかなかありませんでした。

松浦さんの『StanとRでベイズ統計モデリング (Wonderful R)』の241pに書かれている、モデル式12-8や12-9が今回のものに適しているなと思いましたが、コードを上げている方は見当たらなかったです。よしそれならば作ろうかと思った矢先、logics-of-blueさんのStan Advent Calendarの投稿、「Stanで推定する多変量時系列モデル」がかなりどんぴしゃな内容でしたので、コードを拝借してこの一発屋データの推定をしてみようと思います。

まずは、stanのコード

data {
  int T;                       // データ取得期間の長さ
  int performer_num;           // 芸人さんの数
  matrix[T, performer_num] y;  // 観測値
}

parameters {
  vector[T] x;        // 状態の推定値
  vector[performer_num] r; // 芸人さん毎のランダム効果
  real<lower=0> s_w;  // 過程誤差の標準偏差
  real<lower=0> s_v;  // 観測誤差の標準偏差
  real<lower=0> s_r;  // ランダム効果の標準偏差
  //vector[performer_num] sigma; //芸人さん毎の標準偏差
}

model {
  // 状態方程式に従い、状態が遷移する
  for(i in 2:T) {
    x[i] ~ normal(x[i-1], s_w);
  }
  
  // ランダム効果
  r ~ normal(0, s_r);
  
  // 観測方程式に従い、観測値が得られる
  for(i in 1:T) {
    for(j in 1:performer_num) {
      y[i, j] ~ normal(x[i] + r[j], s_v);
    }
  }
}

data {

int T; // データ取得期間の長さ

int performer_num; // 芸人さんの数

matrix[T, performer_num] y; // 観測値

}

parameters {

vector[T] x; // 状態の推定値

vector[performer_num] r; // 芸人さん毎のランダム効果

real<lower=0> s_w; // 過程誤差の標準偏差

real<lower=0> s_v; // 観測誤差の標準偏差

real<lower=0> s_r; // ランダム効果の標準偏差

//vector[performer_num] sigma; //芸人さん毎の標準偏差

}

model {

// 状態方程式に従い、状態が遷移する

for(i in 2:T) {

x[i] ~ normal(x[i-1], s_w);

}

// ランダム効果

r ~ normal(0, s_r);

// 観測方程式に従い、観測値が得られる

for(i in 1:T) {

for(j in 1:performer_num) {

y[i, j] ~ normal(x[i] + r[j], s_v);

}

そしてキックして結果を可視化するためのRコード

library(rstan)
library(bayesplot)

T <- trend_dataset_spread %>% select(-period) %>% nrow()
performer_num <- trend_dataset_spread %>% select(-period) %>% ncol()

data <- list(T = T,
             performer_num = performer_num ,
             y = trend_dataset_spread %>% select(-period))

fit <- stan(file = 'multivariate_time_series.stan',
            data = data,
            seed = 1, 
            iter = 30000,
            warmup = 10000,
            thin = 10
            )

mcmc_rhat(rhat(fit))

# データの整形
stan_df_1 <- fit %>% 
             rstan::extract() %$% x %>% 
             apply(2, quantile, probs = c(0.025, 0.5, 0.975)) %>% 
             t() %>% 
             cbind(1:nrow(trend_dataset_spread)) %>% 
             data.frame
# 列名の変更
colnames(stan_df_1) <- c("lwr", "fit", "upr", "time")
# 結果
head(stan_df_1, n = 3)

ggplot(data = trend_dataset_withcluster) + 
  ggtitle("推定結果（ピコ太郎）") +
  geom_line(aes(x = period, y = trend, color = keyword_en)) + 
  gghighlight(keyword_en == "Pico Taro", use_group_by = FALSE) + 
  geom_line(data = stan_df_1, 
            aes(x = time, y = fit), size = 1.2) +
  geom_ribbon(data = stan_df_1, 
              aes(x = time, ymin = lwr, ymax = upr), alpha = 0.3)

library(rstan)

library(bayesplot)

T <- trend_dataset_spread %>% select(-period) %>% nrow()

performer_num <- trend_dataset_spread %>% select(-period) %>% ncol()

data <- list(T = T,

performer_num = performer_num ,

y = trend_dataset_spread %>% select(-period))

fit <- stan(file = 'multivariate_time_series.stan',

data = data,

seed = 1,

iter = 30000,

warmup = 10000,

thin = 10

)

mcmc_rhat(rhat(fit))

# データの整形

stan_df_1 <- fit %>%

rstan::extract() %$% x %>%

apply(2, quantile, probs = c(0.025, 0.5, 0.975)) %>%

t() %>%

cbind(1:nrow(trend_dataset_spread)) %>%

data.frame

# 列名の変更

colnames(stan_df_1) <- c("lwr", "fit", "upr", "time")

# 結果

head(stan_df_1, n = 3)

ggplot(data = trend_dataset_withcluster) +

ggtitle("推定結果（ピコ太郎）") +

geom_line(aes(x = period, y = trend, color = keyword_en)) +

gghighlight(keyword_en == "Pico Taro", use_group_by = FALSE) +

geom_line(data = stan_df_1,

aes(x = time, y = fit), size = 1.2) +

geom_ribbon(data = stan_df_1,

aes(x = time, ymin = lwr, ymax = upr), alpha = 0.3)

そのまんま実行して、一発屋の時系列の中央値を可視化したらこんな感じになりました。一発屋のトレンドをうまく抽出できているのかなと思います。

今後の改良としては、階層性を持たせ、芸人さんごとのハイパーパラメータを持たせるとかなのですが、正月にでも取り組みたいと思います。（芸人さん以外のデータでやりたい。）

一方で、他にも多変量時系列で何かないか漁っていたのですが、Applied Time Series Analysis for Fisheries and Environmental Sciences : Dynamic factor analysisで紹介されている、Dynamic Factor Analysisというものが面白そうだなと思いました。
bayesdfaというパッケージを用いて、多変量時系列データに存在するであろうトレンドをStanを用いて推定することができるようです。元となった論文には各エリアごとのノルウェーロブスターの個体数のトレンドを推定し、3つのトレンドが発見されたとしています。ただ、同時点間のデータではないという点から今回のデータへの適用は不適切です。

同時点間に観測されていないデータであるという問題を認識した上で、このパッケージを使ってどんなトレンドを抽出できるのか試してみようと思います。

mod_3 = bayesdfa::fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)

rot = bayesdfa::rotate_trends(mod_3)
names(rot)

matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

mod_3 = bayesdfa::fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)

rot = bayesdfa::rotate_trends(mod_3)

names(rot)

matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

徐々に増えてから一気に落ちるトレンドや、一気に増えてから徐々に落ちるトレンドなどがうまく捉えれている気がします。
さらなる試行として、AICのような情報量基準である、Leave One Out Information Criterion (LOOIC)が最も低くなるトレンドの数を探索してみます。

mod_1 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 1)
mod_2 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 2)
mod_3 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)
mod_4 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 4)
mod_5 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 5)

mod_1 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 1)

mod_2 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 2)

mod_3 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 3)

mod_4 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 4)

mod_5 = fit_dfa(y = trend_dataset_spread %>% select(-period) %>% t(), num_trends = 5)

トレンド数を1から5まで指定して実行した結果、5の時が一番LOOICが低くなりました。

rot = bayesdfa::rotate_trends(mod_5)
names(rot)
matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

rot = bayesdfa::rotate_trends(mod_5)

names(rot)

matplot(t(rot$trends_mean), type = "l", lwd = 2, ylab = "mean trend")

まぁ、適切な使い方ではないのですが、徐々に増えてから一気に落ちるトレンドや、一気に増えてから徐々に落ちるトレンドなどが引き続き捉えれているようです。

今後の課題

・Stanによる多変量時系列のモデリングをしてみる。（Dynamic Panel分析とかもできると良い。少なくともStanのドキュメントにはない。）
・Dynamic Factor Analysisの適切な事例での適用をしてみる。

それでは、どうか良い年末をお過ごし下さい！
メリークリスマス！

参考情報

Introduction to gghighlight: Highlight ggplot’s Lines and Points with Predicates
{TSclust} ではじめる時系列クラスタリング
 Applied Time Series Analysis for Fisheries and Environmental Sciences 9.7 Dynamic factor analysis
読了：Montero & Vilar (2014) RのTSclustパッケージで時系列クラスタリング

rstanarmパッケージを使って簡単にベイズモデリングを実行する

はじめに

今回は、rstanarmというパッケージを用いて赤ワインデータを色々といじってみようと思います。
マーケティングの意思決定のための分析などでベイズ統計を使う場面が多々あるのですが、似たような属性のデータがあるのであれば、
一つ一つstanコードを書くのではなく、Rの関数でサクッと実行して試行錯誤していくという形に持っていけたらいいなぁと感じていました。
本気を出すところではstanを、ルーティンワーク的なタスクではrstanarmをみたいな形で使い分けれると良いのではないでしょうか。

rstanarmとは

バックエンドの計算をStanに実行させて、統計モデルの推定を行うためのパッケージ。R上でlm関数のように簡単にベイズ推定を行うことができる。対象ユーザーはベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー。
詳しくはこちら。

インストールする

まずはrstanarmのインストールするのですが、コケまくりました。そのため、バージョンを下げてみることにします。

devtools::install_version("rstanarm", version = "2.17.3", repos = "http://cran.us.r-project.org")

1	devtools::install_version("rstanarm", version = "2.17.3", repos = "http://cran.us.r-project.org")

ここに過去のバージョンがありますが、2.17.4だと動かなかったものの、2.17.3なら動きました。

rstanarmのサンプルを回してみる

今回は、以下の文献を参考にして、大人のIrisとも言える、ワインデータを扱い、質の高いワインかどうかを決める要素を探ります。
How to Use the rstanarm Package | Jonah Gabry and Ben Goodrich

こちらの文献には、ベイズ分析の4つのステップとして以下があげられています。

1.同時分布の特定（同時分布は事前分布と条件付きの尤度をかけ合わせたもの。）
2.MCMCで事後分布を描く
3.モデルがフィットしているか評価する
4.事後予測分布を描き、結果に影響を与える予測項を確認する。

これらのステップをできるだけ素早くできると良いですね。

まずはデータを読み込んで、スケーリングしておきます。（可視化結果は前回と同じなので、載せません。）
加えて、6点以上の評価であれば1を取る二項変数を作成しておきます。

library(tidyverse)
library(GGally)
library(rstanarm)
library(shinystan)
library(loo)

wine_dataset <- read.csv("http://ieor.berkeley.edu/~ieor265/homeworks/winequality-red.csv", sep=";" )

#可視化
wine_dataset_vis <- wine_dataset
ggpairs(wine_dataset_vis)

wine_dataset <- wine_dataset %>% mutate( y = if_else(condition = quality > 6, 1, 0 ))
wine_dataset <- wine_dataset %>% mutate_at(funs(scale(.)), .vars = c(1:11))
wine_dataset <- wine_dataset %>% select(-quality)

library(tidyverse)

library(GGally)

library(rstanarm)

library(shinystan)

library(loo)

wine_dataset <- read.csv("http://ieor.berkeley.edu/~ieor265/homeworks/winequality-red.csv", sep=";" )

#可視化

wine_dataset_vis <- wine_dataset

ggpairs(wine_dataset_vis)

wine_dataset <- wine_dataset %>% mutate( y = if_else(condition = quality > 6, 1, 0 ))

wine_dataset <- wine_dataset %>% mutate_at(funs(scale(.)), .vars = c(1:11))

wine_dataset <- wine_dataset %>% select(-quality)

GLMでロジスティックモデルを推定し、rstanarmで推定した結果と比較します。rstanarmでは傾きや切片の事前分布にスチューデントのt分布を、尤度にロジスティック分布を設定しています。

# GLM Estimation ----------------------------------------------------------
wine_glm_1 <- glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                    total.sulfur.dioxide + density + pH + sulphates + alcohol,
                        data = wine_dataset, family = binomial(link = "logit"))


# Bayesian Estimation with rstanarm ---------------------------------------
wine_bglm_1 <- stan_glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                          total.sulfur.dioxide + density + pH + sulphates + alcohol,
                              data = wine_dataset,
                              family = binomial(link = "logit"), 
                              prior = student_t(df = 7), 
                              prior_intercept = student_t(df = 7),
                              chains = 4, cores = 4, seed = 123)

#glmの結果との確認
cbind("model_glm" = coef(wine_glm_1),
      "model_bayesianglm" = coef(wine_bglm_1 ))

                       model_glm model_bayesianglm
(Intercept)          -2.81452789       -2.85848089
fixed.acidity         0.47871946        0.48115910
volatile.acidity     -0.46215347       -0.46939630
citric.acid           0.11060698        0.10547044
residual.sugar        0.33762729        0.32948671
chlorides            -0.41494490       -0.43932815
free.sulfur.dioxide   0.11318519        0.11910296
total.sulfur.dioxide -0.54377987       -0.56297397
density              -0.48655012       -0.47686448
pH                    0.03461116        0.03067029
sulphates             0.63563064        0.64074793
alcohol               0.80280901        0.82144931

# GLM Estimation ----------------------------------------------------------

wine_glm_1 <- glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol,

data = wine_dataset, family = binomial(link = "logit"))

# Bayesian Estimation with rstanarm ---------------------------------------

wine_bglm_1 <- stan_glm(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol,

data = wine_dataset,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4, cores = 4, seed = 123)

#glmの結果との確認

cbind("model_glm" = coef(wine_glm_1),

"model_bayesianglm" = coef(wine_bglm_1 ))

model_glm model_bayesianglm

(Intercept) -2.81452789 -2.85848089

fixed.acidity 0.47871946 0.48115910

volatile.acidity -0.46215347 -0.46939630

citric.acid 0.11060698 0.10547044

residual.sugar 0.33762729 0.32948671

chlorides -0.41494490 -0.43932815

free.sulfur.dioxide 0.11318519 0.11910296

total.sulfur.dioxide -0.54377987 -0.56297397

density -0.48655012 -0.47686448

pH 0.03461116 0.03067029

sulphates 0.63563064 0.64074793

alcohol 0.80280901 0.82144931

ほとんど係数の大きさが同じであることが確認できます。

ベイズ推定の良いところは事後分布から関心のある係数に関しての取りうる値などをシミュレーションできるところですが、
posterior_interval関数で簡単に計算することができます。

> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "fixed.acidity"), 2)
              2.5% 97.5%
fixed.acidity 0.03   0.9
> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "volatile.acidity"), 2)
                  2.5% 97.5%
volatile.acidity -0.76  -0.2

> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "fixed.acidity"), 2)

2.5% 97.5%

fixed.acidity 0.03 0.9

> round(posterior_interval(wine_bglm_1, prob = 0.95, pars = "volatile.acidity"), 2)

2.5% 97.5%

volatile.acidity -0.76 -0.2

肝心のMCMCの収束診断ですが、shinystanを使います。

やや余談ですが、他のデータセットでshinystanを用いた際に、予測結果にNAsが含まれている場合に、
shinystanが起動しないという問題があり、以下のようなエラー文が吐かれます。

Error in validate_y(y) : NAs not allowed in 'y'.

1	Error in validate_y(y) : NAs not allowed in 'y'.

調べたところ、こちらのgithubにあるように、

launch_shinystan(womensrole_bglm_1,ppd=FALSE)

1	launch_shinystan(womensrole_bglm_1,ppd=FALSE)

のように引数でppd=FALSEのように設定することで、立ち上げることができました。

3つの基準をクリアしているため、収束しています。

係数の分布についても可視化します。

rstanarmの良い点の一つとして、モデルのアップデートが容易に行える点があげられると思いますが、実際、以下のように先程のモデルに変数を追加して推定することができます。
今回は、alcoholを二乗したものを新しい変数として加えます。

#モデルのアップデート
(wine_bglm_2 <- update(wine_bglm_1, formula. = . ~ . + I(alcohol^2)))

1 2	#モデルのアップデート (wine_bglm_2 <- update(wine_bglm_1, formula. = . ~ . + I(alcohol^2)))

次に、looパッケージを用いて、更新したモデルと元のモデルの性能の比較を行います。
looパッケージは統計モデルの予測精度の指標として扱われる、WAIC(Widely Applicable Information Criterion)を計算するためのパッケージで、WICは事後分布から得られる対数尤度の平均や分散からなる値として表されます。looはleave-one-out cross-validationのleave-one-outの頭文字。

さっそく入れようと思ったところ、

Error: is.data.frame(data) || is.matrix(data) is not TRUE

1	Error: is.data.frame(data) \|\| is.matrix(data) is not TRUE

というエラーが出ました。
こちらでも議論されていましたが、

remove.packages("loo")
devtools::install_github("stan-dev/loo", ref = "v1.1.0")

1 2	remove.packages("loo") devtools::install_github("stan-dev/loo", ref = "v1.1.0")

でバージョンを2.0.0から1.1.0に落としたら動きました。

ここで、事後分布が特定のサンプルデータに対して敏感であるかどうかをlooパッケージを用いて可視化します。

loo_bglm_1 <- loo(wine_bglm_1)
loo_bglm_2 <- loo(wine_bglm_2)

par(mfrow = 1:2, mar = c(5,3.8,1,0) + 0.1, las = 3)
plot(loo_bglm_1, label_points = TRUE)
plot(loo_bglm_2, label_points = TRUE)

loo_bglm_1 <- loo(wine_bglm_1)

loo_bglm_2 <- loo(wine_bglm_2)

par(mfrow = 1:2, mar = c(5,3.8,1,0) + 0.1, las = 3)

plot(loo_bglm_1, label_points = TRUE)

plot(loo_bglm_2, label_points = TRUE)

縦軸のshape parameter kは推定の信頼性の指標とされ、大きければ大きいほど信頼できないと見なし、横軸は今回推定したワインデータのデータの番号で、左が元のモデル、右が変数を追加したモデルのものです。
どうやらどちらも0.4未満のkに収まっているようです。参考情報の事例では0.5を超えていましたが、moderate outliersと説明されていたので、今回の推定は問題ないと思われます。

続いてモデルの比較を行います。

> #モデルの比較
> compare_models(loo_bglm_1, loo_bglm_2)
elpd_diff        se 
      5.4       2.8

> #モデルの比較

> compare_models(loo_bglm_1, loo_bglm_2)

elpd_diff se

5.4 2.8

elpd_diffに関しては右のモデルの精度が高ければ正の値を、低ければ負の値を取るようになっています。標準誤差も返されます。
どうやら変数を追加したモデルの方が、ちょっとだけ良さそうです。

続いて、事後予測分布から、どの変数がどのように予測に影響を与えるのかを確かめます。
比較のためにデータを２つほど作成し、両者において一つだけ変数が違うという状況下での、予測される確率の比較を行います。

#パラメータを比較して変数が与える影響を確認する。
newdata <- data.frame(fixed.acidity = c(8.319637,8.319637) ,
                      volatile.acidity = c(0.5278205,0.5278205) ,
                      citric.acid = c(0.2709756,0.2709756),
                      residual.sugar = c(2.538806,2.538806),
                      chlorides = c(0.08746654,0.08746654) ,
                      free.sulfur.dioxide = c(15.87492,15.87492),
                      total.sulfur.dioxide = c(46.46779,46.46779),
                      density = c(0.9967467,0.9967467),
                      pH = c(3.311113,3.311113),
                      sulphates = c(0.6581488,0.6581488),
                      alcohol = c(15.42298,14.42298),
                      'I(alcohol^2)' = c(15.42298^2,14.42298^2) )

y_rep <- posterior_predict(wine_bglm_2, newdata)
summary(y_rep)
summary(apply(y_rep, 1, diff))


> summary(y_rep)
       1                2         
 Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.0000   1st Qu.:0.0000  
 Median :0.0000   Median :0.0000  
 Mean   :0.2402   Mean   :0.2717  
 3rd Qu.:0.0000   3rd Qu.:1.0000  
 Max.   :1.0000   Max.   :1.0000  
> 
> summary(apply(y_rep, 1, diff))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-1.0000  0.0000  0.0000  0.0315  0.0000  1.0000

#パラメータを比較して変数が与える影響を確認する。

newdata <- data.frame(fixed.acidity = c(8.319637,8.319637) ,

volatile.acidity = c(0.5278205,0.5278205) ,

citric.acid = c(0.2709756,0.2709756),

residual.sugar = c(2.538806,2.538806),

chlorides = c(0.08746654,0.08746654) ,

free.sulfur.dioxide = c(15.87492,15.87492),

total.sulfur.dioxide = c(46.46779,46.46779),

density = c(0.9967467,0.9967467),

pH = c(3.311113,3.311113),

sulphates = c(0.6581488,0.6581488),

alcohol = c(15.42298,14.42298),

'I(alcohol^2)' = c(15.42298^2,14.42298^2) )

y_rep <- posterior_predict(wine_bglm_2, newdata)

summary(y_rep)

summary(apply(y_rep, 1, diff))

> summary(y_rep)

1 2

Min. :0.0000 Min. :0.0000

1st Qu.:0.0000 1st Qu.:0.0000

Median :0.0000 Median :0.0000

Mean :0.2402 Mean :0.2717

3rd Qu.:0.0000 3rd Qu.:1.0000

Max. :1.0000 Max. :1.0000

> summary(apply(y_rep, 1, diff))

Min. 1st Qu. Median Mean 3rd Qu. Max.

-1.0000 0.0000 0.0000 0.0315 0.0000 1.0000

他の要素をある一定水準で保った際に、alcoholだけ1度下げることで、平均3%ほど高い評価が得られる確率が高まるという考察となります。

以上で、rstanarmの一連の使い方となるのですが、
一部の関数においては、階層ベイズモデルも行えるので、試してみようと思います。

ただ、階層ベイズにするにも、赤ワインのデータしかないので、グループ変数をどうにか作らないといけません。
あまりやりたくはありませんが、データがないので、説明変数を元にk-means(K=3)によるクラスタリングを行い、それをグループ変数とします。

stan_glmer関数を使えば、以下のような簡単な記述で定数項や係数がグループごとに異なるパラメータの分布に従うとする階層ベイズモデルを推定できます。

stan_glmer(非説明変数 ~ 変数 + ( 変数 - 1 | グループ変数 ),
                         data = wine_dataset_hc,
                         family = binomial(link = "logit"),
                         prior = student_t(df = 7), 
                         prior_intercept = student_t(df = 7),
                         chains = 4,cores = 4, seed = 123,iter = 500)

stan_glmer(非説明変数 ~ 変数 + ( 変数 - 1 | グループ変数 ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

今回は、切片だけがグループごとに異なるモデル、傾きだけがグループごとに異なるモデル、切片も傾きも異なるモデルを作成します。

# Hierachical_Models ------------------------------------------------------
wine_dataset_hc <- wine_dataset
wine_dataset_hc$wine_cluster <- as.numeric(kmeans(wine_dataset_hc %>% select(1:11) ,3)$cluster)
wine_dataset_hc %>% group_by(wine_cluster) %>% summarize(mean(as.numeric(y)),n())

# Hierarchical Models(Varying Intercept) -------------------------------------------------------
wine_bglm_3 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                           total.sulfur.dioxide + density + pH + sulphates + alcohol + ( 1 | wine_cluster ),
                         data = wine_dataset_hc,
                   family = binomial(link = "logit"),
                   prior = student_t(df = 7), 
                   prior_intercept = student_t(df = 7),
                   chains = 4,cores = 4, seed = 123,iter = 500)
coef(wine_bglm_3)

# Hierarchical Models(Varying Slope) ---------------------------------------------------------
wine_bglm_4 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                           total.sulfur.dioxide + density + pH + sulphates + alcohol +
                           ( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                               total.sulfur.dioxide + density + pH + sulphates + alcohol - 1 | wine_cluster ),
                         data = wine_dataset_hc,
                         family = binomial(link = "logit"),
                         prior = student_t(df = 7), 
                         prior_intercept = student_t(df = 7),
                         chains = 4,cores = 4, seed = 123,iter = 500)
coef(wine_bglm_4)

# Hierarchical Models(Varying Intercept and Slope) -------------------------------------------
wine_bglm_5 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                           total.sulfur.dioxide + density + pH + sulphates + alcohol +
                           ( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + 
                               total.sulfur.dioxide + density + pH + sulphates + alcohol | wine_cluster ),
                         data = wine_dataset_hc,
                         family = binomial(link = "logit"),
                         prior = student_t(df = 7), 
                         prior_intercept = student_t(df = 7),
                         chains = 4,cores = 4, seed = 123,iter = 500)
coef(wine_bglm_5)

loo_bglm_1 <- loo(wine_bglm_1)
loo_bglm_3 <- loo(wine_bglm_3)
loo_bglm_4 <- loo(wine_bglm_4)
loo_bglm_5 <- loo(wine_bglm_5)

# Hierachical_Models ------------------------------------------------------

wine_dataset_hc <- wine_dataset

wine_dataset_hc$wine_cluster <- as.numeric(kmeans(wine_dataset_hc %>% select(1:11) ,3)$cluster)

wine_dataset_hc %>% group_by(wine_cluster) %>% summarize(mean(as.numeric(y)),n())

# Hierarchical Models(Varying Intercept) -------------------------------------------------------

wine_bglm_3 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol + ( 1 | wine_cluster ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

coef(wine_bglm_3)

# Hierarchical Models(Varying Slope) ---------------------------------------------------------

wine_bglm_4 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol +

( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol - 1 | wine_cluster ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

coef(wine_bglm_4)

# Hierarchical Models(Varying Intercept and Slope) -------------------------------------------

wine_bglm_5 <- stan_glmer(y ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol +

( fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide +

total.sulfur.dioxide + density + pH + sulphates + alcohol | wine_cluster ),

data = wine_dataset_hc,

family = binomial(link = "logit"),

prior = student_t(df = 7),

prior_intercept = student_t(df = 7),

chains = 4,cores = 4, seed = 123,iter = 500)

coef(wine_bglm_5)

loo_bglm_1 <- loo(wine_bglm_1)

loo_bglm_3 <- loo(wine_bglm_3)

loo_bglm_4 <- loo(wine_bglm_4)

loo_bglm_5 <- loo(wine_bglm_5)

先程紹介した、looパッケージを使って、ベースとなるモデルとの比較を行います。

> #モデルの比較
> compare_models(loo_bglm_1, loo_bglm_3)
elpd_diff        se 
     -0.8       0.5 
> compare_models(loo_bglm_1, loo_bglm_4)
elpd_diff        se 
      1.2       2.5 
> compare_models(loo_bglm_1, loo_bglm_5)
elpd_diff        se 
      1.0       2.5

> #モデルの比較

> compare_models(loo_bglm_1, loo_bglm_3)

elpd_diff se

-0.8 0.5

> compare_models(loo_bglm_1, loo_bglm_4)

elpd_diff se

1.2 2.5

> compare_models(loo_bglm_1, loo_bglm_5)

elpd_diff se

1.0 2.5

うーん、残念ながらどのモデルもベースモデルよりも圧倒的に強いものは無さそうです。

感想

まだまだrstanarmの関数やら機能やら定義を全て把握しきれていないですが、そこらへんがクリアーになれば、これまでのstanでの推定業務において生産性が高まる可能性を感じました。
簡単な階層ベイズモデルくらいなら、非常に直感的に書ける点や、変数の追加によるモデルのアップデートが容易な点などはポイント高いなぁと思います。
とはいえ、実務としてマニュアルでstanコードを作成していくのは必須なので、このパッケージを使うことによって、stanコードの改善に時間をより一層割けるようになるなら、それが一番だと思いました。
あと、「ベイズ初めてです！」という新入りの方とかには慣れ親しんでもらうには良さそうですね。lm関数レベルで実行できてしまうので。
今回、mc-stan.orgの配下にあるページなどを漁る過程で、ベイズ推定結果の可視化などで知らないことにも色々と出会えたので、今後も読み進めていきます。

追記

2018-09-10: Stanを使って変数選択したいにprojpredというパッケージが紹介されており、これを使えば、情報量基準に従った変数選択を簡単に行えるそうです。こうなると、「ベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー」に限らず多くの人が幸せになれるパッケージなのかもしれませんね。

参考情報

Using the loo package (version >= 2.0.0)
Leave-one-outクロスバリデーションの２つのデメリット、からの解決方法
 stan_glmer | Bayesian Generalized Linear Models With Group-Specific Terms Via Stan
WAICを計算してみる
 Package ‘bayesm’
Hierarchical Partial Pooling for Repeated Binary Trials
Leave-one-out cross-validation for non-factorizable models
priors | Prior Distributions And Options
StanとRでベイズ統計モデリング (Wonderful R)

Stanで順序プロビット(Ordered Probit)の推定のためのメモ書き

最近はBayesian Statistics and Marketingという本に関心があって、そこで取り上げられているモデルをStanに落とし込めないか模索しています。
そこで順序プロビット（Ordered Probit）の推定が必要であることがわかったため、Stanでの適用事例を漁っていました。まだマーケティング事例への適用はうまくいってないですが、いったん順序プロビットを簡単にまとめて今後の備忘録としておきます。

順序プロビットとは

被説明変数yが連続潜在変数y∗に対応していると考えるとする。
潜在変数は観察できないが、被説明変数yは観察でき、これらの2つ変数の関係は次のように表される。
（今回扱うデータは3から8までの順序データのため、以下のような表記になる。）

この対応関係は閾値メカニズムと呼ばれている。
各被説明変数をとる確率は以下のように記され、プロビットでは正規分布を扱うため以下のようになる。

これらの選択確率からなる尤度関数を最大にしたものが順序プロビットの推定となる。（c0=−∞でc6=∞とする。σは1とする。）

このように、潜在的な順序関係を想定し、それを満たすように閾値とパラメータを推定する点において、潜在変数を用いたモデルの柔軟性の高さが感じられる。

なぜ順序プロビットを使うのか

マーケティング業務において扱うデータにおいて、NPSやアンケートなど順序尺度の質的変数が多いので、それらのデータを二値データに落とし込んだり、そのまま基数データとして扱うのではなく、適切に扱いたいというモチベーションがあります。加えて、順序尺度の質的変数をもとに予測する際は普通のOLSだと、今回のケースで3を下回ったり、8を超えたりする可能性があり、予測結果として使いにくいです。
アンケートの点数をそのまま被説明変数として回帰しているケースは、データアナリティクスにこだわりの無いメンバーとかであればままあることなので、順序プロビットの民主化というか、布教していきたいと思います。

今回扱うデータ

勝手ながら大人のirisだと思っているワインデータです。今回は赤ワインに絞って、品質に関する順序変数を被説明変数として、各変数との相関を見ていきます。
まずはGGallyパッケージのggpairs()関数を適用して傾向を掴みます。見にくいので是非コードを回して確かめてください。

データに関する説明はワインの味（美味しさのグレード）は予測できるか？（１）で丁寧になされていますので、ご確認ください。

モデル

データセットに含まれる全部を含めて順序プロビットで回帰してみようと思います。
つまり、「酒石酸濃度、酢酸濃度、クエン酸濃度
残留糖分濃度、塩化ナトリウム濃度、遊離亜硫酸濃度
総亜硫酸濃度、密度、ｐＨ、硫酸カリウム濃度、アルコール度数」
の全てを使って赤ワインの質への影響を見ていきます。

Stanコード

最初に、Stanのユーザーガイド2.17の138ページにあるOrdered Probitのサンプルコードを使ってみたのですが、
収束しなかったので、初期値を設定するか弱情報事前分布を導入するかの判断が必要となります。
そこで、jabranhamさんが係数が平均0で分散10の正規分布に従うとするサンプルコードを書かれていたので、そちらを使って推定します。
書き換えているところはデータの制約くらいです。

data{
  int<lower=1> N; // number of obs
  int<lower=3> J; // number of categories
  int<lower=2> K; // num of predictors
  int<lower=0,upper=10> y[N]; // outcome var 
  matrix[N, K] x; // predictor vars 
}

parameters{
  ordered[J-1] tau; // thresholds
  vector[K] beta; // beta coefficients 
}

model{
  vector[J] theta;
  vector[N] xB;
  beta ~ normal(0, 10);
  xB <- x*beta;
  for(n in 1:N){
    theta[1] <- 1 - Phi(xB[n]-tau[1]);
    for(j in 2:J-1)
      theta[j] <- Phi(xB[n]-tau[j-1]) - Phi(xB[n]-tau[j]);
    theta[J] <- Phi(xB[n] - tau[J-1]);
    y[n] ~ categorical(theta);
  }
}

data{

int<lower=1> N; // number of obs

int<lower=3> J; // number of categories

int<lower=2> K; // num of predictors

int<lower=0,upper=10> y[N]; // outcome var

matrix[N, K] x; // predictor vars

}

parameters{

ordered[J-1] tau; // thresholds

vector[K] beta; // beta coefficients

}

model{

vector[J] theta;

vector[N] xB;

beta ~ normal(0, 10);

xB <- x*beta;

for(n in 1:N){

theta[1] <- 1 - Phi(xB[n]-tau[1]);

for(j in 2:J-1)

theta[j] <- Phi(xB[n]-tau[j-1]) - Phi(xB[n]-tau[j]);

theta[J] <- Phi(xB[n] - tau[J-1]);

y[n] ~ categorical(theta);

}

StanをキックするためのRコード

推定結果の可視化を行うためのcommon.Rは松浦さんのGitHubにあるものになります。

library(tidyverse)
library(rstan)
library(GGally)
library(shinystan)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("http://ieor.berkeley.edu/~ieor265/homeworks/winequality-red.csv", sep=";" )

y <- wine_dataset$quality - 2
x <- as.matrix(wine_dataset %>% select(-quality))
x <- scale(x)

# Visualization -----------------------------------------------------------
ggpairs(wine_dataset)


# Estimation --------------------------------------------------------------
stanmodel <- stan_model(file = "orderedprobit.stan")


N <- nrow(x)
J <- 6L
K <- ncol(x)

data_customer_list_test <- list(N=N,
                                J=J,
                                K=K,
                                y=y,
                                x=x)

fit <- stan(file = "orderedprobit.stan",
            data = data_customer_list_test,
            iter = 1000,
            chains = 4)

summary(fit)
traceplot(fit)

# Convergence Check -------------------------------------------------------
launch_shinystan(fit)

# Result Plot -------------------------------------------------------------
source('common.R')

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)


param_names <- c('mcmc', colnames(wine_dataset %>% select(-quality)))

d_est <- data.frame(1:N_mcmc, ms$b)
colnames(d_est) <- param_names
d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])
d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')
d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()
p <- p + theme_bw(base_size=18)
p <- p + coord_flip()
p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')
p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)
p <- p + labs(x='parameter', y='value')
p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))
p

library(tidyverse)

library(rstan)

library(GGally)

library(shinystan)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("http://ieor.berkeley.edu/~ieor265/homeworks/winequality-red.csv", sep=";" )

y <- wine_dataset$quality - 2

x <- as.matrix(wine_dataset %>% select(-quality))

x <- scale(x)

# Visualization -----------------------------------------------------------

ggpairs(wine_dataset)

# Estimation --------------------------------------------------------------

stanmodel <- stan_model(file = "orderedprobit.stan")

N <- nrow(x)

J <- 6L

K <- ncol(x)

data_customer_list_test <- list(N=N,

J=J,

K=K,

y=y,

x=x)

fit <- stan(file = "orderedprobit.stan",

data = data_customer_list_test,

iter = 1000,

chains = 4)

summary(fit)

traceplot(fit)

# Convergence Check -------------------------------------------------------

launch_shinystan(fit)

# Result Plot -------------------------------------------------------------

source('common.R')

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', colnames(wine_dataset %>% select(-quality)))

d_est <- data.frame(1:N_mcmc, ms$b)

colnames(d_est) <- param_names

d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])

d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')

d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()

p <- p + theme_bw(base_size=18)

p <- p + coord_flip()

p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')

p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)

p <- p + labs(x='parameter', y='value')

p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))

結果

まず、MCMCが収束したかどうかの判断ですが、ShinyStanに従うものとします。

ShinyStanによる収束診断をクリアできています。

続いて、推定したパラメータです。

係数の符号がはっきりと分かれている、赤ワインの品質に影響を与えそうな変数としては、volatile acidity（酢酸濃度）、chlorides（塩化ナトリウム濃度）、total.sulfur.dioxide（総亜硫酸濃度）、sulphates（硫酸カリウム濃度）、alcohol（アルコール度数）のようです。

最後に、推定した閾値です。

1と2の閾値が近く、2と3の開きが大きく、あとは比較的均等のようです。

比較

質的変数をそのまま重回帰した際の結果ですが、符号やその大小はあまり変わらないです。やはり予測の際にどの順序尺度の値に対応するかがわかるのが使う利点だと思います。

> dataset <- data.frame(cbind(y, x))
> multiple_regresion_model <- lm(formula = y ~ .,data = dataset)
> summary(multiple_regresion_model)

Call:
lm(formula = y ~ ., data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.68911 -0.36652 -0.04699  0.45202  2.02498 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           3.63602    0.01621 224.372  < 2e-16 ***
fixed.acidity         0.04351    0.04518   0.963   0.3357    
volatile.acidity     -0.19403    0.02168  -8.948  < 2e-16 ***
citric.acid          -0.03556    0.02867  -1.240   0.2150    
residual.sugar        0.02303    0.02115   1.089   0.2765    
chlorides            -0.08821    0.01973  -4.470 8.37e-06 ***
free.sulfur.dioxide   0.04562    0.02271   2.009   0.0447 *  
total.sulfur.dioxide -0.10739    0.02397  -4.480 8.00e-06 ***
density              -0.03375    0.04083  -0.827   0.4086    
pH                   -0.06386    0.02958  -2.159   0.0310 *  
sulphates             0.15533    0.01938   8.014 2.13e-15 ***
alcohol               0.29433    0.02822  10.429  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.648 on 1587 degrees of freedom
Multiple R-squared:  0.3606,    Adjusted R-squared:  0.3561 
F-statistic: 81.35 on 11 and 1587 DF,  p-value: < 2.2e-16

> dataset <- data.frame(cbind(y, x))

> multiple_regresion_model <- lm(formula = y ~ .,data = dataset)

> summary(multiple_regresion_model)

Call:

lm(formula = y ~ ., data = dataset)

Residuals:

Min 1Q Median 3Q Max

-2.68911 -0.36652 -0.04699 0.45202 2.02498

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 3.63602 0.01621 224.372 < 2e-16 ***

fixed.acidity 0.04351 0.04518 0.963 0.3357

volatile.acidity -0.19403 0.02168 -8.948 < 2e-16 ***

citric.acid -0.03556 0.02867 -1.240 0.2150

residual.sugar 0.02303 0.02115 1.089 0.2765

chlorides -0.08821 0.01973 -4.470 8.37e-06 ***

free.sulfur.dioxide 0.04562 0.02271 2.009 0.0447 *

total.sulfur.dioxide -0.10739 0.02397 -4.480 8.00e-06 ***

density -0.03375 0.04083 -0.827 0.4086

pH -0.06386 0.02958 -2.159 0.0310 *

sulphates 0.15533 0.01938 8.014 2.13e-15 ***

alcohol 0.29433 0.02822 10.429 < 2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.648 on 1587 degrees of freedom

Multiple R-squared: 0.3606, Adjusted R-squared: 0.3561

F-statistic: 81.35 on 11 and 1587 DF, p-value: < 2.2e-16

考察

マーケティングにおいて、順序尺度の質的変数を扱う際に順序プロビットを積極的に使っていきたいと思いますが、アンケート分析を行う際は、ユーザーごとの評点の癖が点数に影響を与えている可能性があります。
そのため、点数の付き方がユーザーごとに違うとする階層モデルへの拡張を今後行っていくのが面白いと思いますし、実際に研究されている論文があります。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
stan-examples/limited-dv/oprobit.stan
stan-dev/stan users-guide-2.18.0.pdf
Stanによる順序ロジット回帰
 第9章順序選択モデル：年金投資選択問題
 Wine Quality Data Set
ワインの味（美味しさのグレード）は予測できるか？（１）
世界一簡単な収束[シナイ]Stanコード
 RStanとShinyStanによるベイズ統計モデリング入門

おまけ

数式をブログに載せる際は、こちら
Online LaTeX Equation Editor – create, integrate and download
でインタラクティブに数式を作成し、その結果を
QuickLaTex Publish Math on the Web without compromising quality
に貼り付けて画像を出力しています。

Bayesian Statistics and Marketing – 混合ガウス×階層モデルのマーガリン購買データへの適用

前回の分析では、価格への反応係数の事前分布が正規分布を仮定したモデルを用いていましたが、事後分布から多峰性が観察されました。そこで今回は、各個人の価格への反応係数の事前分布が混合ガウス分布に従うとした場合の事例を扱いたいと思います。

データのおさらい

データ自体は前回のブログと同じですが、先日のTokyo.Rで松浦さんがオススメしていたGGallyパッケージのggpairs関数を用いて、今回扱うデータを可視化してみます。

まず、購買したマーガリンのブランド選択(6ブランド分)と、購入した価格（ドル）からなるデータセットの可視化をすると以下のようになります。

1つ目のブランドが最も多く選択されているようです。購入価格の分布は、ブランドによって多峰性がありそうです。

続いて、家計ごとの属性（家族構成、学歴、職位、退職の有無など）からなるデータセットの可視化をすると以下のようになります。

ほとんどダミー変数なので面白みには欠いていますね。ホワイトカラーの家計が多く、退職していない家計が多く、学歴が低い人が多いようです。年収に関しては対数正規分布に従ってそうです。

モデル

前回と同様に、「Bayesian Statistics and Marketing」の5章に載っている混合ガウスを想定した階層ベイズモデルを扱います。

yi ∼ MNL(Xi, βi) は購買レコードごとの意思決定が多項ロジスティック回帰モデルに従うということを意味し、βiは(説明変数の数×1)のベクトルとなります。
βi = Z∆[i,] + ui は購買レコードごとの価格にかかってくる係数で、その係数が家計の属性データに係数Δ(家計の数×属性データからなる説明変数の数)をかけ合わせたものと潜在的な項の和となります。なお、ここでの属性データには定数項を含めていません。係数Δは平均deltabar、分散A_delta^-1の多変量正規分布に従います。一方で、潜在的な項は平均µind、分散Σindの多変量正規分布に従います。この平均や分散に振られているindが、混合正規分布のパラメータの通し番号となり、多項分布に従います。この多項分布の割当確率がディリクレ分布に従います。最後に混合正規分布の各パラメータは、平均は正規分布に分散は逆ウィシャート分布に従います。

前回との大きな違いは、価格の係数の一部である、潜在的な項において、混合正規分布が仮定されているところになります。
ちなみに、モデルに関しての詳細はbayesmパッケージのマニュアルの61ページ目に記載されていました。

今回のモデルのDAG

Pythonのdaftというモジュールを使うことで、非常に簡単に今回のモデルのDAG(有向非巡回グラフ)を描くことができます。

今回はこちらのPythonコードで描けました。

import daft
from matplotlib import rc

pgm = daft.PGM(shape=[6,6])

# Nodes
pgm.add_node(daft.Node("pvec", r"pvec",1, 2)) # 名前　ラベル　座標(横、縦)
pgm.add_node(daft.Node("ind", r"ind",2, 2))
pgm.add_node(daft.Node("sigma", r"$\Sigma_k$",3, 1))
pgm.add_node(daft.Node("mu", r"$\mu_k$",3, 3))
pgm.add_node(daft.Node("theta", r"$\theta_i$",4, 2))
pgm.add_node(daft.Node("delta", r"$\Delta$",4, 0.5))
pgm.add_node(daft.Node("y", r"$y_i$",5, 2,observed=True))

# Edges
pgm.add_edge("pvec", "ind")
pgm.add_edge("ind", "sigma")
pgm.add_edge("ind", "mu")
pgm.add_edge("sigma", "mu")
pgm.add_edge("sigma", "theta")
pgm.add_edge("mu", "theta")
pgm.add_edge("delta", "theta")
pgm.add_edge("theta", "y")

pgm.render()
pgm.figure.savefig("Hierachical_Models_ForMixtureNormal.png")

import daft

from matplotlib import rc

pgm = daft.PGM(shape=[6,6])

# Nodes

pgm.add_node(daft.Node("pvec", r"pvec",1, 2)) # 名前　ラベル　座標(横、縦)

pgm.add_node(daft.Node("ind", r"ind",2, 2))

pgm.add_node(daft.Node("sigma", r"$\Sigma_k$",3, 1))

pgm.add_node(daft.Node("mu", r"$\mu_k$",3, 3))

pgm.add_node(daft.Node("theta", r"$\theta_i$",4, 2))

pgm.add_node(daft.Node("delta", r"$\Delta$",4, 0.5))

pgm.add_node(daft.Node("y", r"$y_i$",5, 2,observed=True))

# Edges

pgm.add_edge("pvec", "ind")

pgm.add_edge("ind", "sigma")

pgm.add_edge("ind", "mu")

pgm.add_edge("sigma", "mu")

pgm.add_edge("sigma", "theta")

pgm.add_edge("mu", "theta")

pgm.add_edge("delta", "theta")

pgm.add_edge("theta", "y")

pgm.render()

pgm.figure.savefig("Hierachical_Models_ForMixtureNormal.png")

stanコード

今回扱ったstanコードとなります。誤りがある場合はお知らせしていただけると幸いです。

data{
  int<lower=0> N_x; // 購買レコードの数
  int<lower=0> N_z; // 家計の数
  int<lower=0> p_x; // 購買レコードの項目数
  int<lower=0> p_z; // 家計の属性データの項目数
  
  int y[N_x]; // 選択肢
  matrix[N_x, p_x] X; // 説明変数
  matrix[N_z, p_z] Z; // 家計の属性データ
  int<lower=0> hhid[N_x];  // 家計ID

  int K; // 混合ガウス分布の要素数
}

transformed data{
  real nu;
  matrix[p_x, p_x] I; // 購買レコードの説明変数の数の正方行列
  matrix[p_x, p_x] J; // 属性データの説明変数の数の正方行列
  
  nu = p_x + 3; // 説明変数の項に3を足す
  I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成
  J = diag_matrix(rep_vector(1, p_x));
}

parameters{
  vector[p_x] theta_ast[N_z]; // 説明変数の数だけある、購買ごとのパラメータ
  matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ
  vector[p_x] u[N_z]; // 購買レコードごとのパラメータ
  vector[p_x] mu[K]; // 混合分布を構成する平均値
  cov_matrix[p_x] Sigma[K]; // 共分散行列
  simplex[K] pi[N_z]; // シンプレックス（各要素が[0,1]の範囲で合計が1という条件を満たす。）
}

transformed parameters{
  vector[p_x] theta[N_z]; #家計の数だけの係数ベクトル

  for(i in 1:N_z){ // 家計の数だけ繰り返す
    theta[i] = theta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる 
  }
}

model{
  real ps[K]; // 混合正規分布の対数

  for(i in 1:N_x){ // 購買レコードの数だけ繰り返す
    y[i] ~ categorical(softmax(theta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う

    pi[hhid[i]] ~ dirichlet(rep_vector(2, K));
    for(k in 1:K){ // 混合分布の構成要素の数だけ繰り返す
      ps[k] = log(pi[hhid[i]][k]) + multi_normal_lpdf(theta_ast[hhid[i]] | mu[k], Sigma[k]);
    }
    target += log_sum_exp(ps); // 離散パラメータを消去した形で対数尤度を表現する際に必要な計算（周辺化消去）。  
  }

  for(k in 1:K){ // 混合パラメータの数だけ繰り返す
    mu[k] ~ multi_normal(rep_vector(0, p_x), 100*Sigma[k]);
    Sigma[k] ~ inv_wishart(nu, nu*I); 
  }

  for(i in 1:p_z){ // 家計の属性データの数だけ繰り返す
    Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), 100*J); 
  }
}

data{

int<lower=0> N_x; // 購買レコードの数

int<lower=0> N_z; // 家計の数

int<lower=0> p_x; // 購買レコードの項目数

int<lower=0> p_z; // 家計の属性データの項目数

int y[N_x]; // 選択肢

matrix[N_x, p_x] X; // 説明変数

matrix[N_z, p_z] Z; // 家計の属性データ

int<lower=0> hhid[N_x]; // 家計ID

int K; // 混合ガウス分布の要素数

}

transformed data{

real nu;

matrix[p_x, p_x] I; // 購買レコードの説明変数の数の正方行列

matrix[p_x, p_x] J; // 属性データの説明変数の数の正方行列

nu = p_x + 3; // 説明変数の項に3を足す

I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成

J = diag_matrix(rep_vector(1, p_x));

}

parameters{

vector[p_x] theta_ast[N_z]; // 説明変数の数だけある、購買ごとのパラメータ

matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ

vector[p_x] u[N_z]; // 購買レコードごとのパラメータ

vector[p_x] mu[K]; // 混合分布を構成する平均値

cov_matrix[p_x] Sigma[K]; // 共分散行列

simplex[K] pi[N_z]; // シンプレックス（各要素が[0,1]の範囲で合計が1という条件を満たす。）

}

transformed parameters{

vector[p_x] theta[N_z]; #家計の数だけの係数ベクトル

for(i in 1:N_z){ // 家計の数だけ繰り返す

theta[i] = theta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる

}

model{

real ps[K]; // 混合正規分布の対数

for(i in 1:N_x){ // 購買レコードの数だけ繰り返す

y[i] ~ categorical(softmax(theta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う

pi[hhid[i]] ~ dirichlet(rep_vector(2, K));

for(k in 1:K){ // 混合分布の構成要素の数だけ繰り返す

ps[k] = log(pi[hhid[i]][k]) + multi_normal_lpdf(theta_ast[hhid[i]] | mu[k], Sigma[k]);

}

target += log_sum_exp(ps); // 離散パラメータを消去した形で対数尤度を表現する際に必要な計算（周辺化消去）。

}

for(k in 1:K){ // 混合パラメータの数だけ繰り返す

mu[k] ~ multi_normal(rep_vector(0, p_x), 100*Sigma[k]);

Sigma[k] ~ inv_wishart(nu, nu*I);

}

for(i in 1:p_z){ // 家計の属性データの数だけ繰り返す

Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), 100*J);

}

stanをキックするためのコードです。先人が書かれた混合ガウスのスクリプトをHMCで実行した際に26時間ほどかかったので、今回はより複雑なモデルであることから、変分ベイズ法による推定を行ってみることにしました。松浦さんの教科書にあるように、vb関数を用いて変分ベイズ推論を行っています。

library(bayesm)
library(tidyverse)
library(rstan)
library(GGally)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。
hhid_selected <- margarine$choicePrice %>% 
                                  filter(choice %in% c(1,2,3,4,5,7)) %>% 
                                  group_by(hhid) %>% 
                                  summarise(purc_cnt = n()) %>% 
                                  filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。
choicePrice.selected <- margarine$choicePrice %>% 
                                             filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)
#並べにくいので7を6に置き換える。
choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出
demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ
N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）
p <- n_distinct(choicePrice.selected$choice)

#被説明変数
y <- choicePrice.selected$choice

#説明変数
X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く
Z <- demos.selected %>% 
                                   select(-hhid)

#可視化
ggpairs(choicePrice.selected %>% select(2,3,4,5,6,7,9))
ggpairs(Z)


#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。
hhid_index <- demos.selected %>%
                                                   select(hhid) %>% 
                                                   mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする
hhid_x <- choicePrice.selected %>% 
                                                    select(hhid) %>% 
                                                    left_join(hhid_index)

#stanで扱うデータリストの作成
d.dat <- list(N_x=nrow(X), N_z=nrow(Z), 
              p_x=ncol(X), p_z=ncol(Z),
              y=y, X=X, Z=Z,
              hhid = hhid_x$ind,
              K = 3)

#推定
stanmodel <- stan_model(file = "Hierarchical_Mixture.stan")
d.fit_vb <- vb(stanmodel, data=d.dat, seed=123)

library(bayesm)

library(tidyverse)

library(rstan)

library(GGally)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。

hhid_selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7)) %>%

group_by(hhid) %>%

summarise(purc_cnt = n()) %>%

filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。

choicePrice.selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)

#並べにくいので7を6に置き換える。

choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出

demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ

N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）

p <- n_distinct(choicePrice.selected$choice)

#被説明変数

y <- choicePrice.selected$choice

#説明変数

X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く

Z <- demos.selected %>%

select(-hhid)

#可視化

ggpairs(choicePrice.selected %>% select(2,3,4,5,6,7,9))

ggpairs(Z)

#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。

hhid_index <- demos.selected %>%

select(hhid) %>%

mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする

hhid_x <- choicePrice.selected %>%

select(hhid) %>%

left_join(hhid_index)

#stanで扱うデータリストの作成

d.dat <- list(N_x=nrow(X), N_z=nrow(Z),

p_x=ncol(X), p_z=ncol(Z),

y=y, X=X, Z=Z,

hhid = hhid_x$ind,

K = 3)

#推定

stanmodel <- stan_model(file = "Hierarchical_Mixture.stan")

d.fit_vb <- vb(stanmodel, data=d.dat, seed=123)

実行結果

今回はK={1,3,5,10}の混合要素数で推定を行いました。

K=1のケース

K=3のケース

K=5のケース

K=10のケース

Kが小さいと散らばりが比較的小さそうに見えます。逆に、Kが大きくなると散らばりが出てくるようです。

おわりに

　２回に渡って、マーケティングデータを用いたベイズ統計モデリングを学んできましたが、数式からstanコードに落とし込む作業の際に、stanの関数をある程度知らないとやりにくいという至極当然なことを実感しました。これまでに使ってきたstanの関数は限られたものしか扱っていなかったと言えます。特に今回は離散パラメータをstanで扱うパーツがあったので、先行研究や松浦さんの本を読みながらの手探りが多かったです。
　あと、複雑な階層ベイズモデルを扱う際に、頭の中を整理しないと手が止まってしまう感じがあったので、数式と対応するコードを横に並べながら進めました。
　マーケティングにおいては、顧客の属性ごとに多峰性のあるような事例を扱うことが多く、かつ各々のサンプルサイズも期待できないことが多いので、今回の学びを分析業務で試してみたいと思います。

参考情報

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)
Package ‘bayesm’
Multivariate Gaussian Mixture Model done properly
StanとRでベイズ統計モデリング (Wonderful R)
機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)
daftでグラフィカルモデル

目次

はじめに

1.『ベイズモデリングの世界』（岩波書店）

2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』（コロナ社）

3.『構造的因果モデルの基礎』（共立出版）

4.『現場で使える！PyTorch開発入門 深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』（翔泳社）

5.『作ってわかる! アンサンブル学習アルゴリズム入門』（シーアンドアール研究所）

6.『数理統計学―基礎から学ぶデータ解析』（内田老鶴圃）

7.『44の例題で学ぶ統計的検定と推定の解き方』（オーム社）

8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』（技術評論社）

9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)』

10.『scikit-learnとTensorFlowによる実践機械学習』（オライリージャパン）

11.『AIアルゴリズムマーケティング 自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)』

12.『入門 統計的因果推論』（朝倉書店）

13.『実践 ベイズモデリング -解析技法と認知モデル-』（朝倉書店）

14.『機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)』

15.『その問題、数理モデルが解決します』（ベレ出版）

16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』（CCCメディアハウス）

17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』（北大路書房）

18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』（共立出版）

19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』（SBクリエイティブ）

20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』（オライリージャパン）

[+α]『プログラマのためのGoogle Cloud Platform入門 サービスの全体像からクラウドネイティブアプリケーション構築まで』（翔泳社）

終わりに

はじめに

データについて

分析方針

分析概要

データ整形

訓練データとテストデータの分割

クロスバリデーション用のデータの作成

stanの実行

クロスバリデーション結果の出力

テストデータでの予測

非階層モデルとの比較

おわりに

参考情報

はじめに

生存時間分析とは

生存時間分析をマーケティングで使う際の用途

先行研究

生存時間分析で使えるデータセット

Stanでの実行例

おわりに

参考文献

はじめに

項目反応理論（Item Response Theory：IRT）とは

回答が2値変数のモデル

回答が順序変数のモデル（まずい < まぁまぁ < おいしい）

ビジネスでの適用可能性について

データ

モデルの推定

結果の解釈

おわりに

参考情報

はじめに

分析工程

データの収集

データの整形

可視化

時系列クラスタリングの適用

トレンドの推定

今後の課題

参考情報

はじめに

rstanarmとは

インストールする

rstanarmのサンプルを回してみる

感想

追記

参考情報

順序プロビットとは

なぜ順序プロビットを使うのか

今回扱うデータ

モデル

Stanコード

StanをキックするためのRコード

結果

比較

4.『現場で使える！PyTorch開発入門深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』（翔泳社）

11.『AIアルゴリズムマーケティング自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)』

12.『入門統計的因果推論』（朝倉書店）

13.『実践ベイズモデリング -解析技法と認知モデル-』（朝倉書店）

14.『機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)』

[+α]『プログラマのためのGoogle Cloud Platform入門サービスの全体像からクラウドネイティブアプリケーション構築まで』（翔泳社）