Bayesian Statistics and Marketingの5章 – 家計の異質性を考慮した階層ベイズモデル

はじめに

ゴールデンウィークで実家に持ち込む本としてチョイスしたのが、2005年出版の「Bayesian Statistics and Marketing」です。大学院のときに購入して、ちょっとしか読んでませんでした。

この本は、字面の通りマーケティング関連の分析に関してベイズ統計を使ってアプローチするというもので、この書籍のために作られた、Rのbayesmというパッケージの紹介もあり、理論だけでなくRで実践することもできます。1章から7章までの全ての分析事例に対して実行可能な関数が用意されています。（CRANにあるdocumentも120p程度と割と大きめのパッケージです。）

和書で言うと、東北大学の照井先生の「ベイズモデリングによるマーケティング分析」などがありますが、その82pでもBayesian Statistics and Marketingとbayesmパッケージが紹介されています。

今回は、5章に載っている階層ベイズモデルを用いた、家計の異質性を考慮したブランド選択モデルの分析を紹介します。加えて、GitHubでstanによる再現を試みている方がいらっしゃったので、その方のコードの紹介も行います。

最近はこれまで以上にベイズ統計が流行ってきていますが、マーケティング×ベイズの書籍は限られている印象なので、少しでもリサーチのお役に立てれば幸いです。

目的

マーガリンの購買データから、ブランドごと、家計ごとのマーガリン価格に対しての反応の違いを明らかにしたい。

データ

bayesmパッケージにある、margarineデータ。data(margarine)で呼び出せ、詳細はcranのドキュメントに載っています。

Household Panel Data on Margarine Purchasesには、516家計の購買データと、家計ごとのデモグラフィック情報が収められています。1991年の論文のデータとなるので、かなり昔のデータです。

購買データは価格（USドル）と選択したブランドのID（10種類）
デモグラフィック情報はfamily size（家族構成）、学歴、職位、退職の有無などのダミー変数

今回の事例では、5回以上購買した家計に限定して分析しているため、
313家計・3405の購買レコードからなるデータセットとなります。

モデル

家計ごとに異なる、マーガリン価格に対する反応を想定。各マーガリンのブランドの価格に対するパラメータの数は家計の数だけある。
価格に対する反応は家計の属性によっても決まる。
という前提に立ち、以下のセッティングで推論していきます。
6つのブランド選択に関する多項ロジスティックモデル（カテゴリカル分布とsoftmax関数の適用）
1階層目はブランドごとの価格を説明変数とし、価格に対する反応係数をかけ合わせたものを多項ロジスティックモデルの入力とする。
2階層目はブランドの価格に対する反応係数が家計ごとの定数項と属性データに属性ごとの係数をかけ合わせたものからなる。
家計ごとの定数項は平均0、分散V_betaの正規分布に従う。
属性ごとの係数は平均vec(delta_bar)、分散V_betaクロネッカーのデルタA^(-1)の正規分布に従う。
分散V_betaは平均υ、分散Vの逆ウィシャート分布に従う。
A = 0.01、υ = 6 + 3 = 9、V = υI(Iは単位行列)

$\textit{ l } $ \beta_i | y_i , X_i $ [ Multinomial Logit ]$

$B = Z \varDelta + U , u_i \sim N(0,V_\beta)$

$vec(\varDelta | V_\beta ) \sim N(vec(\bar{\varDelta}), V_\beta \otimes A^{-1} )$

$V_\beta \sim IW(\upsilon ,V )$

コード

kefitsさんがいくつかの章に登場するbayesmでの実践例をstanに書き直されているようですので、そちらのコードで学ばせていただこうと思います。
https://github.com/kefits/Bayesian-Statistics-and-Marketing

以下が、stanのコードとなっています。ここでは、Hierarchical_MNL.stanとして保存します。

data{
  int<lower=0> N_x; // 購買レコードの数
  int<lower=0> N_z; // 家計の数
  int<lower=0> p_x; // 購買レコードの項目数
  int<lower=0> p_z; // 家計の属性データの項目数
  
  int y[N_x]; // 選択肢
  matrix[N_x, p_x] X; // 説明変数
  matrix[N_z, p_z] Z; // 家計の属性データ
  int<lower=0> hhid[N_x];  // 家計ID
}

transformed data{
  real nu;
  matrix[p_x, p_x] I; // 説明変数の数の正方行列
  
  nu = p_x + 3; // 説明変数の項に3を足す
  I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成
}

parameters{
  vector[p_x] beta_ast[N_z]; // 説明変数の数だけある、家計ごとのパラメータ
  matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ
  cov_matrix[p_x] V_b; // 共分散行列
}

transformed parameters{
  vector[p_x] beta[N_z]; #家計の数だけの係数ベクトル
  matrix[p_x, p_x] L_b; #共分散行列（beta(家計ごとの係数の共分散)）
  matrix[p_x, p_x] L_d; #共分散行列（delta(属性データの係数の共分散)）
  
  L_b = cholesky_decompose(V_b); // 共分散行列のコレスキー因子をもとめる
  L_d = cholesky_decompose(100*V_b); // 共分散行列に0.01で割ったもののコレスキー因子をもとめる
  for(i in 1:N_z){
    beta[i] = beta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる 
  }
}

model{
  for(i in 1:N_x){
    y[i] ~ categorical(softmax(beta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う
  }
  for(i in 1:p_z){
    Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), L_d); // コレスキー因子（L_d）を引数にとる多変量正規分布(推定の高速化のために用いることがある。)
  }
  beta_ast ~ multi_normal_cholesky(rep_vector(0, p_x), L_b); // コレスキー因子（L_b）を引数にとる多変量正規分布
  V_b ~ inv_wishart(nu, nu*I); // 正規分布の共分散行列の共役事前分布として逆ウィシャート分布を利用
}

data{

int<lower=0> N_x; // 購買レコードの数

int<lower=0> N_z; // 家計の数

int<lower=0> p_x; // 購買レコードの項目数

int<lower=0> p_z; // 家計の属性データの項目数

int y[N_x]; // 選択肢

matrix[N_x, p_x] X; // 説明変数

matrix[N_z, p_z] Z; // 家計の属性データ

int<lower=0> hhid[N_x]; // 家計ID

}

transformed data{

real nu;

matrix[p_x, p_x] I; // 説明変数の数の正方行列

nu = p_x + 3; // 説明変数の項に3を足す

I = diag_matrix(rep_vector(1, p_x)); // 1を繰り返しp_x個並べた対角行列を作成

}

parameters{

vector[p_x] beta_ast[N_z]; // 説明変数の数だけある、家計ごとのパラメータ

matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ

cov_matrix[p_x] V_b; // 共分散行列

}

transformed parameters{

vector[p_x] beta[N_z]; #家計の数だけの係数ベクトル

matrix[p_x, p_x] L_b; #共分散行列（beta(家計ごとの係数の共分散)）

matrix[p_x, p_x] L_d; #共分散行列（delta(属性データの係数の共分散)）

L_b = cholesky_decompose(V_b); // 共分散行列のコレスキー因子をもとめる

L_d = cholesky_decompose(100*V_b); // 共分散行列に0.01で割ったもののコレスキー因子をもとめる

for(i in 1:N_z){

beta[i] = beta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる

}

model{

for(i in 1:N_x){

y[i] ~ categorical(softmax(beta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う

}

for(i in 1:p_z){

Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), L_d); // コレスキー因子（L_d）を引数にとる多変量正規分布(推定の高速化のために用いることがある。)

}

beta_ast ~ multi_normal_cholesky(rep_vector(0, p_x), L_b); // コレスキー因子（L_b）を引数にとる多変量正規分布

V_b ~ inv_wishart(nu, nu*I); // 正規分布の共分散行列の共役事前分布として逆ウィシャート分布を利用

}

以下はstanをキックするためのRコードです。

library(bayesm)
library(dplyr)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。
hhid_selected <- margarine$choicePrice %>% 
                  filter(choice %in% c(1,2,3,4,5,7)) %>% 
                  group_by(hhid) %>% 
                  summarise(purc_cnt = n()) %>% 
                  filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。
choicePrice.selected <- margarine$choicePrice %>% 
                          filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)
#並べにくいので7を6に置き換える。
choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出
demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ
N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）
p <- n_distinct(choicePrice.selected$choice)

#被説明変数
y <- choicePrice.selected$choice

#説明変数
X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く
Z <- demos.selected %>% 
        select(-hhid)

#定数項を1列目に追加する
Z <- data.frame(intercept = rep(1, nrow(Z))) %>% 
        bind_cols(Z)

#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。
hhid_index <- demos.selected %>%
                select(hhid) %>% 
                mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする
hhid_x <- choicePrice.selected %>% 
            select(hhid) %>% 
            left_join(hhid_index)

#stanで扱うデータリストの作成
d.dat <- list(N_x=nrow(X), N_z=nrow(Z), 
              p_x=ncol(X), p_z=ncol(Z),
              y=y, X=X, Z=Z,
              hhid = hhid_x$ind)

#推定
d.fit <- stan("../Chapter5/Hierarchical_MNL.stan",
              data = d.dat,
              iter = 500,
              chains = 4)

library(bayesm)

library(dplyr)

library(rstan)

rstan_options(auto_write = TRUE)

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

data("margarine")

#1,2,3,4,5,7の商品に関してデータを抽出し、家計IDごとにカウントし、5件以上のものに絞る。

hhid_selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7)) %>%

group_by(hhid) %>%

summarise(purc_cnt = n()) %>%

filter(purc_cnt >= 5)

#今回扱う商品のカラムだけを抽出し、先ほど絞ったユーザーのリストに合致するデータでフィルターする。

choicePrice.selected <- margarine$choicePrice %>%

filter(choice %in% c(1,2,3,4,5,7) & hhid %in% hhid_selected$hhid)

#並べにくいので7を6に置き換える。

choicePrice.selected$choice[choicePrice.selected$choice == 7] <- 6

#家計ごとに関する属性データの抽出

demos.selected <- margarine$demos %>% filter(hhid %in% hhid_selected$hhid)

#データサイズ

N <- nrow(choicePrice.selected)

#選択肢の数（特に使っているデータではない。）

p <- n_distinct(choicePrice.selected$choice)

#被説明変数

y <- choicePrice.selected$choice

#説明変数

X <- choicePrice.selected %>% select(3,4,5,6,7,9)

#家計の属性データから家計IDを除く

Z <- demos.selected %>%

select(-hhid)

#定数項を1列目に追加する

Z <- data.frame(intercept = rep(1, nrow(Z))) %>%

bind_cols(Z)

#家計の属性データから家計IDを抽出し、1から行数までのインデックスを付与する。

hhid_index <- demos.selected %>%

select(hhid) %>%

mutate(ind = seq(1,nrow(demos.selected)))

#購買データの家計IDを抽出し、先ほどのインデックスとjoinする

hhid_x <- choicePrice.selected %>%

select(hhid) %>%

left_join(hhid_index)

#stanで扱うデータリストの作成

