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でモデルを識別した結果、キャリーオーバー効果や形状効果を想定した設定が選ばれた。
残差に自己相関があったため、今回のモデルはまだ不完全で改善の余地がある。

[R]ボージョレ・ヌーボーのコメントに対してLDATSパッケージを使って時系列トピックモデルを扱う

はじめに

先日、某勉強会でLTをしました。その際に10秒だけ紹介したRのパッケージについて記事を書いてみようと思います。

LDATSパッケージについて

時系列でのトピックモデルを推定することができるパッケージです。
やっていることとしてはLDAでトピックを推定して次元を減らし、そのトピックの多変量時系列に関してベイズ手法による変化点検知のためのパラメータ推定を行っているようです。GitHubの該当しそうなソースコードに多変量のデータに対するsoftmax関数での回帰をやっているとの記述がある。（multinomial Bayesian Time Series analysis）

元となっている論文を見る限り、BoW(Bag of Words)を想定して作っておらず、20~30程度のグループからなるデータに対して適用するのがちょうど良いです。アクセスログのページカテゴリや、マーケティングの顧客セグメントであればそんなに数は多くないので扱いやすいと思います。

データ

Webサイトから集めてきたボージョレ・ヌーボーのキャッチコピー14年分を今回は扱います。実は販売店側のキャッチコピーとワイン委員会が決めた評価が存在します。私の知っている世界は販売店側のキャッチコピーだけでした。

試してみた

今回はとにかく動くことだけを考えて、汚いコードとなっております。やっていることとしては、キャッチコピーを販売側とワイン委員会側のものを一つにつないで、数字を正規表現で「数字」に変換し、RMeCabで形態素解析をし、LDATS向けの形式のデータを作成していきます。
途中で、日本語の文字化け問題を回避するためにGoogle翻訳を使って単語名を置き換えています。
1時系列につき1文書となるようにデータを作っていく必要があるのですが、今回はボージョレ・ヌーボーのキャッチコピーなので最初から1時系列につき1文書となっているため都合が良いです。
データとソースコードはこちら。

library(tidyverse)
library(RMeCab)
library(LDATS)

wine_rating <- read_csv(file = "TimeseriesLDA/dataset.csv")

wine_rating <- wine_rating %>% filter(!is.na(sales_catch_copy_text),
                                      !is.na(ratings),
                                      year > 1999,
                                      !is.na(commission_text))

wine_rating$sales_catch_copy_text <- gsub(pattern = "[0-9]",
                                          replacement = "数字",
                                          x = wine_rating$sales_catch_copy_text)
wine_rating$commission_text <- gsub(pattern = "[0-9]",
                                    replacement = "数字",
                                    x = wine_rating$commission_text)

wine_rating$bind_text <- paste0(wine_rating$sales_catch_copy_text ,
                                wine_rating$commission_text) 


# Bag of wordsの生成
res <- docMatrixDF(wine_rating$bind_text,minFreq=3)
res <- data.frame(res)
# View(rownames(res))

# Google Spread Sheetの=GOOGLETRANSLATE(C18,"ja","en")で変換した英語のデータを読み込む
translate_df <- read_csv(file = "TimeseriesLDA/translate.csv",col_names = FALSE)
colnames(translate_df) <- c("word_ja", "word_en")

word_translate <- data.frame(word_ja=rownames(res))
word_translate <- word_translate %>% left_join(translate_df, by="word_ja")

rownames(res) <- word_translate$word_en

word_vector <- row.names(res)
colnames(res) <- wine_rating$year

# LDATSで扱えるデータ構造を作成
for (i in 1:nrow(res)) {
  nam <- paste( word_vector[i], sep = "")
  assign(nam, as.integer(res[i, ]))
}

# ここでの変数名が可視化の際に表示される
document_term_table <- data.frame(list(sa=sa,
                                       fruits=fruits,
                                       Greatness=Greatness,
                                       Can=Can,
                                       workmanship=workmanship,
                                       taste=taste,
                                       quality=quality,
                                       Year=Year,
                                       Thenumbers=Thenumbers,
                                       Highest=Highest,
                                       fruit=fruit,
                                       Great=Great,
                                       delicate=delicate,
                                       complexity=complexity,
                                       rich=rich,
                                       past=past,
                                       fragrance=fragrance
                                       ))

