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’

蒙古タンメン中本コーパスに対してのLDAの適用とトピック数の探索

モチベーション

前回の記事では、Webスクレイピングにより入手した、蒙古タンメン中本の口コミデータに関して、Word2Vecを適用した特徴量エンジニアリングの事例を紹介しました。
今回はせっかく興味深いデータがあるので、どのようなトピックがあるのかをLDAを適用したいと思います。加えて、これまで記事で扱ってきたLDAの事例では評価指標であるPerplexityやCoherenceを扱ってこなかったことから、トピック数がどれくらいであるべきなのか、考察も含めて行いたいと思います。以前扱った階層ディリクレ過程であれば、トピック数を事前に決める必要が無いのですが、今回は扱わないものとします。

環境

・MacBook Pro
・Python3.5
・R version 3.4.4

Gensimで行うLDA

今回もPythonのGensimライブラリを用いて行います。

パープレキシティ
- テストデータに対して計算
- 負の対数尤度で、低いほどよい。
  - パープレキシティが低いと、高い精度で予測できるよい確率モデルと見なされる。汎化能力を表す指標。
  - トピックの数をいくらでも増やせばパープレキシティは下がる傾向が出ている。
  - 教科書でのパープレキシティの事例に関しては、トピック数を増やせば低くなるという傾向が出ている。

以下のコードでパープレキシティを計算します。

import pandas as pd
from gensim import corpora, models
import  gensim as  gensim

import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import FloatProgress
from IPython.display import display, clear_output

nakamoto_corpus = pd.read_pickle("nakamoto_corpus.pickle")

output_list = pd.DataFrame(columns=["y","x"])
plt.ion()
fig = plt.figure()
axe = fig.add_subplot(111)

for iterations in range(2, 100):
                        clear_output(wait = True)

                        nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)
                        nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

                        texts = [ ]
                        for line in nakamoto_train.text_wakati:
                            texts.append(line.split())

                        texts_test = [ ]
                        for line in nakamoto_test.text_wakati:
                            texts_test.append(line.split())

                        # 辞書作成
                        dictionary = corpora.Dictionary(texts)
                        dictionary.filter_extremes(no_below=20, no_above=0.3)

                        # コーパスを作成
                        corpus = [dictionary.doc2bow(text) for text in texts]
                        corpus_test = [dictionary.doc2bow(text) for text in texts_test]

                        # LDA の計算
                        topic_N = 1 + iterations
                        parameters = topic_N
                        lda = gensim.models.ldamodel.LdaModel(
                            corpus=corpus,
                            alpha='auto',
                            num_topics=topic_N,
                            id2word=dictionary
                        )

                        d = {'y':lda.log_perplexity(chunk=corpus_test),'x':parameters}
                        df = pd.DataFrame(data=d,index=[0])
                        output_list = output_list.append(df)
                        
                        axe.plot(output_list.x,output_list.y)
                        fig.set_size_inches(8, 8)
                        display(fig)
                        axe.cla()
                        
                        for i in range(topic_N):
                            print('TOPIC:', i, '__', lda.print_topic(i))

import pandas as pd

from gensim import corpora, models

import gensim as gensim

import matplotlib.pyplot as plt

import numpy as np

from ipywidgets import FloatProgress

from IPython.display import display, clear_output

nakamoto_corpus = pd.read_pickle("nakamoto_corpus.pickle")

output_list = pd.DataFrame(columns=["y","x"])

plt.ion()

fig = plt.figure()

axe = fig.add_subplot(111)

for iterations in range(2, 100):

clear_output(wait = True)

nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)

nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

texts = [ ]

for line in nakamoto_train.text_wakati:

texts.append(line.split())

texts_test = [ ]

for line in nakamoto_test.text_wakati:

texts_test.append(line.split())

# 辞書作成

dictionary = corpora.Dictionary(texts)

dictionary.filter_extremes(no_below=20, no_above=0.3)

# コーパスを作成

corpus = [dictionary.doc2bow(text) for text in texts]

corpus_test = [dictionary.doc2bow(text) for text in texts_test]

# LDA の計算

topic_N = 1 + iterations

parameters = topic_N

lda = gensim.models.ldamodel.LdaModel(

corpus=corpus,

alpha='auto',

num_topics=topic_N,

id2word=dictionary

)

d = {'y':lda.log_perplexity(chunk=corpus_test),'x':parameters}

df = pd.DataFrame(data=d,index=[0])

output_list = output_list.append(df)

axe.plot(output_list.x,output_list.y)

fig.set_size_inches(8, 8)

display(fig)

axe.cla()

for i in range(topic_N):

print('TOPIC:', i, '__', lda.print_topic(i))

実際に、中本コーパスで計算したトピック数に対してのパープレキシティは以下のように推移しました。

Ldaのモデル選択におけるperplexityの評価によると、
”複数のトピック数で比べて、Perplexityが最も低いものを選択する。」という手法は人間にとって有益なモデルを選択するのに全く役に立たない可能性がある。”と記されています。

『トピックモデルによる統計的潜在意味解析』には、”識別問題の特徴量として使う場合は識別問題の評価方法で決定すればよい”とあるので、目的によってはパープレキシティにこだわらなくても良いと思われます。

今回のケースだと、パープレキシティだけだと、決めかねてしまいますね。

コヒーレンス
- トピックごとの単語間類似度の平均
- トピック全体のコヒーレンスが高ければ、良い学習アルゴリズムとみなす。

以下のコードでコヒーレンスを計算します。

output_list = pd.DataFrame(columns=["y","x"])
plt.ion()
fig = plt.figure()
axe = fig.add_subplot(111)

for iterations in range(2, 100):
                        clear_output(wait = True)

                        nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)
                        nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

                        texts = [ ]
                        for line in nakamoto_train.text_wakati:
                            texts.append(line.split())

                        texts_test = [ ]
                        for line in nakamoto_test.text_wakati:
                            texts_test.append(line.split())

                        # 辞書作成
                        dictionary = corpora.Dictionary(texts)
                        dictionary.filter_extremes(no_below=20, no_above=0.3)

                        # コーパスを作成
                        corpus = [dictionary.doc2bow(text) for text in texts]
                        corpus_test = [dictionary.doc2bow(text) for text in texts_test]

                        # LDA の計算
                        topic_N = 1 + iterations
                        parameters = topic_N
                        lda = gensim.models.ldamodel.LdaModel(
                            corpus=corpus,
                            alpha='auto',
                            num_topics=topic_N,
                            id2word=dictionary
                        )

                        
                        cm = models.coherencemodel.CoherenceModel(model=lda, corpus=corpus, coherence='u_mass')  # tm is the trained topic model                        
                        d = {'y':cm.get_coherence(),'x':parameters}
                        df = pd.DataFrame(data=d,index=[0])
                        output_list = output_list.append(df)
                        
                        axe.plot(output_list.x,output_list.y)
                        fig.set_size_inches(8, 8)
                        display(fig)
                        axe.cla()
                        
                        for i in range(topic_N):
                            print('TOPIC:', i, '__', lda.print_topic(i))

output_list = pd.DataFrame(columns=["y","x"])

plt.ion()

fig = plt.figure()

axe = fig.add_subplot(111)

for iterations in range(2, 100):

clear_output(wait = True)

nakamoto_test = nakamoto_corpus.sample(frac=0.1, replace=True)

nakamoto_train= nakamoto_corpus[~nakamoto_corpus.index.isin(nakamoto_test.index)]

texts = [ ]

for line in nakamoto_train.text_wakati:

texts.append(line.split())

texts_test = [ ]

for line in nakamoto_test.text_wakati:

texts_test.append(line.split())

# 辞書作成

dictionary = corpora.Dictionary(texts)

dictionary.filter_extremes(no_below=20, no_above=0.3)

# コーパスを作成

corpus = [dictionary.doc2bow(text) for text in texts]

corpus_test = [dictionary.doc2bow(text) for text in texts_test]

# LDA の計算

topic_N = 1 + iterations

parameters = topic_N

lda = gensim.models.ldamodel.LdaModel(

corpus=corpus,

alpha='auto',

num_topics=topic_N,

id2word=dictionary

)

cm = models.coherencemodel.CoherenceModel(model=lda, corpus=corpus, coherence='u_mass') # tm is the trained topic model

d = {'y':cm.get_coherence(),'x':parameters}

df = pd.DataFrame(data=d,index=[0])

output_list = output_list.append(df)

axe.plot(output_list.x,output_list.y)

fig.set_size_inches(8, 8)

display(fig)

axe.cla()

for i in range(topic_N):

print('TOPIC:', i, '__', lda.print_topic(i))

実際に推定してみたところ、トピック数が20を超えたあたりからコヒーレンスが下がる傾向があるので、
それ以上のトピック数は追い求めない方が良いのかもしれません。

Rでもやってみる

Rでトピック数を決める良い方法がないか調べてみたところ、ldatuningとかいうパッケージがあることがわかりました。複数の論文（Griffiths2004, CaoJuan2009, Arun2010,Deveaud2014）で扱われている手法を元に、適切なトピック数を探れるようです。このパッケージを紹介しているブログの事例では、90から140の範囲で最適なトピック数となることが示されています。詳しくはこちらを見てください。
Select number of topics for LDA model

以下のコードで実行しました。一部、驚異のアニヲタさんのコードを拝借しております。なお、ldaパッケージのlexicalize関数を用いることで、ldatuningに入力するデータを作成することができます。

library(ldatuning)
library(topicmodels)
library(tidyverse)

nakamoto_data <- read_csv(file = "nakamoto_dataset.csv")

TFIDF <- function(corpus, progress=FALSE){ 
  res <- matrix(0, nr=length(corpus$vocab), nc=4)
  dimnames(res) <- list(corpus$vocab, c("documents", "count", "freq", "score"))
  res[, "documents"] <- length(corpus$documents)
  wordset <- mapply(function(x) x[1,], corpus$documents) 
  allfreq <- matrix(unlist(corpus$documents), nr=2)
  wordfreq <- tapply(allfreq[2,], allfreq[1,], sum) 
  for(v in seq(corpus$vocab)){ 
    count_docs <- sum(sapply(lapply(wordset, "==", v-1), any)) #
    res[v, "freq"] <- count_docs
    if(progress){ 
      pb <- txtProgressBar(min=1, max=length(corpus$vocab), style=3)
      setTxtProgressBar(pb, v)
    }
  }
  res[, "count"] <- wordfreq
  res[, "score"] <- log(res[, "count"]) * log(res[, "documents"]/res[, "freq"])
  return(as.data.frame(res))
}

lex1 <- lda::lexicalize(nakamoto_data$text_wakati_remove_freq)

s0 <- TFIDF(lex1, TRUE)
term1 <- rownames(s0)[s0$score > 0]
lex2 <- list(documents=lda::lexicalize(nakamoto_data$text_wakati_remove_freq, vocab=term1), vocab=term1)
dtm2 <- ldaformat2dtm(lex2$documents, lex2$vocab)

result <- FindTopicsNumber(
  dtm2,
  topics = seq(from = 2, to = 100, by = 1),
  metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),
  method = "Gibbs",
  control = list(seed = 77),
  mc.cores = 2L,
  verbose = TRUE
)