d.dat <- list(N_x=nrow(X), N_z=nrow(Z),

p_x=ncol(X), p_z=ncol(Z),

y=y, X=X, Z=Z,

hhid = hhid_x$ind)

#推定

d.fit <- stan("../Chapter5/Hierarchical_MNL.stan",

data = d.dat,

iter = 500,

chains = 4)

実行結果

Core i5、8GBメモリのMacBook Proで40分ほどかかりました。

traceplot(d.fit)で以下のように4回の試行結果が描かれますが、収束しているようです。

summary関数を使えばわかりますが、3913行ものパラメータたちのサマリーが得られます。

313家計の家計ごとのブランドに対するパラメータ(1878個)
313家計の家計ごとのブランドに対する潜在パラメータ(1878個)
6ブランドの係数の共分散行列(36個)
6ブランドの係数の分散のハイパーパラメータの行列(36個)
6ブランドの属性データ(8つ)に対する係数(48個)
6ブランドの属性データに対する係数の共分散行列(36個)
lp(log posterior(確率密度の和でモデル比較で扱う。))(1個)

64番目の家計の各ブランドの価格に対する係数の分布を確認すると、4番目・5番目のブランドの係数が他のブランドに比べて小さいことがわかります。

続いて、家計ごとの係数に関して集計し、係数ごとの相関係数を見てみると、各ブランドごとに正の相関、負の相関がありそうです。

#トレースプロット
traceplot(d.fit)

#係数のサマリー
summary_table <- summary(d.fit)$summary

draws <- extract(d.fit)
beta <- as.data.frame(draws$beta)
Delta <- as.data.frame(draws$Delta)
V_b <- as.data.frame(draws$V_b)

hhid_info <- inner_join(hhid_index, hhid_selected)

# 1000行*313列のデータを313000行*1列のデータにしたい。
for (i in 1:6) {
  nam <- paste("beta", i, sep = "")
  assign(nam, beta[,(1+313*(i-1)):(313*(i))] %>% tidyr::gather(key, value))
}

beta_matrix <- beta1 %>% bind_cols(beta2,beta3,beta4,beta5,beta6)
beta_matrix <- beta_matrix %>% select(-starts_with("key"))

#相関係数
cor(beta_matrix)

            value     value1     value2     value3     value4    value5
value   1.0000000  0.5902734  0.4864998 -0.1798877 -0.4781558 0.3188025
value1  0.5902734  1.0000000  0.6134343 -0.2484286 -0.4441728 0.2002954
value2  0.4864998  0.6134343  1.0000000  0.1336322 -0.4512474 0.3663121
value3 -0.1798877 -0.2484286  0.1336322  1.0000000  0.6149186 0.2671819
value4 -0.4781558 -0.4441728 -0.4512474  0.6149186  1.0000000 0.1591574
value5  0.3188025  0.2002954  0.3663121  0.2671819  0.1591574 1.0000000

#トレースプロット

traceplot(d.fit)

#係数のサマリー

summary_table <- summary(d.fit)$summary

draws <- extract(d.fit)

beta <- as.data.frame(draws$beta)

Delta <- as.data.frame(draws$Delta)

V_b <- as.data.frame(draws$V_b)

hhid_info <- inner_join(hhid_index, hhid_selected)

# 1000行*313列のデータを313000行*1列のデータにしたい。

for (i in 1:6) {

nam <- paste("beta", i, sep = "")

assign(nam, beta[,(1+313*(i-1)):(313*(i))] %>% tidyr::gather(key, value))

}

beta_matrix <- beta1 %>% bind_cols(beta2,beta3,beta4,beta5,beta6)

beta_matrix <- beta_matrix %>% select(-starts_with("key"))

#相関係数

cor(beta_matrix)

value value1 value2 value3 value4 value5

value 1.0000000 0.5902734 0.4864998 -0.1798877 -0.4781558 0.3188025

value1 0.5902734 1.0000000 0.6134343 -0.2484286 -0.4441728 0.2002954

value2 0.4864998 0.6134343 1.0000000 0.1336322 -0.4512474 0.3663121

value3 -0.1798877 -0.2484286 0.1336322 1.0000000 0.6149186 0.2671819

value4 -0.4781558 -0.4441728 -0.4512474 0.6149186 1.0000000 0.1591574

value5 0.3188025 0.2002954 0.3663121 0.2671819 0.1591574 1.0000000

最後に、家計ごとに集計した、ブランドに対する価格反応係数の事後分布を描きます。

~~多峰性などはなく、正規分布に従っているようです。他のブランドと比較して、5番目の係数が小さいようです。~~

というのは誤りで、一週間後に気づいたのですが、家計ごとのブランドごとの係数の事後分布の平均値をプロットするべきでした。
正しくはこちらです。

事前情報として正規分布を仮定していましたが、係数に関して正規分布に従っていません。
そのため、事前情報として対称性のあるような正規分布を扱うのは適切ではなさそうです。

おわりに

2005年の本とは言え、十分に使いみちのある本だと思いました。まだまだ扱いきれていないですが、引き続き勉強していきます。
この本にはケーススタディが5つほどあるのですが、それのstanコード化などをしていけばかなり力がつくような気がします。

マーケティングの部署で働くデータアナリストにとって、マーケティング×ベイズの話は非常にモチベーションの上がるところなので、こういう文献を今後も見つけていきたい。

参考文献

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)
Bayesian Statistics and Marketingのサポートサイト
 ベイズモデリングによるマーケティング分析
 StanとRでベイズ統計モデリング (Wonderful R)
RStanのおさらいをしながら読む岩波DS 1 Shinya Uryu
Stanのlp__とは何なのかうなどん
 ‘LP__’ IN STAN OUTPUT
Package ‘bayesm’

R Advent Calendar 2017 rvestを用いてポケモンデータをスクレイピング&分析してみた

R Advent Calendar 2017の11日目を担当するMr_Sakaueです。
今回はrvestパッケージを用いて、友人がハマっているポケモンの情報を集めてみようと思います。
もっとも、業務でWebスクレイピングする際はPythonでBeautifulSoupやSeleniumを使うことがほとんどなのですが、たまにはRでやってみようと思います。

目次
・やりたいこと
・rvestについて
・データの取得と集計と可視化と分析
・まとめ
・参考情報

やりたいこと

今回はポケモンたちのデータを集めた上で、以下の内容を行いたいと思います。

ポケモンのサイトから種族値を取得
ポケモンの種族値を標準化して再度ランキング
ポケモンのレア度や経験値に関する情報を取得
レア度や経験値と相関しそうな種族値を探る

今回扱った全てのコードはこちらに載せております。
https://github.com/KamonohashiPerry/r_advent_calendar_2017/tree/master

※種族値はゲームにおける隠しパラメータとして設定されている、ポケモンの能力値とされている。

rvestについて

rvestはRでWebスクレイピングを簡単に行えるパッケージです。ここでの説明は不要に思われますが、今回はread_html()、html_nodes()、html_text()、html_attr()の4つ関数を用いました。

基本的に以下の3ステップでWebの情報を取得することができます。

STEP1
read_html()でHTMLからソースコードを取得する。(Pythonでいう、requestとBeautifulSoup)
STEP2
html_nodes()でソースコードから指定した要素を抽出する。(PythonでいうところのfindAll)
STEP3
html_text()やhtml_attr()で抽出した要素からテキストやリンクを抽出する。(Pythonでいうところのget(‘href’)など)

データの取得と集計と可視化

検索エンジンで検索してだいたい1位のサイトがあったので、そちらのWebサイトに載っているポケモンの種族値の一覧をスクレイピング対象とさせていただきます。

ポケモンのサイトから種族値を取得

library(rvest)
library(tidyverse)
library(magrittr)
library(reshape2)

# htmlソースコードを読み込む
pokemon_ranking <- read_html("https://yakkun.com/sm/status_list.htm")

# class属性がtdタグのノードを抽出
node_extracted <- html_nodes(pokemon_ranking, "td")

# ノードからテキストを抽出して行列にして、名前以外を数値に型変換して、変数名を変更する。
pokemon_data <- data.frame(matrix(html_text(node_extracted),
                       ncol = 9,byrow = TRUE),stringsAsFactors = FALSE) %>% 
                set_colnames(c('id', 'name', 'Hit_Points', 'Attack',
                               'Defense', 'Special_Attack',
                               'Special_Defense', 'Speed', 'Total')) %>%
                mutate_at(vars(-name), as.numeric)


# ポケモン別のページを取得するためのURLの取得
pokemon_link <- pokemon_ranking %>% html_nodes("td") %>% html_nodes('a') %>% html_attr('href')
pokemon_link <- gsub(x = pokemon_link,
                     pattern = './zukan',
                     replacement = "https://yakkun.com/sm/zukan")
# ポケモン別のURLを先ほどのデータに加える
pokemon_data <- pokemon_data %>% mutate(url = pokemon_link)

library(rvest)

library(tidyverse)

library(magrittr)

library(reshape2)

# htmlソースコードを読み込む

pokemon_ranking <- read_html("https://yakkun.com/sm/status_list.htm")

# class属性がtdタグのノードを抽出

node_extracted <- html_nodes(pokemon_ranking, "td")

# ノードからテキストを抽出して行列にして、名前以外を数値に型変換して、変数名を変更する。

pokemon_data <- data.frame(matrix(html_text(node_extracted),

ncol = 9,byrow = TRUE),stringsAsFactors = FALSE) %>%

set_colnames(c('id', 'name', 'Hit_Points', 'Attack',

'Defense', 'Special_Attack',

'Special_Defense', 'Speed', 'Total')) %>%

mutate_at(vars(-name), as.numeric)

# ポケモン別のページを取得するためのURLの取得

pokemon_link <- pokemon_ranking %>% html_nodes("td") %>% html_nodes('a') %>% html_attr('href')

pokemon_link <- gsub(x = pokemon_link,

pattern = './zukan',

replacement = "https://yakkun.com/sm/zukan")

# ポケモン別のURLを先ほどのデータに加える

pokemon_data <- pokemon_data %>% mutate(url = pokemon_link)

以上のコードを実行すれば、こんな感じでポケモンの種族値一覧を得る事ができます。

とりあえず、種族値合計（Total Tribal Value 以下、TTV）のランキングの上位を確認してみます。知らないんですが、メガミュウツーとかいうイカつそうなポケモンが上位にいるようです。昭和の世代には縁のなさそうなポケモンばかりですねぇ。

■TTVランキング

取得した種族値を項目別に集計したり、Boxプロットを描いてみます。どうやら、攻撃の平均が高く、ヒットポイントや素早さの平均は低いようです。

# 集計
pokemon_data_melt <- melt(pokemon_data %>% select(-url), id.vars = 'name')
pokemon_data_melt %>% 
    group_by(variable) %>% 
    summarise(mean = mean(value),
              median = median(value),
              sd = sd(value),
              max = max(value),
              min = min(value),
              cv = sd/mean)

# 集計

pokemon_data_melt <- melt(pokemon_data %>% select(-url), id.vars = 'name')