# 共変量データセット
document_covariate_table <- data.frame(list(year=as.integer(colnames(res))),
                                       list(rating=as.integer(wine_rating$ratings)))

test_set <- list(document_term_table=document_term_table,
                 document_covariate_table=document_covariate_table)

# 時系列トピックモデルの実行
r_LDATS <- LDA_TS(test_set,
                  topics = 3:6, 
                  nseeds = 2,
                  formulas = ~1,  
                  nchangepoints = 1:2,
                  timename = "year")
# 対数尤度などの出力
print(r_LDATS)

# 時系列トピックモデルの可視化（先行研究に準拠）
plot(r_LDATS)

library(tidyverse)

library(RMeCab)

library(LDATS)

wine_rating <- read_csv(file = "TimeseriesLDA/dataset.csv")

wine_rating <- wine_rating %>% filter(!is.na(sales_catch_copy_text),

!is.na(ratings),

year > 1999,

!is.na(commission_text))

wine_rating$sales_catch_copy_text <- gsub(pattern = "[0-9]",

replacement = "数字",

x = wine_rating$sales_catch_copy_text)

wine_rating$commission_text <- gsub(pattern = "[0-9]",

replacement = "数字",

x = wine_rating$commission_text)

wine_rating$bind_text <- paste0(wine_rating$sales_catch_copy_text ,

wine_rating$commission_text)

# Bag of wordsの生成

res <- docMatrixDF(wine_rating$bind_text,minFreq=3)

res <- data.frame(res)

# View(rownames(res))

# Google Spread Sheetの=GOOGLETRANSLATE(C18,"ja","en")で変換した英語のデータを読み込む

translate_df <- read_csv(file = "TimeseriesLDA/translate.csv",col_names = FALSE)

colnames(translate_df) <- c("word_ja", "word_en")

word_translate <- data.frame(word_ja=rownames(res))

word_translate <- word_translate %>% left_join(translate_df, by="word_ja")

rownames(res) <- word_translate$word_en

word_vector <- row.names(res)

colnames(res) <- wine_rating$year

# LDATSで扱えるデータ構造を作成

for (i in 1:nrow(res)) {

nam <- paste( word_vector[i], sep = "")

assign(nam, as.integer(res[i, ]))

}

# ここでの変数名が可視化の際に表示される

document_term_table <- data.frame(list(sa=sa,

fruits=fruits,

Greatness=Greatness,

Can=Can,

workmanship=workmanship,

taste=taste,

quality=quality,

Year=Year,

Thenumbers=Thenumbers,

Highest=Highest,

fruit=fruit,

Great=Great,

delicate=delicate,

complexity=complexity,

rich=rich,

past=past,

fragrance=fragrance

))

# 共変量データセット

document_covariate_table <- data.frame(list(year=as.integer(colnames(res))),

list(rating=as.integer(wine_rating$ratings)))

test_set <- list(document_term_table=document_term_table,

document_covariate_table=document_covariate_table)

# 時系列トピックモデルの実行

r_LDATS <- LDA_TS(test_set,

topics = 3:6,

nseeds = 2,

formulas = ~1,

nchangepoints = 1:2,

timename = "year")

# 対数尤度などの出力

print(r_LDATS)

# 時系列トピックモデルの可視化（先行研究に準拠）

plot(r_LDATS)

こちらは論文の図と同じものだとドキュメントの説明にあったので、論文の説明を見る限り、表すものとしては以下のようです。

一番上の積み上げグラフはトピックごとの単語の割合を表しています。
二番目の折れ線グラフはLDAによって推定されたトピックの時系列推移です。
三番目のヒストグラムは二番目の時系列における変化点を集計したものです。
四番目の折れ線グラフはモデルが推定したトピック割合の変化点の前後での推移です。

今回の図では文字が潰れていて見にくいですが、

トピック1はボキャブラリーが比較的リッチなコメント（「フルーティー」「フレグランス」「複雑」）
トピック2は数字を用いたコメント（「何年に一度の！」みたいな）
トピック3はボキャブラリーが貧相なコメント（「すごい！」みたいな）

のようです。
二番目の折れ線グラフを見る限り、周期的に数字を用いたコメントが現れているように思われます。四番目の折れ線グラフの変化点を見る限り、近年は数字を用いたコメントが相対的に減ってきて、リッチなボキャブラリーになってきているようです。