FindTopicsNumber_plot(values = result)

library(ldatuning)

library(topicmodels)

library(tidyverse)

nakamoto_data <- read_csv(file = "nakamoto_dataset.csv")

TFIDF <- function(corpus, progress=FALSE){

res <- matrix(0, nr=length(corpus$vocab), nc=4)

dimnames(res) <- list(corpus$vocab, c("documents", "count", "freq", "score"))

res[, "documents"] <- length(corpus$documents)

wordset <- mapply(function(x) x[1,], corpus$documents)

allfreq <- matrix(unlist(corpus$documents), nr=2)

wordfreq <- tapply(allfreq[2,], allfreq[1,], sum)

for(v in seq(corpus$vocab)){

count_docs <- sum(sapply(lapply(wordset, "==", v-1), any)) #

res[v, "freq"] <- count_docs

if(progress){

pb <- txtProgressBar(min=1, max=length(corpus$vocab), style=3)

setTxtProgressBar(pb, v)

}

res[, "count"] <- wordfreq

res[, "score"] <- log(res[, "count"]) * log(res[, "documents"]/res[, "freq"])

return(as.data.frame(res))

}

lex1 <- lda::lexicalize(nakamoto_data$text_wakati_remove_freq)

s0 <- TFIDF(lex1, TRUE)

term1 <- rownames(s0)[s0$score > 0]

lex2 <- list(documents=lda::lexicalize(nakamoto_data$text_wakati_remove_freq, vocab=term1), vocab=term1)

dtm2 <- ldaformat2dtm(lex2$documents, lex2$vocab)

result <- FindTopicsNumber(

dtm2,

topics = seq(from = 2, to = 100, by = 1),

metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),

method = "Gibbs",

control = list(seed = 77),

mc.cores = 2L,

verbose = TRUE

)

FindTopicsNumber_plot(values = result)

これを見る限りは、60〜70個の辺りに落ち着くのでしょうか。

トピックの吐き出し

Rでの結果から、60個程度のトピックで推定し、各記事に割り当てが最大のトピックを付与して、トピック別の口コミ評価をみてみようと思います。

以下のコードではトピック別の口コミ評価のしやすさからtopicmodelsパッケージを用いた推定となっています。

nbo_topics <- 60
lda_estimate <- topicmodels::LDA(dtm2,control=list(verbose=1, alpha = 0.1),
                                 k = nbo_topics,
                                 method = "Gibbs")

#トピックの上位5単語を確認する
terms_each_topics <- data.frame(terms(lda_estimate,5))
topic_keywords <- data.frame(topic_keywords_5 = apply(t(terms_each_topics),1,paste,collapse=","))
topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

#割り振られた最大の確率のトピックを抽出し、口コミデータと統合する
topics_each_document <- data.frame(topic_id=topics(lda_estimate,1))
topics_each_document <- nakamoto_data %>% bind_cols(topics_each_document)

#トピックごとの口コミ評価を計算する
topic_rating_summary <- topics_each_document %>% 
                            group_by(topic_id) %>% 
                            summarise(average_rating = mean(rating),
                                      count=n()) %>% 
                            left_join(topic_keywords, by = "topic_id") %>% 
                            arrange(desc(average_rating))

nbo_topics <- 60

lda_estimate <- topicmodels::LDA(dtm2,control=list(verbose=1, alpha = 0.1),

k = nbo_topics,

method = "Gibbs")

#トピックの上位5単語を確認する

terms_each_topics <- data.frame(terms(lda_estimate,5))

topic_keywords <- data.frame(topic_keywords_5 = apply(t(terms_each_topics),1,paste,collapse=","))

topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

#割り振られた最大の確率のトピックを抽出し、口コミデータと統合する

topics_each_document <- data.frame(topic_id=topics(lda_estimate,1))

topics_each_document <- nakamoto_data %>% bind_cols(topics_each_document)

#トピックごとの口コミ評価を計算する

topic_rating_summary <- topics_each_document %>%

group_by(topic_id) %>%

summarise(average_rating = mean(rating),

count=n()) %>%

left_join(topic_keywords, by = "topic_id") %>%

arrange(desc(average_rating))

口コミ評価の点数が上位のトピックはこんな感じです。

口コミ評価の点数が下位のトピックはこんな感じです。

中本は社会人2〜3年目で新規メディアの立ち上げのストレス解消で数回行きましたが、北極の赤さは異常だと思います。北極を食べたり、トッピングする余裕のある人、ましてや辛さを倍にするという時点で口コミ評価も高くなると考えるのは自然なのかもしれません。

参考情報

トピックモデル (機械学習プロフェッショナルシリーズ)
トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)
models.ldamodel – Latent Dirichlet Allocation
Ldaのモデル選択におけるperplexityの評価
 pythonでgensimを使ってトピックモデル(LDA)を行う
 gensim0.8.6のチュートリアルをやってみた【コーパスとベクトル空間】
LDA 実装の比較
 Jupyter notebookにMatplotlibでリアルタイムにチャートを書く
 Inferring the number of topics for gensim’s LDA – perplexity, CM, AIC, and BIC
Select number of topics for LDA model
47の心得シリーズをトピックモデルで分類する。 – 驚異のアニヲタ社会復帰への道

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版

RのFactoRizationMachinesパッケージを用いたFMのページビューデータへの適用

Googleタグマネージャーで集めたアクセスログデータを用いて、前回と同様に記事のレコメンドにチャレンジしてみようと思います。FactoRizationMachinesパッケージという便利そうなパッケージの存在も知れたことから、今回は以前から気になっていたFactorization Machineを扱います。

【目次】
・Factorization Machine(FM)の概要
・パッケージ紹介とインストール
・サンプルデータの構造把握と前処理
・FMの実行
・結果
・参考文献

Factorization Machine(FM)の概要

組み合わせ特徴量を扱う教師あり学習モデル。行列分解とSVMを合体させた手法。
スパースになりやすいデータの予測問題で扱う。
1ユーザーのある商品に対しての評価を、1評価1行として表して、ユーザーとアイテムの交互作用の特徴ベクトルを扱う。
相互作用項に関して、時間や文脈などを自由に入れられる。
相互作用項を次元圧縮する際の要素数を事前に決める必要がある。
Matrix Factorizationよりも精度が良いとされている。特徴量エンジニアリングなどで使われているようです。（Click-Through Rate Prediction）

パッケージ紹介とインストール

FactoRizationMachinesパッケージは線形SVMと2次のFMと高次のFMを実行することができ、引数で正則化項も加えることができます。現段階においては回帰のみで分類問題への適用は今後の開発となるようです。CRANから普通に
install.packages(‘FactoRizationMachines’)
でインストールします。libFMexeパッケージの場合は、libFMをインストールしてパスを指定しておく必要がありますが、このパッケージに関しては不要となります。

サンプルデータの構造把握と前処理

FactoRizationMachinesパッケージのサンプルコードにおいては、MovieLensのデータがサンプルデータとして載せられていました。ユーザーのID(整数)、映画のID(整数)、評価(整数、5段階)、日時(整数)からなるデータに対して、sparseMatrixに変換していました。

今回は、前回の投稿非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみると同じデータを使って、アクセスログデータに適用しようと思います。FactoRizationMachinesの形式に合わせるために、このブログのアクセスログも、クッキーのIDを整数に、記事のIDを整数に、閲覧回数を5段階(5以上を5に変換)に、日時を整数に変更しています。

library(RGA)
library(tidyverse)
library(reshape2)

authorize()
prof <-list_profiles()

start_date <- "2017-07-04"
end_date   <- "2017-10-07"

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

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_")

#ユーザーIDの整数化
user <- unique(accesslogdata$dimension1)
user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化
articles <- unique(accesslogdata$pagePath)
articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合
accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))
accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")
accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出
last_date <- accesslogdata %>% select(user_number,article_number,date) %>% 
                arrange(user_number,desc(date))
last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計
accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>% 
                        select(user_number,article_number,pageview)
accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

library(RGA)

library(tidyverse)

library(reshape2)

authorize()

prof <-list_profiles()

start_date <- "2017-07-04"

end_date <- "2017-10-07"

#demention until 7

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

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dimension1,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",

fetch.by = "day")

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_")

#ユーザーIDの整数化

user <- unique(accesslogdata$dimension1)

user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化

articles <- unique(accesslogdata$pagePath)

articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合

accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))

accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")

accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出

last_date <- accesslogdata %>% select(user_number,article_number,date) %>%

arrange(user_number,desc(date))

last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計

accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>%

select(user_number,article_number,pageview)

accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

FMの実行

デフォルトの設定c(1, 10)では線形のウェイトが有効で、2次の項の要素数が10で正則化項なしのFMを実行することになります。引数に関する詳しい情報はPackage ‘FactoRizationMachines’に書かれています。今回はサンプルを参考に正則化項ありでモデルを実行します。まず、アクセスログデータに対して、ユーザーのIDからなる整数ベクトル、記事のIDからなる整数ベクトル、セッションのあった日時のデータからなる整数ベクトルを作成し、sparseMatrix関数を用いて元データを変形し、80%のデータをトレーニングに、20%のデータをテストに割り当てます。さらに、テストデータに関して、予測値との平均二乗誤差を計算します。

library(FactoRizationMachines)
library(Matrix)

user = accesslogdata_sum$user_number
items = accesslogdata_sum$article_number + max(user)
wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form
data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)
target = accesslogdata_sum$pageview

# Subset data to training and test data
set.seed(123)
subset = sample.int(nrow(data),nrow(data)*.8)
data.train = data[subset,]
data.test = data[-subset,]
target.train = target[subset]
target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel
model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine
# with second-order 10 factors (default) and regularization
model = FM.train(data.train,target.train,regular=0.1)
# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

library(FactoRizationMachines)

library(Matrix)

user = accesslogdata_sum$user_number

items = accesslogdata_sum$article_number + max(user)

wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form

data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)

target = accesslogdata_sum$pageview

# Subset data to training and test data

set.seed(123)

subset = sample.int(nrow(data),nrow(data)*.8)

data.train = data[subset,]

data.test = data[-subset,]

target.train = target[subset]

target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel

model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine

# with second-order 10 factors (default) and regularization

model = FM.train(data.train,target.train,regular=0.1)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

結果

各モデルについての平均二乗誤差を計算しています。
線形モデルや高次元モデルよりも、2次の項を持つFMが精度が高いようです。

