Bayesian Statistics and Marketing – 混合ガウス×階層モデルのマーガリン購買データへの適用

前回の分析では、価格への反応係数の事前分布が正規分布を仮定したモデルを用いていましたが、事後分布から多峰性が観察されました。そこで今回は、各個人の価格への反応係数の事前分布が混合ガウス分布に従うとした場合の事例を扱いたいと思います。

データのおさらい

データ自体は前回のブログと同じですが、先日のTokyo.Rで松浦さんがオススメしていたGGallyパッケージのggpairs関数を用いて、今回扱うデータを可視化してみます。

まず、購買したマーガリンのブランド選択(6ブランド分)と、購入した価格（ドル）からなるデータセットの可視化をすると以下のようになります。

1つ目のブランドが最も多く選択されているようです。購入価格の分布は、ブランドによって多峰性がありそうです。

続いて、家計ごとの属性（家族構成、学歴、職位、退職の有無など）からなるデータセットの可視化をすると以下のようになります。

ほとんどダミー変数なので面白みには欠いていますね。ホワイトカラーの家計が多く、退職していない家計が多く、学歴が低い人が多いようです。年収に関しては対数正規分布に従ってそうです。

モデル

前回と同様に、「Bayesian Statistics and Marketing」の5章に載っている混合ガウスを想定した階層ベイズモデルを扱います。

yi ∼ MNL(Xi, βi) は購買レコードごとの意思決定が多項ロジスティック回帰モデルに従うということを意味し、βiは(説明変数の数×1)のベクトルとなります。
βi = Z∆[i,] + ui は購買レコードごとの価格にかかってくる係数で、その係数が家計の属性データに係数Δ(家計の数×属性データからなる説明変数の数)をかけ合わせたものと潜在的な項の和となります。なお、ここでの属性データには定数項を含めていません。係数Δは平均deltabar、分散A_delta^-1の多変量正規分布に従います。一方で、潜在的な項は平均µind、分散Σindの多変量正規分布に従います。この平均や分散に振られているindが、混合正規分布のパラメータの通し番号となり、多項分布に従います。この多項分布の割当確率がディリクレ分布に従います。最後に混合正規分布の各パラメータは、平均は正規分布に分散は逆ウィシャート分布に従います。

前回との大きな違いは、価格の係数の一部である、潜在的な項において、混合正規分布が仮定されているところになります。
ちなみに、モデルに関しての詳細はbayesmパッケージのマニュアルの61ページ目に記載されていました。

今回のモデルのDAG

Pythonのdaftというモジュールを使うことで、非常に簡単に今回のモデルのDAG(有向非巡回グラフ)を描くことができます。

今回はこちらのPythonコードで描けました。

import daft
from matplotlib import rc

pgm = daft.PGM(shape=[6,6])

# Nodes
pgm.add_node(daft.Node("pvec", r"pvec",1, 2)) # 名前　ラベル　座標(横、縦)
pgm.add_node(daft.Node("ind", r"ind",2, 2))
pgm.add_node(daft.Node("sigma", r"$\Sigma_k$",3, 1))
pgm.add_node(daft.Node("mu", r"$\mu_k$",3, 3))
pgm.add_node(daft.Node("theta", r"$\theta_i$",4, 2))
pgm.add_node(daft.Node("delta", r"$\Delta$",4, 0.5))
pgm.add_node(daft.Node("y", r"$y_i$",5, 2,observed=True))

# Edges
pgm.add_edge("pvec", "ind")
pgm.add_edge("ind", "sigma")
pgm.add_edge("ind", "mu")
pgm.add_edge("sigma", "mu")
pgm.add_edge("sigma", "theta")
pgm.add_edge("mu", "theta")
pgm.add_edge("delta", "theta")
pgm.add_edge("theta", "y")

pgm.render()
pgm.figure.savefig("Hierachical_Models_ForMixtureNormal.png")

import daft

from matplotlib import rc