おわりに

時系列トピックモデルをカジュアルに試せる面白そうなパッケージだなと思い、LDATSパッケージを触ってみましたが、そもそもBoWなどを想定して作られているパッケージではないので、単語数が多いような分析ではそもそも可視化ができず使いにくいだろうなと思いました。マーケティングなどでユーザーのセグメントの推移を分析したい場合などにちょうど良いのだろうと思われます。

参考情報

[1] Long‐term community change through multiple rapid transitions in a desert rodent community
[2] Latent Dirichlet Allocation coupled with Bayesian Time Series analyses
[3] Package ‘LDATS’

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ではじめるベイズ統計モデリングによるデータ分析入門』, 講談社

ベイジアン線形回帰モデルの式変形とRでのギブスサンプリングの適用

今回は特に目新しい手法というわけでもなく、線形回帰モデルのギブスサンプリングについて忘備録として残しておきたいと思います。
ベイジアン線形回帰モデルはプログラミング言語で言う、Hello World!的なものなので、あえてブログで取り上げる必要があると考えていないのですが、導出をしては忘れの繰り返しが嫌なので自分のために残しておこうと考えました。加えて、Stanのありがたみを感じられ、Stanへのコミットメントが増すのではないかとも期待しています。

・モデル
・数式の展開
・Rのコードの紹介
・おわりに
・参考情報

モデル

東北大学の照井教授の『ベイズモデリングによるマーケティング分析』に載せられている表記に従い、以下のように記します。

説明変数の数がk個の正規線形モデル

を考える。その場合、尤度関数は

となる。

係数パラメータの事前分布や条件付きの誤差分散の事前分布は以下のように設定する。（βは正規分布に従い、σ2|βは逆ガンマ分布に従う。）

数式の展開

私が大学院生だった時に、数式の展開をどう進めるかを手っ取り早く知る方法としては、「ネットに上がっている海外の大学院の講義資料を漁る」という作戦を取っていました。こうすることで数学のセンスがそれほど高くなくても、理解し進めることができました。今回に関してもおそらく、わかりやすく解説している海外の研究者がいるはずだと思い漁ってみたところ、コロンビア大学の機械学習の講義資料を見つけることができました。

資料はこちらのPDF（Course Notes for Bayesian Models for Machine Learning）で126ページにもなっていますが、導出のステップなどが非常に丁寧に書かれています。

それでは、今回の講義ノートを参考にしながら、線形モデルにおいて、ギブスサンプリングを行うところまでの導出を行いたいと思います。

まず、同時事後分布を以下の左辺のように置き、ベイズの定理を用いて右辺のように表記する。

次に、条件付き確率の定義と先程の尤度関数から以下のようになる。

yが与えられたもとでのp(y)は一定のため、比例している分子だけを残すと以下のようになる。

同時事後分布に事前分布の関数を代入していくと、

となる。両辺について対数を取ると、

となる。ここでβやσ2についての事前分布の形状から、同時事後分布におけるβやσ2について整理するための目標となる形状を確かめる。
まず、βはp(β)の定義より、対数を取りβについて整理すると、

となる。つまり、1/B0や1/B0・β0に該当する表現を先程の対数を取った同時事後分布から得ることを目標とする。
他方、σ2についても同様に、p(β|σ2)の定義より対数を取りσ2について整理すると、

となる。つまり、ν0やδ0に該当する表現を、同じく対数を取った同時事後分布から得ることを目標とする。

以上のパラメータごとの目標とする形状になるように各々のパラメータについて、対数を取った同時事後分布を整理する。

まずはβについてまとめ、関係のない項をconst.にする。

先程もとめた目標の形状を当てはめると以下のようになる。

よって、βの事後分布は以下のようになる。

他方、σ2についても同様に、関係のない項をconst.にし、目標の形状にまとめると以下のようになる。

目標の形状と比較すると以下のようになる。

よって、σ2の事後分布は以下のようになる。

Rのコードの紹介

条件付き事後分布からβやσ2の従う分布の形状がわかったので、それらを使ってRでギブスサンプリングを行います。先日、たまたま見つけた線形回帰モデルのギブスサンプリングのRのソースコードを拝借しようと思います。

ギブスサンプリングでは、先程導出した条件付き分布からβ→σ2と交互にサンプリングしていきます。それを記述したRコードは以下の通りです。