> model = SVM.train(data.train,target.train)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.514266
> model = FM.train(data.train,target.train,regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.369983
> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.388247

> model = SVM.train(data.train,target.train)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.514266

> model = FM.train(data.train,target.train,regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.369983

> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.388247

こちらは、この中で性能の良かった2次の項を持つFMの予測結果とテストデータの結果をプロットしたものです。4点を超える値をあまり予測できていないようです。今回はサンプルを回しただけなので、本来であれば次元の数kや正則化のセッティングをいろいろいじったり、相互作用項を新しく追加するなどして精度を高めることが必要です。

結果の比較だけでは仕事で使えないので、実際に予測した結果を取り出したいと思います。
実際に運用するとなると、ページIDを所与として、ページビュー数を0とおいて（型をそろえるため。NULLだとエラーになった）、任意のタイミング(date)を想定して、モデルにデータを適用し、評価の高いものをサジェストするスタイルになるのではないでしょうか。

この結果だと、ユーザー98に記事1を見せることに対して4.02点が与えられています。

参考文献

Factorization Machinesを今更読みました
 Factorization Machines
High-order factorization machines with R #tokyor 61
Factorization Machinesのおはなし。
libFMexeを動かすまで (R Wrapper for the libFM Executable参照記事)
一歩Matrix Factorization、二歩Factorization Machines、三歩Field-aware Factorization Machines…『分解、三段突き！！』
[論文] Factorization Machines (ICDM 2010) 読んだ 22:41
Factorization machines with r
Factorization Machinesについて調べてみた

非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみる

概要

2017年7月より当ブログにおいて、Google AnalyticsのカスタムディメンジョンにUser IDを取得して残せるようにしました。これにより、ユーザーの一人一人が閲覧した記事の分析を行うことが可能となっています。そこで今回は、記事の数も50近くなってきたことから、レコメンド手法の適用を行ってみたいと思います。Exampleデータでレコメンドをしても、ドメイン知識がないと、どれくらい良さそうな推薦なのかの判断が付きにくいと考えていたので、自分で運営しているサイトのデータであれば評価をつけやすいのではないかと思われます。
手法としてはNMF（non-Negative Matrix Factorization）を適用します。

【目次】
・前処理
・データ確認
・行列因子分解（Matrix Factorization）
・レコメンド結果の確認
・感想
・参考文献

前処理

ユーザー単位の閲覧した記事データを作るためにGAのAPIで取得したデータを整形し、集計します。このようなデータを手に入れるには自身のブログにおいて、GTM(Google Tag Maneger)を設置して、GTMのコンテナにカスタムJS変数などを格納する必要があります。詳しくはこちらのブログを参照ください。（ 2016年の新定番！ユーザーエクスプローラーをもっと活用するための簡単な方法）

#パッケージの読み込み
library(RGA)
library(tidyverse)

#GAの認証
authorize()
prof <-list_profiles()

#取得期日の指定
start_date <- "2017-07-04"
end_date   <- "2017-09-16"


#GAのAPIからデータ取得（ディメンションの数は7つまで）
accesslogdata <- get_ga(profileId = prof$id[2],
              start.date = start_date,
              end.date = end_date,
              dimensions = "ga:pagePath,
                            ga:dimension1,
                            ga:dateHourMinute,
                            ga:deviceCategory,
                            ga:userType,
                            ga:fullReferrer",
              sort = "-ga:sessions",
              metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",
              fetch.by = "day")

#記事のURLだけに絞り込む
accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

#ユーザーごとに重複した記事の情報を削除する
df_dups <- accesslogdata[c("pagePath","dimension1")]
accesslogdata <- accesslogdata[!duplicated(df_dups),]
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

articles <- unique(accesslogdata$pagePath)

#記事ごとの閲覧有無に関するカラムを生成する
for(i in 1:length(articles)){
  accesslogdata <- accesslogdata %>% mutate(variable=if_else(pagePath==articles[i],1,0))
  accesslogdata <- eval(parse(text=paste0("accesslogdata %>% rename( ",articles[i],"=variable)")))
}

#クッキー単位で接触した記事ごとに集計する
accesslogdata <- accesslogdata %>% select(-pagePath,-dateHourMinute,-deviceCategory,-userType,-fullReferrer,-sessions,-goal1Completions,-dateHour,-pageviews) %>% 
                    group_by(dimension1) %>% summarise_all(funs(sum))

#パッケージの読み込み

library(RGA)

library(tidyverse)

#GAの認証

authorize()

prof <-list_profiles()

#取得期日の指定

start_date <- "2017-07-04"

end_date <- "2017-09-16"

#GAのAPIからデータ取得（ディメンションの数は7つまで）

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

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dimension1,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",

fetch.by = "day")

#記事のURLだけに絞り込む

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

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

#ユーザーごとに重複した記事の情報を削除する

df_dups <- accesslogdata[c("pagePath","dimension1")]

accesslogdata <- accesslogdata[!duplicated(df_dups),]

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

articles <- unique(accesslogdata$pagePath)

#記事ごとの閲覧有無に関するカラムを生成する

for(i in 1:length(articles)){

accesslogdata <- accesslogdata %>% mutate(variable=if_else(pagePath==articles[i],1,0))

accesslogdata <- eval(parse(text=paste0("accesslogdata %>% rename( ",articles[i],"=variable)")))

}

#クッキー単位で接触した記事ごとに集計する

accesslogdata <- accesslogdata %>% select(-pagePath,-dateHourMinute,-deviceCategory,-userType,-fullReferrer,-sessions,-goal1Completions,-dateHour,-pageviews) %>%

group_by(dimension1) %>% summarise_all(funs(sum))

このような形で行にユーザーが、列に記事が、各要素にアクセスの有無が入るようになっています。本来は評価データを用いますが、評価データがないので、アクセスしたかどうかの1-0データとして扱います。

データ確認

#ユーザーごとの接触ページ数を集計し、ヒストグラムで描く
session <- data.frame(PageView=rowSums(accesslogdata[,-1]))
ggplot(data = session, aes(x = PageView)) + geom_histogram(binwidth = 1) +
    ylab("User Count") + ggtitle("Article Page View Each Users") + 
    theme(plot.title = element_text(hjust = 0.5))

#ユーザーごとの接触ページ数を集計し、ヒストグラムで描く

session <- data.frame(PageView=rowSums(accesslogdata[,-1]))

ggplot(data = session, aes(x = PageView)) + geom_histogram(binwidth = 1) +

ylab("User Count") + ggtitle("Article Page View Each Users") +

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

以下より、行列の要素のうち、値が入っている要素は全体の2.9%であることがわかります。参考文献は2.6%をスパースと言っていたことから、今回のデータもスパースであると考えることができます。

> X <- as(accesslogdata[,-1],"matrix")
> table(X)
X
     0      1 
122741   3777 
> table(X)[2]/(table(X)[1]+table(X)[2])
         1 
0.02985346

> X <- as(accesslogdata[,-1],"matrix")

> table(X)

0 1

122741 3777

> table(X)[2]/(table(X)[1]+table(X)[2])

0.02985346

行列因子分解（Matrix Factorization）

行列因子分解はユーザーごとのアイテム購買履歴などのデータからなる行列の因子分解を行うもので、レコメンドにおいては因子分解の結果として、ゼロ以上の要素を持つ場合に推薦対象とするなどの活用がされています。非負値行列因子分解は要素が非負となるように制約を付けて行列を分解するもので、正の値を取ることから推薦に使いやすいとされています。

今回はNMFパッケージを用いて、ほぼデフォルトな設定で結果を眺めてみますが。本来であれば以下のような考慮が必要となります。最適化などをがっつり進めようとするとクロスバリデーションなどを行ったり、計算リソースが結構かかったりすると思われます。

アルゴリズムの選択
- brunet,KL,lee,Frobenius,offset,nsNMF,ls-nmf,pe-nmf,siNMF,snmf/r,snmf/lがある。
- デフォルトではbrunetとなっている。
シーディングメソッドの選択（初期値の決め方）
- 固定値(none),ランダム(random),ica,nndsvdがある。
- デフォルトではランダムとなっている。
実行回数の設定
- １００〜200回を選ぶことが多い。
並列処理
- 処理量が多いので最適化などを考えると並列処理も重要。
因子分解ランクの推定
- ランクの数が多いと残差は減るが、過学習になりやすい。クロスバリデーションを行う必要がある。

library(NMF)

#NMFの実行（6次元以上の次元を削除）
res <- nmf(accesslogdata[,-1], 5, seed = 1234, .options = "t")

#試行回数と損失関数としてのカルバックライブラー情報量の可視化
plot(res)

#基底行列の取得
w <- basis(res)

#係数行列の取得
h <- coef(res)

#レコメンデーション評価値を得る
accesslogdata.nmf <- w %*% h

#確認するユーザーの行番号
userid <- c(377,987,2584)

#RAWデータとレコメンデーション結果の比較
recommend_result <- data.frame(round(t(accesslogdata.nmf[userid,]), 2),
           t(accesslogdata[userid,-1]))

#記事タイトルの読み込み
article_title <- read_csv(file = "article_titles.csv")

recommend_result <- left_join(rownames_to_column(recommend_result), article_title, by = c("rowname" = "article_id")) %>% 
                        mutate(rowname=title) %>% select(-title)

library(NMF)

#NMFの実行（6次元以上の次元を削除）

res <- nmf(accesslogdata[,-1], 5, seed = 1234, .options = "t")

#試行回数と損失関数としてのカルバックライブラー情報量の可視化

plot(res)

#基底行列の取得

w <- basis(res)

#係数行列の取得

h <- coef(res)

#レコメンデーション評価値を得る

accesslogdata.nmf <- w %*% h

#確認するユーザーの行番号

userid <- c(377,987,2584)

#RAWデータとレコメンデーション結果の比較

recommend_result <- data.frame(round(t(accesslogdata.nmf[userid,]), 2),

t(accesslogdata[userid,-1]))

#記事タイトルの読み込み

article_title <- read_csv(file = "article_titles.csv")

recommend_result <- left_join(rownames_to_column(recommend_result), article_title, by = c("rowname" = "article_id")) %>%

mutate(rowname=title) %>% select(-title)

レコメンド結果の確認

各ユーザーに対して推薦されたアイテムのうち、比較的値の大きなものに関してコメントをしたいと思います。

ユーザー1
「XGBoostやパラメータチューニングの仕方に関する調査」を閲覧したユーザー1に対して、「XGBoostのパラメータチューニング実践 with Python」を推薦できている一方で、「OpenCV&Pythonで画像の類似度を計算させる〜イケメンの顔比較」など関係性の薄そうな記事も推薦してしまっています。

ユーザー2
「人工知能学会全国大会2017のWebマーケティングで参考になりそうな研究9選」を閲覧したユーザー2に対して、「Kaggleで使われた特徴量エンジニアリングとアルゴリズムまとめ」を推薦できている。

ユーザー3
「洋楽の歌詞データでDoc2vecを実行してみる」を閲覧したユーザー3に対して、「LDA（潜在的ディリクレ配分法）まとめ手法の概要と試行まで」を推薦できている一方で、「ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる」と関係性の薄そうな記事も推薦している。本来、近いと思われる、「Word2Vecでクラシックの楽曲情報をコーパスとして類似度を出してみる」を推薦できていない。

感想

結果に関して、納得のいくものといかないもので半々な印象でデフォルトの設定では満足できない。
NMFはクロスバリデーションなどにより過学習を避けた形での最適なランク数やアルゴリズムの選定、初期値の設定などが必要と思われるので、そこらへんの取り組みも今後行いたい。
FM(Factorization Machine)は回帰項を追加できるので、今後挑戦してみたい。
1-0のデータではなく評価のデータが欲しい。

参考文献

岩波データサイエンス Vol.5
非負値行列因子分解（NMF）によるレコメンドのちょっとした例
 An introduction to NMF package
非負ッ値行列・ファクトリゼーション～半群じゃないから難しくないもん！

ベイジアンネットワークを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章

Billboard100位以内の楽曲の歌詞情報にLDAを適用してみた

目次

・はじめに
・データ収集
・Rによる分析
・LDAの結果
・参考文献

はじめに

前回の投稿でBillboardの週次洋楽ランキングデータをWebスクレイピングで取得し、楽曲の消費サイクルのような順位の挙動を確かめることができました。（某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた）今回は、歌詞の情報を用いて順位データとつなぐことにより、どのような単語の入っている洋楽がBillboardにおいてTop10に入る傾向があるのかをLDAを行うことで確かめたいと思います。

データ収集

残念なことに、Billboardのサイトに歌詞の情報は載っていません。そこで、洋楽の歌詞が取り上げられている某サイトをPython(3系)でWebスクレイピングし、名寄せを頑張って順位データと歌詞データを繋ぎます。

幸いなことに某サイトのURLに規則性があったので、アーティスト名からなるURLを生成し、そのURLをWebスクレイピングして楽曲のリストを集め、今回のBillboardのランキングに入った楽曲のみに絞ります。

#アーティストの楽曲一覧の取得

import urllib
from bs4 import BeautifulSoup
from urllib.request import urlopen
import requests
from urllib.error import HTTPError, URLError
import csv, re, time
from requests.exceptions import ConnectionError

f = open('artist_url_list.csv', 'r')
dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。
data01 =[] #URL
data02 =[] #曲名
data03 =[] #link

for row in dataReader:
       for url in row:
            time.sleep(10.0) #sleep(秒指定)
            try:
                    r = requests.get(url)
                    soup =  BeautifulSoup(r.content, 'html.parser')
                    
                    for body in soup.findAll("td",{'class':'colfirst'}):
                        for link in body.findAll("a"):
                                data01.append(url)
                                data02.append(''.join(link.findAll(text=True)))
                                data03.append(link.get("href"))
                                data = zip(data01,data02,data03)
                                #CSV出力
                                with open('artistpage_result.csv','wt',errors='backslashreplace') as fout:
                                    writecsv = csv.writer(fout,lineterminator='\n')
                                    writecsv.writerows(data)                                   
                                
                                    
            except HTTPError as e:
                print(e.code)
                
            except URLError:
                print("URLError")

#アーティストの楽曲一覧の取得

import urllib

from bs4 import BeautifulSoup

from urllib.request import urlopen

import requests

from urllib.error import HTTPError, URLError

import csv, re, time

from requests.exceptions import ConnectionError

f = open('artist_url_list.csv', 'r')

dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。

data01 =[] #URL

data02 =[] #曲名

data03 =[] #link

for row in dataReader:

for url in row:

time.sleep(10.0) #sleep(秒指定)

try:

r = requests.get(url)

soup = BeautifulSoup(r.content, 'html.parser')

for body in soup.findAll("td",{'class':'colfirst'}):

for link in body.findAll("a"):

data01.append(url)

data02.append(''.join(link.findAll(text=True)))

data03.append(link.get("href"))

data = zip(data01,data02,data03)

#CSV出力

with open('artistpage_result.csv','wt',errors='backslashreplace') as fout:

writecsv = csv.writer(fout,lineterminator='\n')

writecsv.writerows(data)

except HTTPError as e:

print(e.code)

except URLError:

print("URLError")

楽曲をランキングに含まれるもののみに絞ったら、歌詞詳細ページを取得します。

#歌詞詳細の歌詞該当部分のみ抽出

import urllib
from bs4 import BeautifulSoup
from urllib.request import urlopen
import requests
from urllib.error import HTTPError, URLError
import csv, re, time
from requests.exceptions import ConnectionError

f = open('song_detail_url.csv', 'r')
dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。
data01 =[] #URL
data02 =[] #歌詞

for row in dataReader:
       for url in row:
            time.sleep(10.0) #sleep(秒指定)
            try:
                    r = requests.get(url)
                    soup =  BeautifulSoup(r.content, 'html.parser')
                    
                    for body in soup.findAll("div",{'id':'content_h'}):
                                data01.append(url)
                                data02.append(''.join(body.findAll(text=True)))
                                data = zip(data01,data02)
                                #CSV出力
                                with open('lyrics_result.csv','wt',errors='backslashreplace') as fout:
                                    writecsv = csv.writer(fout,lineterminator='\n')
                                    writecsv.writerows(data)                                   
                                
                                    
            except HTTPError as e:
                print(e.code)
                           
            except URLError:
                print("URLError")

#歌詞詳細の歌詞該当部分のみ抽出

import urllib

from bs4 import BeautifulSoup

from urllib.request import urlopen

import requests

from urllib.error import HTTPError, URLError

import csv, re, time

from requests.exceptions import ConnectionError

f = open('song_detail_url.csv', 'r')

dataReader = csv.reader(f)

#結果の出力用のデータフレームを作る。

data01 =[] #URL

data02 =[] #歌詞

for row in dataReader:

for url in row:

time.sleep(10.0) #sleep(秒指定)

try:

r = requests.get(url)

soup = BeautifulSoup(r.content, 'html.parser')

for body in soup.findAll("div",{'id':'content_h'}):

data01.append(url)

data02.append(''.join(body.findAll(text=True)))

data = zip(data01,data02)

#CSV出力

with open('lyrics_result.csv','wt',errors='backslashreplace') as fout:

writecsv = csv.writer(fout,lineterminator='\n')

writecsv.writerows(data)

except HTTPError as e:

print(e.code)

except URLError:

print("URLError")

うまいこと歌詞情報を手に入れることができました。ざっと947曲です。

	url	lyrics
0	http://www.lyricsfreak.com/e/eminem/love+the+w...	Just gonna stand there and watch me burnBut th...
1	http://www.lyricsfreak.com/t/taio+cruz/dynamit...	I came to dance-dance-dance-dance (Yeah)I hate...
2	http://www.lyricsfreak.com/t/taylor+swift/mine...	Oh, oh, ohOh, oh, ohYou were in college, worki...
3	http://www.lyricsfreak.com/e/enrique+iglesias/...	One life, one loveEnrique Iglesias, PitbullY'a...
4	http://www.lyricsfreak.com/b/bob/airplanes_208...	Can we pretend that airplanesIn the night sky ...
5	http://www.lyricsfreak.com/m/mike+posner/coole...	If I could write you a song,And make you fall ...
6	http://www.lyricsfreak.com/j/jason+derulo/ridi...	Yea yeah, yeah, yeah, yeah,I'm feeling like a ...
7	http://www.lyricsfreak.com/t/travie+mccoy/bill...	I wanna be a billionaire so freakin' badBuy al...
8	http://www.lyricsfreak.com/d/drake/find+your+l...	I'm more than just an option (hey, hey, hey) R...
9	http://www.lyricsfreak.com/u/usher/omg_2087748...	Oh my goshBaby let meDid it again, so Imma let...
10	http://www.lyricsfreak.com/b/bob/magic_2087969...	I got the magic in meEvery time I touch that t...
11	http://www.lyricsfreak.com/n/nicki+minaj/your+...	[Chorus]Shawty I'm a only tell you this once, ...
12	http://www.lyricsfreak.com/m/maroon+5/misery_2...	Oh yeahOh yeahSo scared of breaking itThat you...
13	http://www.lyricsfreak.com/t/train/hey+soul+si...	Hey, hey, heyYour lipstick stains on the front...
14	http://www.lyricsfreak.com/b/bruno+mars/just+t...	Oh her eyes, her eyesMake the stars look like ...
15	http://www.lyricsfreak.com/l/lady+gaga/alejand...	I know that we are young,And I know that you m...
16	http://www.lyricsfreak.com/l/la+roux/bulletpro...	Been there, done that, messed aroundI'm having...
17	http://www.lyricsfreak.com/f/flo+rida/club+can...	You know I know howTo make 'em stop and stare ...
18	http://www.lyricsfreak.com/s/shontelle/impossi...	I remember years agoSomeone told me I should t...
19	http://www.lyricsfreak.com/p/paramore/the+only...	When I was youngerI saw my daddy cryAnd curse ...
20	http://www.lyricsfreak.com/u/usher/there+goes+...	Yeah, Right,Usher baby, OKYeah man, rightThere...

url lyrics

0 http://www.lyricsfreak.com/e/eminem/love+the+w... Just gonna stand there and watch me burnBut th...

1 http://www.lyricsfreak.com/t/taio+cruz/dynamit... I came to dance-dance-dance-dance (Yeah)I hate...

2 http://www.lyricsfreak.com/t/taylor+swift/mine... Oh, oh, ohOh, oh, ohYou were in college, worki...

3 http://www.lyricsfreak.com/e/enrique+iglesias/... One life, one loveEnrique Iglesias, PitbullY'a...

4 http://www.lyricsfreak.com/b/bob/airplanes_208... Can we pretend that airplanesIn the night sky ...

5 http://www.lyricsfreak.com/m/mike+posner/coole... If I could write you a song,And make you fall ...

6 http://www.lyricsfreak.com/j/jason+derulo/ridi... Yea yeah, yeah, yeah, yeah,I'm feeling like a ...

7 http://www.lyricsfreak.com/t/travie+mccoy/bill... I wanna be a billionaire so freakin' badBuy al...

8 http://www.lyricsfreak.com/d/drake/find+your+l... I'm more than just an option (hey, hey, hey) R...

9 http://www.lyricsfreak.com/u/usher/omg_2087748... Oh my goshBaby let meDid it again, so Imma let...

10 http://www.lyricsfreak.com/b/bob/magic_2087969... I got the magic in meEvery time I touch that t...

11 http://www.lyricsfreak.com/n/nicki+minaj/your+... [Chorus]Shawty I'm a only tell you this once, ...

12 http://www.lyricsfreak.com/m/maroon+5/misery_2... Oh yeahOh yeahSo scared of breaking itThat you...

13 http://www.lyricsfreak.com/t/train/hey+soul+si... Hey, hey, heyYour lipstick stains on the front...

14 http://www.lyricsfreak.com/b/bruno+mars/just+t... Oh her eyes, her eyesMake the stars look like ...

15 http://www.lyricsfreak.com/l/lady+gaga/alejand... I know that we are young,And I know that you m...

16 http://www.lyricsfreak.com/l/la+roux/bulletpro... Been there, done that, messed aroundI'm having...

17 http://www.lyricsfreak.com/f/flo+rida/club+can... You know I know howTo make 'em stop and stare ...

18 http://www.lyricsfreak.com/s/shontelle/impossi... I remember years agoSomeone told me I should t...

19 http://www.lyricsfreak.com/p/paramore/the+only... When I was youngerI saw my daddy cryAnd curse ...

20 http://www.lyricsfreak.com/u/usher/there+goes+... Yeah, Right,Usher baby, OKYeah man, rightThere...

Rによる分析

ここから、Rにてテキストマイニングを行いたいと思います。まず、tmパッケージを用いて、不要語（stop word）を除去します。具体的にはtheとかyouとかを除外しています。

#歌詞データの読み込み
lyrics_dataset <- read.csv(file = "lyrics_result.csv",as.is = TRUE,header = FALSE)
colnames(lyrics_dataset) <- c("link","lyrics")

library(tm)

#歌詞を小文字にする
document_dataset$lyrics <- tolower(document_dataset$lyrics)

#stop wordを除去する
stopwords_regex = paste(stopwords('en'), collapse = '\\b|\\b')
stopwords_regex = paste0('\\b', stopwords_regex, '\\b')
document_dataset$lyrics <- stringr::str_replace_all(document_dataset$lyrics, stopwords_regex, '')

#歌詞データの読み込み

lyrics_dataset <- read.csv(file = "lyrics_result.csv",as.is = TRUE,header = FALSE)

colnames(lyrics_dataset) <- c("link","lyrics")

library(tm)

#歌詞を小文字にする

document_dataset$lyrics <- tolower(document_dataset$lyrics)

#stop wordを除去する

stopwords_regex = paste(stopwords('en'), collapse = '\\b|\\b')

stopwords_regex = paste0('\\b', stopwords_regex, '\\b')

document_dataset$lyrics <- stringr::str_replace_all(document_dataset$lyrics, stopwords_regex, '')

続いて、LDAを実行できるtopicmodelsパッケージで扱えるようにするために、テキストデータに以下の処理を施します。

#文書単語行列の作成

corpus <- Corpus(VectorSource(document_dataset$lyrics))
inspect(corpus)
dtm <- DocumentTermMatrix(corpus)
findFreqTerms(dtm)

#文書単語行列の作成

corpus <- Corpus(VectorSource(document_dataset$lyrics))

inspect(corpus)

dtm <- DocumentTermMatrix(corpus)

findFreqTerms(dtm)

あとは以下のコードでLDAを実行するだけです。トピック数はアドホックに20としています。研究者の方、いい加減ですみません。

#LDAの実行

library(topicmodels)

nbo_topics <- 20
lda <- LDA(dtm,control=list(verbose=1), k = nbo_topics,method = "Gibbs")

#LDAの実行

library(topicmodels)

nbo_topics <- 20

lda <- LDA(dtm,control=list(verbose=1), k = nbo_topics,method = "Gibbs")

LDAの結果

まずは推定されたトピックごとの上位10単語をみてみます。トピック1はラブソングとかでしょうか。トピック17にパリピっぽい単語が、トピック18にスラングが含まれていますね。

#トピックの上位10単語を確認する
terms_each_topics <- data.frame(terms(lda,10))

1 2	#トピックの上位10単語を確認する terms_each_topics <- data.frame(terms(lda,10))

> terms_each_topics
   Topic.1 Topic.2 Topic.3 Topic.4 Topic.5 Topic.6 Topic.7 Topic.8 Topic.9 Topic.10 Topic.11   Topic.12 Topic.13 Topic.14 Topic.15 Topic.16 Topic.17 Topic.18 Topic.19　Topic.20
1     love    know    feel  better   wanna     hey     one  chorus    back     yeah    gonna      never     like      got      ’m     keep     stop      ain      let      get
2     like     now   heart   world    want    said     ooh     way    time     baby      low       will      new     like   don’t     even     just     like      can     good
3     make    just    life    whoa    take     old    call   verse    know     girl     hear      still     city     back     ’re  getting    hands     shit      say    night
4    touch    need    away     run    rock   every   cause     can    come   little     mean       ever     high    right    young      one      put     fuck      believe      ain
5     know   think    just   light     see     woo  gettin     pre    long     like    sound       eyes     bout     wish      ’s     give    party      got      fly     kind
6   nobody   cause   break  things    kiss    left    born    tell    like     just    shaky     always     ride     hold     ’ll     lose    crazy     hook      fall   really
7     baby    much    real    find    come    told     day     got    best    right     just everything     yeah     boom      que   please     live      gon      made    sleep
8    cause    give  enough    show    make nothing   makes    like    home      get     solo      leave      fun     know   ain’t   youand   lights    nigga      first champion
9     name    mind   every     see    body   daddy    came    used     til      can  tonight       hope     know      one     para     just     play   niggas      words    catch
10  loving  really    find waiting tonight   sweet   stand    made alright     look   wicked       stay      get     come   can’t  without      see    money      lonely      til

> terms_each_topics

Topic.1 Topic.2 Topic.3 Topic.4 Topic.5 Topic.6 Topic.7 Topic.8 Topic.9 Topic.10 Topic.11 Topic.12 Topic.13 Topic.14 Topic.15 Topic.16 Topic.17 Topic.18 Topic.19　Topic.20

1 love know feel better wanna hey one chorus back yeah gonna never like got ’m keep stop ain let get

2 like now heart world want said ooh way time baby low will new like don’t even just like can good

3 make just life whoa take old call verse know girl hear still city back ’re getting hands shit say night

4 touch need away run rock every cause can come little mean ever high right young one put fuck believe ain

5 know think just light see woo gettin pre long like sound eyes bout wish ’s give party got fly kind

6 nobody cause break things kiss left born tell like just shaky always ride hold ’ll lose crazy hook fall really

7 baby much real find come told day got best right just everything yeah boom que please live gon made sleep

8 cause give enough show make nothing makes like home get solo leave fun know ain’t youand lights nigga first champion

9 name mind every see body daddy came used til can tonight hope know one para just play niggas words catch

10 loving really find waiting tonight sweet stand made alright look wicked stay get come can’t without see money lonely til

見ずらいので、行を一つにまとめて、トピックにidを割り振ります。

#トピックの上位10単語をまとめて、トピックにidをふる
topic_keywords <- data.frame(topic_keywords_10=apply(t(terms_each_topics),1,paste,collapse=","))
topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

#トピックの上位10単語をまとめて、トピックにidをふる

topic_keywords <- data.frame(topic_keywords_10=apply(t(terms_each_topics),1,paste,collapse=","))

topic_keywords <- topic_keywords %>% mutate(topic_id=1:n())

> topic_keywords
                                              topic_keywords_10 topic_id
1       love,like,make,touch,know,nobody,baby,cause,name,loving        1
2          know,now,just,need,think,cause,much,give,mind,really        2
3        feel,heart,life,away,just,break,real,enough,every,find        3
4      better,world,whoa,run,light,things,find,show,see,waiting        4
5          wanna,want,take,rock,see,kiss,come,make,body,tonight        5
6          hey,said,old,every,woo,left,told,nothing,daddy,sweet        6
7           one,ooh,call,cause,gettin,born,day,makes,came,stand        7
8              chorus,way,verse,can,pre,tell,got,like,used,made        8
9           back,time,know,come,long,like,best,home,til,alright        9
10           yeah,baby,girl,little,like,just,right,get,can,look       10
11     gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked       11
12 never,will,still,ever,eyes,always,everything,leave,hope,stay       12
13               like,new,city,high,bout,ride,yeah,fun,know,get       13
14             got,like,back,right,wish,hold,boom,know,one,come       14
15        ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t       15
16   keep,even,getting,one,give,lose,please,youand,just,without       16
17         stop,just,hands,put,party,crazy,live,lights,play,see       17
18           ain,like,shit,fuck,got,hook,gon,nigga,niggas,money       18
19         let,can,say,believe,fly,fall,made,first,words,lonely       19
20      get,good,night,ain,kind,really,sleep,champion,catch,til       20

> topic_keywords

topic_keywords_10 topic_id

1 love,like,make,touch,know,nobody,baby,cause,name,loving 1

2 know,now,just,need,think,cause,much,give,mind,really 2

3 feel,heart,life,away,just,break,real,enough,every,find 3

4 better,world,whoa,run,light,things,find,show,see,waiting 4

5 wanna,want,take,rock,see,kiss,come,make,body,tonight 5

6 hey,said,old,every,woo,left,told,nothing,daddy,sweet 6

7 one,ooh,call,cause,gettin,born,day,makes,came,stand 7

8 chorus,way,verse,can,pre,tell,got,like,used,made 8

9 back,time,know,come,long,like,best,home,til,alright 9

10 yeah,baby,girl,little,like,just,right,get,can,look 10

11 gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked 11

12 never,will,still,ever,eyes,always,everything,leave,hope,stay 12

13 like,new,city,high,bout,ride,yeah,fun,know,get 13

14 got,like,back,right,wish,hold,boom,know,one,come 14

15 ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t 15

16 keep,even,getting,one,give,lose,please,youand,just,without 16

17 stop,just,hands,put,party,crazy,live,lights,play,see 17

18 ain,like,shit,fuck,got,hook,gon,nigga,niggas,money 18

19 let,can,say,believe,fly,fall,made,first,words,lonely 19

20 get,good,night,ain,kind,really,sleep,champion,catch,til 20

最後に、BillboardでTop10に入ったかどうかのデータを作っておき、そのデータと各歌詞を繋ぎ、各歌詞ごとに割りふられた確率が最大のトピックで集計をします。

#上位10位に入ったことがあるかどうかのダミーを作成
top_10_songs <- merge_dataset %>% group_by(link.y) %>% summarise(top_10=sum(top_10))
top_10_songs <- top_10_songs %>% mutate(top_10_dummy=ifelse(top_10>0,1,0))

#割り振られた最大の確率のトピックを抽出し、歌詞データと統合する
topics_each_document <- data.frame(topic_id=topics(lda,1))
topics_each_document <- cbind(link.y=document_dataset$link.y,topics_each_document)

#上位10位ダミーを繋ぎ、トピックの上位10単語の表現も繋ぐ
topics_each_document <- topics_each_document %>% left_join(top_10_songs,by="link.y")
topics_each_document <- topics_each_document %>% left_join(topic_keywords,by="topic_id")

#トピックごとのBillboardのTop10ランクイン割合をもとめる
topics_each_document %>% group_by(topic_keywords_10) %>% summarise(top_10_dummy=mean(top_10_dummy),count=n())

# A tibble: 20 × 3
                                              topic_keywords_10 `mean(top_10_dummy)` count
                                                         <fctr>                <dbl> <int>
1         ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t           0.12195122    41
2            ain,like,shit,fuck,got,hook,gon,nigga,niggas,money           0.14705882    68
3           back,time,know,come,long,like,best,home,til,alright           0.21666667    60
4      better,world,whoa,run,light,things,find,show,see,waiting           0.25000000    40
5              chorus,way,verse,can,pre,tell,got,like,used,made           0.14814815    54
6        feel,heart,life,away,just,break,real,enough,every,find           0.10416667    48
7       get,good,night,ain,kind,really,sleep,champion,catch,til           0.22857143    35
8      gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked           0.22857143    35
9              got,like,back,right,wish,hold,boom,know,one,come           0.13888889    36
10         hey,said,old,every,woo,left,told,nothing,daddy,sweet           0.09756098    41
11   keep,even,getting,one,give,lose,please,youand,just,without           0.18918919    37
12         know,now,just,need,think,cause,much,give,mind,really           0.21333333    75
13         let,can,say,believe,fly,fall,made,first,words,lonely           0.15217391    46
14               like,new,city,high,bout,ride,yeah,fun,know,get           0.14705882    34
15      love,like,make,touch,know,nobody,baby,cause,name,loving           0.17021277    47
16 never,will,still,ever,eyes,always,everything,leave,hope,stay           0.22222222    72
17          one,ooh,call,cause,gettin,born,day,makes,came,stand           0.28205128    39
18         stop,just,hands,put,party,crazy,live,lights,play,see           0.25000000    48
19         wanna,want,take,rock,see,kiss,come,make,body,tonight           0.10810811    37
20           yeah,baby,girl,little,like,just,right,get,can,look           0.14814815    54

#上位10位に入ったことがあるかどうかのダミーを作成

top_10_songs <- merge_dataset %>% group_by(link.y) %>% summarise(top_10=sum(top_10))

top_10_songs <- top_10_songs %>% mutate(top_10_dummy=ifelse(top_10>0,1,0))

#割り振られた最大の確率のトピックを抽出し、歌詞データと統合する

topics_each_document <- data.frame(topic_id=topics(lda,1))

topics_each_document <- cbind(link.y=document_dataset$link.y,topics_each_document)

#上位10位ダミーを繋ぎ、トピックの上位10単語の表現も繋ぐ

topics_each_document <- topics_each_document %>% left_join(top_10_songs,by="link.y")

topics_each_document <- topics_each_document %>% left_join(topic_keywords,by="topic_id")

#トピックごとのBillboardのTop10ランクイン割合をもとめる

topics_each_document %>% group_by(topic_keywords_10) %>% summarise(top_10_dummy=mean(top_10_dummy),count=n())

# A tibble: 20 × 3

topic_keywords_10 `mean(top_10_dummy)` count

1 ’m,don’t,’re,young,’s,’ll,que,ain’t,para,can’t 0.12195122 41

2 ain,like,shit,fuck,got,hook,gon,nigga,niggas,money 0.14705882 68

3 back,time,know,come,long,like,best,home,til,alright 0.21666667 60

4 better,world,whoa,run,light,things,find,show,see,waiting 0.25000000 40

5 chorus,way,verse,can,pre,tell,got,like,used,made 0.14814815 54

6 feel,heart,life,away,just,break,real,enough,every,find 0.10416667 48

7 get,good,night,ain,kind,really,sleep,champion,catch,til 0.22857143 35

8 gonna,low,hear,mean,sound,shaky,just,solo,tonight,wicked 0.22857143 35

9 got,like,back,right,wish,hold,boom,know,one,come 0.13888889 36

10 hey,said,old,every,woo,left,told,nothing,daddy,sweet 0.09756098 41

11 keep,even,getting,one,give,lose,please,youand,just,without 0.18918919 37

12 know,now,just,need,think,cause,much,give,mind,really 0.21333333 75

13 let,can,say,believe,fly,fall,made,first,words,lonely 0.15217391 46

14 like,new,city,high,bout,ride,yeah,fun,know,get 0.14705882 34

15 love,like,make,touch,know,nobody,baby,cause,name,loving 0.17021277 47

16 never,will,still,ever,eyes,always,everything,leave,hope,stay 0.22222222 72

17 one,ooh,call,cause,gettin,born,day,makes,came,stand 0.28205128 39

18 stop,just,hands,put,party,crazy,live,lights,play,see 0.25000000 48

19 wanna,want,take,rock,see,kiss,come,make,body,tonight 0.10810811 37

20 yeah,baby,girl,little,like,just,right,get,can,look 0.14814815 54

#トピックごとのランクイン割合を高い順に棒グラフで描写

library(ggplot2)
ggplot(topics_and_top10, aes(x=reorder(topic_keywords_10,top_10_dummy), y=top_10_dummy)) +
  geom_bar(stat='identity') + coord_flip() + xlab("topics")

#トピックごとのランクイン割合を高い順に棒グラフで描写

library(ggplot2)

ggplot(topics_and_top10, aes(x=reorder(topic_keywords_10,top_10_dummy), y=top_10_dummy)) +

geom_bar(stat='identity') + coord_flip() + xlab("topics")

BillboardのTop10ランクイン割合の高いトピックTop3
「one,ooh,call,cause,gettin,born,day,makes,came,stand」
「better,world,whoa,run,light,things,find,show,see,waiting」・・・明るい感じ？
「stop,just,hands,put,party,crazy,live,lights,play,see」・・・パリピぽい

BillboardのTop10ランクイン割合の低いトピックTop3
「wanna,want,take,rock,see,kiss,come,make,body,tonight」・・・欲求系？
「feel,heart,life,away,just,break,real,enough,every,find」・・・癒し系？
「hey,said,old,every,woo,left,told,nothing,daddy,sweet」

あまり洋楽を聴かないので、得られたトピックの解釈が中々できないのがもどかしいです。ただ、スラングの歌詞を含む歌詞はそんなにランクイン割合が悪いわけではなさそうですね。洋楽をもっと聴いて、前処理などもう少し工夫してリベンジしたいですね。

参考文献

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

Pythonクローリング&スクレイピング -データ収集・解析のための実践開発ガイド-

モダンなRによるテキスト解析 topicmodels: An R Package for Fitting Topic Models

某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた

はじめに

知人より、洋楽の流行りに疎いのでキャッチアップしたいという要望があり、某洋楽ヒットチャートの週次ランキングとTop100のデータを大量に集めてみようと思うに至りました。今回は深い考察を行うには至っていませんが、簡単にRにて集計・可視化を行います。

データ収集

Webスクレイピング対象の某洋楽ヒットチャートの週次ランキングは今週の順位・先週の順位・アーティスト名・曲名・詳細ページへのリンクなどが載せられおり、毎週土曜日更新されています。サイト内から導線はありませんが、URLのパラメータに法則があるため、うまく収集できます。今回は2010年8月〜2017年6月の約7年分のデータを集めます。

URLのリストをCSVで読み込み、BeautifulSoupでタグを指定して抽出します。

#ライブラリの読み込み
import urllib
from bs4 import BeautifulSoup
from urllib.request import urlopen
import requests
from urllib.error import HTTPError
import csv, time

#URLのリストを読み込む
f = open('url_list.csv', 'r')
dataReader = csv.reader(f)

#結果の出力用のリストを作る
data01 =[] #URL
data02 =[] #順位
data03 =[] #title&artist
data04 =[] #link
data05 =[] #artist

for row in dataReader:
       for url in row:
            time.sleep(10.0) #sleep(秒指定)
            try:
                    r = requests.get(url)
                    soup =  BeautifulSoup(r.content, 'html.parser')
                    
                    for body in soup.findAll("tbody"):
                        for detail in body.findAll("tr"):
                            for ranking in detail.findAll("td",{'class':'rank_td'}):
                                for content in detail.findAll("div",{'class':'name_detail'}):
                                    for artist_name in content.findAll("span"):
                                        for link in artist_name.findAll("a"):
                                            data01.append(url)
                                            data02.append(''.join(ranking.findAll(text=True)))
                                            data03.append(''.join(content.findAll(text=True)))
                                            data04.append(link.get("href"))
                                            data05.append(''.join(artist_name.findAll(text=True)))
                                            data = zip(data01,data02,data03,data04,data05)
                                            #CSV出力
                                            with open('result.csv','wt',errors='backslashreplace') as fout:
                                                writecsv = csv.writer(fout,lineterminator='\n')
                                                writecsv.writerows(data)                                   
                                
                                    
            except HTTPError as e:
                print(e.code)

#ライブラリの読み込み

import urllib

from bs4 import BeautifulSoup

from urllib.request import urlopen

import requests

from urllib.error import HTTPError

import csv, time

#URLのリストを読み込む

f = open('url_list.csv', 'r')

dataReader = csv.reader(f)

#結果の出力用のリストを作る

data01 =[] #URL

data02 =[] #順位

data03 =[] #title&artist

data04 =[] #link

data05 =[] #artist

for row in dataReader:

for url in row:

time.sleep(10.0) #sleep(秒指定)

try:

r = requests.get(url)

soup = BeautifulSoup(r.content, 'html.parser')

for body in soup.findAll("tbody"):

for detail in body.findAll("tr"):

for ranking in detail.findAll("td",{'class':'rank_td'}):

for content in detail.findAll("div",{'class':'name_detail'}):

for artist_name in content.findAll("span"):

for link in artist_name.findAll("a"):

data01.append(url)

data02.append(''.join(ranking.findAll(text=True)))

data03.append(''.join(content.findAll(text=True)))

data04.append(link.get("href"))

data05.append(''.join(artist_name.findAll(text=True)))

data = zip(data01,data02,data03,data04,data05)

#CSV出力

with open('result.csv','wt',errors='backslashreplace') as fout:

writecsv = csv.writer(fout,lineterminator='\n')

writecsv.writerows(data)

except HTTPError as e:

print(e.code)

データ取得後は簡単にpandasのstr.replaceで整形すると、以下のような結果になります。今週の順位と先週の順位が引っ付いてしまっています。

#Webスクレイピングしたデータをpandasで整形
import pandas as pd

#データ読み込み
data_list = pd.read_csv("result.csv", header = None, delimiter=",", encoding='utf-8')
data_list.columns = ['url', 'ranking','title&artist','link','artist']

#文字列の置換
data_list['url'] = data_list['url'].str.replace('置換前','置換後')
data_list['title&artist'] = data_list['title&artist'].str.replace('置換前','置換後')
data_list['ranking'] = data_list['ranking'].str.replace('置換前','置換後')

#CSVに保存
data_list.to_csv("dataset.csv",index=False)

#Webスクレイピングしたデータをpandasで整形

import pandas as pd

#データ読み込み

data_list = pd.read_csv("result.csv", header = None, delimiter=",", encoding='utf-8')

data_list.columns = ['url', 'ranking','title&artist','link','artist']

#文字列の置換

data_list['url'] = data_list['url'].str.replace('置換前','置換後')

data_list['title&artist'] = data_list['title&artist'].str.replace('置換前','置換後')

data_list['ranking'] = data_list['ranking'].str.replace('置換前','置換後')

#CSVに保存

data_list.to_csv("dataset.csv",index=False)

ここから横着してRで整形し、各週の順位データを作成しました。

#引っ付いたランキングに関するデータを分離して、データテーブルに加える。
dataset <- strsplit(ranking_dataset$ranking, "from")
rankings <- data.frame(current_ranking=as.integer(),
                       previous_ranking=as.integer())
for( i in 1:nrow(ranking_dataset)){
  rankings[i,1] <- as.integer(dataset[[i]][1])
  rankings[i,2] <- as.integer(dataset[[i]][2])
}

ranking_dataset <- cbind(ranking_dataset,rankings)

ranking_dataset <- ranking_dataset %>% select(url,current_ranking,previous_ranking,title.artist,artist,link)
colnames(ranking_dataset) <- c("date","current_ranking","previous_ranking","title",
                               "artist","link")

#引っ付いたランキングに関するデータを分離して、データテーブルに加える。

dataset <- strsplit(ranking_dataset$ranking, "from")

rankings <- data.frame(current_ranking=as.integer(),

previous_ranking=as.integer())

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

rankings[i,1] <- as.integer(dataset[[i]][1])

rankings[i,2] <- as.integer(dataset[[i]][2])

}