pgm = daft.PGM(shape=[6,6])

# Nodes

pgm.add_node(daft.Node("pvec", r"pvec",1, 2)) # 名前　ラベル　座標(横、縦)

pgm.add_node(daft.Node("ind", r"ind",2, 2))

pgm.add_node(daft.Node("sigma", r"$\Sigma_k$",3, 1))

pgm.add_node(daft.Node("mu", r"$\mu_k$",3, 3))

pgm.add_node(daft.Node("theta", r"$\theta_i$",4, 2))

pgm.add_node(daft.Node("delta", r"$\Delta$",4, 0.5))

pgm.add_node(daft.Node("y", r"$y_i$",5, 2,observed=True))

# Edges

pgm.add_edge("pvec", "ind")

pgm.add_edge("ind", "sigma")

pgm.add_edge("ind", "mu")

pgm.add_edge("sigma", "mu")

pgm.add_edge("sigma", "theta")

pgm.add_edge("mu", "theta")

pgm.add_edge("delta", "theta")

pgm.add_edge("theta", "y")

pgm.render()

pgm.figure.savefig("Hierachical_Models_ForMixtureNormal.png")

stanコード

今回扱ったstanコードとなります。誤りがある場合はお知らせしていただけると幸いです。

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

  int K; // 混合ガウス分布の要素数
}

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

parameters{
  vector[p_x] theta_ast[N_z]; // 説明変数の数だけある、購買ごとのパラメータ
  matrix[p_z, p_x] Delta; // 属性データの説明変数の数×購買データの説明変数の数だけのパラメータ
  vector[p_x] u[N_z]; // 購買レコードごとのパラメータ
  vector[p_x] mu[K]; // 混合分布を構成する平均値
  cov_matrix[p_x] Sigma[K]; // 共分散行列
  simplex[K] pi[N_z]; // シンプレックス（各要素が[0,1]の範囲で合計が1という条件を満たす。）
}

transformed parameters{
  vector[p_x] theta[N_z]; #家計の数だけの係数ベクトル

  for(i in 1:N_z){ // 家計の数だけ繰り返す
    theta[i] = theta_ast[i] + Delta' * Z[i]'; // 係数は家計属性ごとの特徴に異質なDeltaとbeta_astの和で決まる 
  }
}

model{
  real ps[K]; // 混合正規分布の対数

  for(i in 1:N_x){ // 購買レコードの数だけ繰り返す
    y[i] ~ categorical(softmax(theta[hhid[i]] .* to_vector(X[i]))); //カテゴリカル分布にsoftmaxを組み合わせて多項ロジスティック回帰を行う

    pi[hhid[i]] ~ dirichlet(rep_vector(2, K));
    for(k in 1:K){ // 混合分布の構成要素の数だけ繰り返す
      ps[k] = log(pi[hhid[i]][k]) + multi_normal_lpdf(theta_ast[hhid[i]] | mu[k], Sigma[k]);
    }
    target += log_sum_exp(ps); // 離散パラメータを消去した形で対数尤度を表現する際に必要な計算（周辺化消去）。  
  }

  for(k in 1:K){ // 混合パラメータの数だけ繰り返す
    mu[k] ~ multi_normal(rep_vector(0, p_x), 100*Sigma[k]);
    Sigma[k] ~ inv_wishart(nu, nu*I); 
  }

  for(i in 1:p_z){ // 家計の属性データの数だけ繰り返す
    Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), 100*J); 
  }
}

data{

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

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

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

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

int y[N_x]; // 選択肢

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

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

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

int K; // 混合ガウス分布の要素数

}

transformed data{

real nu;

matrix[p_x, p_x] I; // 購買レコードの説明変数の数の正方行列

matrix[p_x, p_x] J; // 属性データの説明変数の数の正方行列

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

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

J = diag_matrix(rep_vector(1, p_x));

}