# function for blocked gibbs sampler
block_gibbs <- function(y, x, iter, burnin, trim){
  # initialize gibbs
  xprimex_inv <- solve(t(x)%*%x) # calculate once for repeated use in sampler
  s <- numeric(iter) # shell for phi
  b <- matrix(nrow=iter, ncol = 4) # shell for betas
  s[1] <- 6 # initial phi value to start sampler
  
  # phi hyperparameters
  a <- .5
  g <- 10000
  
  # gibbs sampling
  for(i in 2:iter ){
    b[i,] <- rmvnorm(n = 1, 
                   mean = ((xprimex_inv%*%t(x))%*%y), 
                   sigma = s[i-1]*xprimex_inv )
    
    s[i] <- rinvgamma(n = 1, 
                      shape = (n/2 + a), 
                      rate = .5*( t((y - x%*%t(t(b[i,])) ))%*%(y - x%*%t(t(b[i,])) ) ) + g)
  }
  
  # apply burnin and trimming  
  keep_draws <- seq(burnin,iter,trim)
  s <- s[keep_draws]
  b <- b[keep_draws,]
  
  # format and output
  joint_post<-data.frame(b=b,s=s)
  colnames(joint_post)[1:(ncol(x))]<-paste0('B',0:(ncol(x)-1) )
  
  joint_post_long<-gather(joint_post,keep_draws) %>%
    rename(param=keep_draws, draw=value) %>%
    mutate(iter=rep(keep_draws,ncol(joint_post)))
  
  return(joint_post_long)
}

# function for blocked gibbs sampler

block_gibbs <- function(y, x, iter, burnin, trim){

# initialize gibbs

xprimex_inv <- solve(t(x)%*%x) # calculate once for repeated use in sampler

s <- numeric(iter) # shell for phi

b <- matrix(nrow=iter, ncol = 4) # shell for betas

s[1] <- 6 # initial phi value to start sampler

# phi hyperparameters

a <- .5

g <- 10000

# gibbs sampling

for(i in 2:iter ){

b[i,] <- rmvnorm(n = 1,

mean = ((xprimex_inv%*%t(x))%*%y),

sigma = s[i-1]*xprimex_inv )

s[i] <- rinvgamma(n = 1,

shape = (n/2 + a),

rate = .5*( t((y - x%*%t(t(b[i,])) ))%*%(y - x%*%t(t(b[i,])) ) ) + g)

}

# apply burnin and trimming

keep_draws <- seq(burnin,iter,trim)

s <- s[keep_draws]

b <- b[keep_draws,]

# format and output

joint_post<-data.frame(b=b,s=s)

colnames(joint_post)[1:(ncol(x))]<-paste0('B',0:(ncol(x)-1) )

joint_post_long<-gather(joint_post,keep_draws) %>%

rename(param=keep_draws, draw=value) %>%

mutate(iter=rep(keep_draws,ncol(joint_post)))

return(joint_post_long)

}

先程導出したβの事後分布である正規分布からのサンプリングの後（15~17行目）、そのサンプリングしたβを用いて、同じく先程導出したσ2の事後分布である逆ガンマ分布からサンプリングし（19~21行目）、それを指定した回数だけ繰り返し、所定の数まではバーンインとして除外します。（25~27行目）こうして導出した数式と、ギブスサンプリングのコードを見比べると理解が捗ると思いました。

実際に、先程のGitHubのソースコードを回してみると、以下のようにギブスサンプリングのイタレーションのプロットや、パラメータの事後分布を確認できます。

全体のコードはこちらです。

################################################################################
###### 0 - Packages and Simulate Data
################################################################################
library(mvtnorm)
library(invgamma)
library(ggplot2)
library(dplyr)
library(tidyr)
library(xtable)
set.seed(200)

n <- 50 # number of observation
# simulate model matrix
x <- cbind(1, rnorm(n, 0, 1), rnorm(n, 5,10),rnorm(n, 100,10))

# true beta coefficients
tb <- c(1000, 50, -50, 10)

# true phi
ts <- 10000
I <- diag(1,n,n) # identity matrix used for covariance matrix

# simulate outcome for regression 
y <- t(rmvnorm(1, x%*%tb, ts*I))

# simulate many outcomes...used later for asymptotic evaluations
y_list <- replicate(1000, t(rmvnorm(1, x%*%tb, ts*I)),simplify = FALSE)