ranking_dataset <- cbind(ranking_dataset,rankings)

ranking_dataset <- ranking_dataset %>% select(url,current_ranking,previous_ranking,title.artist,artist,link)

colnames(ranking_dataset) <- c("date","current_ranking","previous_ranking","title",

"artist","link")

データ確認

データ構造はこのような形です。

#データ構造を確認する
> str(ranking_dataset)
'data.frame':	22254 obs. of  6 variables:
 $ date            : chr  "2010/08/21" "2010/08/21" "2010/08/21" "2010/08/21" ...
 $ current_ranking : int  1 2 3 4 5 6 7 8 10 11 ...
 $ previous_ranking: int  1 3 NA 2 5 4 6 9 11 8 ...
 $ title           : chr  "Love The Way You Lie" "Dynamite" "Mine" "California Gurls" ...
 $ artist          : chr  "Eminem Featuring Rihanna" "Taio Cruz" "Taylor Swift" "Katy Perry Featuring Snoop Dogg" ...
 $ link            : chr  "/artists/detail/306903" "/artists/detail/464754" "/artists/detail/319433" "/artists/detail/450934" ...

#データ構造を確認する

> str(ranking_dataset)

'data.frame': 22254 obs. of 6 variables:

$ date : chr "2010/08/21" "2010/08/21" "2010/08/21" "2010/08/21" ...