parameters{

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

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

vector[p_x] u[N_z]; // 購買レコードごとのパラメータ

vector[p_x] mu[K]; // 混合分布を構成する平均値

cov_matrix[p_x] Sigma[K]; // 共分散行列

simplex[K] pi[N_z]; // シンプレックス（各要素が[0,1]の範囲で合計が1という条件を満たす。）

}

transformed parameters{

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

for(i in 1:N_z){ // 家計の数だけ繰り返す

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

}

model{

real ps[K]; // 混合正規分布の対数

for(i in 1:N_x){ // 購買レコードの数だけ繰り返す

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

pi[hhid[i]] ~ dirichlet(rep_vector(2, K));

for(k in 1:K){ // 混合分布の構成要素の数だけ繰り返す

ps[k] = log(pi[hhid[i]][k]) + multi_normal_lpdf(theta_ast[hhid[i]] | mu[k], Sigma[k]);

}

target += log_sum_exp(ps); // 離散パラメータを消去した形で対数尤度を表現する際に必要な計算（周辺化消去）。

}

for(k in 1:K){ // 混合パラメータの数だけ繰り返す

mu[k] ~ multi_normal(rep_vector(0, p_x), 100*Sigma[k]);

Sigma[k] ~ inv_wishart(nu, nu*I);

}

for(i in 1:p_z){ // 家計の属性データの数だけ繰り返す

Delta[i] ~ multi_normal_cholesky(rep_vector(0, p_x), 100*J);

}

stanをキックするためのコードです。先人が書かれた混合ガウスのスクリプトをHMCで実行した際に26時間ほどかかったので、今回はより複雑なモデルであることから、変分ベイズ法による推定を行ってみることにしました。松浦さんの教科書にあるように、vb関数を用いて変分ベイズ推論を行っています。

library(bayesm)
library(tidyverse)
library(rstan)
library(GGally)

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

data("margarine")

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

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

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

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

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

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

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

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

#可視化
ggpairs(choicePrice.selected %>% select(2,3,4,5,6,7,9))
ggpairs(Z)


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

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

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

#推定
stanmodel <- stan_model(file = "Hierarchical_Mixture.stan")
d.fit_vb <- vb(stanmodel, data=d.dat, seed=123)

library(bayesm)

library(tidyverse)

library(rstan)

library(GGally)

rstan_options(auto_write = TRUE)

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

data("margarine")

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

hhid_selected <- margarine$choicePrice %>%

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

group_by(hhid) %>%

summarise(purc_cnt = n()) %>%

filter(purc_cnt >= 5)

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

choicePrice.selected <- margarine$choicePrice %>%

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

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

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

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

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

#データサイズ

N <- nrow(choicePrice.selected)

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

p <- n_distinct(choicePrice.selected$choice)

#被説明変数

y <- choicePrice.selected$choice

#説明変数

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

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

Z <- demos.selected %>%

select(-hhid)

#可視化

ggpairs(choicePrice.selected %>% select(2,3,4,5,6,7,9))

ggpairs(Z)

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

hhid_index <- demos.selected %>%

select(hhid) %>%

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

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

hhid_x <- choicePrice.selected %>%

select(hhid) %>%

left_join(hhid_index)

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

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

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

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

hhid = hhid_x$ind,

K = 3)

#推定

stanmodel <- stan_model(file = "Hierarchical_Mixture.stan")

d.fit_vb <- vb(stanmodel, data=d.dat, seed=123)

実行結果

今回はK={1,3,5,10}の混合要素数で推定を行いました。

K=1のケース

K=3のケース

K=5のケース

K=10のケース

Kが小さいと散らばりが比較的小さそうに見えます。逆に、Kが大きくなると散らばりが出てくるようです。

おわりに

　２回に渡って、マーケティングデータを用いたベイズ統計モデリングを学んできましたが、数式からstanコードに落とし込む作業の際に、stanの関数をある程度知らないとやりにくいという至極当然なことを実感しました。これまでに使ってきたstanの関数は限られたものしか扱っていなかったと言えます。特に今回は離散パラメータをstanで扱うパーツがあったので、先行研究や松浦さんの本を読みながらの手探りが多かったです。
　あと、複雑な階層ベイズモデルを扱う際に、頭の中を整理しないと手が止まってしまう感じがあったので、数式と対応するコードを横に並べながら進めました。
　マーケティングにおいては、顧客の属性ごとに多峰性のあるような事例を扱うことが多く、かつ各々のサンプルサイズも期待できないことが多いので、今回の学びを分析業務で試してみたいと思います。

参考情報

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)
Package ‘bayesm’
Multivariate Gaussian Mixture Model done properly
StanとRでベイズ統計モデリング (Wonderful R)
機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)
daftでグラフィカルモデル

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版