################################################################################
###### 1 - Run Blocked Gibbs Sampler
################################################################################

# function for blocked gibbs sampler
block_gibbs <- function(y, x, iter, burnin, trim){
  # initialize gibbs
  xprimex_inv <- solve(t(x)%*%x) # calculate once for repeated use in sampler
  s <- numeric(iter) # shell for phi
  b <- matrix(nrow=iter, ncol = 4) # shell for betas
  s[1] <- 6 # initial phi value to start sampler
  
  # phi hyperparameters
  a <- .5
  g <- 10000
  
  # gibbs sampling
  for(i in 2:iter ){
    b[i,] <- rmvnorm(n = 1, 
                   mean = ((xprimex_inv%*%t(x))%*%y), 
                   sigma = s[i-1]*xprimex_inv )
    
    s[i] <- rinvgamma(n = 1, 
                      shape = (n/2 + a), 
                      rate = .5*( t((y - x%*%t(t(b[i,])) ))%*%(y - x%*%t(t(b[i,])) ) ) + g)
  }
  
  # apply burnin and trimming  
  keep_draws <- seq(burnin,iter,trim)
  s <- s[keep_draws]
  b <- b[keep_draws,]
  
  # format and output
  joint_post<-data.frame(b=b,s=s)
  colnames(joint_post)[1:(ncol(x))]<-paste0('B',0:(ncol(x)-1) )
  
  joint_post_long<-gather(joint_post,keep_draws) %>%
    rename(param=keep_draws, draw=value) %>%
    mutate(iter=rep(keep_draws,ncol(joint_post)))
  
  return(joint_post_long)
}

# run gibbs sampler with specified parameters
post_dist <- block_gibbs(y = y, x = x, iter = 500000, burnin = 100000, trim = 50)

################################################################################
###### 2 - Summarize and Visualize Posterior Distributions 
################################################################################

# calculate posterior summary statistics (stats not used in rest of code)
post_sum_stats<-post_dist %>%
  group_by(param) %>%
  summarise(median=median(draw),
            lwr=quantile(draw,.025),
            upr=quantile(draw,.975)) %>%
  mutate(true_vals=c(tb,ts))

# merge on summary statistics
post_dist <- post_dist %>%
  left_join(post_sum_stats, by='param')

# plot MCMC Chains
ggplot(post_dist,aes(x=iter,y=draw)) +
  geom_line() +
  geom_hline(aes(yintercept=true_vals, col='red'), show.legend=FALSE)+
  facet_grid(param ~ .,scale='free_y',switch = 'y') +
  theme_bw() + 
  xlab('Gibbs Sample Iteration') + ylab('MCMC Chains') + 
  ggtitle('Gibbs Sampler MCMC Chains by Parameter')

# plot Posterior Distributions
ggplot(post_dist,aes(x=draw)) +
  geom_histogram(aes(x=draw),bins=50) +
  geom_vline(aes(xintercept = true_vals,col='red'), show.legend = FALSE) +
  facet_grid(. ~ param, scale='free_x',switch = 'y') +
  theme_bw() + 
  xlab('Posterior Distributions') + ylab('Count') + 
  ggtitle('Posterior Distributions of Parameters (true values in red)')

100

101

102

103

104

105

106

107

################################################################################

###### 0 - Packages and Simulate Data

################################################################################

library(mvtnorm)

library(invgamma)

library(ggplot2)

library(dplyr)

library(tidyr)

library(xtable)

set.seed(200)

n <- 50 # number of observation

# simulate model matrix

x <- cbind(1, rnorm(n, 0, 1), rnorm(n, 5,10),rnorm(n, 100,10))

# true beta coefficients

tb <- c(1000, 50, -50, 10)

# true phi

ts <- 10000

I <- diag(1,n,n) # identity matrix used for covariance matrix

# simulate outcome for regression

y <- t(rmvnorm(1, x%*%tb, ts*I))

# simulate many outcomes...used later for asymptotic evaluations

y_list <- replicate(1000, t(rmvnorm(1, x%*%tb, ts*I)),simplify = FALSE)

################################################################################

###### 1 - Run Blocked Gibbs Sampler

################################################################################

# function for blocked gibbs sampler