$ current_ranking : int 1 2 3 4 5 6 7 8 10 11 ...

$ previous_ranking: int 1 3 NA 2 5 4 6 9 11 8 ...

$ title : chr "Love The Way You Lie" "Dynamite" "Mine" "California Gurls" ...

$ artist : chr "Eminem Featuring Rihanna" "Taio Cruz" "Taylor Swift" "Katy Perry Featuring Snoop Dogg" ...

$ link : chr "/artists/detail/306903" "/artists/detail/464754" "/artists/detail/319433" "/artists/detail/450934" ...

まず、どんな楽曲やアーティストがランキングに入っているのかを簡単に確認してみます。ほとんど聞いたことない人の名前、曲名ですが。。

#楽曲ごとのランクイン回数の上位20人
> ranking_dataset %>% count(title) %>% arrange(desc(n)) %>% print(n=20)
# A tibble: 1,785 × 2
                          title     n
                          <chr> <int>
1                   Radioactive    91
2                          Sail    78
3             Party Rock Anthem    68
4                Counting Stars    67
5                       Animals    66
6           Rolling In The Deep    65
7                        Ho Hey    62
8                         Sorry    62
9  Somebody That I Used To Know    61
10                       Demons    60
11                    All Of Me    58
12                       Lights    57
13 Dark Horse Featuring Juicy J    56
14                  Some Nights    56
15                 Stay With Me    55
16                        Hello    54
17      Can't Stop The Feeling!    52
18     Cheap Thrills Featuring     52
19            Don't Let Me Down    52
20                      Pompeii    52
# ... with 1,765 more rows