pokemon_data_melt %>%

group_by(variable) %>%

summarise(mean = mean(value),

median = median(value),

sd = sd(value),

max = max(value),

min = min(value),

cv = sd/mean)

■種族値のサマリー

# Box-Plotを描く
ggplot(data = pokemon_data_melt %>% filter(!(variable %in% c('id','Total'))),
       aes(x = variable, y = value)) +
   geom_boxplot() + ggtitle("Tribal Value") + 
   theme(plot.title = element_text(hjust = 0.5)) + coord_flip()

# Box-Plotを描く

ggplot(data = pokemon_data_melt %>% filter(!(variable %in% c('id','Total'))),

aes(x = variable, y = value)) +

geom_boxplot() + ggtitle("Tribal Value") +

theme(plot.title = element_text(hjust = 0.5)) + coord_flip()

■種族値のBoxプロット

ポケモンの種族値を標準化して再度ランキング

さて、攻撃の平均が高かったり、ヒットポイントと素早さの平均が低かったりしたので、各々の項目を標準化した上で、再度ランキングを作ってみたいと思います。

pokemon_data_standardized <- pokemon_data
pokemon_data_standardized <- pokemon_data_standardized %>% 
                              mutate_at(vars(Hit_Points,
                                             Attack,
                                             Defense,
                                             Special_Attack,
                                             Special_Defense,
                                             Speed),funs(scale(.) %>% as.vector))

pokemon_data_standardized <- pokemon_data_standardized %>% 
                              mutate(Total = rowSums(select(.,c(3:8))))

pokemon_data_standardized_ranking <- pokemon_data_standardized %>% arrange(desc(Total))
pokemon_data_standardized_ranking <- pokemon_data_standardized_ranking %>% mutate(standardized_ranking = 1:n())


pokemon_data_standardized_melt <- melt(pokemon_data_standardized %>% select(-url), id.vars = 'name')
pokemon_data_standardized_melt %>% 
  group_by(variable) %>% 
  summarise(mean = mean(value),
            median = median(value),
            sd = sd(value),
            max = max(value),
            min = min(value),
            cv = sd/mean)

pokemon_data_standardized <- pokemon_data

pokemon_data_standardized <- pokemon_data_standardized %>%

mutate_at(vars(Hit_Points,

Attack,

Defense,

Special_Attack,

Special_Defense,

Speed),funs(scale(.) %>% as.vector))

pokemon_data_standardized <- pokemon_data_standardized %>%

mutate(Total = rowSums(select(.,c(3:8))))

pokemon_data_standardized_ranking <- pokemon_data_standardized %>% arrange(desc(Total))

pokemon_data_standardized_ranking <- pokemon_data_standardized_ranking %>% mutate(standardized_ranking = 1:n())

pokemon_data_standardized_melt <- melt(pokemon_data_standardized %>% select(-url), id.vars = 'name')

pokemon_data_standardized_melt %>%

group_by(variable) %>%

summarise(mean = mean(value),

median = median(value),

sd = sd(value),

max = max(value),

min = min(value),

cv = sd/mean)

■標準化した種族値のサマリー

平均0、分散1にできているようです。

# Box-Plotを描く
ggplot(data = pokemon_data_standardized_melt %>% filter(!(variable %in% c('id','Total','Total_standardized'))),
       aes(x = variable, y = value)) +
  geom_boxplot() + ggtitle("Tribal Value") + 
  theme(plot.title = element_text(hjust = 0.5)) + coord_flip()

# Box-Plotを描く

ggplot(data = pokemon_data_standardized_melt %>% filter(!(variable %in% c('id','Total','Total_standardized'))),

aes(x = variable, y = value)) +

geom_boxplot() + ggtitle("Tribal Value") +

theme(plot.title = element_text(hjust = 0.5)) + coord_flip()

■標準化した種族値のBoxプロット

他よりも低かったヒットポイントと、高かった攻撃がならされていることが確認できます。

■標準化前後でのTTVランキングのギャップが大きかったものをピックアップ

ラッキーが144位ほど出世しています。攻撃が低く、ヒットポイントの高いラッキーが標準化により優遇されるようになったと考える事ができます。ポケモン大会の上位ランカーである後輩社員もラッキーは手強いですと言っていたのでまんざらでもないのでしょう。

ポケモンのレア度や経験値に関する情報を取得

今回のサイトには、個別にポケモン別のページが用意されており、そちらから、ゲットしやすさや経験値に関する情報を抽出します。

# ポケモンの個別ページの情報を格納するデータフレームの作成
pokemon_detail_database <- data.frame(url = as.character(),
                                      name = as.character(),
                                      rarity = as.integer(),
                                      experience = as.integer())

# ポケモン別のURLからゲットしやすさなどを抽出するための関数
Pokemon_Detail_Get <- function(pokemon_url){
  pokemon_detail <- read_html(pokemon_url)
  
  # XPathで名前とゲットしやすさと経験値タイプを取得
  node_extracted_pokemon_name <- pokemon_detail %>% html_nodes(xpath="//tr[1]") %>% html_text()
  node_extracted_pokemon_name <- node_extracted_pokemon_name[1]
  
  node_extracted_pokemon_get <- pokemon_detail %>% html_nodes(xpath="//tr[24]/td[2]") %>% html_text()
  node_extracted_pokemon_get <- as.integer(gsub(x = node_extracted_pokemon_get[1], pattern = "\u00A0", replacement = ""))
  
  node_extracted_pokemon_exp <- pokemon_detail %>% html_nodes(xpath="//tr[26]/td[2]") %>% html_text()
  node_extracted_pokemon_exp <- as.integer(gsub(x = node_extracted_pokemon_exp[1], pattern = "万", replacement = "0000"))
  
  pokemon_detail_data <- data.frame(url = pokemon_url,
                                    name = node_extracted_pokemon_name,
                                    rarity = node_extracted_pokemon_get,
                                    experience = node_extracted_pokemon_exp)
  return(pokemon_detail_data)
  
  Sys.sleep(30)
}

# ポケモン別のページをスクレイピングする
pokemon_detail_database <- map_dfr(pokemon_link ,
                               ~Pokemon_Detail_Get(.))

# 重複したURLを削除する
pokemon_detail_database <- pokemon_detail_database %>% distinct(url, .keep_all = TRUE)

# 種族値のデータとゲットしやすさなどのデータを繋ぎこむ
pokemon_data_standardized <- pokemon_data_standardized %>% left_join(pokemon_detail_database %>% select(-name), by ="url")

# ポケモンの個別ページの情報を格納するデータフレームの作成

pokemon_detail_database <- data.frame(url = as.character(),

name = as.character(),

rarity = as.integer(),

experience = as.integer())

# ポケモン別のURLからゲットしやすさなどを抽出するための関数

Pokemon_Detail_Get <- function(pokemon_url){

pokemon_detail <- read_html(pokemon_url)

# XPathで名前とゲットしやすさと経験値タイプを取得

node_extracted_pokemon_name <- pokemon_detail %>% html_nodes(xpath="//tr[1]") %>% html_text()

node_extracted_pokemon_name <- node_extracted_pokemon_name[1]

node_extracted_pokemon_get <- pokemon_detail %>% html_nodes(xpath="//tr[24]/td[2]") %>% html_text()

node_extracted_pokemon_get <- as.integer(gsub(x = node_extracted_pokemon_get[1], pattern = "\u00A0", replacement = ""))

node_extracted_pokemon_exp <- pokemon_detail %>% html_nodes(xpath="//tr[26]/td[2]") %>% html_text()

node_extracted_pokemon_exp <- as.integer(gsub(x = node_extracted_pokemon_exp[1], pattern = "万", replacement = "0000"))

pokemon_detail_data <- data.frame(url = pokemon_url,

name = node_extracted_pokemon_name,

rarity = node_extracted_pokemon_get,

experience = node_extracted_pokemon_exp)

return(pokemon_detail_data)

Sys.sleep(30)

}

# ポケモン別のページをスクレイピングする

pokemon_detail_database <- map_dfr(pokemon_link ,

~Pokemon_Detail_Get(.))

# 重複したURLを削除する

pokemon_detail_database <- pokemon_detail_database %>% distinct(url, .keep_all = TRUE)

# 種族値のデータとゲットしやすさなどのデータを繋ぎこむ

pokemon_data_standardized <- pokemon_data_standardized %>% left_join(pokemon_detail_database %>% select(-name), by ="url")

以上のコードを実行すれば、やや時間がかかりますが、全ポケモンのゲットしやすさや経験値のデータを抽出する事ができます。それらの情報がゲットできたら、まずは可視化します。

# ゲットしやすさのヒストグラム
ggplot(data = pokemon_data_standardized, aes(x = rarity)) + geom_histogram() 

# 経験値のヒストグラム
ggplot(data = pokemon_data_standardized, aes(x = experience)) + geom_histogram()

# ゲットしやすさのヒストグラム

ggplot(data = pokemon_data_standardized, aes(x = rarity)) + geom_histogram()

# 経験値のヒストグラム

ggplot(data = pokemon_data_standardized, aes(x = experience)) + geom_histogram()

■ゲットしやすさのヒストグラム

ゲットのしやすさは、小さいほど捕まえる難易度が高くなっています。難易度の高いポケモンである0が多過ぎるので、このデータは欠損値が0になっているのではないかと疑われます。

■経験値のヒストグラム

経験値は、レベル100になるまでに要する経験値をさしています。ほとんどが100万程度となっているようです。

■ゲットしやすさと標準化TTVの散布図

やはり、ゲットしやすさに関してはデータに不備があるようで、コラッタ（アローラの姿）のような雑魚ポケのゲットのしやすさが0だったり、伝説のポケモンであるネクロズマが255だったりします。ただ、上限と下限のデータを間引けば右下がりの傾向が見られそうです。

■経験値と標準化TTVの散布図

経験値が多く必要にも関わらず、TTVが低い集団があります。どうやらこの集団に属するのは、「キノガッサ」・「マクノシタ」・「イルミーゼ」・「ゴクリン」・「シザリガー」などで、一回しか進化しないポケモンのようです。これらのポケモンは育てにくく、TTVの低い、コスパの悪そうなポケモンと考えることができるのではないでしょうか。（技や特性によってはバリューあるかもしれませんが。）

レア度や経験値と相関しそうな種族値を探る

先ほどのレア度に関しては、データがおかしそうだったので、レア度0と255に関しては除外してみます。

# おかしそうなレア度0と255のデータを除外する。
pokemon_data_standardized_filtered <- pokemon_data_standardized %>% filter(rarity > 0, rarity < 255)

# ゲットのしやすさと標準化TTV
ggplot(data = pokemon_data_standardized_filtered, aes(x = rarity, y = Total)) + 
  geom_point() + ylab('Total Tribal Value')

# おかしそうなレア度0と255のデータを除外する。

pokemon_data_standardized_filtered <- pokemon_data_standardized %>% filter(rarity > 0, rarity < 255)

# ゲットのしやすさと標準化TTV