確率的プログラミングライブラリ「Edward」まとめ

Edwardで何ができるのか知らなかったので、忘備録として残しておきます。

目次
・Edwardとは
・Edwardでできること
・参考スライド
・参考文献

Edwardとは

・LDAで有名なコロンビア大のBlei先生の研究室で、2016年より開発されている確率的プログラミング((　プログラミング言語の変数をモデルの構成要素として使うプログラミング))のPythonライブラリ。
・積み木のように明快な形で確率的モデリングを行うことができる。（モデル→推論→評価　を一括でできる。）
・ベイズ統計と機械学習、深層学習、確率的プログラミングを融合させている。
・計算の際にTensorFlowを用いている。TensorBoardを可視化の際に用いることもできる。
・計算速度がStanやPyMC3よりも速い。GPUを用いた高速化も可能。((　pip install tensorflow-gpuでGPU版のTensolFlowを入れておく必要がある。))
・統計学者のGeorge Edward Pelham Boxから名前を取っている。

Edwardでできること

一般的なベイズ推定は当然ながら、深層学習向けのベイズ適用系の事例が豊富なようです。

・ベイズ線形回帰( Supervised Learning (Regression) )
・バッチトレーニング（巨大なデータセットにおける学習で用いる）( Batch Training )
・Tensorboardを用いた可視化（ Tensorboard ）
・Automated Transformations( Automated Transformations )
・線形混合効果モデル( Linear Mixed Effects Models )
・教師あり学習による分類( Supervised Learning (Classification) )
・教師なし学習( Unsupervised Learning )
・ニューラルネットワークの潜在空間モデル( Latent Space Models for Neural Data )
・混合密度ネットワーク( Mixture Density Networks )
・GAN( Generative Adversarial Networks )
・確率的デコーダー( Probabilistic Decoder )
・ネットワークの推論( Inference Networks )
・ベイジアンニューラルネットワーク( Bayesian Neural Network )
・確率的PCA(主成分分析)( Probabilistic PCA )

jupyterのコードたちはblei-lab/edwardのnotebookに載っています。

2層のニューラルネットワークへのベイズ推定の適用(
Bayesian Deep Learning with Edward (and a trick using Dropout) – Andrew Rowan – PyData London 2017)
Gounosyの方のブログによると、CTR予測などで扱うことができるようです。

参考スライド

確率的プログラミングライブラリEdward from Yuta Kashino

ベイズ統計の数理と深層学習 @db analytics showcase Sapporo 2017 from Takeshi Yoneda

（DL hacks輪読）Bayesian Neural Network from Masahiro Suzuki

参考文献

Edward: A library for probabilistic modeling, inference, and criticism

DEEP PROBABILISTIC PROGRAMMING

EdwardでBayesian DNN＋Variational Inferenceをやってみた話

DEEP PROBABILISTIC PROGRAMMING —”深層学習＋ベイズ”のライブラリ— Edwardの紹介

【Edward】MCMCの数学的基礎からStochastic Gradient Langevin Dynamicsの実装まで

Hello, world! Stan, PyMC3, and Edward
stanの開発者の方がstanとPyMC3とEdwardを比較しています。