#楽曲ごとのランクイン回数の上位20人

> ranking_dataset %>% count(title) %>% arrange(desc(n)) %>% print(n=20)

# A tibble: 1,785 × 2

title n

1 Radioactive 91

2 Sail 78

3 Party Rock Anthem 68

4 Counting Stars 67

5 Animals 66

6 Rolling In The Deep 65

7 Ho Hey 62

8 Sorry 62

9 Somebody That I Used To Know 61

10 Demons 60

11 All Of Me 58

12 Lights 57

13 Dark Horse Featuring Juicy J 56

14 Some Nights 56

15 Stay With Me 55

16 Hello 54

17 Can't Stop The Feeling! 52

18 Cheap Thrills Featuring 52

19 Don't Let Me Down 52

20 Pompeii 52

# ... with 1,765 more rows

#アーティストごとのランクイン回数の上位20位
> ranking_dataset %>% count(artist) %>% arrange(desc(n)) %>% print(n=20)
# A tibble: 876 × 2
             artist     n
              <chr> <int>
1             Drake   350
2        Bruno Mars   325
3           Rihanna   301
4      Taylor Swift   299
5             Adele   294
6   Imagine Dragons   258
7     One Direction   254
8        Katy Perry   237
9     Justin Bieber   204
10      Chris Brown   195
11 Carrie Underwood   192
12  Lady Antebellum   184
13      Nicki Minaj   183
14      OneRepublic   183
15   Ellie Goulding   182
16         Maroon 5   182
17       The Weeknd   182
18       Ed Sheeran   169
19    Blake Shelton   167
20          Beyonce   152
# ... with 856 more rows

