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/

ディリクレNBDモデルのマーケティング分野での適用に関して色々調べてみた

どの企業のマーケティング担当者も、自分が属する市場や自分の会社の商品の需要予測ないし購買行動の予測に関心があると思います。それらを予測できる手法として、ディリクレNBDモデルが候補にあがると思います。
以前、『確率思考の戦略論 USJでも実証された数学マーケティングの力』を読んだ際に、需要のシミュレーションにディリクレNBDモデルが使えるということを知りましたが、繰り返し利用するサービスや消費財向けのアプローチなのかなと思って、自分が扱うサービスへの適用可能性の低さから距離を置いてきました。
ところが、文献によると様々な分野でディリクレNBDモデルを使って生命保険やクレジットカードなどのサービス市場需要の予測がうまくいっているケースがあるらしく、真剣に向き合ってみようと思い調べてみることにしました。

歴史

マーケティングの分野では、消費者の繰り返し購買行動やダブルジョパディの法則など様々な現象を説明するための研究と実践がなされてきました。
そのために確率的なモデルが色々と適用されていきました。

“The Dirichlet model in marketing”という論文の図を引用しますが、1952年ごろからブランドロイヤルティーの研究が始まり、1958年には線形モデルが用いられ、1959年にはNBD(Negative Binomial Distribution：負の二項分布)を用いた購買予測などがなされるようになりました。
それらのアプローチを様々な対象に適用しつつ、1975年にはBBD(Beta Binomial Distribution：ベータ二項分布)が、ついに1984年にはNBD-DirichletモデルがGoodhardtらによって考案されるに至りました。

理論の概要

本家の論文（The Dirichlet: A Comprehensive Model of Buying Behaviour）や“Calculation of Theoretical Brand Performance Measures from the Parameters of the Dirichlet Model”を読まれることをお勧めしますが、この理論が想定しているものは以下の通りです。

・各々の消費者がポワソン過程に従って、負の二項分布に従い購買行動をする。
・消費者ごとの購買率はガンマ分布に従う。
・消費者の利用可能なブランドの選択は多項分布に従う。
・ブランドの選択確率自体は様々な消費者ごとに、多変量ベータ分布、あるいはディリクレ分布に従う。

なお、ディリクレモデルは負の二項分布とディリクレ多変量分布の二つの確率密度関数の結合からなり、二つの分布はそれぞれ独立であることを想定しています。
マーケットにおける購買自体の行動をつかさどる分布（負の二項分布）と、何を選ぶかをつかさどる分布（ディリクレ多変量分布）が独立しているということになります。

まず、負の二項分布についてですが、確率密度関数は以下の通りとなります。
$$ f_{\gamma, \beta}(k) = \frac{\Gamma (\gamma + k)}{\Gamma(\gamma)k!}\frac{\beta^k}{(1+\beta)^{\gamma + k}} \\ \textrm{ for } k = 0, 1, 2,…$$

負の二項分布は二つのパラメータを持ちます。
・形状（shape）パラメータ（$\gamma$）
・尺度（scale）パラメータ（$\beta$）

続いて、ディリクレ多変量分布ですが、確率密度関数は以下の通りとなります。

$$ f_{\alpha_1, \alpha_2, \dots, \alpha_h}(r_1, r_2, \dots, r_h | r_1 + r_2 + \dots + r_h = k ) \\
\frac{\Gamma \left ( \sum_{j=1}^{h} \alpha_j \right ) k! }{ \Gamma \left ( \sum_{j=1}^{h} \alpha_j + k \right ) } \prod_{j=1}^{h} \frac{\Gamma \left ( \alpha_j + r_j \right ) }{r_j! \Gamma \left ( \alpha_j \right )}$$

$r_h$などはそのマーケットにおける特定ブランドの購買率で、ブランド数がhだけあるとしています。
なお、
$$ r_1 + r_2 + \dots + r_h = k $$
で、足し合わせるとkがマーケット自体の購買率になるようになっています。ディリクレ多変量分布はブランドの数だけ、h個の正のパラメータ$ \alpha_1,\alpha_2, \dots, \alpha_h $を持ちます。

ディリクレモデルは、以上の負の二項分布とディリクレ多変量分布を結合したもので、確率密度関数は以下の通りになります。

$$ f_{\gamma, \beta, \alpha_1, \alpha_2, \dots, \alpha_h}(r_1, r_2, \dots, r_h ) = \\ f_{\gamma, \beta}(k) f_{\alpha_1, \alpha_2, \dots, \alpha_h}(r_1, r_2, \dots, r_h | r_1 + r_2 + \dots + r_h = k )$$

この確率密度関数に色々かけて期待値を計算することで以下のようなものを算出することができます。
・あるブランドの購買者ごとの、あるブランドの購買の数の理論値
・プロダクトクラスでの購買の数の平均値
・あるブランドのみを買う母集団の割合
・購買者の平均購買頻度

諸々の式の導出は「[確率思考の戦略論] 1.確率理論の導入とプレファレンスの数学的説明」に詳しく書かれているので、見てみることをお勧めします。
数理統計学の知識とか、ベータ関数・ガンマ関数などの知識があれば式変形などに関してスムーズに理解できると思います。

ただ、この手法自体への批判などもないわけではなく、先ほどあげた論文、“The Dirichlet model in marketing”では、消費者行動を理解するのにつながるモデルではないとか、意思決定を扱っているモデルではないとする批判などがあるようです。ただ、幅広いマーケットで予測されてきた実績もあることから、モデルの限界を知りつつ適用するマーケットを増やしていくことが実務家に求められているのかなと思われます。

適用例

ディリクレモデルが様々な市場で適用されている事例が、『ディリクレモデルの境界条件 ― サービスへの適用可能性と限界 ―』という論文に記されていました。

当初は、インスタントコーヒー市場、衣料用洗剤市場、シャンプー市場などの最寄品（計画的に購入されることが少なく、単価は低く、何度も繰り返し購入される製品・サービス）のみで近似できるとされていましたが、
近年では銀行市場、クレジットカード市場、ガソリンスタンド市場、スーパーマーケット市場、テレビ番組市場、ファーストフード市場、オーケストラ市場、フットボールリーグ市場、野球市場などのサービス業や、スポーツウェアや自動車などの比較的高額な市場においても近似できるとされています。

データについて

ディリクレモデルの理論値を計算するには、4つのデータが必要とされています。
・1.そのカテゴリを購買した人々の全体に占める割合
・2.カテゴリに占める、いずれかの製品を購買した人々に対して記録された当該製品カテゴリの購買回数の平均
・3.各ブランドを一度でも購入した人々の割合
・4.各ブランドを購買した人々による各ブランドの購買回数の平均値