ggplot(data = pokemon_data_standardized_filtered, aes(x = rarity, y = Total)) +

geom_point() + ylab('Total Tribal Value')

■ゲットしやすさと標準化TTVの散布図

やはり除外する事で、理想的な右下がりの傾向を示す散布図が得られたと思います。
さて、各種族値がレア度にどれだけ相関しているのかを分析したいのですが、その前にレア度を表す二項変数を作成します。

■ゲットしやすさが50以下であれば1、それ以外を0にする変数を作成

pokemon_data_standardized_filtered <- pokemon_data_standardized_filtered %>% mutate(y = ifelse(rarity <= 50, 1, 0))

1	pokemon_data_standardized_filtered <- pokemon_data_standardized_filtered %>% mutate(y = ifelse(rarity <= 50, 1, 0))

続いて、各種族値を説明変数として、レア度を目的変数としたロジスティック回帰モデルの推定をrstanで実行させます。

■stanコード

data {
  int N;
  real Hit_Points[N];
  real Attack[N];
  real Defense[N];
  real Special_Attack[N];
  real Special_Defense[N];
  real Speed[N];
  int<lower=0, upper=1> Y[N];
}


parameters {
  real b[7];
}


model {
  for (n in 1:N)
    Y[n] ~ bernoulli_logit(b[1] + b[2]*Hit_Points[n] + b[3]*Attack[n] + b[4]*Defense[n] + b[5]*Special_Attack[n] + b[6]*Special_Defense[n] + b[7]*Speed[n]);
}

data {

int N;

real Hit_Points[N];

real Attack[N];

real Defense[N];

real Special_Attack[N];

real Special_Defense[N];

real Speed[N];

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

}

parameters {

real b[7];

}

model {

for (n in 1:N)

Y[n] ~ bernoulli_logit(b[1] + b[2]*Hit_Points[n] + b[3]*Attack[n] + b[4]*Defense[n] + b[5]*Special_Attack[n] + b[6]*Special_Defense[n] + b[7]*Speed[n]);

}

■rstanでロジスティック回帰を行い、推定結果を可視化するコード

library(rstan)

N <- nrow(pokemon_data_standardized_filtered)

data <- list(N = N,
             Hit_Points = pokemon_data_standardized_filtered$Hit_Points,
             Attack = pokemon_data_standardized_filtered$Attack,
             Defense = pokemon_data_standardized_filtered$Defense,
             Special_Attack = pokemon_data_standardized_filtered$Special_Attack,
             Special_Defense = pokemon_data_standardized_filtered$Special_Defense,
             Speed = pokemon_data_standardized_filtered$Speed,
             Y = pokemon_data_standardized_filtered$y)

fit <- stan(file = 'logistic_regression.stan',
            data = data,
            seed = 1234)

summary(fit)

traceplot(fit)

source('common.R')

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7))
d_est <- data.frame(1:N_mcmc, ms$b)
colnames(d_est) <- param_names
d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])
d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')
d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()
p <- p + theme_bw(base_size=18)
p <- p + coord_flip()
p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')
p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)
p <- p + labs(x='parameter', y='value')
p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))
p

library(rstan)

N <- nrow(pokemon_data_standardized_filtered)

data <- list(N = N,

Hit_Points = pokemon_data_standardized_filtered$Hit_Points,

Attack = pokemon_data_standardized_filtered$Attack,

Defense = pokemon_data_standardized_filtered$Defense,

Special_Attack = pokemon_data_standardized_filtered$Special_Attack,

Special_Defense = pokemon_data_standardized_filtered$Special_Defense,

Speed = pokemon_data_standardized_filtered$Speed,

Y = pokemon_data_standardized_filtered$y)

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

data = data,

seed = 1234)

summary(fit)

traceplot(fit)

source('common.R')

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7))

d_est <- data.frame(1:N_mcmc, ms$b)

colnames(d_est) <- param_names

d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])

d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')

d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()

p <- p + theme_bw(base_size=18)

p <- p + coord_flip()

p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')

p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)

p <- p + labs(x='parameter', y='value')

p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))

■MCMCのシミュレーション結果のトレースプロット

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

■ロジスティック回帰の推定結果

見にくいので、推定結果を松浦さんの「StanとRでベイズ統計モデリング」にあるコードを用いて可視化します。

■推定結果の可視化

どうやら、0を含まない係数について見てみると、b3（攻撃）、b5（特殊攻撃）、b6（特殊防御）が高いほど、レア度が増す傾向があるようです。珍しいポケモンは攻撃が強いという傾向があると言えるのではないでしょうか。

まとめ

rvestは簡単にスクレイピングできて便利。
ポケモンデータは色々整備されてそうで今後も分析したら面白そう。
珍しいポケモンは「攻撃」、「特殊攻撃」、「特殊防御」が高い傾向がある。
経験値が必要なのにTTVの低い、コスパの悪そうなポケモンたちがいる。

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

参考情報

階層ベイズモデルの直帰率分析への適用 with rstan

松浦さんの『StanとRでベイズ統計モデリング』の8章の階層ベイズがすごくわかりやすいなぁと思いつつも、自分の持っているデータで試していなかったので、これを機に実践してみようと思います。
やや変数を追加しているくらいで大した変更点はありませんが、題材としては当ブログのアクセスログにおける直帰率に関するデータで、どのような要素が直帰率に影響を与えるのかを分析します。

目次
・モデル概要
・前処理
・推定
・結果（非階層モデルとの比較）
・参考文献

モデル概要

モデルは8章のロジスティック回帰の階層モデルに一部変数を追加していますが、ほぼそのままです。記事ごとのパラメータやリファラーごとのパラメータを想定しています。
Nは記事数でnはそのインデックス、Cはリファラーの数でcはそのインデックス、Iはログとして残っているセッションの数でiはそのインデックスとなっています。hatebuは記事のはてぶ数、stringlineは記事の行数、holidayは休日or祝日ダミー変数、daytimeは12:00~18:00なら1をとるダミー変数、revisitedは再訪問ユーザーなら1を取るダミー変数となっています。記事ごと・リファラーごとに直帰のしやすさが違う（パラメータが従う正規分布のパラメータがそれぞれ異なる）という仮定のもとに立ったモデルとなります。

$$x[i] = b_{1} + x_{記事}[記事ID[i]] \\ + x_{リファラー}[リファラーID[i]] + x_{セッション}[i] $$

$$q[i] = inverselogit(x[i]) $$

$$Y[i] \sim Bernoulli(q[i]) $$

$$x_{記事}[n] = b_{2}hatebu + b_{3}stringline[n] \\ + b_{記事間の差}[n] $$

$$b_{記事間の差}[n] \sim Normal(0, \sigma_{記事番号}) $$

$$x_{リファラー}[c] = b_{リファラー間の差}[c] $$

$$b_{リファラー間の差}[c] \sim Normal(0, \sigma_{リファラー番号}) $$

$$x_{セッション}[i] = b_{4}holiday[i] 　+ b_{5}divice[i] \\ + b_{6}daytime[i] + b_{7}revisited[i] $$

前処理

GAのAPIからデータを取得して1セッション1記事になるようにデータを作成しています。数ヶ月で25000件ほどデータがあったのですが、計算に時間がかかるので、データ数を2400件くらいにサンプリングしています。

library(RGA)
library(tidyverse)
library(Nippon)

authorize()
prof <-list_profiles()

start_date <- "2017-04-01"
end_date   <- "2017-10-20"

accesslogdata <- get_ga(profileId = prof$id[2],
                        start.date = start_date,
                        end.date = end_date,
                        dimensions = "ga:pagePath,
                        ga:dateHourMinute,
                        ga:deviceCategory,
                        ga:userType,
                        ga:referralPath,
                        ga:fullReferrer",
                        sort = "-ga:sessions",
                        metrics = "ga:sessions,ga:bounces",
                        fetch.by = "day")

#reshape_url
accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

accesslogdata_filtered <- accesslogdata %>% filter(sessions == 1)
accesslogdata_filtered <- accesslogdata_filtered %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'t.co/'),'twitter', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'b.hatena.ne.jp/'),'b.hatena.ne.jp', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'d.hatena.ne.jp/'),'d.hatena.ne.jp', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'facebook.com/'),'facebook.com', fullReferrer))

#refferer list
refferer_cat <- c("google", "(direct)", "twitter", 
                  "yahoo", "b.hatena.ne.jp", "facebook.com",
                  "bing")

daytime_cat <- c("12","13","14","15","16","17","18")

#making weekday data and daytime data
accesslogdata_filtered <- accesslogdata_filtered %>% 
                            filter(fullReferrer %in% refferer_cat) %>% 
                            mutate(date=as.Date(format(substr(accesslogdata_filtered$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d")) %>% 
                            mutate(hourminutes = substr(accesslogdata_filtered$dateHourMinute,start = 9,stop = 10 )) %>% 
                            mutate(holiday = ifelse(is.jholiday(date), 1, 0)) %>% 
                            mutate(daytime = ifelse(hourminutes %in% daytime_cat, 1, 0)) %>% 
                            mutate(device = ifelse(deviceCategory == "desktop", 1, 0)) %>% 
                            mutate(revisited = ifelse(userType == "Returning Visitor", 1, 0))

selected_dataset <- accesslogdata_filtered %>% 
                      select(pagePath,fullReferrer,
                             device,revisited,holiday,daytime,bounces) %>% 
                      filter(!(pagePath %in% c("article_10", "article_53") )) %>% 
                      bind_cols(rand =runif(nrow(selected_dataset), min = 0, max = 1)) %>% 
                      filter(rand <= 2400/nrow(selected_dataset))

referer_cat <- selected_dataset %>% 
                  select(fullReferrer) %>% 
                  distinct() %>% 
                  mutate(referer_id = 1:n())

article_cat <- selected_dataset %>% 
                  select(pagePath) %>% 
                  distinct() %>% 
                  mutate(article_id = 1:n())

selected_dataset <- selected_dataset %>% 
                      left_join(article_cat, by="pagePath") %>% 
                      left_join(referer_cat, by="fullReferrer") %>% 
                      select(-pagePath, -fullReferrer)

article_data <- read_csv(file = "kamonohashiperry_text.csv")
article_data$url <- gsub(article_data$url,
                         pattern = "http://kamonohashiperry.com/archives/",replacement = "article_")

article_data <- article_data %>% left_join(article_cat, by=c("url"="pagePath"))
article_data <- article_data[!(is.na(article_data$article_id)),]

library(RGA)

library(tidyverse)

library(Nippon)

authorize()

prof <-list_profiles()

start_date <- "2017-04-01"

end_date <- "2017-10-20"

accesslogdata <- get_ga(profileId = prof$id[2],

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:referralPath,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:bounces",

fetch.by = "day")

#reshape_url

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))

accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

accesslogdata_filtered <- accesslogdata %>% filter(sessions == 1)

accesslogdata_filtered <- accesslogdata_filtered %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'t.co/'),'twitter', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'b.hatena.ne.jp/'),'b.hatena.ne.jp', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'d.hatena.ne.jp/'),'d.hatena.ne.jp', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'facebook.com/'),'facebook.com', fullReferrer))