#アーティストごとのランクイン回数の上位20位

> ranking_dataset %>% count(artist) %>% arrange(desc(n)) %>% print(n=20)

# A tibble: 876 × 2

artist n

1 Drake 350

2 Bruno Mars 325

3 Rihanna 301

4 Taylor Swift 299

5 Adele 294

6 Imagine Dragons 258

7 One Direction 254

8 Katy Perry 237

9 Justin Bieber 204

10 Chris Brown 195

11 Carrie Underwood 192

12 Lady Antebellum 184

13 Nicki Minaj 183

14 OneRepublic 183

15 Ellie Goulding 182

16 Maroon 5 182

17 The Weeknd 182

18 Ed Sheeran 169

19 Blake Shelton 167

20 Beyonce 152

# ... with 856 more rows

続いて、2010年8月〜2017年6月の間に100位以内に入った数を楽曲ごとにヒストグラムにしてみます。べき乗分布な形かと思いきや、20回前後で盛り上がっているのが気になりますね。

#タイトルごとのランクイン回数の集計値をヒストグラムに示す
title_ranking <- ranking_dataset %>% count(title) %>% arrange(desc(n))
ggplot(data = title_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title")

#タイトルごとのランクイン回数の集計値をヒストグラムに示す

title_ranking <- ranking_dataset %>% count(title) %>% arrange(desc(n))

ggplot(data = title_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title")

#楽曲のランクイン回数のサマリー
> summary(title_ranking$n)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    1.00    9.00   12.47   20.00   91.00

#楽曲のランクイン回数のサマリー

> summary(title_ranking$n)

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

1.00 1.00 9.00 12.47 20.00 91.00

中央値が9週間なので、意外と長い期間Top100には入っています。

続いて、100位以内に入った数をアーティストごとにヒストグラムにしてみます。こちらはべき乗分布のような形になっています。

#アーティストごとのランクイン回数の集計値をヒストグラムに示す
artist_ranking <- ranking_dataset %>% count(artist) %>% arrange(desc(n))
ggplot(data = artist_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Artist")

#アーティストごとのランクイン回数の集計値をヒストグラムに示す

artist_ranking <- ranking_dataset %>% count(artist) %>% arrange(desc(n))

ggplot(data = artist_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Artist")

Top10入りの楽曲の実態

Top10に入っている楽曲のみに絞って、ヒストグラムを描いてみます。

#Top10入りのデータに絞ってヒストグラムを描く
top_10 <- ranking_dataset %>% filter(current_ranking<=10) %>% count(title) %>% arrange(desc(n))
ggplot(data = top_10,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title in Top10")

#Top10入りのデータに絞ってヒストグラムを描く

top_10 <- ranking_dataset %>% filter(current_ranking<=10) %>% count(title) %>% arrange(desc(n))

ggplot(data = top_10,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title in Top10")

#Top10に入っている楽曲の10位以内ランクイン回数サマリー
> summary(top_10$n)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   3.000  10.000   9.779  15.000  32.000

#Top10に入っている楽曲の10位以内ランクイン回数サマリー

> summary(top_10$n)

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

1.000 3.000 10.000 9.779 15.000 32.000

Top10に入ったら、10週近くは10位以内に含まれるようです。上位はすぐに取って代わられるのかと思いきや、人気が人気を呼ぶとかなのでしょうか。確か、某ECサイトの方が、生キャラメルは売れるから売れたんだとか言っていた気がします。

順位の推移

100位以内にランクインした回数が最も多かった楽曲のTop10に関して、時系列プロットをしてみます。

#100位内ランクイン回数が上位の10曲に関して時系列プロットを行う。
title_ranking_top10 <- title_ranking[1:10,]
title_ranking_top10 <- ranking_dataset %>% filter(title %in% title_ranking_top10$title)
title_ranking_top10$date <- as.Date(title_ranking_top10$date)
ggplot(data = title_ranking_top10,aes(x = date,
                                      y = current_ranking,
                                      colour=title)) + geom_line() + 
  　　　　　　　　　　　　　　　　　　ylab("current ranking") + 
  　　　　　　　　　　　　　　　　　　ggtitle("Time Series Plot of Top10")

#100位内ランクイン回数が上位の10曲に関して時系列プロットを行う。

title_ranking_top10 <- title_ranking[1:10,]

title_ranking_top10 <- ranking_dataset %>% filter(title %in% title_ranking_top10$title)

title_ranking_top10$date <- as.Date(title_ranking_top10$date)

ggplot(data = title_ranking_top10,aes(x = date,

y = current_ranking,

colour=title)) + geom_line() +

　　　　　　　　　　　　　　　　　　ylab("current ranking") +

　　　　　　　　　　　　　　　　　　ggtitle("Time Series Plot of Top10")

初回に上位にランクインして、後は下がるだけの楽曲や、じわじわとランキングを上げていく楽曲などが観察されています。楽曲の消費のサイクルみたいなものがあるのでしょう。

今後について

せっかく面白そうなデータが手に入ったので、リンク先も辿って、どのような楽曲やアーティストの特徴が人気に繋がりうるのか見てみるのも良いですね。あと、もう少し洋楽聴いてみようと思います。私はクラシック音楽とジャズしかウォークマンに入っていないので。

参考文献

Pythonクローリング&スクレイピング -データ収集・解析のための実践開発ガイド-

ラーメン二郎の某飲食店レビューサイトデータに対して共分散構造分析をしてみる

データ

ラーメン二郎に関して、某飲食店レビューサイトのデータをWebスクレイピングしたもので、「料理・味」・「サービス」・「雰囲気」・「CP」・「酒・ドリンク」の評価項目に関して、1~5の実数値が割り振られています。ラーメン二郎の店舗数40店のうち、欠損のなかった37店舗の評価データとなります。

記述統計量・ヒストグラム

まずは記述統計量とヒストグラムを見てみます。

> describe(jirou_dataset[,-1])
                 vars  n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
taste               1 37 3.58 0.09   3.59    3.59 0.01 3.10 3.73  0.63 -3.77    17.97 0.02
service             2 37 3.51 0.11   3.55    3.53 0.04 3.10 3.59  0.49 -2.64     6.41 0.02
atmosphere          3 37 3.47 0.16   3.53    3.50 0.03 3.00 3.59  0.59 -2.11     3.06 0.03
cost_performance    4 37 3.61 0.12   3.59    3.61 0.01 3.07 3.95  0.88 -1.46     9.16 0.02
drink               5 37 3.08 0.09   3.07    3.07 0.10 2.99 3.33  0.34  1.00     0.39 0.01

> describe(jirou_dataset[,-1])

vars n mean sd median trimmed mad min max range skew kurtosis se

taste 1 37 3.58 0.09 3.59 3.59 0.01 3.10 3.73 0.63 -3.77 17.97 0.02

service 2 37 3.51 0.11 3.55 3.53 0.04 3.10 3.59 0.49 -2.64 6.41 0.02

atmosphere 3 37 3.47 0.16 3.53 3.50 0.03 3.00 3.59 0.59 -2.11 3.06 0.03

cost_performance 4 37 3.61 0.12 3.59 3.61 0.01 3.07 3.95 0.88 -1.46 9.16 0.02

drink 5 37 3.08 0.09 3.07 3.07 0.10 2.99 3.33 0.34 1.00 0.39 0.01

#各評価項目に関してヒストグラムを描く
hist_jirou_dataset <- melt(jirou_dataset)
g <- ggplot(data = hist_jirou_dataset,
            aes(x = value,
                y = ..density..)) +
  geom_histogram(alpha = 0.5,position = "identity") +
  geom_density(alpha = 0)
g + facet_wrap(~variable,nrow=5)

#各評価項目に関してヒストグラムを描く

hist_jirou_dataset <- melt(jirou_dataset)

g <- ggplot(data = hist_jirou_dataset,

aes(x = value,

y = ..density..)) +

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

geom_density(alpha = 0)

g + facet_wrap(~variable,nrow=5)

まぁ、お酒を飲むところではないので、drinkは低いですよね。しかしながら、総じて3点台なのは驚きです。昔行ったことがあって、雰囲気は決して良くはないはずなので。

共分散構造分析

今回、検証したい仮説は「ラーメン店としての質が二郎愛につながるかどうか」です。

ラーメン店としての質に繋がりそうな評価項目
「料理・味」
「CP」
「酒・ドリンク」

二郎愛に繋がりそうな評価項目
「料理・味」
「サービス」
「雰囲気」

加えて、料理・味とCPは関係していそうな項目なので、その点もパスにおいて考慮しておきます。

以下の図は仮説のイメージです。

パスを描写する際の細かい指定は参考文献を参照されると良いと思います。
さっそく、モデルを作って実行します。

#lavaanパッケージのためにモデルを描写する。
library(lavaan)

model_jirou <- 'f1 =~ cost_performance + drink + taste
                f2 =~ taste + atmosphere + service
                f2 ~ f1
                cost_performance ~~ taste'

fit <- sem(model = model_jirou,
           data = jirou_dataset[,-1],
           estimator="MLM")

summary(fit)

#lavaanパッケージのためにモデルを描写する。

library(lavaan)

model_jirou <- 'f1 =~ cost_performance + drink + taste

f2 =~ taste + atmosphere + service

f2 ~ f1

cost_performance ~~ taste'

fit <- sem(model = model_jirou,

data = jirou_dataset[,-1],

estimator="MLM")

summary(fit)

こちらが推定結果です。

> summary(fit)
lavaan (0.5-23.1097) converged normally after 100 iterations

  Number of observations                            37

  Estimator                                         ML      Robust
  Minimum Function Test Statistic                2.057       1.965
  Degrees of freedom                                 2           2
  P-value (Chi-square)                           0.358       0.374
  Scaling correction factor                                  1.047
    for the Satorra-Bentler correction

Parameter Estimates:

  Information                                 Expected
  Standard Errors                           Robust.sem

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  f1 =~                                               
    cost_performnc    1.000                           
    drink            -0.153    0.134   -1.139    0.255
    taste             0.381    0.194    1.958    0.050
  f2 =~                                               
    taste             1.000                           
    atmosphere        4.008    1.972    2.032    0.042
    service           3.147    1.423    2.211    0.027

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  f2 ~                                                
    f1                0.095    0.130    0.731    0.465

Covariances:
                      Estimate  Std.Err  z-value  P(>|z|)
 .cost_performance ~~                                    
   .taste               -0.001    0.007   -0.115    0.908

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .cost_performnc    3.606    0.020  179.760    0.000
   .drink             3.081    0.014  214.293    0.000
   .taste             3.582    0.015  237.571    0.000
   .atmosphere        3.471    0.027  129.316    0.000
   .service           3.509    0.019  187.760    0.000
    f1                0.000                           
   .f2                0.000                           

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .cost_performnc   -0.007    0.017   -0.390    0.696
   .drink             0.007    0.002    3.170    0.002
   .taste             0.002    0.003    0.712    0.477
   .atmosphere        0.005    0.004    1.273    0.203
   .service          -0.000    0.001   -0.329    0.742
    f1                0.021    0.017    1.242    0.214
   .f2                0.001    0.001    1.282    0.200

> summary(fit)

lavaan (0.5-23.1097) converged normally after 100 iterations

Number of observations 37

Estimator ML Robust

Minimum Function Test Statistic 2.057 1.965

Degrees of freedom 2 2

P-value (Chi-square) 0.358 0.374

Scaling correction factor 1.047

for the Satorra-Bentler correction

Parameter Estimates:

Information Expected

Standard Errors Robust.sem

Latent Variables:

Estimate Std.Err z-value P(>|z|)

f1 =~

cost_performnc 1.000

drink -0.153 0.134 -1.139 0.255

taste 0.381 0.194 1.958 0.050

f2 =~

taste 1.000

atmosphere 4.008 1.972 2.032 0.042

service 3.147 1.423 2.211 0.027

Regressions:

Estimate Std.Err z-value P(>|z|)

f2 ~

f1 0.095 0.130 0.731 0.465

Covariances:

Estimate Std.Err z-value P(>|z|)

.cost_performance ~~

.taste -0.001 0.007 -0.115 0.908

Intercepts:

Estimate Std.Err z-value P(>|z|)

.cost_performnc 3.606 0.020 179.760 0.000

.drink 3.081 0.014 214.293 0.000

.taste 3.582 0.015 237.571 0.000

.atmosphere 3.471 0.027 129.316 0.000

.service 3.509 0.019 187.760 0.000

f1 0.000

.f2 0.000

Variances:

Estimate Std.Err z-value P(>|z|)

.cost_performnc -0.007 0.017 -0.390 0.696

.drink 0.007 0.002 3.170 0.002

.taste 0.002 0.003 0.712 0.477

.atmosphere 0.005 0.004 1.273 0.203

.service -0.000 0.001 -0.329 0.742

f1 0.021 0.017 1.242 0.214

.f2 0.001 0.001 1.282 0.200

lavaan (0.5-23.1097) converged normally after 100 iterationsとあるので、適切に推定されたようです。自由度が2とデータ数が少ないためかなりギリギリな推定となっています。ただ、推定すべき母数の数よりもデータ数が一応多い状態ではあります。二郎インスパイア系の店のデータも集めた方がいいかもしれないですね。

モデルの評価

モデルの評価として適合度と母数の推定に関して見ていきます。

適合度
・適合度指標であるCFI(Comparative Fit Index)が1なので、適合度に関しては良さそうです。
・同じく適合度指標であるTLI(Tucker-Lewis Index)が0.998なので、適合度に関しては良さそうです。
・0.05以下であれば当てはまりが良いとされるRMSEAは0.028なので、当てはまりは良さそうです。
・0に近いほど良いとされるSRMRは0.024となっています。

母数の推定

f2 ~ f1はラーメン店としての質と二郎愛の関係を想定したものですが、p値が0.47と全然だめでした。ラーメン店としての質が二郎愛に繋がるという仮説は正しいとは言えないです。

参考文献にもあるよう、係数の解釈を行いやすくするために標準化推定値を求めます。

> standardizedsolution(fit)
                lhs op              rhs est.std     se      z pvalue
1                f1 =~ cost_performance   1.208  0.524  2.303  0.021
2                f1 =~            drink  -0.257  0.178 -1.448  0.148
3                f1 =~            taste   0.612  0.325  1.883  0.060
4                f2 =~            taste   0.400  0.120  3.327  0.001
5                f2 =~       atmosphere   0.901  0.066 13.723  0.000
6                f2 =~          service   1.016  0.049 20.694  0.000
7                f2  ~               f1   0.382  0.338  1.132  0.258
8  cost_performance ~~            taste  -0.201  1.651 -0.122  0.903
9  cost_performance ~~ cost_performance  -0.459  1.267 -0.362  0.717
10            drink ~~            drink   0.934  0.091 10.226  0.000
11            taste ~~            taste   0.278  0.399  0.698  0.485
12       atmosphere ~~       atmosphere   0.187  0.118  1.582  0.114
13          service ~~          service  -0.033  0.100 -0.333  0.740
14               f1 ~~               f1   1.000  0.000     NA     NA
15               f2 ~~               f2   0.854  0.258  3.309  0.001
16 cost_performance ~1                   29.960  8.666  3.457  0.001
17            drink ~1                   35.715  4.676  7.638  0.000
18            taste ~1                   39.595 15.318  2.585  0.010
19       atmosphere ~1                   21.553  4.333  4.974  0.000
20          service ~1                   31.293  7.949  3.937  0.000
21               f1 ~1                    0.000  0.000     NA     NA
22               f2 ~1                    0.000  0.000     NA     NA

> standardizedsolution(fit)

lhs op rhs est.std se z pvalue

1 f1 =~ cost_performance 1.208 0.524 2.303 0.021

2 f1 =~ drink -0.257 0.178 -1.448 0.148

3 f1 =~ taste 0.612 0.325 1.883 0.060

4 f2 =~ taste 0.400 0.120 3.327 0.001

5 f2 =~ atmosphere 0.901 0.066 13.723 0.000

6 f2 =~ service 1.016 0.049 20.694 0.000

7 f2 ~ f1 0.382 0.338 1.132 0.258

8 cost_performance ~~ taste -0.201 1.651 -0.122 0.903

9 cost_performance ~~ cost_performance -0.459 1.267 -0.362 0.717

10 drink ~~ drink 0.934 0.091 10.226 0.000

11 taste ~~ taste 0.278 0.399 0.698 0.485

12 atmosphere ~~ atmosphere 0.187 0.118 1.582 0.114

13 service ~~ service -0.033 0.100 -0.333 0.740

14 f1 ~~ f1 1.000 0.000 NA NA

15 f2 ~~ f2 0.854 0.258 3.309 0.001

16 cost_performance ~1 29.960 8.666 3.457 0.001

17 drink ~1 35.715 4.676 7.638 0.000

18 taste ~1 39.595 15.318 2.585 0.010

19 atmosphere ~1 21.553 4.333 4.974 0.000

20 service ~1 31.293 7.949 3.937 0.000

21 f1 ~1 0.000 0.000 NA NA

22 f2 ~1 0.000 0.000 NA NA

これによると、f1(ラーメン店としての質)に対してはコストパフォーマンスが最も影響を与えるようです。味よりもコストパフォーマンスが勝っているという考察になりますが、それはそれで面白いですね。他方、f2(二郎愛)に対しては味よりも雰囲気・サービスが影響を与えるようです。

semPlotパッケージを用いてパス図を出力してみます。

#パス図の出力
library(semPlot)

semPaths(fit,whatLabels = "stand",optimizeLatRes = T)

#パス図の出力

library(semPlot)

semPaths(fit,whatLabels = "stand",optimizeLatRes = T)

変数名が3文字に省略されているのですが、修正方法がパッと見つからなかったので、そのまま載せています。

初歩の初歩ですが、一通りの進め方がわかったので、今後も共分散構造分析にチャレンジしてみたいと思います。

参考文献

共分散構造分析 R編―構造方程式モデリング