これらのデータはアンケートなどで収集する必要があるため、ディリクレモデルは総じて時間もお金もかかる手法であると言うこともできると思います。先ほどあげた、『ディリクレモデルの境界条件 ― サービスへの適用可能性と限界 ―』にはWEB調査でどんなアンケートを集めるのかの詳細まで書かれているので、自社の属するマーケットでこのモデルを適用する際は、それらを参考にすると良いと思います。

実践するコード

コードに関しては、こちらのColab（2_nb.ipynb）を書かれている方がおられたので、それを試すのも良いと思います。

ライブラリに関して、2016年と決して新しくないですが、R言語ではNBDdirichletというライブラリがあります。手法が1984年のものですから別に問題はないと思います。

# install.packages("NBDdirichlet")

library(NBDdirichlet)

cat.pen <- 0.56 # 市場浸透度
cat.buyrate <- 2.6 # 購買頻度
brand.share <- c(0.25, 0.19, 0.1, 0.1, 0.09, 0.08, 0.03, 0.02) # ブランドのマーケットシェア
brand.pen.obs <- c(0.2,0.17,0.09,0.08,0.08,0.07,0.03,0.02) # ブランド浸透度
brand.name <- c("Colgate DC", "Macleans","Close Up","Signal","ultrabrite",
                "Gibbs SR","Boots Priv. Label","Sainsbury Priv. Lab.") # ブランド名

dobj <- dirichlet(cat.pen, cat.buyrate, brand.share, brand.pen.obs, brand.name)
print(dobj)
summary(dobj)

# install.packages("NBDdirichlet")

library(NBDdirichlet)

cat.pen <- 0.56 # 市場浸透度

cat.buyrate <- 2.6 # 購買頻度

brand.share <- c(0.25, 0.19, 0.1, 0.1, 0.09, 0.08, 0.03, 0.02) # ブランドのマーケットシェア

brand.pen.obs <- c(0.2,0.17,0.09,0.08,0.08,0.07,0.03,0.02) # ブランド浸透度

brand.name <- c("Colgate DC", "Macleans","Close Up","Signal","ultrabrite",

"Gibbs SR","Boots Priv. Label","Sainsbury Priv. Lab.") # ブランド名

dobj <- dirichlet(cat.pen, cat.buyrate, brand.share, brand.pen.obs, brand.name)

print(dobj)

summary(dobj)

こちらを実行した結果は、以下の通りです。

> print(dobj)
Number of Brands in the Category = 8 
Brand List : Colgate DC : Macleans : Close Up : Signal : ultrabrite : Gibbs SR : Boots Priv. Label : Sainsbury Priv. Lab.
Brands' Market Shares: 0.25 0.19 0.1 0.1 0.09 0.08 0.03 0.02 
Brands' Penetration:   0.2 0.17 0.09 0.08 0.08 0.07 0.03 0.02 
Multiple of Base Time Period: 1 , Current M = 1.456 

Category Penetration = 0.56 , with Buying Rate = 2.6 
Estimated Dirichlet Model Parameters:
NBD: M = 1.46 ,  K = 0.78 ;  Dirichlet: S = 1.3

> print(dobj)

Number of Brands in the Category = 8

Brand List : Colgate DC : Macleans : Close Up : Signal : ultrabrite : Gibbs SR : Boots Priv. Label : Sainsbury Priv. Lab.

Brands' Market Shares: 0.25 0.19 0.1 0.1 0.09 0.08 0.03 0.02

Brands' Penetration: 0.2 0.17 0.09 0.08 0.08 0.07 0.03 0.02

Multiple of Base Time Period: 1 , Current M = 1.456

Category Penetration = 0.56 , with Buying Rate = 2.6

Estimated Dirichlet Model Parameters:

NBD: M = 1.46 , K = 0.78 ; Dirichlet: S = 1.3

このdirichlet関数を使うことで、データを与えるだけで、ディリクレモデルのパラメータMやKやSが簡単に求まりました。

summaryで推定した諸々の値を確認できます。引数で色々と調整するので、実際に使う際はリファレンス（Package ‘NBDdirichlet’）を見てください。頻度のカットオフ値やヘビーユーザーの域値や、ブランド重複を見る際の基準ブランドなどを指定できます。

・buy：理論的なブランド浸透度、購買率、そのブランドの買い手によるマーケット内での購買率
・freq：そのブランドの購買の数の分布（以下の場合、6以上はまとめられている）
・heavy：特定の頻度の買い手におけるマーケットでの理論的なブランド浸透度、購買率
・dup：特定のブランドの購買者が他のブランドも購買している割合

> summary(dobj)
$buy
                     pen.brand pur.brand pur.cat
Colgate DC                0.20      1.82    3.16
Macleans                  0.16      1.76    3.22
Close Up                  0.09      1.68    3.30
Signal                    0.09      1.68    3.30
ultrabrite                0.08      1.67    3.31
Gibbs SR                  0.07      1.66    3.32
Boots Priv. Label         0.03      1.62    3.37
Sainsbury Priv. Lab.      0.02      1.61    3.38

$freq
                        0    1    2    3    4 5   6+
Colgate DC           0.80 0.12 0.04 0.02 0.01 0 0.01
Macleans             0.84 0.10 0.03 0.01 0.01 0 0.00
Close Up             0.91 0.06 0.02 0.01 0.00 0 0.00
Signal               0.91 0.06 0.02 0.01 0.00 0 0.00
ultrabrite           0.92 0.05 0.02 0.01 0.00 0 0.00
Gibbs SR             0.93 0.05 0.01 0.01 0.00 0 0.00
Boots Priv. Label    0.97 0.02 0.01 0.00 0.00 0 0.00
Sainsbury Priv. Lab. 0.98 0.01 0.00 0.00 0.00 0 0.00

$heavy
                     Penetration Avg Purchase Freq
Colgate DC                  0.34              1.61
Macleans                    0.27              1.57
Close Up                    0.15              1.51
Signal                      0.15              1.51
ultrabrite                  0.13              1.50
Gibbs SR                    0.12              1.49
Boots Priv. Label           0.05              1.46
Sainsbury Priv. Lab.        0.03              1.45

$dup
          Colgate DC             Macleans             Close Up               Signal 
                1.00                 0.19                 0.10                 0.10 
          ultrabrite             Gibbs SR    Boots Priv. Label Sainsbury Priv. Lab. 
                0.09                 0.08                 0.03                 0.02