#refferer list

refferer_cat <- c("google", "(direct)", "twitter",

"yahoo", "b.hatena.ne.jp", "facebook.com",

"bing")

daytime_cat <- c("12","13","14","15","16","17","18")

#making weekday data and daytime data

accesslogdata_filtered <- accesslogdata_filtered %>%

filter(fullReferrer %in% refferer_cat) %>%

mutate(date=as.Date(format(substr(accesslogdata_filtered$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d")) %>%

mutate(hourminutes = substr(accesslogdata_filtered$dateHourMinute,start = 9,stop = 10 )) %>%

mutate(holiday = ifelse(is.jholiday(date), 1, 0)) %>%

mutate(daytime = ifelse(hourminutes %in% daytime_cat, 1, 0)) %>%

mutate(device = ifelse(deviceCategory == "desktop", 1, 0)) %>%

mutate(revisited = ifelse(userType == "Returning Visitor", 1, 0))

selected_dataset <- accesslogdata_filtered %>%

select(pagePath,fullReferrer,

device,revisited,holiday,daytime,bounces) %>%

filter(!(pagePath %in% c("article_10", "article_53") )) %>%

bind_cols(rand =runif(nrow(selected_dataset), min = 0, max = 1)) %>%

filter(rand <= 2400/nrow(selected_dataset))

referer_cat <- selected_dataset %>%

select(fullReferrer) %>%

distinct() %>%

mutate(referer_id = 1:n())

article_cat <- selected_dataset %>%

select(pagePath) %>%

distinct() %>%

mutate(article_id = 1:n())

selected_dataset <- selected_dataset %>%

left_join(article_cat, by="pagePath") %>%

left_join(referer_cat, by="fullReferrer") %>%

select(-pagePath, -fullReferrer)

article_data <- read_csv(file = "kamonohashiperry_text.csv")

article_data$url <- gsub(article_data$url,

pattern = "http://kamonohashiperry.com/archives/",replacement = "article_")

article_data <- article_data %>% left_join(article_cat, by=c("url"="pagePath"))

article_data <- article_data[!(is.na(article_data$article_id)),]

推定

stanコードはこちらになります。

data {
  int N; #the number of article
  int C; #the number of referer
  int I; #the number of log
  int<lower=0> hatebu[N]; 
  real<lower=0> stringline[N];
  int<lower=1, upper=N> article_id[I];
  int<lower=1, upper=C> referer_id[I];
  real<lower=0, upper=1> holiday[I];
  real<lower=0, upper=1> device[I];
  real<lower=0, upper=1> daytime[I];
  real<lower=0, upper=1> revisited[I];
  int<lower=0, upper=1> Y[I];
}

parameters {
  real b[7];
  real b_P[N];
  real b_C[C];
  real<lower=0> s_P;
  real<lower=0> s_C;
}

transformed parameters {
  real x_P[N];
  real x_C[C];
  real x_J[I];
  real x[I];
  real q[I];
  for (n in 1:N)
    x_P[n] = b[2]*stringline[n] + b[3]*hatebu[n] + b_P[n];
  for (c in 1:C)
    x_C[c] = b_C[c];
  for (i in 1:I) {
    x_J[i] = b[4]*holiday[i] + b[5]*device[i] + b[6]*daytime[i] + b[7]*revisited[i];
    x[i] = b[1] + x_P[article_id[i]] + x_C[referer_id[i]] + x_J[i];
    q[i] = inv_logit(x[i]);
  }
}

model {
  for (n in 1:N)
    b_P[n] ~ normal(0, s_P);
  for (c in 1:C)
    b_C[c] ~ normal(0, s_C);
  for (i in 1:I)
    Y[i] ~ bernoulli(q[i]);
}

data {

int N; #the number of article

int C; #the number of referer

int I; #the number of log

int<lower=0> hatebu[N];

real<lower=0> stringline[N];

int<lower=1, upper=N> article_id[I];

int<lower=1, upper=C> referer_id[I];

real<lower=0, upper=1> holiday[I];

real<lower=0, upper=1> device[I];

real<lower=0, upper=1> daytime[I];

real<lower=0, upper=1> revisited[I];

int<lower=0, upper=1> Y[I];

}

parameters {

real b[7];

real b_P[N];

real b_C[C];

real<lower=0> s_P;

real<lower=0> s_C;

}

transformed parameters {

real x_P[N];

real x_C[C];

real x_J[I];

real x[I];

real q[I];

for (n in 1:N)

x_P[n] = b[2]*stringline[n] + b[3]*hatebu[n] + b_P[n];

for (c in 1:C)

x_C[c] = b_C[c];

for (i in 1:I) {

x_J[i] = b[4]*holiday[i] + b[5]*device[i] + b[6]*daytime[i] + b[7]*revisited[i];

x[i] = b[1] + x_P[article_id[i]] + x_C[referer_id[i]] + x_J[i];

q[i] = inv_logit(x[i]);

}

model {

for (n in 1:N)

b_P[n] ~ normal(0, s_P);

for (c in 1:C)

b_C[c] ~ normal(0, s_C);

for (i in 1:I)

Y[i] ~ bernoulli(q[i]);

}

rstanを用いたstan実行用のRコードです。ヴァイオリンプロットで主要な係数の分布を見る処理も書かれています。

library(rstan)

N <- nrow(article_cat)
C <- nrow(referer_cat)
I <- nrow(selected_dataset)

data <- list(N = N,
             C = C,
             hatebu = article_data$hatebu,
             stringline = article_data$str_lines,
             article_id = selected_dataset$article_id,
             referer_id = selected_dataset$referer_id,
             holiday = selected_dataset$holiday,
             device = selected_dataset$device,
             daytime = selected_dataset$daytime,
             revisited = selected_dataset$revisited,
             Y = selected_dataset$bounces)

fit <- stan(file = "model/model8-8_access_analysis.stan",
            data = data,
            pars = c("b", "b_P", "b_C", "s_P", "s_C", "q"),
            seed = 1234)

source('../common.R')

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7), 's_P', 's_C')
d_est <- data.frame(1:N_mcmc, ms$b, ms$s_P, ms$s_C)
colnames(d_est) <- param_names
d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])
d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')
d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()
p <- p + theme_bw(base_size=18)
p <- p + coord_flip()
p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')
p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)
p <- p + labs(x='parameter', y='value')
p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))
p

library(rstan)

N <- nrow(article_cat)

C <- nrow(referer_cat)

I <- nrow(selected_dataset)

data <- list(N = N,

C = C,

hatebu = article_data$hatebu,

stringline = article_data$str_lines,

article_id = selected_dataset$article_id,

referer_id = selected_dataset$referer_id,

holiday = selected_dataset$holiday,

device = selected_dataset$device,

daytime = selected_dataset$daytime,

revisited = selected_dataset$revisited,

Y = selected_dataset$bounces)

fit <- stan(file = "model/model8-8_access_analysis.stan",

data = data,

pars = c("b", "b_P", "b_C", "s_P", "s_C", "q"),

seed = 1234)

source('../common.R')

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7), 's_P', 's_C')

d_est <- data.frame(1:N_mcmc, ms$b, ms$s_P, ms$s_C)

colnames(d_est) <- param_names

d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])

d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')

d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()

p <- p + theme_bw(base_size=18)

p <- p + coord_flip()

p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')

p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)

p <- p + labs(x='parameter', y='value')

p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))

結果

係数を見る限りは、符号の向きが確かなのはb5(PCダミー)とb7(再訪問ユーザーダミー)なので、PCの方が直帰しにくく、再訪問ユーザーの方が直帰しにくいという傾向があると考えることができます。

教科書ではAUCを非階層モデルと比較していましたので、比較してみようと思います。
AUCの計算を行うためのコードもGithubに載っていましたのでそちらを使います。

#calculate auc
library(pROC)

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)
spec <- seq(from=0, to=1, len=201)
probs <- c(0.1, 0.5, 0.9)

auces <- numeric(N_mcmc)
m_roc <- matrix(nrow=length(spec), ncol=N_mcmc)
for (i in 1:N_mcmc) {
  roc_res <- roc(selected_dataset$bounces, ms$q[i,])
  auces[i] <- as.numeric(roc_res$auc)
  m_roc[,i] <- coords(roc_res, x=spec, input='specificity', ret='sensitivity')
}
quantile(auces, prob=probs)

#calculate auc

library(pROC)

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

spec <- seq(from=0, to=1, len=201)

probs <- c(0.1, 0.5, 0.9)

auces <- numeric(N_mcmc)

m_roc <- matrix(nrow=length(spec), ncol=N_mcmc)

for (i in 1:N_mcmc) {

roc_res <- roc(selected_dataset$bounces, ms$q[i,])

auces[i] <- as.numeric(roc_res$auc)

m_roc[,i] <- coords(roc_res, x=spec, input='specificity', ret='sensitivity')

}

quantile(auces, prob=probs)

＃ロジスティック回帰の階層ベイズ推定でのAUC
> quantile(auces, prob=probs)
      10%       50%       90% 
0.6683084 0.6770603 0.6846594


#ロジスティック回帰のベイズ推定でのAUC
> quantile(auces_non_hiral, prob=probs_non_hiral)
      10%       50%       90% 
0.5463287 0.5623697 0.5776898

＃ロジスティック回帰の階層ベイズ推定でのAUC

> quantile(auces, prob=probs)

10% 50% 90%

0.6683084 0.6770603 0.6846594

#ロジスティック回帰のベイズ推定でのAUC

> quantile(auces_non_hiral, prob=probs_non_hiral)

10% 50% 90%

0.5463287 0.5623697 0.5776898

80%が良いとされているAUCには程遠いですが、記事やリファラーごとの差を考慮しない非階層のものよりもAUCが高いと言えます。
ちなみに、教科書の例のAUCは80%ほどでした。
Webマーケのデータ分析においてロジットは汎用性が高いで、今回のコードを土台に色々と業務で試していこうと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル原著第2版

RStanで学部時代の研究を振り返ってみる

研究概要

大学時代に実験経済学で行った実験結果のデータがUSBに入っていたので、振り返って分析などをしてみたいと思います。

研究目的
　ピア効果に関して、競争相手が自分よりも秀でた人がいいのか劣った人がいいのかを確かめる。

実験方法
・1分間で100マス計算を2セット解いてもらう。（めちゃ速い人には3枚目も渡した）
・実験開始後、実験対象のクラスによって、途中で「平均的なクラスは○○マスまで進んでいます！」とアナウンスします。アナウンスすることで、競争相手のレベルを知り、焦るなり余裕を感じるなりしてもらおうという計画です。
なお、対照群はアナウンスをしていません。アナウンス内容は「平均告知（１８秒）、上告知（１５秒）、超上告知（１２秒）、下告知（２０秒）」と4パターンとなります。
・計算が間違っているものは加点しません。