Pythonで体験するベイズ推論 PyMCによるMCMC入門

ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる

ベイジアンネットワークの知見が無かったので、調べた情報をまとめています。一応、載せているスクリプトでRを用いて予測するということができます。

【目次】
・ベイジアンネットワークとは
・ベイジアンネットワークの用途
・ベイジアンネットワークの推定のステップ
・Rでの実行例
・おまけ：Webサービスへの応用例
・参考文献

ベイジアンネットワークとは

・複数の確率変数の間の定性的な依存関係をグラフ構造によって表し、個々の変数の間の定量的な関係を条件付確率で表した確率モデル。前提として、有向非循環（グルグルと回らないグラフ）となっているグラフの構造を持つものに限定している。
・入力となる変数と出力となる変数はモデルの中では区別されない。
・時間という明確な因果関係などをモデルに組み込みやすいので、系列データなどを扱うケースが多い。モデル設計者がデータが生成されるプロセスを考慮しながらモデルを組んでいける。
・非循環性とd分離の仮定のみによって導かれる現在考えられる最も自然な離散モデルであり、現在の様々なモデルの中でも、最も表現力と予測力を持つモデルとされる。

ベイジアンネットワークの用途

・物体追跡
・ジェスチャ認識
・Webサイトなどのレコメンデーションサービス
・広告配信
・メルマガなどの配信最適化
など画像処理問題から消費者行動問題までいろいろな分野で活用されているようです。

ベイジアンネットワークの推定のステップ

モデルの選択や推定は学習 (learning)と呼ばれ、以下の2段階のステップを踏みます。

Step1：構造学習(structure learning)…データからネットワーク構造を学習する。
Step2：パラメータ学習(parameter learning)…Step1で学習した構造によって意味付けした分布のパラメータを学習する。

ベイジアンネットワークの推定は条件付き確率やMAP(Maximum a posteriori)推定量などを用います。構造学習の際に、事前知識を導入することができます。

Rでの実行例

・データ
何千もの主要な免疫系細胞から取り出した11種類のリン酸タンパク質とリン脂質の測定値からなるデータです。以下のサイトよりダウンロードできます。
Supporting Online Material Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data

・パッケージ
ベイジアンネットワークを実行できるパッケージはいろいろありますが、参考文献に従いbnlearnを用います。加えて、ネットワーク図の可視化のためにRgraphvizを用います。Rgraphvizがcranでは対応したものがなかったので、こちらを参考にして、インストールしました。
Provides plotting capabilities for R graph objects

・ネットワークの図示

#ネットワークの図示
library(Rgraphviz)
library(bnlearn)

spec = paste("[PKC][PKA|PKC][praf|PKC:PKA][pmek|PKC:PKA:praf]",
             "[p44.42|pmek:PKA][pakts473|p44.42:PKA][P38|PKC:PKA]",
             "[pjnk|PKC:PKA][plcg][PIP3|plcg][PIP2|plcg:PIP3]")

net = model2network(spec)
class(net)

graphviz.plot(net, shape = "ellipse")

#ネットワークの図示

library(Rgraphviz)

library(bnlearn)

spec = paste("[PKC][PKA|PKC][praf|PKC:PKA][pmek|PKC:PKA:praf]",

"[p44.42|pmek:PKA][pakts473|p44.42:PKA][P38|PKC:PKA]",

"[pjnk|PKC:PKA][plcg][PIP3|plcg][PIP2|plcg:PIP3]")

net = model2network(spec)

class(net)

graphviz.plot(net, shape = "ellipse")

・データ確認

#データ確認
> protein_dataset <-  read.csv(file = "protain_dataset.csv",as.is = TRUE)
> head(protein_dataset,5)
   praf pmek  plcg PIP2  PIP3 p44.42 pakts473  PKA   PKC  P38 pjnk