> summary(dobj)

$buy

pen.brand pur.brand pur.cat

Colgate DC 0.20 1.82 3.16

Macleans 0.16 1.76 3.22

Close Up 0.09 1.68 3.30

Signal 0.09 1.68 3.30

ultrabrite 0.08 1.67 3.31

Gibbs SR 0.07 1.66 3.32

Boots Priv. Label 0.03 1.62 3.37

Sainsbury Priv. Lab. 0.02 1.61 3.38

$freq

0 1 2 3 4 5 6+

Colgate DC 0.80 0.12 0.04 0.02 0.01 0 0.01

Macleans 0.84 0.10 0.03 0.01 0.01 0 0.00

Close Up 0.91 0.06 0.02 0.01 0.00 0 0.00

Signal 0.91 0.06 0.02 0.01 0.00 0 0.00

ultrabrite 0.92 0.05 0.02 0.01 0.00 0 0.00

Gibbs SR 0.93 0.05 0.01 0.01 0.00 0 0.00

Boots Priv. Label 0.97 0.02 0.01 0.00 0.00 0 0.00

Sainsbury Priv. Lab. 0.98 0.01 0.00 0.00 0.00 0 0.00

$heavy

Penetration Avg Purchase Freq

Colgate DC 0.34 1.61

Macleans 0.27 1.57

Close Up 0.15 1.51

Signal 0.15 1.51

ultrabrite 0.13 1.50

Gibbs SR 0.12 1.49

Boots Priv. Label 0.05 1.46

Sainsbury Priv. Lab. 0.03 1.45

$dup

Colgate DC Macleans Close Up Signal

1.00 0.19 0.10 0.10

ultrabrite Gibbs SR Boots Priv. Label Sainsbury Priv. Lab.

0.09 0.08 0.03 0.02

この推定結果を使って、売上がどれくらいになるかが高精度にわかるのであれば、ブランドの浸透度が上がるような施策を打つことで、あるいは競合がそれらを打つことで売上がどう変わっていくかをシミュレーションすることもできるかもしれません。
その背景にはダブルジョパディの法則などがあると思うと結構面白そうですね。

参考文献

・The Dirichlet: A Comprehensive Model of Buying Behaviour
・Calculation of Theoretical Brand Performance Measures from the Parameters of the Dirichlet Model
・The Dirichlet model in marketing
・NBDdirichlet-package: NBD-Dirichlet model of consumer buying behavior
・Package ‘NBDdirichlet’
・ディリクレモデルの境界条件 ― サービスへの適用可能性と限界 ―
・[確率思考の戦略論] 1.確率理論の導入とプレファレンスの数学的説明
・Negative Binomial distribution
・Dirichlet distribution

[PRの巻]RユーザのためのRStudio[実践]入門（改訂2版）

今回はご献本PRの巻です。お題はこちらの本、『改訂2版 RユーザのためのRStudio[実践]入門〜tidyverseによるモダンな分析フローの世界』です。

リアルでもバーチャルでも大変お世話になっているy__mattuさんより、ありがたく献本いただけました。どうやら、今回の改訂では56ページほどページ数が増えているようです。

目次やら詳細は技術評論者のページをご覧ください。

端的に言うと、データ分析でR言語を使ったことがあるユーザーで、Rの比較的新しいライブラリの使い方を知りたい方、RStudioを使って生産性高く分析をしたい方にピッタリの本です。データ収集・前処理・可視化・レポーティングを学ぶことができます。

元よりこの本のユーザーであったことから、これまで色々なマーケター・エンジニアにこの本を勧めてきました。

社内にデータ分析に関心のある人が一定数いて、そのような人にR言語を勧める際にこの本はちょうど良いと思っています。
この本は、統計学や機械学習などについて触れてはいないですが、データ分析のほとんどの工程は前処理や可視化ですし、それを行えずして統計学も機械学習も成し遂げることはできないので、現実を知ってもらう上でも良い本です。
この本で自由自在に前処理・EDAができるようになったら、統計学や機械学習について書かれたR言語の本をやるように伝えています。

ここでは、初めてこの本を手に取る方、1版をすでにお持ちの方のそれぞれのお気持ちになって見どころを述べていきます。

初めてこの本を手に取る方
- 見どころ1
  RStudioの画面の説明やショートカットの説明などが丁寧で、最初の段階で生産性の高いコーディングのための準備が捗る。
- 見どころ2
  Webスクレイピングにしっかりと章が割かれている。
  データ分析を学ぶに際して、多くのケースで手元にデータがなかったり、関心のあるデータが表形式で提供されていなかったりします。そのため、Webスクレイピングをしっかりと学べるのはデータ分析を自分でやっていくために不可欠な技術だろうと思います。
- 見どころ3
1版をすでにお持ちの方
- 見どころ1
  stringrやlubridateなどのライブラリの使い方についてAppendixで40ページ以上が割かれています。RやPythonを行ったり来たりする生活をしていると忘れがちなので、充実のAppendixは嬉しいと思います。stringrは何が出来たっけ？てなりますし、lubridateも非常にしばしば忘れる気がします。
- 見どころ2
  tidyrの記述が増えていました。Pythonに浮気をしていると、Rのモダンな書き方に付いていけなくなるかもしれないので、こういった本でキャッチアップできるのは良いことに思います。

1版を会社のデスクに置いたまんまにしているので、今度しっかり比較して追記しようと思います。