block_gibbs <- function(y, x, iter, burnin, trim){

# initialize gibbs

xprimex_inv <- solve(t(x)%*%x) # calculate once for repeated use in sampler

s <- numeric(iter) # shell for phi

b <- matrix(nrow=iter, ncol = 4) # shell for betas

s[1] <- 6 # initial phi value to start sampler

# phi hyperparameters

a <- .5

g <- 10000

# gibbs sampling

for(i in 2:iter ){

b[i,] <- rmvnorm(n = 1,

mean = ((xprimex_inv%*%t(x))%*%y),

sigma = s[i-1]*xprimex_inv )

s[i] <- rinvgamma(n = 1,

shape = (n/2 + a),

rate = .5*( t((y - x%*%t(t(b[i,])) ))%*%(y - x%*%t(t(b[i,])) ) ) + g)

}

# apply burnin and trimming

keep_draws <- seq(burnin,iter,trim)

s <- s[keep_draws]

b <- b[keep_draws,]

# format and output

joint_post<-data.frame(b=b,s=s)

colnames(joint_post)[1:(ncol(x))]<-paste0('B',0:(ncol(x)-1) )

joint_post_long<-gather(joint_post,keep_draws) %>%

rename(param=keep_draws, draw=value) %>%

mutate(iter=rep(keep_draws,ncol(joint_post)))

return(joint_post_long)

}

# run gibbs sampler with specified parameters

post_dist <- block_gibbs(y = y, x = x, iter = 500000, burnin = 100000, trim = 50)

################################################################################

###### 2 - Summarize and Visualize Posterior Distributions

################################################################################

# calculate posterior summary statistics (stats not used in rest of code)

post_sum_stats<-post_dist %>%

group_by(param) %>%

summarise(median=median(draw),

lwr=quantile(draw,.025),

upr=quantile(draw,.975)) %>%

mutate(true_vals=c(tb,ts))

# merge on summary statistics

post_dist <- post_dist %>%

left_join(post_sum_stats, by='param')

# plot MCMC Chains

ggplot(post_dist,aes(x=iter,y=draw)) +

geom_line() +

geom_hline(aes(yintercept=true_vals, col='red'), show.legend=FALSE)+

facet_grid(param ~ .,scale='free_y',switch = 'y') +

theme_bw() +

xlab('Gibbs Sample Iteration') + ylab('MCMC Chains') +

ggtitle('Gibbs Sampler MCMC Chains by Parameter')

# plot Posterior Distributions

ggplot(post_dist,aes(x=draw)) +

geom_histogram(aes(x=draw),bins=50) +

geom_vline(aes(xintercept = true_vals,col='red'), show.legend = FALSE) +

facet_grid(. ~ param, scale='free_x',switch = 'y') +

theme_bw() +

xlab('Posterior Distributions') + ylab('Count') +

ggtitle('Posterior Distributions of Parameters (true values in red)')

おわりに

シンプルなモデルですらこれだけ導出に手間がかかるということからも、Stanなどの確率的プログラミング言語のありがたみは非常に大きいなと思いました。こうして残すことで今後忘れたとしてもすぐに思い出せる気がします。
しかしながら、Stanでは事前分布と尤度を指定してしまえば、事後分布を計算し、知りたいパラメータについて解いた条件付き分布からサンプリングしてくれるわけですから、研究者の寿命を伸ばしたと言っても過言ではないかもしれません。

参考情報

[1]John Paisley (2016), “Course Notes for Bayesian Models for Machine Learning”, Columbia University
[2]照井伸彦 (2008), 『ベイズモデリングによるマーケティング分析』, 東京電機大学出版局
[3]須山敦志 (2017), 『機械学習スタートアップシリーズベイズ推論による機械学習入門』, 講談社
[4]stablemarkets,BayesianTutorials/MultipleLinearReg/multiplelinearreg.r

[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］項目反応理論の応用でフリースタイルダンジョン登場ラッパーの強さをランキングしてみた」

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
に貼り付けて画像を出力しています。

目次

はじめに

LDATSパッケージについて

データ

試してみた

おわりに

参考情報

はじめに

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の実行

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

テストデータでの予測

非階層モデルとの比較

おわりに

参考情報

モデル

数式の展開

Rのコードの紹介

おわりに

参考情報

はじめに

生存時間分析とは

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

先行研究

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

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入門サービスの全体像からクラウドネイティブアプリケーション構築まで』（翔泳社）