実験対象
　某国立大学の経済学部生の1~2年の必修科目履修者217名（先生に交渉して授業の開始5分を頂いて実験を行いました。）
　内部進学やスポ専などがない分、計算能力的にある程度近い集団ではないかと思われます。

検証方法
　アナウンスごとに100マス計算の点数の水準が変わりうるのかを回帰分析などで判断。

データ可視化

以下、実験カテゴリごとの略記です。
下告知（２０秒）・・・slow_20
上告知（１５秒）・・・fast_15
超上告知（１２秒）・・・fastest_12
平均告知（１８秒）・・・average_18
対照群・・・baseline

データ構造の確認です。

> str(dataset)
'data.frame':	217 obs. of  4 variables:
 $ categories    : chr  "baseline" "baseline" "baseline" "baseline" ...
 $ points        : int  72 79 81 81 98 99 100 101 102 104 ...
 $ errors        : int  0 0 0 0 2 1 0 0 1 0 ...
 $ genuine_points: int  72 79 81 81 96 98 100 101 101 104 ...

> str(dataset)

'data.frame': 217 obs. of 4 variables:

$ categories : chr "baseline" "baseline" "baseline" "baseline" ...

$ points : int 72 79 81 81 98 99 100 101 102 104 ...

$ errors : int 0 0 0 0 2 1 0 0 1 0 ...

$ genuine_points: int 72 79 81 81 96 98 100 101 101 104 ...

平均値、中央値、標準偏差、サンプルサイズを出してみます。

> dataset %>% group_by(categories) %>% summarise(average=mean(genuine_points),
+                                                median=median(genuine_points),
+                                                stdev=sd(genuine_points),
+                                                sample=n())
# A tibble: 5 × 5
  categories  average median    stdev sample
       <chr>    <dbl>  <dbl>    <dbl>  <int>
1 average_18 120.8605    117 28.08455     43
2   baseline 121.0000    121 22.19109     46
3    fast_15 123.0357    126 24.32355     56
4 fastest_12 123.6774    125 24.09756     31
5    slow_20 126.3902    126 23.20547     41

> dataset %>% group_by(categories) %>% summarise(average=mean(genuine_points),

+ median=median(genuine_points),

+ stdev=sd(genuine_points),

+ sample=n())

# A tibble: 5 × 5

categories average median stdev sample

1 average_18 120.8605 117 28.08455 43

2 baseline 121.0000 121 22.19109 46

3 fast_15 123.0357 126 24.32355 56

4 fastest_12 123.6774 125 24.09756 31

5 slow_20 126.3902 126 23.20547 41

中央値で見てみると、baselineに対してわずかですが点数に違いがありそうに見えます。

実験種別で点数に関するヒストグラムと確率密度関数を確認してみます。

library(ggplot2)

g <- ggplot(data = dataset,
            aes(x = genuine_points,
                y = ..density..)) +
            geom_histogram(alpha = 0.5,position = "identity") +
            geom_density(alpha = 0)
g + facet_wrap(~categories,nrow=5)

library(ggplot2)

g <- ggplot(data = dataset,

aes(x = genuine_points,

y = ..density..)) +

geom_histogram(alpha = 0.5,position = "identity") +

geom_density(alpha = 0)

g + facet_wrap(~categories,nrow=5)

baselineが多峰性がありそうなのが気になります。average_18は低そうに見えますね。

RStanで重回帰

『StanとRでベイズ統計モデリング』にあるコードを参考にしています。正規分布を事前分布にした線形回帰モデルです。
被説明変数が点数、説明変数が実験種別のダミー変数だけからなります。

library(rstan)
library(dummies)

dummies <- dummy.data.frame(dataset, sep = "_", names = c("categories"))

analytical_dataset <- dummies %>% select(categories_average_18,
                                         categories_fast_15,
                                         categories_fastest_12,
                                         categories_slow_20,
                                         genuine_points)


data <- list(N=nrow(analytical_dataset),
             average_18=analytical_dataset$categories_average_18,
             fast_15=analytical_dataset$categories_fast_15,
             fastest_12=analytical_dataset$categories_fastest_12,
             slow_20=analytical_dataset$categories_slow_20,
             genuine_points=analytical_dataset$genuine_points)


stan_code <- "
data{
int N; //the number of student
int<lower=0> genuine_points[N];
real<lower=0, upper=1> average_18[N];
real<lower=0, upper=1> fast_15[N];
real<lower=0, upper=1> fastest_12[N];
real<lower=0, upper=1> slow_20[N];
}

parameters{
real b1;
real b2;
real b3;
real b4;
real b5;
real<lower=0> sigma;
}

transformed parameters{
real mu[N];
for(n in 1:N)
mu[n] = b1 + b2*average_18[n] + b3*fast_15[n] + b4*fastest_12[n] + b5*slow_20[n];
}

model{
for(n in 1:N)
genuine_points[n] ~ normal(mu[n], sigma);
}
"

fit <- stan(model_code =stan_code, data=data, seed=1234)
fit.summary <-data.frame(summary(fit)$summary)
head(fit.summary,6)

library(rstan)

library(dummies)

dummies <- dummy.data.frame(dataset, sep = "_", names = c("categories"))

analytical_dataset <- dummies %>% select(categories_average_18,

categories_fast_15,

categories_fastest_12,

categories_slow_20,

genuine_points)

data <- list(N=nrow(analytical_dataset),

average_18=analytical_dataset$categories_average_18,

fast_15=analytical_dataset$categories_fast_15,

fastest_12=analytical_dataset$categories_fastest_12,

slow_20=analytical_dataset$categories_slow_20,

genuine_points=analytical_dataset$genuine_points)

stan_code <- "

data{

int N; //the number of student

int<lower=0> genuine_points[N];

real<lower=0, upper=1> average_18[N];

real<lower=0, upper=1> fast_15[N];

real<lower=0, upper=1> fastest_12[N];

real<lower=0, upper=1> slow_20[N];

}

parameters{

real b1;

real b2;

real b3;

real b4;

real b5;

real<lower=0> sigma;

}

transformed parameters{

real mu[N];

for(n in 1:N)

mu[n] = b1 + b2*average_18[n] + b3*fast_15[n] + b4*fastest_12[n] + b5*slow_20[n];

}

model{

for(n in 1:N)

genuine_points[n] ~ normal(mu[n], sigma);

}

fit <- stan(model_code =stan_code, data=data, seed=1234)

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

head(fit.summary,6)

結果

traceplot(fit)でMCMCのサンプリング結果を確認する。

収束しているように見えます。

以下推定結果ですが、残念ながらベイズ予測区間において符号の逆転が起きていないものはなかったので、アナウンスによる効果があるとは言えないようです。ただ、slow_20の係数がおしいですね。少なくとも他の実験種別よりも、アナウンス効果があるかもしれないという考察に止まりそうです。

> head(fit.summary,6)
             mean   se_mean       sd      X2.5.       X25.        X50.       X75.     X97.5.    n_eff      Rhat
b1    121.0253094 0.1034836 3.460792 114.214104 118.793988 121.0140383 123.213054 127.902759 1118.429 0.9998152
b2     -0.1195318 0.1326887 4.957153 -10.009827  -3.372795  -0.1304577   3.117195   9.694184 1395.715 1.0004034
b3      1.9440496 0.1248433 4.661571  -7.494857  -1.043060   1.9756929   5.007503  10.827640 1394.228 1.0000830
b4      2.5460694 0.1436941 5.447930  -8.026616  -1.128587   2.4997953   6.220891  13.044467 1437.426 1.0008688
b5      5.3029965 0.1351662 4.988323  -4.479865   2.004280   5.3108476   8.543292  15.458379 1361.988 1.0004238
sigma  24.6158717 0.0211252 1.150481  22.501612  23.811858  24.5900106  25.368057  26.970197 2965.905 1.0017414

> head(fit.summary,6)

mean se_mean sd X2.5. X25. X50. X75. X97.5. n_eff Rhat

b1 121.0253094 0.1034836 3.460792 114.214104 118.793988 121.0140383 123.213054 127.902759 1118.429 0.9998152

b2 -0.1195318 0.1326887 4.957153 -10.009827 -3.372795 -0.1304577 3.117195 9.694184 1395.715 1.0004034

b3 1.9440496 0.1248433 4.661571 -7.494857 -1.043060 1.9756929 5.007503 10.827640 1394.228 1.0000830

b4 2.5460694 0.1436941 5.447930 -8.026616 -1.128587 2.4997953 6.220891 13.044467 1437.426 1.0008688

b5 5.3029965 0.1351662 4.988323 -4.479865 2.004280 5.3108476 8.543292 15.458379 1361.988 1.0004238

sigma 24.6158717 0.0211252 1.150481 22.501612 23.811858 24.5900106 25.368057 26.970197 2965.905 1.0017414

分布でも見てみます。アナウンス効果が0を確実に超えているとは言えないですね。

library (reshape)
library(dplyr)
library(ggplot2)
post   <- extract (fit, permuted = F)
m.post <- melt (post)
m.post <- m.post %>% filter(parameters %in% c("b1","b2","b3","b4","b5"))
graph  <- ggplot (m.post, aes(x = value))
graph  <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw() 
plot (graph)

library (reshape)

library(dplyr)

library(ggplot2)

post <- extract (fit, permuted = F)

m.post <- melt (post)

m.post <- m.post %>% filter(parameters %in% c("b1","b2","b3","b4","b5"))

graph <- ggplot (m.post, aes(x = value))

graph <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw()

plot (graph)

結局、学部時代のレポートと結論は変わらないのですが、係数が0よりも大きい確率という観点で結果に向き合えたのは良かったと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)

大学におけるStanの講義資料などを探してみた

ブログよりは大学の講義の方が体系立てて学べるのではないか、効率的に勉強できるのではないかと思い、各大学が公開しているStanに関するサイトを調べてみました。調べ方は非常に簡単で、Google検索で「stan site:大学のドメイン名」でヒットした上位を基本的に見ています。

東京大学、京都大学、東北大学、大阪大学、慶應大学、早稲田大学、名古屋大学、同志社大学、etc…などを見ましたがなかなかweb上で公開されている資料が見つかりませんでした。非公開か大学ドメイン下での公開がされていないのだと思われます。そのため、検索対象を海外にまで広げてみました。（次回は教員のwebサイトを漁ってみようと思います。）

結論として、学ぶのにちょうどよいと思えるのは、神戸大学の資料とStanford大学の資料でした。今後はこの二つの資料も使って学習を進めていこうと思います。

神戸大学

政治学方法論 II (Research Methods in Political Science II)
ベイズ統計学の授業の内容が公開されています。教科書は「Bayesian Data Analysis, 3rd Edition. CRC Press.」です。ちなみに、こちらは無料のPDFが公開されています。（Bayesian Data Analysis, Third Edition(PDF)）

階層モデルとStan によるベイズ推定
階層ベイズモデルの説明とstanのコードが記されており、学習が捗ります。

講義のスライドはこちらにあります。（ yukiyanai/rm2-Bayes ）

東京工業大学