[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’

RStudioをdockerで使える、Rockerを触ってみた

はじめに

dockerは会社のエンジニア向けの勉強会で紹介されていて、存在は知ってはいたが使っていませんでした。先日参加したTokyo.RでRockerというものが紹介されており、よし使ってみようと思うに至りました。
今回は遅ればせながら、Rockerを使ってみて、何かを回して、保存するという一連のチュートリアルをやってみます。

用語のざっくり理解

専門家の方に怒られたらあれですが、私の理解はこんな感じです。

Dockerfile
- OS、言語、ライブラリ、バージョンなどをレシピの食材っぽい感じでテキスト形式で書かれたもの
イメージ
- コンテナの元。Dockerfileをもとにdockerで生成できる。あるいはdocker hubとか言うのサイトから入手できる。ここからコンテナ（RStudioとかを回したい環境のこと）を起動できる。
コンテナ
- RStudioとかを起動させておく実行環境。
永続化
- ファイルの内容が消えないように残すこと。立ち上げたコンテナを永続化しないと都度データやインストールしたパッケージが消える。

Rockerはdocker hubで提供されているイメージです。RockerのGitHubにはDockerfileも当然あげられていました。自分であれやこれや設定を追加するにはこれを参考にいじることになるのだと思います。

とにかくやってみる

まずは、dockerをインストールします。

Install Docker Desktop on Mac
docker自体はだいぶ前にインストールしたので、入れ方を忘れていますが、ここを見ればわかるはず。

ターミナルで以下を実行

docker run -e PASSWORD=perio --rm -p 8787:8787 rocker/tidyverse

1	docker run -e PASSWORD=perio --rm -p 8787:8787 rocker/tidyverse

ブラウザなどで以下を表示

http://localhost:8787/

1	http://localhost:8787/

DALEXのサンプルコードを回してみる。DALEXはTokyo.Rでも紹介されている機械学習の解釈可能性に関する手法をあらかた実践できるパッケージです。
moDel Agnostic Language for Exploration and eXplanation

以下の図はDALEXなどを使ってGLMで推定した結果の各特徴量の影響度をプロットしたものです。

ターミナルでコンテナのIDを確認する。

docker ps

docker ps

分析結果やソースコードやデータを保存したいので、コンテナの永続化を行います。

docker commit f9317be953b5 perio/rocker_hadley

1	docker commit f9317be953b5 perio/rocker_hadley

ここではユーザー名をperio、イメージ名をrocker_hadleyにしています。

永続化したコンテナのイメージがきちんと存在するか確認してみます。

$ docker images

REPOSITORY                     TAG                 IMAGE ID            CREATED             SIZE
perio/rocker_hadley            latest              f3403d6a9d85        4 hours ago         2.21GB

$ docker images

REPOSITORY TAG IMAGE ID CREATED SIZE

perio/rocker_hadley latest f3403d6a9d85 4 hours ago 2.21GB

ありますね。

さて、先程のコンテナを止めてみましょう。コンテナIDを指定するだけです。

docker stop f9317be953b5

1	docker stop f9317be953b5

永続化したコンテナを読み込んでRstudioを起動してみます。

docker run -e PASSWORD=perio --rm -p 8787:8787 perio/rocker_hadley

1	docker run -e PASSWORD=perio --rm -p 8787:8787 perio/rocker_hadley

実は、誤って最新の状態で保存せずに再度立ち上げてしまいましたｗDALEXのサンプルコードを回すために色々パッケージインストールしたので、それが消えてしまった。残念。
元のimageにはないRDataがあるので、一応、永続化はできているので安心です。うむ、次からはきちんと保存をしよう。

これでどのPCでも同じ条件で分析していけるので、良さそうです。
~~ただ、一部で入らないパッケージ（sfパッケージとか）があったので、より良いDockerfileを探し求めたいとも思いますね。~~

sfを含む空間解析系のパッケージを扱うものとして、rocker/rspatialがあります。
記事中に登場するtidyverseイメージをベースにしているので同様に使えます。tidyverseイメージで入れるよりは大分楽になるかと。 https://t.co/YBUnZEjEKV

— Uryu Shinya (@u_ribo) January 26, 2020

u_riboさんがsfパッケージが動くdocker imageを教えてくれました。これまで、一つの環境にギュッとパッケージを詰め込んでやってきたんですが、目的に応じてdocker imageを選ぶのが良いのでしょう。

docker run -e PASSWORD=perio --rm -p 8787:8787 rocker/geospatial

1	docker run -e PASSWORD=perio --rm -p 8787:8787 rocker/geospatial

せっかくなので先日作成した事故物件分析用のコードを回してみましょう。

mapviewもインストールできて、問題なく使えました。u_riboさんありがとう。

参考情報

[1] The Rocker Project Docker Containers for the R Environment
[2] The Rocker Images: choosing a container
[3] rocker/tidyverse
[4] How to save data
[5] DALEXverse and fraud detection

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

R advent calendar 2019 RSelenium、jpmesh、sfパッケージで東京23区の事故物件を分析してみよう！

はじめに

今回で3回目となるR advent Calendarですが、前回は「一発屋芸人の検索トレンドのデータ」を扱い、前々回は「ポケモンのデータ」を扱いました。今回は人の命に関わるようなデータを扱ってみたいと思い、某サイトから東京都の23区内における事故物件の住所と詳細を集めてきました。どのようなエリアが事故が起きやすいのかの分析を行います。(以下では、事故物件をAP(Accident Property)と呼びます。)

R advent calendar用に向けたデータ収集コードができたが、データがけっこう汚いのでアノテーション工数を覚悟したほうが良さそう。
今回は人の命に関わるデータですね。

— SKUE (@Mr_Sakaue) December 8, 2019

分析工程

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

データの収集

APに関する情報を某サイトより集める必要があります。そこで必要なライブラリとしては、RSeleniumやtwitteRがあげられます。
twitteRが必要な理由は、APに関するサイトにAPの一覧ページがなく、公式アカウントがAPに関するページのリンクを投稿しているところにあります。ただ、私が以前使っていた頃とはTwitterAPIの仕様が変わり、3ヶ月よりも前の情報にアクセスできなくなっていました。そのため、今後のデータに関してはTwitterAPIでいいのですが、過去のものに関しては別アプローチが必要となります。
また、APに関するサイトはJavaScriptで地図が表示されているだけなので、RSeleniumを使って地図をクリックさせ、表示された情報をスクレイピングするなどの処理が必要となります。
当初の想定ではTwitterのデータ収集、リンク先の情報をRSeleniumでスクレイピングするだけの簡単な仕事だと思っていたのですが、過去のデータにアクセスできないので、地図上で一つ一つ見つけていくためにRSeleniumだけで頑張ることにしました。(私の過去のアドカレ史上、一番面倒なデータとなりました。)

誰もやらないと思うのですが、一応手順を記しておきます。

RSeleniumだけでMapからAPの情報を抽出するための手順
1.都内の住所一覧を収集
2.検索窓に住所を入力
3.検索結果一覧の上位5件をクリック
4.一度地図を引くことでAPを広い範囲で捉えれるようにする
5.APのマークの要素を取得し、1件ずつクリックし、表示されたAPの情報をデータフレームに格納する

こちらが取得できたデータです。

RSeleniumでAPの
・住所
・発生時期（フリーテキスト）
・AP詳細
・投稿日
を集めることができるので、その住所データに対して、Yahoo!のジオコードに関するAPIを利用します。(利用申請が必要なはずです。4年前くらいに申請していたのでそのまま使えました。)
Yahoo!のAPIを使えば、住所から緯度経度の情報を取得することができます。

APの緯度経度がわかれば、jpmeshパッケージを用いて1kmメッシュやら10kmメッシュやらのメッシュデータに変換することができます。
jpmeshを用いてメッシュデータに変換し、メッシュ単位でAPの発生件数を集計します。

データ収集用のソースコードは思いのほか長くなってしまったので、GitHubにあげておきました。
https://github.com/KamonohashiPerry/r_advent_calendar_2019

ここで再度、手順を整理しておきましょう。

Twitterに出てきたものだけを取得（直近3ヶ月〜）する場合、
run_tweet_collect.RでTweetを収集
↓
run_selenium.RでAPの情報をスクレイピング
↓
run_map_api.Rで住所から緯度経度の取得
↓
making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

直近3ヶ月以前のものを取得する場合、
run_selenium_from_map.Rで地図上から直接APの情報を取得する
↓
making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

データの整形

APの1kmメッシュデータを手に入れたら、kokudosuuchiパッケージを使って国土地理院の収集したデータをつなぎこみます。手順としては、以下のとおりです。

まずは推計人口というそのエリアの人口の予測値です。今回は2010年のものを抽出しました。こちらは1kmメッシュのデータなので、変換することなく使えて都合が良いです。

library(kokudosuuchi)

# 推計人口
estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",
                                      cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,
                                         city_code=estimated_population$Mesh3_POP_00$CITY_CODE,
                                         population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>% 
                                mutate(one_k_mesh=as.integer(one_k_mesh))

library(kokudosuuchi)

# 推計人口

estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",

cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,

city_code=estimated_population$Mesh3_POP_00$CITY_CODE,

population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>%

mutate(one_k_mesh=as.integer(one_k_mesh))

続いて、2015年の公示地価を抽出しました。

# 公示地価
published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",
                                    cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),
                                      address=as.character(published_land_price$`L01-15`$L01_019),
                                      area_code=as.integer(published_land_price$`L01-15`$L01_017))