1 39.20  8.9  4.26 9.22  1.88   25.9     46.6  223  3.02 31.6 65.5
2 42.60 20.2 12.10 1.45 40.00   30.2     51.9  264  5.94 14.1 16.1
3  5.88 17.9  3.43 2.19 18.10  105.0    300.0 1928  6.98 20.5 12.9
4 34.30 18.8  3.65 1.70 16.50  352.0    519.0 1263 18.30 59.4 35.9
5 24.60 16.0  6.92 9.39  8.06   15.7     53.8  346  9.31 18.1 12.5

#データ確認

> protein_dataset <- read.csv(file = "protain_dataset.csv",as.is = TRUE)

> head(protein_dataset,5)

praf pmek plcg PIP2 PIP3 p44.42 pakts473 PKA PKC P38 pjnk

1 39.20 8.9 4.26 9.22 1.88 25.9 46.6 223 3.02 31.6 65.5

2 42.60 20.2 12.10 1.45 40.00 30.2 51.9 264 5.94 14.1 16.1

3 5.88 17.9 3.43 2.19 18.10 105.0 300.0 1928 6.98 20.5 12.9

4 34.30 18.8 3.65 1.70 16.50 352.0 519.0 1263 18.30 59.4 35.9

5 24.60 16.0 6.92 9.39 8.06 15.7 53.8 346 9.31 18.1 12.5

・構造学習
データセットから、ネットワーク構造を推定します。推定するための手法は様々あるようですが、今回はヒルクライムアルゴリズムを用います。

#構造学習の実行
> print(hc(protein_dataset))

  Bayesian network learned via Score-based methods

  model:
   [PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]
   [PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]
  nodes:                                 11 
  arcs:                                  11 
    undirected arcs:                     0 
    directed arcs:                       11 
  average markov blanket size:           2.18 
  average neighbourhood size:            2.00 
  average branching factor:              1.00 

  learning algorithm:                    Hill-Climbing 
  score:                                 BIC (Gauss.) 
  penalization coefficient:              3.348517 
  tests used in the learning procedure:  185 
  optimized:                             TRUE

#構造学習の実行

> print(hc(protein_dataset))

Bayesian network learned via Score-based methods

model:

[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]

[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]

nodes: 11

arcs: 11

undirected arcs: 0

directed arcs: 11

average markov blanket size: 2.18

average neighbourhood size: 2.00

average branching factor: 1.00

learning algorithm: Hill-Climbing

score: BIC (Gauss.)

penalization coefficient: 3.348517

tests used in the learning procedure: 185

optimized: TRUE

構造学習より得られたネットワーク構造を可視化します。因果関係がダメな感じのネットワークが出来ていないかチェックする必要がありますね。機械的に作るだけでなく、背景的知識も考慮するという進め方が推奨されています。

#推定したネットワークを可視化
spec.estimated = paste("[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]",
                        "[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]")
net.estimated = model2network(spec.estimated)
class(net.estimated)
graphviz.plot(net.estimated, shape = "ellipse")

#推定したネットワークを可視化

spec.estimated = paste("[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]",

"[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]")

net.estimated = model2network(spec.estimated)

class(net.estimated)

graphviz.plot(net.estimated, shape = "ellipse")

・パラメータ学習
先ほど、データセットから求めた構造について、以下のコードでパラメータ推定します。パラメータの推定や構造推定は試行錯誤するところなので、こんな簡単に済む訳では無いようです。

#パラメータ学習
bn.fit(学習した構造, データセット,method = "mle")

1 2	#パラメータ学習 bn.fit(学習した構造, データセット,method = "mle")

・予測
テストデータを用意して、データセットから求めた構造・パラメータを用いて、任意の変数の予測を行います。予測したい変数の値は欠損していたらエラーで予測してくれないですが、全部0にしておけば回ります。運用上は予測したい変数の値はわからないはずなので。