勉強用(STAN)
stanのコードが載っていました。
正規分布、線形回帰モデル、混合正規分布、ニューラルネットワーク、多種粒子Totally Asymmetric Simple Exclusion Process、混合正規分布でのクラスタリング、ロジスティック回帰などのコードがあるようです。

こちらはstanの説明用の資料です。（Stanによるハミルトニアンモンテカルロ法を用いたサンプリングについて）

Stanford University

Statistical Rethinking A Bayesian Course with Examples in R and Stan
youtubeで2015年の講義が見れるようです。( Statistical Rethinking Winter 2015 )
講義のスライドも公開されています。( Talks by Richard McElreath )
ゴーレムをモデルの引き合いに出して紹介しているのを見て、ユーモアセンスあるなぁと思いました。2016年版の資料も今後アップされると思うので、見逃せないですね。

Colombia University

Home page for the book, “Bayesian Data Analysis”
stanの開発チームの方がコロンビア大学の研究者なので、絶対にあるだろうと思いましたが、スライドとかは特にありませんでした。学生の講義ノートは筆記体で画像になっているので、あまり読むことはお勧めはしません。
私として嬉しいのは。Rstanで教科書のコードを実行するためのスクリプトがGitHubで公開されていることでした。（ avehtari/BDA_R_demos/demos_rstan/ ）

RstanでCVRの前後比較をするためのコード

目的

データサイエンス界隈の方がP値での意思決定に警鐘を鳴らしている昨今、施策実施に関するCVRの前後比較をχ2乗検定のP値を用いるのではなく、ベイズ統計学によるアプローチにチャレンジしてみたいと思いました。『基礎からのベイズ統計学』の8章で取り上げられていた比率データに対してのベイズ統計学的アプローチをもとに、stanを用いて事後分布から意思決定をするための進め方を紹介します。

進め方

・データの整形
・stanコード作成
・rstanでの引数の指定
・rでの可視化

データの特徴

Webマーケティング界隈では大変に多用するデータだと思いますが、実験を行ったユーザーに対しての開封・非開封、これまで通りのユーザーの開封・非開封の自然数からなるデータです。

stanコード

stanコードは
・dataブロック
・parametersブロック
・transformed parametersブロック（今回は不使用）
・modelブロック
・generated quantitiesブロック
からなります。
今回は自然数のデータであることから、ディリクレ分布を事前分布に設定するために、parametersブロックにおいてsimplexを指定しています。（教科書の比率データのものをそのまま使っています。）
modelは二項データしか出てこないので、二項分布を用いています。generated quantitiesブロックでは各々の比率、比率の差、比率の差が0を超える確率・0.01を超える確率、リスク比、リスク比が1を超える確率、オッズ比などを出力するようにしています。

data{
	int<lower=0> N[2]; #自然数からなるベクトルNの指定
	int n[2,2]; #整数からなる行列nの指定
}

parameters{
	simplex[2] p[2]; #ディリクレ分布を事前分布に設定したもとでの確率pを指定
}

model{
	for(i in 1:2){
		for(j in 1:2){
			n[i,j] ~ binomial(N[j], p[j][i]); #二項分布
		}
	}
}
generated quantities{
	real p11;
	real p10;
	real p01;
	real p00;
	real d;
	real delta_over;
	real delta_over_onep;
	real RR;
	real RRover;
	real OR;
	p11 <- p[1][1];
	p10 <- p[1][2];
	p01 <- p[2][1];
	p00 <- p[2][2];
	d <- p11 - p01; #比率の差
	delta_over <- if_else(d > 0,1,0);
	delta_over_onep <- if_else(d > 0.01,1,0);
	RR <- p11/p01; #リスク比
	RRover <- if_else(RR > 2,1,0);
	OR <- (p11/p10) / (p01/p00); #オッズ比
}

data{

int<lower=0> N[2]; #自然数からなるベクトルNの指定

int n[2,2]; #整数からなる行列nの指定

}

parameters{

simplex[2] p[2]; #ディリクレ分布を事前分布に設定したもとでの確率pを指定

}

model{

for(i in 1:2){

for(j in 1:2){

n[i,j] ~ binomial(N[j], p[j][i]); #二項分布

}

generated quantities{

real p11;

real p10;

real p01;

real p00;

real d;

real delta_over;

real delta_over_onep;

real RR;

real RRover;

real OR;

p11 <- p[1][1];

p10 <- p[1][2];

p01 <- p[2][1];

p00 <- p[2][2];

d <- p11 - p01; #比率の差

delta_over <- if_else(d > 0,1,0);

delta_over_onep <- if_else(d > 0.01,1,0);

RR <- p11/p01; #リスク比

RRover <- if_else(RR > 2,1,0);

OR <- (p11/p10) / (p01/p00); #オッズ比

}

rコード

以下は、stanをrで実行し、ggplot2などで可視化するためのコードが記されています。

library(rstan)
library (reshape)
library (ggplot2)

scr <- "model871.stan" #stanコード名

#カウントデータを用意します。
N <- c(1123, 1200)
n <- structure(.Data = c(106, 1017, 46, 1154), .Dim = c(2, 2))

#stanで用いるデータの型に変換
data <-list(N=N, n=n)

#パラメータの設定
par<-c("p","d","delta_over","delta_over_onep","RR","RRover","OR")
#pは各々の比率
#dは比率の差
#delta_overは比率の差が0を超える確率
#delta_over_onepは比率の差が0.1を超える確率
#RRはリスク比
#RRoverはRRが1を超える確率
#ORはオッズ比

war<-1000               #バーンインの期間を指定しています。
ite<-11000              #試行回数をしていしています。
see<-12345              #乱数の種
dig<-3                  #有効数字
cha<-1                  #連鎖構成数

#stanの実行
fit <- stan(file = scr, data = data, warm=war, iter=ite, seed=see,
	      pars=par,chains=cha)

#結果の出力
print(fit,pars=par,digits_summary=dig)

#事後分布の可視化
post   <- extract (fit, permuted = F)
m.post <- melt (post)
graph  <- ggplot (m.post, aes(x = value))
graph  <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw()
plot (graph)

library(rstan)

library (reshape)

library (ggplot2)

scr <- "model871.stan" #stanコード名

#カウントデータを用意します。

N <- c(1123, 1200)

n <- structure(.Data = c(106, 1017, 46, 1154), .Dim = c(2, 2))

#stanで用いるデータの型に変換

data <-list(N=N, n=n)

#パラメータの設定

par<-c("p","d","delta_over","delta_over_onep","RR","RRover","OR")

#pは各々の比率

#dは比率の差

#delta_overは比率の差が0を超える確率

#delta_over_onepは比率の差が0.1を超える確率

#RRはリスク比

#RRoverはRRが1を超える確率

#ORはオッズ比

war<-1000 #バーンインの期間を指定しています。

ite<-11000 #試行回数をしていしています。

see<-12345 #乱数の種

dig<-3 #有効数字

cha<-1 #連鎖構成数

#stanの実行

fit <- stan(file = scr, data = data, warm=war, iter=ite, seed=see,

pars=par,chains=cha)

#結果の出力

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

#事後分布の可視化

post <- extract (fit, permuted = F)

m.post <- melt (post)

graph <- ggplot (m.post, aes(x = value))

graph <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw()

plot (graph)

推定結果&可視化

今回の例では、実験を行ったユーザーのCVRの差が0以上の確率（delta_over）が1.0なので、ほぼ確実に差があると言えそうです。0.01以上差がある確率も1.0なので1%以上は差があると言えそうです。リスク比（RR）に関しては2.47と実験しない場合と比べて2.47倍程度CVを高めています。オッズ比（OR）は2.63とあるので、実験によるCV増大効果が2.63倍あると考えることができます。χ2乗検定では、二つの集団が独立かどうかを検定していますが、ベイズ統計学に従えば、「1%を超える確率」を算出することが容易なので、ディレクターなどに説明する際は圧倒的に理解を得られそうな気がします。

> print(fit,pars=par,digits_summary=dig)
Inference for Stan model: model871.
1 chains, each with iter=11000; warmup=1000; thin=1; 
post-warmup draws per chain=10000, total post-warmup draws=10000.

                 mean se_mean    sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
p[1,1]          0.095   0.000 0.006 0.083 0.091 0.095 0.099 0.107  5645    1
p[1,2]          0.905   0.000 0.006 0.893 0.901 0.905 0.909 0.917  5645    1
p[2,1]          0.039   0.000 0.004 0.031 0.036 0.039 0.041 0.047  5488    1
p[2,2]          0.961   0.000 0.004 0.953 0.959 0.961 0.964 0.969  5488    1
d               0.056   0.000 0.007 0.042 0.051 0.056 0.061 0.071  5662    1
delta_over      1.000   0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000  NaN
delta_over_onep 1.000   0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000  NaN
RR              2.477   0.004 0.303 1.933 2.267 2.453 2.662 3.136  5454    1
RRover          0.955   0.003 0.207 0.000 1.000 1.000 1.000 1.000  5749    1
OR              2.633   0.005 0.341 2.026 2.396 2.606 2.842 3.377  5463    1

Samples were drawn using NUTS(diag_e) at Sun Mar 13 23:19:36 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

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

Inference for Stan model: model871.

1 chains, each with iter=11000; warmup=1000; thin=1;

post-warmup draws per chain=10000, total post-warmup draws=10000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

p[1,1] 0.095 0.000 0.006 0.083 0.091 0.095 0.099 0.107 5645 1

p[1,2] 0.905 0.000 0.006 0.893 0.901 0.905 0.909 0.917 5645 1

p[2,1] 0.039 0.000 0.004 0.031 0.036 0.039 0.041 0.047 5488 1

p[2,2] 0.961 0.000 0.004 0.953 0.959 0.961 0.964 0.969 5488 1

d 0.056 0.000 0.007 0.042 0.051 0.056 0.061 0.071 5662 1

delta_over 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000 NaN

delta_over_onep 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000 NaN

RR 2.477 0.004 0.303 1.933 2.267 2.453 2.662 3.136 5454 1

RRover 0.955 0.003 0.207 0.000 1.000 1.000 1.000 1.000 5749 1

OR 2.633 0.005 0.341 2.026 2.396 2.606 2.842 3.377 5463 1

Samples were drawn using NUTS(diag_e) at Sun Mar 13 23:19:36 2016.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

参考文献

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門
 rstanでちょこちょこ

RstanでGoogleトレンドの推定

rstanを用いて、Googleトレンドデータの予測モデルを推定してみます。

ほとんど岩波データサイエンスのものですが、Googleトレンドのデータを月ごとの季節性を加味した状態空間モデルを用いて予測してみました。

今回の分析では、
・modelのstanコード(stan)
・Rでstanを動かすためのコード(R)
（・可視化のためのコード(R)）←必須ではない
を用意します。

データですが、GoogleTrendのサイトで任意のキーワードで検索して、
その時系列データをCSVでダウンロードすれば手に入ります。（ちょっと見つけにくい）

データの形式はシンプルで、
先頭にY
とおいて後はトレンドの値を行ごとに置いていけばいけます。

要はN行1列データをテキストファイルに保存すればOKです。（1行目はY）

まずstanのコードですが、岩波データサイエンスのサンプルコードの季節を4から12に変えています。（たったこれだけ）
Googleトレンドのデータは月単位でも結構値がふれることがあるので、月ごとに応じた潜在的な変数が必要だと思いました。

data {
   int<lower=1> T;
   int<lower=1> T_next;
   real Y[T];
}

parameters {
   real mu[T];
   real s[T];
   real<lower=0> s_mu;
   real<lower=0> s_s;
   real<lower=0> s_r;
}

model {
   for(t in 2:T)
      mu[t] ~ normal(mu[t-1], s_mu);
   for(t in 12:T)
      s[t] ~ normal(-s[t-1]-s[t-2]-s[t-3]-s[t-4]-s[t-5]-s[t-6]-s[t-7]-s[t-8]-s[t-9]-s[t-10]-s[t-11], s_s);
   for(t in 1:T)
      Y[t] ~ normal(mu[t]+s[t], s_r);
}

generated quantities {
   real mu_all[T+T_next];
   real s_all[T+T_next];
   real y_next[T_next];

   for (t in 1:T){
      mu_all[t] <- mu[t];
      s_all[t] <- s[t];
   }
   for (t in (T+1):(T+T_next)){
      mu_all[t] <- normal_rng(mu_all[t-1], s_mu);
      s_all[t] <- normal_rng(-s_all[t-1]-s_all[t-2]-s_all[t-3]-s_all[t-4]-s_all[t-5]-s_all[t-6]-s_all[t-7]-s_all[t-8]-s_all[t-9]-s_all[t-10]-s_all[t-11], s_s);
   }
   for (t in 1:T_next)
      y_next[t] <- normal_rng(mu_all[T+t]+s_all[T+t], s_r);
}

data {

int<lower=1> T;

int<lower=1> T_next;

real Y[T];

}

parameters {

real mu[T];

real s[T];

real<lower=0> s_mu;

real<lower=0> s_s;

real<lower=0> s_r;

}

model {

for(t in 2:T)

mu[t] ~ normal(mu[t-1], s_mu);

for(t in 12:T)

s[t] ~ normal(-s[t-1]-s[t-2]-s[t-3]-s[t-4]-s[t-5]-s[t-6]-s[t-7]-s[t-8]-s[t-9]-s[t-10]-s[t-11], s_s);

for(t in 1:T)

Y[t] ~ normal(mu[t]+s[t], s_r);

}

generated quantities {

real mu_all[T+T_next];

real s_all[T+T_next];

real y_next[T_next];

for (t in 1:T){

mu_all[t] <- mu[t];

s_all[t] <- s[t];

}

for (t in (T+1):(T+T_next)){

mu_all[t] <- normal_rng(mu_all[t-1], s_mu);

s_all[t] <- normal_rng(-s_all[t-1]-s_all[t-2]-s_all[t-3]-s_all[t-4]-s_all[t-5]-s_all[t-6]-s_all[t-7]-s_all[t-8]-s_all[t-9]-s_all[t-10]-s_all[t-11], s_s);

}

for (t in 1:T_next)

y_next[t] <- normal_rng(mu_all[T+t]+s_all[T+t], s_r);

}

Rでstanを動かすためのコードですが、ここはサンプルコードとほぼ一緒です。

library(rstan)
setwd("任意")

d <- read.csv('data-trend4.txt', header=TRUE)
T <- nrow(d)
T_next <- 8
data <- list(T=T, T_next=T_next, Y=d$Y)

stanmodel <- stan_model(file='model_custom.stan')
fit <- sampling(
   stanmodel, data=data, pars=c('mu_all','s_all','y_next','s_mu','s_s','s_r'),
   iter=10200, warmup=200, thin=10, chains=3,
   seed=123
)

library(rstan)

setwd("任意")

d <- read.csv('data-trend4.txt', header=TRUE)

T <- nrow(d)

T_next <- 8

data <- list(T=T, T_next=T_next, Y=d$Y)

stanmodel <- stan_model(file='model_custom.stan')

fit <- sampling(

stanmodel, data=data, pars=c('mu_all','s_all','y_next','s_mu','s_s','s_r'),

iter=10200, warmup=200, thin=10, chains=3,

seed=123

)

可視化のためのコードについてもサンプルコードとほぼ一緒です。

library(ggplot2)

# after estimation

d_obs <- data.frame(X=1:T, Y=d$Y)
p <- ggplot()
p <- p + theme_bw() + theme(text=element_text(size=18))
p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', alpha=0.8, size=2)
p <- p + labs(x='Time [month]', y='trend')
p <- p + coord_cartesian(xlim=c(0.9, 152.1))
ggsave(file='fig2-top-left.png', plot=p, dpi=300, width=6, height=4)


makeDataFrameQuantile <- function(x, y_smp){
  qua <- apply(y_smp, 2, quantile, prob=c(0.1, 0.25, 0.5, 0.75, 0.9))
  d_est <- data.frame(X=x, t(qua))
  colnames(d_est) <- c('X', 'p10', 'p25', 'p50', 'p75', 'p90')
  return(d_est)
}

plotTimecourse <- function(file, d_est, d_obs){
  p <- ggplot()
  p <- p + theme_bw() + theme(text=element_text(size=18))
  p <- p + geom_vline(xintercept=T, linetype='dashed')
  p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p10, ymax=p90), fill='black', alpha=0.25)
  p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p25, ymax=p75), fill='black', alpha=0.5)
  p <- p + geom_line(data=d_est, aes(x=X, y=p10), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p90), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p25), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p75), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p50), color='black', size=0.4)
  if (!is.null(d_obs)){
    p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', size=2, alpha=0.9)
  }
  p <- p + labs(x='month', y='trend')
  p <- p + coord_cartesian(xlim=c(0.9, 152.1))
  ggsave(file=file, plot=p, dpi=300, width=6, height=4)
}