for (i in 1:nrow(published_land_price_df)) {
  published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])
  published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])
}

# 東京23区の地価に絞り込む
published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),
                                                                    stringr::str_detect(address,"区"))
published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)
published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

for (i in 1:nrow(published_land_price_df_tokyo)) {
  published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]
}

# 公示地価

published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",

cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),

address=as.character(published_land_price$`L01-15`$L01_019),

area_code=as.integer(published_land_price$`L01-15`$L01_017))

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

published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])

published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])

}

# 東京23区の地価に絞り込む

published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),

stringr::str_detect(address,"区"))

published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)

published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

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

published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]

}

こちらはメッシュデータではないので、緯度経度の情報から1kmメッシュのデータに変換する必要があります。後で行います。

単純集計・可視化

今回のデータセットのデータ数は3919件です。本当は7000件以上はあると思われますが、マップから取ってくるという勝手上、なかなか全てを取り切ることができませんでした。

まずは、1kmメッシュごとのAP発生件数のヒストグラムです。

g <- ggplot(data = estimated_population_df_tokyo %>% filter(accident_count>0),
            aes(accident_count)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))
g <- g +geom_histogram(bins = 50)
g <- g + ggtitle("1kmメッシュごとのAP件数のヒストグラム")
g

g <- ggplot(data = estimated_population_df_tokyo %>% filter(accident_count>0),

aes(accident_count)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))

g <- g +geom_histogram(bins = 50)

g <- g + ggtitle("1kmメッシュごとのAP件数のヒストグラム")

1kmメッシュにおける人口のヒストグラムです。

公示地価のヒストグラムです。

1kmメッシュにおける人口あたりのAP件数のヒストグラムです。

分析

ここでは、色々な軸でAPのデータに向き合ってみようと思います。

APの発生件数の集計

人口が多いところがAPの発生件数が多いところだと思われますが、とりあえず確認します。