#ベイジアンネットワークによる予測
library(forecast)
training.set = protein_dataset[1:405, ] #パラメータ推定用の訓練用データ
test.set = protein_dataset[406:810, ]  #テスト用データ
baysian_structure = hc(training.set)   #訓練用データでベイジアンネットワークの構造を学習
fitted = bn.fit(baysian_structure, training.set,method = "mle")     #パラメータの学習
test.set2 <- test.set
test.set2$PKC <- 0 #本番想定でこの値を知らないということにしておく。
pred = predict(fitted, "PKC", test.set2)  #テストデータが与えられたもとでのPKCの予測
head(cbind(pred, test.set[, "PKC"]))      #予測値と実績値の比較
accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出


> head(cbind(pred, test.set[, "PKC"]))      #予測値と実績値の比較
         pred     
[1,] 19.76000 22.3
[2,] 21.66300 26.4
[3,] 11.82080 27.9
[4,] 23.85982 53.3
[5,] 21.82902 34.3
[6,] 17.29319 16.8
> accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出
              ME     RMSE      MAE       MPE     MAPE
Test set 0.47376 10.02839 7.479115 -90.60837 114.6982

#ベイジアンネットワークによる予測

library(forecast)

training.set = protein_dataset[1:405, ] #パラメータ推定用の訓練用データ

test.set = protein_dataset[406:810, ] #テスト用データ

baysian_structure = hc(training.set) #訓練用データでベイジアンネットワークの構造を学習

fitted = bn.fit(baysian_structure, training.set,method = "mle") #パラメータの学習

test.set2 <- test.set

test.set2$PKC <- 0 #本番想定でこの値を知らないということにしておく。

pred = predict(fitted, "PKC", test.set2) #テストデータが与えられたもとでのPKCの予測

head(cbind(pred, test.set[, "PKC"])) #予測値と実績値の比較

accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出

> head(cbind(pred, test.set[, "PKC"])) #予測値と実績値の比較

pred

[1,] 19.76000 22.3

[2,] 21.66300 26.4

[3,] 11.82080 27.9

[4,] 23.85982 53.3

[5,] 21.82902 34.3

[6,] 17.29319 16.8

> accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出

ME RMSE MAE MPE MAPE

Test set 0.47376 10.02839 7.479115 -90.60837 114.6982

#予測結果の可視化
library(ggplot2)
df <- data.frame(prediction=pred,actual=test.set[, "PKC"])
p <- ggplot(df,aes(x = prediction,y = actual)) + geom_point() 
p <- p + geom_line(data = data.frame(x = c(0,40), y = c(0,40)),aes(x = x, y = y), colour = "red")
p

#予測結果の可視化

library(ggplot2)

df <- data.frame(prediction=pred,actual=test.set[, "PKC"])

p <- ggplot(df,aes(x = prediction,y = actual)) + geom_point()

p <- p + geom_line(data = data.frame(x = c(0,40), y = c(0,40)),aes(x = x, y = y), colour = "red")

全然45度線上に乗っていないので、あまり精度は高くないようです。

cpquery関数と言って、条件付き確率を予測できる関数もあるようです。
Perform conditional probability queries

おまけ：Webサービスへの応用例

参考文献には、Webサイトの閲覧データをもとに、ユーザーの行動を予測するというモデルの紹介がなされていました。ただ、Webサイトの場合、ページ数や商品の数・ユーザーの数も膨大なことから、ネットワーク構造の推定が難しいようです。そこで、ネットワークを作るに際して、変数を確率的潜在意味解析（pLSA）で扱いやすい数に絞るなどの工夫をされていました。自社のサイトに関しても適用する際に、変数の多さは待った無しだと思うので、いざ適用する際は情報圧縮技術を駆使したいですね。

参考文献

ベイジアンネットワーク技術ユーザ・顧客のモデル化と不確実性推論

確率的グラフィカルモデル

Learning Bayesian Networks in R an Example in Systems Biology

グラフィカルモデル入門

ベイジアンネットワークを用いたWeb レコメンデーションシステムの開発

PRML8章

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).