la <- rstan::extract(fit)
d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$mu_all)
plotTimecourse(file='fig2-bottom-left.png', d_est=d_est, d_obs=d_obs)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$s_all)
plotTimecourse(file='fig2-bottom-right.png', d_est=d_est, d_obs=NULL)

d_est <- makeDataFrameQuantile(x=(T+1):(T+T_next), y_smp=la$y_next)
d_est <- rbind(data.frame(X=T, p10=d$Y[T], p25=d$Y[T], p50=d$Y[T], p75=d$Y[T], p90=d$Y[T]), d_est)
plotTimecourse(file='fig2-top-right.png', d_est=d_est, d_obs=d_obs)

library(ggplot2)

# after estimation

d_obs <- data.frame(X=1:T, Y=d$Y)

p <- ggplot()

p <- p + theme_bw() + theme(text=element_text(size=18))

p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', alpha=0.8, size=2)

p <- p + labs(x='Time [month]', y='trend')

p <- p + coord_cartesian(xlim=c(0.9, 152.1))

ggsave(file='fig2-top-left.png', plot=p, dpi=300, width=6, height=4)

makeDataFrameQuantile <- function(x, y_smp){

qua <- apply(y_smp, 2, quantile, prob=c(0.1, 0.25, 0.5, 0.75, 0.9))

d_est <- data.frame(X=x, t(qua))

colnames(d_est) <- c('X', 'p10', 'p25', 'p50', 'p75', 'p90')

return(d_est)

}

plotTimecourse <- function(file, d_est, d_obs){

p <- ggplot()

p <- p + theme_bw() + theme(text=element_text(size=18))

p <- p + geom_vline(xintercept=T, linetype='dashed')

p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p10, ymax=p90), fill='black', alpha=0.25)

p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p25, ymax=p75), fill='black', alpha=0.5)

p <- p + geom_line(data=d_est, aes(x=X, y=p10), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p90), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p25), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p75), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p50), color='black', size=0.4)

if (!is.null(d_obs)){

p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', size=2, alpha=0.9)

}

p <- p + labs(x='month', y='trend')

p <- p + coord_cartesian(xlim=c(0.9, 152.1))

ggsave(file=file, plot=p, dpi=300, width=6, height=4)

}

la <- rstan::extract(fit)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$mu_all)

plotTimecourse(file='fig2-bottom-left.png', d_est=d_est, d_obs=d_obs)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$s_all)

plotTimecourse(file='fig2-bottom-right.png', d_est=d_est, d_obs=NULL)

d_est <- makeDataFrameQuantile(x=(T+1):(T+T_next), y_smp=la$y_next)

d_est <- rbind(data.frame(X=T, p10=d$Y[T], p25=d$Y[T], p50=d$Y[T], p75=d$Y[T], p90=d$Y[T]), d_est)

plotTimecourse(file='fig2-top-right.png', d_est=d_est, d_obs=d_obs)

以上を実行した結果、以下のような図が出てきます。

こちらは実際の時系列データのプロットです。

8期先までの予測です。

8期先までの予測範囲です。信頼区間90%までの範囲となっています。

推定した潜在的な季節性のデータをプロットしています。

ついでに、4月までのデータを用いて、5~9月の予測を行い、その比較を行っています。

5月が大きく外れましたが、その後はある程度当てれているように見えます。
5月も当てれるようなモデルを作りたいものですね。

参考文献

<br />

Rstanの参考文献（インストール・使い方・実践）

Rstanに関する情報を集めたものです。
・インストール
・使い方
・実践
について載せています。
随時更新します。

インストール関連

Windows 7にRStanをインストールする
http://cordea.hatenadiary.com/entry/2013/11/20/201209

Building R for Windows
https://cran.r-project.org/bin/windows/Rtools/index.html

【R】OSXでRStanの導入と簡単な例題【MCMC】
http://www.fisproject.jp/2015/04/rstan/

MCMCの計算にStanを使ってみた（超基礎・導入編）
http://tjo.hatenablog.com/entry/2013/11/06/201735

R stan導入公開版
http://www.slideshare.net/KojiKosugi/r-stan

使い方

RStan Getting Started
https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started

Computation in R and Stan
http://www.stat.columbia.edu/~gelman/book/software.pdf

stan Documentation
http://mc-stan.org/documentation/

Stan: A Probabilistic Programming Language
http://www.stat.columbia.edu/~gelman/research/published/stan-paper-revision-feb2015.pdf

The Stan Modeling Language
http://mlss2014.hiit.fi/mlss_files/2-stan.pdf

Stanで統計モデリングを学ぶ(7): 時系列の「トレンド」を目視ではなくきちんと統計的に推定する
http://tjo.hatenablog.com/entry/2014/12/05/190105

Bayesian linear mixed models using Stan: A
tutorial for psychologists, linguists, and cognitive
scientists
http://arxiv.org/pdf/1506.06201v1.pdf

実践

R で状態空間モデル: 状態空間時系列分析入門を {rstan} で再現したい
http://sinhrks.hatenablog.com/entry/2015/05/28/071124

書籍

StanとRでベイズ統計モデリング (Wonderful R)

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門