> estimated_population_df_tokyo %>% group_by(市区町村) %>% summarize(sum_accident=sum(accident_count), population=sum(population), mean_accident_percapita=mean(accident_percapita_100)) %>% arrange(desc(sum_accident))
Simple feature collection with 49 features and 4 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 139.15 ymin: 35.5 xmax: 139.9125 ymax: 35.84167
epsg (SRID):    NA
proj4string:    NA
# A tibble: 49 x 5
   市区町村 sum_accident population mean_accident_percap…                                                                     poligon
   <chr>           <dbl>      <dbl>                 <dbl>                                                                  <GEOMETRY>
 1 世田谷区          343    952572.                0.0337 POLYGON ((139.675 35.60833, 139.6875 35.60833, 139.6875 35.6, 139.675 35.6…
 2 新宿区            302    398145.                0.0795 POLYGON ((139.7 35.69167, 139.7 35.68333, 139.6875 35.68333, 139.6875 35.6…
 3 大田区            291    803602.                0.03   MULTIPOLYGON (((139.775 35.55, 139.7875 35.55, 139.7875 35.55833, 139.775 …
 4 板橋区            283    471726.                0.0570 POLYGON ((139.7125 35.75, 139.7125 35.74167, 139.7 35.74167, 139.7 35.75, …
 5 練馬区            276    575022.                0.0484 POLYGON ((139.6875 35.74167, 139.6875 35.73333, 139.675 35.73333, 139.6625…
 6 杉並区            268    549477.                0.0477 POLYGON ((139.65 35.675, 139.65 35.66667, 139.6375 35.66667, 139.6375 35.6…
 7 中野区            214    332012.                0.0619 MULTIPOLYGON (((139.675 35.69167, 139.675 35.68333, 139.6625 35.68333, 139…
 8 豊島区            203    278974.                0.0723 POLYGON ((139.725 35.725, 139.725 35.71667, 139.7125 35.71667, 139.7125 35…
 9 港区              159    188168.                0.137  MULTIPOLYGON (((139.7875 35.63333, 139.7875 35.625, 139.775 35.625, 139.77…
10 品川区            159    325433.                0.0382 POLYGON ((139.7625 35.59167, 139.7625 35.58333, 139.75 35.58333, 139.75 35…
# … with 39 more rows

> estimated_population_df_tokyo %>% group_by(市区町村) %>% summarize(sum_accident=sum(accident_count), population=sum(population), mean_accident_percapita=mean(accident_percapita_100)) %>% arrange(desc(sum_accident))

Simple feature collection with 49 features and 4 fields

geometry type: GEOMETRY

dimension: XY

bbox: xmin: 139.15 ymin: 35.5 xmax: 139.9125 ymax: 35.84167

epsg (SRID): NA

proj4string: NA

# A tibble: 49 x 5

市区町村 sum_accident population mean_accident_percap… poligon

1 世田谷区 343 952572. 0.0337 POLYGON ((139.675 35.60833, 139.6875 35.60833, 139.6875 35.6, 139.675 35.6…

2 新宿区 302 398145. 0.0795 POLYGON ((139.7 35.69167, 139.7 35.68333, 139.6875 35.68333, 139.6875 35.6…

3 大田区 291 803602. 0.03 MULTIPOLYGON (((139.775 35.55, 139.7875 35.55, 139.7875 35.55833, 139.775 …

4 板橋区 283 471726. 0.0570 POLYGON ((139.7125 35.75, 139.7125 35.74167, 139.7 35.74167, 139.7 35.75, …

5 練馬区 276 575022. 0.0484 POLYGON ((139.6875 35.74167, 139.6875 35.73333, 139.675 35.73333, 139.6625…

6 杉並区 268 549477. 0.0477 POLYGON ((139.65 35.675, 139.65 35.66667, 139.6375 35.66667, 139.6375 35.6…

7 中野区 214 332012. 0.0619 MULTIPOLYGON (((139.675 35.69167, 139.675 35.68333, 139.6625 35.68333, 139…

8 豊島区 203 278974. 0.0723 POLYGON ((139.725 35.725, 139.725 35.71667, 139.7125 35.71667, 139.7125 35…

9 港区 159 188168. 0.137 MULTIPOLYGON (((139.7875 35.63333, 139.7875 35.625, 139.775 35.625, 139.77…

10 品川区 159 325433. 0.0382 POLYGON ((139.7625 35.59167, 139.7625 35.58333, 139.75 35.58333, 139.75 35…

# … with 39 more rows

世田谷区は最も人口が多いことから、AP発生件数では一番となっています。続いて、歌舞伎町などがある新宿が来ています。しかしながら、人口に占めるAP発生件数で言うと、港区がかなり高く出ているのがわかります。

人口あたりのAP件数

ここでは、メッシュデータをsfパッケージ用のオブジェクトに変換して、1kmにおける人口あたりのAP発生割合を可視化しています。

library(jpmesh)
library(sf)
library(tidyverse)
library(kokudosuuchi)
library(jpndistrict)
library(mapview)

# load dataset
load(file = "accident_df_with_coordinate_from_map.RData")

# メッシュデータの追加
accident_df <- accident_df %>% mutate(k_mesh="")
accident_df <- accident_df %>% filter(latitude != "", !is.na(upload_date))

# 緯度、軽度の情報から任意のメッシュのデータを得る
for (i in 1:nrow(accident_df)) {
  accident_df$k_mesh[i] <- jpmesh::coords_to_mesh(
                                        as.numeric(accident_df$longitude[i]),
                                        as.numeric(accident_df$latitude[i]),
                                        mesh_size = "1km")
}

# メッシュを整数にする
accident_df <- accident_df %>% mutate(k_mesh=as.integer(k_mesh))

# 東京都に絞る
accident_df <- accident_df %>% filter(stringr::str_detect(formatted_address,"東京"))
accident_df %>% group_by(k_mesh) %>% count() %>% arrange(desc(n))

# 公示地価
published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",
                                    cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),
                                      address=as.character(published_land_price$`L01-15`$L01_019),
                                      area_code=as.integer(published_land_price$`L01-15`$L01_017))

for (i in 1:nrow(published_land_price_df)) {
  published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])
  published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])
}

# 東京の地価に絞り込む
published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),
                                                                    stringr::str_detect(address,"区"))
published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)
published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

for (i in 1:nrow(published_land_price_df_tokyo)) {
  published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]
}

# 緯度経度を1kmメッシュにする
for (i in 1:nrow(published_land_price_df)) {
  published_land_price_df$mesh[i] <- jpmesh::coords_to_mesh(
    as.numeric(published_land_price_df$longitude[i]),
    as.numeric(published_land_price_df$latitude[i]),
    mesh_size = "1km")
}

published_land_price_df <- published_land_price_df %>% mutate(mesh=as.integer(mesh))
published_land_price_df_summary <- published_land_price_df %>% group_by(mesh) %>% summarise(mean_price=mean(price))

# 推計人口
estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",
                                   cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,
                                       city_code=estimated_population$Mesh3_POP_00$CITY_CODE,
                                       population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>% 
                                mutate(one_k_mesh=as.integer(one_k_mesh))

accident_summary <- accident_df %>%
                        group_by(k_mesh) %>%
                        summarise(accident_count=n()) %>% 
                        arrange(desc(accident_count))

estimated_population_df <- estimated_population_df %>% 
                                        left_join(accident_summary,
                                                  by = c("one_k_mesh"="k_mesh"))

estimated_population_df$accident_count <- replace_na(estimated_population_df$accident_count, 0)

estimated_population_df <- estimated_population_df %>% mutate(accident_percapita=accident_count/population)

estimated_population_df <- estimated_population_df %>% 
                              left_join(published_land_price_df_summary,
                                        by=c("one_k_mesh"="mesh"))

estimated_population_df <- estimated_population_df %>% 
                                mutate(accident_ratio_class=if_else(estimated_population_df$accident_percapita > 0.05, "0.05以上",
                                                                    if_else(estimated_population_df$accident_percapita > 0.01, "0.01以上","0.01未満")))
estimated_population_df <- estimated_population_df %>% mutate(accident_percapita_100 = round(estimated_population_df$accident_percapita*100,2)  )


estimated_population_df_tokyo <- estimated_population_df %>% filter(city_code>=13000, city_code<13300)
estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(poligon="")

for (i in 1:nrow(estimated_population_df_tokyo)) {
  estimated_population_df_tokyo$poligon[i] <- export_mesh(estimated_population_df_tokyo$one_k_mesh[i])
}
# sfオブジェクトに変換
estimated_population_df_tokyo <- sf::st_as_sf(estimated_population_df_tokyo)

# マップでの可視化
estimated_population_df_tokyo %>% mapview::mapview(zcol = "accident_percapita_100")

100

101

102

103

104

105

106

107

108

library(jpmesh)

library(sf)

library(tidyverse)

library(kokudosuuchi)

library(jpndistrict)

library(mapview)

# load dataset

load(file = "accident_df_with_coordinate_from_map.RData")

# メッシュデータの追加

accident_df <- accident_df %>% mutate(k_mesh="")

accident_df <- accident_df %>% filter(latitude != "", !is.na(upload_date))

# 緯度、軽度の情報から任意のメッシュのデータを得る

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

accident_df$k_mesh[i] <- jpmesh::coords_to_mesh(

as.numeric(accident_df$longitude[i]),

as.numeric(accident_df$latitude[i]),

mesh_size = "1km")

}

# メッシュを整数にする

accident_df <- accident_df %>% mutate(k_mesh=as.integer(k_mesh))

# 東京都に絞る

accident_df <- accident_df %>% filter(stringr::str_detect(formatted_address,"東京"))

accident_df %>% group_by(k_mesh) %>% count() %>% arrange(desc(n))

# 公示地価

published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",

cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),

address=as.character(published_land_price$`L01-15`$L01_019),

area_code=as.integer(published_land_price$`L01-15`$L01_017))

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

published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])

published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])

}

# 東京の地価に絞り込む

published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),

stringr::str_detect(address,"区"))

published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)

published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

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

published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]

}

# 緯度経度を1kmメッシュにする

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

published_land_price_df$mesh[i] <- jpmesh::coords_to_mesh(

as.numeric(published_land_price_df$longitude[i]),

as.numeric(published_land_price_df$latitude[i]),

mesh_size = "1km")

}

published_land_price_df <- published_land_price_df %>% mutate(mesh=as.integer(mesh))

published_land_price_df_summary <- published_land_price_df %>% group_by(mesh) %>% summarise(mean_price=mean(price))

# 推計人口

estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",

cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,

city_code=estimated_population$Mesh3_POP_00$CITY_CODE,

population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>%

mutate(one_k_mesh=as.integer(one_k_mesh))

accident_summary <- accident_df %>%

group_by(k_mesh) %>%

summarise(accident_count=n()) %>%

arrange(desc(accident_count))

estimated_population_df <- estimated_population_df %>%

left_join(accident_summary,

by = c("one_k_mesh"="k_mesh"))

estimated_population_df$accident_count <- replace_na(estimated_population_df$accident_count, 0)

estimated_population_df <- estimated_population_df %>% mutate(accident_percapita=accident_count/population)

estimated_population_df <- estimated_population_df %>%

left_join(published_land_price_df_summary,

by=c("one_k_mesh"="mesh"))

estimated_population_df <- estimated_population_df %>%

mutate(accident_ratio_class=if_else(estimated_population_df$accident_percapita > 0.05, "0.05以上",

if_else(estimated_population_df$accident_percapita > 0.01, "0.01以上","0.01未満")))

estimated_population_df <- estimated_population_df %>% mutate(accident_percapita_100 = round(estimated_population_df$accident_percapita*100,2) )

estimated_population_df_tokyo <- estimated_population_df %>% filter(city_code>=13000, city_code<13300)

estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(poligon="")

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

estimated_population_df_tokyo$poligon[i] <- export_mesh(estimated_population_df_tokyo$one_k_mesh[i])

}

# sfオブジェクトに変換

estimated_population_df_tokyo <- sf::st_as_sf(estimated_population_df_tokyo)

# マップでの可視化

estimated_population_df_tokyo %>% mapview::mapview(zcol = "accident_percapita_100")

こちらのmapviewパッケージで作ったマップはインタラクティブにいじることができます。ぜひ関心のあるエリアでいじってみてください。

1kmメッシュ人口あたりのAP発生件数（×100）の可視化

<br />

1kmメッシュでのAP発生件数の可視化

<br />

比率ベースで、色の明るいメッシュのところを見ると、港区、中央区、新宿区、渋谷区などがAPが発生しやすいようです。件数ベースで言うと新宿が一番多いですね。
一番色が明るい港区はてっきり六本木ではないかと思ったのですが、新橋から日比谷にかけたエリアでした。会社員による自○が多いようです。恐ろしいものです。

APの名前の集計

APの名前を集計してみます。これは別にこの名前だからAPになりやすいというわけではなく、単純に数が多いだけの可能性がありますし、実際にそうだろうと思われます。AP発生率を知るには、APではないものも含めた全物件名に占めるAP発生件名を手に入れないといけませんが、全物件名を収集するのが難しいことから単純に頻度の集計となります。今回は、wordcloud2パッケージを使って、ワードクラウドにしてみます。文字が大きいと頻度が高いものとなります。

library('wordcloud2')
library(RMeCab)

load(file = "accident_df_with_coordinate__from_map.RData")

address_bow <- docMatrixDF(accident_df$address,minFreq=1)
word_list <- data.frame(word=rownames(address_bow), count=rowSums(address_bow))
word_list %>% filter(count<100, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

library('wordcloud2')

library(RMeCab)

load(file = "accident_df_with_coordinate__from_map.RData")

address_bow <- docMatrixDF(accident_df$address,minFreq=1)

word_list <- data.frame(word=rownames(address_bow), count=rowSums(address_bow))

word_list %>% filter(count<100, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

ハイツ、荘、コーポ、マンション、号棟、アパート、ハウスなどが多く出現しているようです。ただ、物件の名前としても頻度が高いとも考えられますね。

地価と人口あたりのAP発生件数の関係

ここでは地価のデータとAP発生の関係性について見てみます。

estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(bin = ntile(mean_price, 10))
vis_df <- estimated_population_df_tokyo %>% filter(accident_count>0) %>% group_by(bin) %>% summarize(mean_accident_percapita=mean(accident_percapita),
                                                                                                     count=n(),
                                                                                                     mean_price=mean(mean_price))

g <- ggplot(data = vis_df, aes(x =  bin, y =mean_accident_percapita)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))
g <- g +geom_point() + stat_smooth(se=T,fullrange = T,level = 0.95)
g <- g + ggtitle("地価の階級値と平均人口あたりのAP発生件数")
g

estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(bin = ntile(mean_price, 10))

vis_df <- estimated_population_df_tokyo %>% filter(accident_count>0) %>% group_by(bin) %>% summarize(mean_accident_percapita=mean(accident_percapita),

count=n(),

mean_price=mean(mean_price))

g <- ggplot(data = vis_df, aes(x = bin, y =mean_accident_percapita)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))

g <- g +geom_point() + stat_smooth(se=T,fullrange = T,level = 0.95)

g <- g + ggtitle("地価の階級値と平均人口あたりのAP発生件数")

地価の階級値（10個のパーセンタイルに分割）を横軸に、縦軸に人口あたりAP発生数をおくと、地価が上がるに従い人口あたりAP発生数が高まる傾向があります。これは、人口密度が高く地価の高いところではAPが発生しやすいということを示しているのではないでしょうか。人口密度が高いと地価があがる、人口密度が高いと治安が悪くなるという可能性が考えられます。

APの詳細の集計

ここではAPになってしまった詳細の内容について先ほどと同様に形態素解析を行いワードクラウドにしてみます。

accident_bow <- docMatrixDF(accident_df$detail,minFreq=1)
word_list <- data.frame(word=rownames(accident_bow), count=rowSums(accident_bow))
word_list %>% filter(count<300, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

accident_bow <- docMatrixDF(accident_df$detail,minFreq=1)

word_list <- data.frame(word=rownames(accident_bow), count=rowSums(accident_bow))

word_list %>% filter(count<300, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

どうやら孤独死が多いようです。高齢者の人口構成比が関係しているのだろうと思われます。

APの発生時期に関するテキストマイニング

ここでは、発生時期に含まれる四桁の数字を集計して、何年くらいのAPが多いのかをざっくりと掴みます。

date_bow <- docMatrixDF(accident_df$date,minFreq=1)
word_list <- data.frame(word=rownames(date_bow), count=rowSums(date_bow))
word_list %>% filter(stringr::str_detect(word, "[0-9][0-9][0-9][0-9]")) %>% wordcloud2()

date_bow <- docMatrixDF(accident_df$date,minFreq=1)

word_list <- data.frame(word=rownames(date_bow), count=rowSums(date_bow))

word_list %>% filter(stringr::str_detect(word, "[0-9][0-9][0-9][0-9]")) %>% wordcloud2()

どうやら昔のデータはあまり登録されていないようです。記憶が確かではないかもしれませんし、古すぎるものは消されているのかもしれませんね。あのサイトはユーザー生成コンテンツ（UGC）なので、投稿する人はそこまで昔のことをわざわざ投稿しないのかもしれないですね。

APの詳細に関するテキストマイニング

ここではトピック数10として、topicmodelsパッケージを使いLDAを行います。

library(RMeCab)
library(topicmodels)

accident_df$detail_fixed <- gsub(x = accident_df$detail,
                                 pattern = "\\(|\\)|（|）|/|\\.|:|_|-|[0-9]|[０-９]",
                                 replacement = " ")

# Bag of wordsの生成
res <- docMatrixDF(accident_df$detail_fixed,minFreq=4)
colnames(res) <- accident_df$k_mesh
res <- t(res)
res <- res[rowSums(res) >=1, ]

# トピック数
k <- 10
# 出力単語数
i <- 5
# 出力トピック数
j <- 5

# LDAの推定
LDA_estimate <- LDA(res, k,method="Gibbs",control=list(verbose=1))

# トピックごとの単語の出力
terms_each_topics <- data.frame(terms(LDA_estimate,i))

library(RMeCab)

library(topicmodels)

accident_df$detail_fixed <- gsub(x = accident_df$detail,

pattern = "\$|\$|（|）|/|\\.|:|_|-|[0-9]|[０-９]",

replacement = " ")

# Bag of wordsの生成

res <- docMatrixDF(accident_df$detail_fixed,minFreq=4)

colnames(res) <- accident_df$k_mesh

res <- t(res)

res <- res[rowSums(res) >=1, ]

# トピック数

k <- 10

# 出力単語数

i <- 5

# 出力トピック数

j <- 5

# LDAの推定

LDA_estimate <- LDA(res, k,method="Gibbs",control=list(verbose=1))

# トピックごとの単語の出力

terms_each_topics <- data.frame(terms(LDA_estimate,i))

なかなか恐ろしいキーワードが多いですが、なんとなくですがうまく分類されているのではないかと思われます。

トピック1は男性の不幸
トピック2は不動産屋に言われた告知事項
トピック3は孤独死
トピック4は病死
トピック5は火災・転落・事故
トピック6は事故のあった建物に関する記載
トピック7は腐乱した事例
トピック8は建物に関して不明であることの記載
トピック9は心理的瑕疵あり
トピック10は自○

となっているように思われます。まさかこのようなデータにトピックモデルを使うことになるとは。

おわりに

今回はR言語のみを用いて、APに関するデータを収集し、地図にプロットしたり他のメッシュデータとつなぎ合わせて分析をするなどしました。APが発生しやすいエリア、APと地価との関係、APのテキストマイニングなど興味深い結果が得られたと思います。
一つ残念なのは、時系列情報がフリーテキストなので、APがどのエリアでどの頻度で発生していくのかの分析のコストが高く、今回は時系列情報を用いた分析にチャレンジできませんでした。
今後はタクシーの需要推定の分析で行われているように、メッシュ単位でのAP発生確率の推定などを機械学習で行えると面白いなと思います。どなたか一緒にアノテーションしましょう！

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

参考情報

数多くの方々の記事を見てどうにか仕上げることができました。感謝します。

[1]【追記あり】sfパッケージでシェープファイルを読み込んでmapviewパッケージで可視化するまで
[2]How to use mesh cord in R
[3]Rを使ってワードクラウドを作ってみました
[4]国土数値情報ダウンロードサービスWeb APIからデータを取得するためのRパッケージです
[5]東京の地価公示データを眺める
[6]Chapter 1 Introduction to spatial data in R
[7][翻訳] RSelenium vignette: RSeleniumの基本
[8]RからYahoo!のジオコーディングを利用する方法
[9]EMBEDDING A LEAFLET MAP ON WORDPRESS
[10]mapview advanced controls
[11]RSeleniumでChromeからファイルをダウンロードするディレクトリを指定する方法
[12]Selenium Serverが立ち上がらないときはportが被っているかも！？
[13]brew install selenium-server-standalone
[14]ナウでヤングなRの環境変数管理方法
[15]タクシードライバー向け需要予測について
[16]LDA with topicmodels package for R, how do I get the topic probability for each term?
[17]dplyr — 高速data.frame処理

[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（ダウンロード注意）

目次

はじめに

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

終わりに

はじめに

分析工程

データの収集

データの整形

単純集計・可視化

分析

APの発生件数の集計

人口あたりのAP件数

APの名前の集計

地価と人口あたりのAP発生件数の関係

APの詳細の集計

APの発生時期に関するテキストマイニング

APの詳細に関するテキストマイニング

おわりに

参考情報

はじめに

データについて

分析方針

分析概要

データ整形

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

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

stanの実行

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

テストデータでの予測

非階層モデルとの比較

おわりに

参考情報

モデル

数式の展開

Rのコードの紹介

おわりに

参考情報

はじめに

生存時間分析とは

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

先行研究

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

Stanでの実行例

おわりに

参考文献

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

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

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

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

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

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