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サイトに載っているポケモンの種族値の一覧をスクレイピング対象とさせていただきます。

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

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

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

■TTVランキング

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

■種族値のサマリー

■種族値のBoxプロット

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

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

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

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

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

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

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

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

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

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

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

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

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

■経験値のヒストグラム

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

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

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

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

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

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

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

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

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

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

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

■stanコード

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

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

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

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

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

■推定結果の可視化

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

まとめ

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

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

参考情報

データサイエンティストのための最新知識と実践 Rではじめよう! [モダン]なデータ分析
StanとRでベイズ統計モデリング (Wonderful R)
【R言語】rvestパッケージによるウェブスクレイピング その2
Receiving NAs when scraping links (href) with rvest

階層ベイズモデルの直帰率分析への適用 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件くらいにサンプリングしています。

推定

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

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

結果

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

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

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

参考文献

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

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

研究概要

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

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

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

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

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

データ可視化

以下、実験カテゴリごとの略記です。
下告知(20秒)・・・slow_20
上告知(15秒)・・・fast_15
超上告知(12秒)・・・fastest_12
平均告知(18秒)・・・average_18
対照群・・・baseline

データ構造の確認です。

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

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

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


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

RStanで重回帰

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

結果

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

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

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

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

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

参考文献

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

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

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

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

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

神戸大学

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

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

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

東京工業大学

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

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

Stanford University

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

Colombia University

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

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

目的

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

進め方

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

データの特徴

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

スクリーンショット 2016-03-13 22.30.16

stanコード

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

rコード

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

推定結果&可視化

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

posterior_distribution_2

参考文献

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

RstanでGoogleトレンドの推定

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

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

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

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

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

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

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

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

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

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

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

fig2-top-right
8期先までの予測です。

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

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

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

prediction_trend

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

参考文献

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

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

  • インストール関連
  • Windows 7にRStanをインストールする

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

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

    MCMCの計算にStanを使ってみた(超基礎・導入編)

    R stan導入公開版

  • 使い方
  • RStan Getting Started
    https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started

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

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

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

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

    Stanで統計モデリングを学ぶ(7): 時系列の「トレンド」を目視ではなくきちんと統計的に推定する

    Bayesian linear mixed models using Stan: A
    tutorial for psychologists, linguists, and cognitive
    scientists

    http://arxiv.org/pdf/1506.06201v1.pdf

  • 実践
  • R で 状態空間モデル: 状態空間時系列分析入門を {rstan} で再現したい

    書籍

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

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