参加できなかった第74回TokyoRのキャッチアップと結婚式について

自身の結婚式があったため、参加できなかった第74回目のTokyoRについてキャッチアップするために、公開資料を読んだ際に感じたメモをここに記します。
そして、世のデータサイエンティストが結婚式の際に苦しまないように参考になりそうな情報を少し書きました。

今回はTogetterも初めて作ってみました。(第74回R勉強会@東京(#TokyoR)のタグが付いたものたち

初心者セッション

初心者セッション1 – Data Import & Export –

資料なし

初心者セッション2 – Data Handling –

https://ymattu.github.io/TokyoR74/slide.html#/

いつもながら良い資料です。まだ社内でアクティブにRを広めてはいないですが、広めるならこの資料が良いですね。
上から下に直感的に書けるのは初学者には大事なので。
ただ、SQLの知識がないメンバーとかには補助教材が必要ですね。

登場するパッケージはdplyrとlubridateとstringrとforcatsとpurrr。

初心者セッション3 – Plot & Visualization –

資料なし

応用セッション

How LINE Corp Use R to Compete in a Data-Driven World

資料なし

LINEでのRの活用最前線の話なのでしょうが、資料がないのが寂しいです。
きっと、タイムラインで共有されていたこの記事をベースに話されているのではないか。

LINE の全社員が必要に応じて担当サービスのデータを分析できる環境の構築

  • 総勢50名の機械学習エンジニア・データサイエンティスト・データプランナー・データエンジニアが活躍できる分析基盤を作っていますよと。
  • Hadoop クラスタのデータの全社公開&活用のために、エンドユーザ向けのWebインターフェース「OASIS」をゼロから新規に開発。
    Spark, Spark SQL, PySpark, SparkRおよびPrestoクエリを叩くことができる。
    約20のサービス・部署で利用され、月間利用者数は約200人(データ関連の人以外で150人は触っていることになる。すごく層が厚そう。)

あと、ブリスベンの写真がTwitterで写り込んでいたので、UseR!2018のお話をされているのだろうと思われます。
linerパッケージについても語られている模様。

UseR!2018に参加し、社内Rパッケージ「liner」の活用事例を紹介しました

これのことでしょう。写真も合致している。
「いらすとや」の画像がワールドワイドに使われているのがシュールでいいですね。

「データの取得、分析、レポーティング、そして結果の共有にいたるまで、様々な便利機能を提供」と書かれています。
図を見る限りは、

  • RStudioでの解析結果をDBにカジュアルに保存したり、通知したり、ドキュメント化したりできる
  • PrestoやらHiveQLなどもRStudio上で実行できる(Tab補完とかもしてくれるんでしょうか?そこは聞いてみないとわからない。)
  • コーポレートカラーを適用したggplot2を利用できる
  • A/Bテストの結果をShinyのアプリで確認できる

などの機能があるようです。50人もいる高単価な人達の時間を節約できるという点でも、すごくインパクトのある取り組みですね。

tidyeval入門以前

(speakerdeckのembedに若干苦戦しましたw)

湯谷さんの考えるtidyevalについて英語で書かれています。
環境に応じた値をRは自動で引っ張ってくるけど、たまに干渉してしまうことがあり、実際dplyrとstatsは共にfilter()関数を持っていると。
どの優先順位をもたせるかのコントロールって難しそうですね。
干渉を避けるためにquosureやunquoteというのが説明されています。私の理解が追いついていないので、正直めちゃわかった感はないのですが、
様々な環境下で動くパッケージを作ろうと考える際は不可避な領域なのではないでしょうか。

LT

「うまい飯が作りたい」

recipesパッケージの紹介です。面白い方なんですね。
tidymodelsパッケージの中に内包されているそうな。

recipesパッケージの使い方を丁寧に説明してくださっています。

目的変数と説明変数を最初に明示的に指定し、前処理の手法をパイプ演算子でつないでいくだけ。
前処理のステップは使い回せるとのことで、似たようなデータを扱う場合は使いまわして楽をできるとのこと。
これは試してみる価値がありそうですね。k-nn法を用いた欠損値補完なども関数として用意されているようです。

マジレスすると、モテるかどうかは容姿や性格によるところが大きいと思うので、
身なりを清潔に保つとか、連れて行くと喜ばれそうな場所や体験を提供するとかが近道な気がしますよね。

不連続回帰とrdrobustパッケージの紹介

計量経済学系のバックグラウンドをお持ちのkagglerの方のLTです。
イスラム教の社会が、女性の高校修了率に与える影響について、統計的因果推論をされています。
rdrobustというパッケージを用いて、不連続回帰(RDD)という手法を実践され、イスラム教政治がランダム割当になりやすい状況を作っています。
全データで推定すると負の影響が推定されていましたが、RDDを使うことで正の効果があるという結果となりました。これは他の学部領域での先行研究とも整合的とのこと。

PCAや対応分析で補完要素を使う

資料なし

CiNii API その2

資料なし

Rで健康体

資料は後ほど公開される模様

データのみは公開されている


https://github.com/weda-654/my_health_log
0 forks.
0 stars.
0 open issues.
Recent commits:

Soccer × Attribution Analysis

u++さんのLTですね!アクセスログや第三者配信データを用いた分析でよく扱われるアトリビューション分析をサッカーの貢献度に利用したという話です。
ChannelAttributionパッケージを使われています。

マーケティングの業務で使う際は、これまで見逃していた意外なページや意外な参照元などがこの分析で見えてきたりします。
ただ、価値はあってもコストがかかっては元も子もないので、コストに関する記述があったのも実務で使われている方の視点だなと思いました。

サッカーはルールくらいしか知らないので誰が意外なのかはわからないですw

地理空間データの交差検証、正しくできていますか?

地理空間データにおける交差検証する際の手法として、Spacial Cross-ValidationとTarget-oriented cross-validationが挙げられています。
空間データ向けのパッケージとしては、sfやCASTが、学習周りではmlrやcaretが扱われていました。
Referenceがあるのが嬉しいですね。

reticulateパッケージとデータサイエンスフロー

資料なし

結婚式について

なぜ開くのか

大事な思い出づくりのため。

工程

  • Willing to Payの決定(全てはここ)
    • 持ち込みし放題の式場
    • 大学OB割引などがある式場
    • 料理が美味しい(国賓を迎えたりしているか?)
    • 荘厳な雰囲気
    • アクセスが良い
  • 教会式か人前式か神前式かの選定(教会式だとめっちゃ高かったりするところもある)
  • 誘う対象の選定
    • 共に勉強を頑張った仲
    • 共に仕事を頑張った仲
    • 親族
  • 住所の聞き出し
    • Googleフォームを活用
  • 大量の切手の購入(送付用と返信用も)
    • 郵便局は21時まで空いていたりするので助かった。(どこもそうなのかな?)
  • 上司や友人へのスピーチの依頼、乾杯の依頼
    • どういうオーディエンスなのかを事前に伝えておく。スピーチ作成者の負担を軽減する。
  • 招待状の作成、発送
  • 料理の試食
    • いろんな種類を食べるのでお腹ぱんぱん。
    • 料理は1000~2000円程度の予算アップはした方が面白い。
  • ウェディングケーキのデザインのすり合わせ
  • テーブルクロス、花、引き出物、ネームプレートなどのすり合わせ
  • 司会者との打ち合わせ
  • 自分のスピーチの作成(新郎新婦ともに)
    • 意外性と感動と笑いを織り交ぜるのが良い。
  • メッセージカードの作成
    • スプレッドシートに書きなぐり、それを手書きでひたすら書ききるのみ。
  • 式場音楽の選定
  • イベントの進行のすり合わせ
  • DVDの作成(オープニング・プロフィール・エンディング)
  • 大量のピン札の調達(交通費は全額支給しろという親の教え)
  • ウェルカムボードの作成(ダイソーで4~500円くらいで材料は手に入る)

工夫した点

  • DVDを3枚自作した
    • お願いすると結構お金がかかるし、理想形に近づけるための試行錯誤の回数も限られる。
    • ハイスペックPCを持っているデータサイエンティストなら動画の編集に耐えられるはず。
    • マックのiMovieを使えばGUIで簡単に動画作成ができる。
  • 十分な大きさのメッセージカード
    • わざわざ来てくれた友人との思い出をとにかく書きまくる。

結果としてよかったこと

  • DVDの自作
    • 3枚とも笑いを提供でき、上映後は拍手している人もいたようです。(裏手にいたので会場の様子は直接見れなかったですが)
  • 定型文を避けたスピーチ
    • 書く前によくある定型文を見たんですが、面白くないと思い、独自のスピーチにしました。結果として大爆笑を提供できました。
  • 普段から写真や動画を撮りまくっていたことで、DVD制作の素材が潤沢にあった。
    • Googleフォトに写真をアップしまくれば容量に悩まずに済むので、とにかく日常的に撮り続けましょう。

やはり、結婚式はある意味でエンターテインメントなので、自分が工夫できるところ(スピーチとムービー)は少し頑張ってみるといいのかな、と思いました。

反省点

  • 結婚式前日は有給休暇を取るべき
    • 仕事を19時であがったけど、その後の準備で疲れ果てた。目にクマが若干できた。

結婚を頑張るエンジニアやデータサイエンティストの皆さんへ

結婚式の準備は色々と時間がかかります。
世の優秀なデータサイエンティストが、そのようなことに時間を割きすぎるのは社会的な損失なので、この参考情報を元に少しでも楽に準備をしていただけると幸いですね。

rstanarmパッケージを使って簡単にベイズモデリングを実行する

はじめに

今回は、rstanarmというパッケージを用いて赤ワインデータを色々といじってみようと思います。
マーケティングの意思決定のための分析などでベイズ統計を使う場面が多々あるのですが、似たような属性のデータがあるのであれば、
一つ一つstanコードを書くのではなく、Rの関数でサクッと実行して試行錯誤していくという形に持っていけたらいいなぁと感じていました。
本気を出すところではstanを、ルーティンワーク的なタスクではrstanarmをみたいな形で使い分けれると良いのではないでしょうか。

rstanarmとは

バックエンドの計算をStanに実行させて、統計モデルの推定を行うためのパッケージ。R上でlm関数のように簡単にベイズ推定を行うことができる。対象ユーザーはベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー。
詳しくはこちら

インストールする

まずはrstanarmのインストールするのですが、コケまくりました。そのため、バージョンを下げてみることにします。

ここに過去のバージョンがありますが、2.17.4だと動かなかったものの、2.17.3なら動きました。

rstanarmのサンプルを回してみる

今回は、以下の文献を参考にして、大人のIrisとも言える、ワインデータを扱い、質の高いワインかどうかを決める要素を探ります。
How to Use the rstanarm Package | Jonah Gabry and Ben Goodrich

こちらの文献には、ベイズ分析の4つのステップとして以下があげられています。

  • 1.同時分布の特定(同時分布は事前分布と条件付きの尤度をかけ合わせたもの。)
  • 2.MCMCで事後分布を描く
  • 3.モデルがフィットしているか評価する
  • 4.事後予測分布を描き、結果に影響を与える予測項を確認する。

これらのステップをできるだけ素早くできると良いですね。

まずはデータを読み込んで、スケーリングしておきます。(可視化結果は前回と同じなので、載せません。)
加えて、6点以上の評価であれば1を取る二項変数を作成しておきます。

 

GLMでロジスティックモデルを推定し、rstanarmで推定した結果と比較します。rstanarmでは傾きや切片の事前分布にスチューデントのt分布を、尤度にロジスティック分布を設定しています。

ほとんど係数の大きさが同じであることが確認できます。

ベイズ推定の良いところは事後分布から関心のある係数に関しての取りうる値などをシミュレーションできるところですが、
posterior_interval関数で簡単に計算することができます。

肝心のMCMCの収束診断ですが、shinystanを使います。

やや余談ですが、他のデータセットでshinystanを用いた際に、予測結果にNAsが含まれている場合に、
shinystanが起動しないという問題があり、以下のようなエラー文が吐かれます。

調べたところ、こちらのgithubにあるように、

のように引数でppd=FALSEのように設定することで、立ち上げることができました。

3つの基準をクリアしているため、収束しています。

係数の分布についても可視化します。

rstanarmの良い点の一つとして、モデルのアップデートが容易に行える点があげられると思いますが、実際、以下のように先程のモデルに変数を追加して推定することができます。
今回は、alcoholを二乗したものを新しい変数として加えます。

次に、looパッケージを用いて、更新したモデルと元のモデルの性能の比較を行います。
looパッケージは統計モデルの予測精度の指標として扱われる、WAIC(Widely Applicable Information Criterion)を計算するためのパッケージで、WICは事後分布から得られる対数尤度の平均や分散からなる値として表されます。looはleave-one-out cross-validationのleave-one-outの頭文字。

さっそく入れようと思ったところ、

というエラーが出ました。
こちらでも議論されていましたが、

でバージョンを2.0.0から1.1.0に落としたら動きました。

ここで、事後分布が特定のサンプルデータに対して敏感であるかどうかをlooパッケージを用いて可視化します。

縦軸のshape parameter kは推定の信頼性の指標とされ、大きければ大きいほど信頼できないと見なし、横軸は今回推定したワインデータのデータの番号で、左が元のモデル、右が変数を追加したモデルのものです。
どうやらどちらも0.4未満のkに収まっているようです。参考情報の事例では0.5を超えていましたが、moderate outliersと説明されていたので、今回の推定は問題ないと思われます。

続いてモデルの比較を行います。

elpd_diffに関しては右のモデルの精度が高ければ正の値を、低ければ負の値を取るようになっています。標準誤差も返されます。
どうやら変数を追加したモデルの方が、ちょっとだけ良さそうです。

続いて、事後予測分布から、どの変数がどのように予測に影響を与えるのかを確かめます。
比較のためにデータを2つほど作成し、両者において一つだけ変数が違うという状況下での、予測される確率の比較を行います。

他の要素をある一定水準で保った際に、alcoholだけ1度下げることで、平均3%ほど高い評価が得られる確率が高まるという考察となります。

以上で、rstanarmの一連の使い方となるのですが、
一部の関数においては、階層ベイズモデルも行えるので、試してみようと思います。

ただ、階層ベイズにするにも、赤ワインのデータしかないので、グループ変数をどうにか作らないといけません。
あまりやりたくはありませんが、データがないので、説明変数を元にk-means(K=3)によるクラスタリングを行い、それをグループ変数とします。

stan_glmer関数を使えば、以下のような簡単な記述で定数項や係数がグループごとに異なるパラメータの分布に従うとする階層ベイズモデルを推定できます。

今回は、切片だけがグループごとに異なるモデル、傾きだけがグループごとに異なるモデル、切片も傾きも異なるモデルを作成します。

先程紹介した、looパッケージを使って、ベースとなるモデルとの比較を行います。

うーん、残念ながらどのモデルもベースモデルよりも圧倒的に強いものは無さそうです。

感想

まだまだrstanarmの関数やら機能やら定義を全て把握しきれていないですが、そこらへんがクリアーになれば、これまでのstanでの推定業務において生産性が高まる可能性を感じました。
簡単な階層ベイズモデルくらいなら、非常に直感的に書ける点や、変数の追加によるモデルのアップデートが容易な点などはポイント高いなぁと思います。
とはいえ、実務としてマニュアルでstanコードを作成していくのは必須なので、このパッケージを使うことによって、stanコードの改善に時間をより一層割けるようになるなら、それが一番だと思いました。
あと、「ベイズ初めてです!」という新入りの方とかには慣れ親しんでもらうには良さそうですね。lm関数レベルで実行できてしまうので。
今回、mc-stan.orgの配下にあるページなどを漁る過程で、ベイズ推定結果の可視化などで知らないことにも色々と出会えたので、今後も読み進めていきます。

追記

2018-09-10: Stanを使って変数選択したいにprojpredというパッケージが紹介されており、これを使えば、情報量基準に従った変数選択を簡単に行えるそうです。こうなると、「ベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー」に限らず多くの人が幸せになれるパッケージなのかもしれませんね。

参考情報

Using the loo package (version >= 2.0.0)
Leave-one-outクロスバリデーションの2つのデメリット、からの解決方法
stan_glmer | Bayesian Generalized Linear Models With Group-Specific Terms Via Stan
WAICを計算してみる
Package ‘bayesm’
Hierarchical Partial Pooling for Repeated Binary Trials
Leave-one-out cross-validation for non-factorizable models
priors | Prior Distributions And Options
StanとRでベイズ統計モデリング (Wonderful R)

Stanで順序プロビット(Ordered Probit)の推定のためのメモ書き

最近はBayesian Statistics and Marketingという本に関心があって、そこで取り上げられているモデルをStanに落とし込めないか模索しています。
そこで順序プロビット(Ordered Probit)の推定が必要であることがわかったため、Stanでの適用事例を漁っていました。まだマーケティング事例への適用はうまくいってないですが、いったん順序プロビットを簡単にまとめて今後の備忘録としておきます。

順序プロビットとは

被説明変数yが連続潜在変数y∗に対応していると考えるとする。
潜在変数は観察できないが、被説明変数yは観察でき、これらの2つ変数の関係は次のように表される。
(今回扱うデータは3から8までの順序データのため、以下のような表記になる。)

この対応関係は閾値メカニズムと呼ばれている。
各被説明変数をとる確率は以下のように記され、プロビットでは正規分布を扱うため以下のようになる。

これらの選択確率からなる尤度関数を最大にしたものが順序プロビットの推定となる。(c0=−∞でc6=∞とする。σは1とする。)

このように、潜在的な順序関係を想定し、それを満たすように閾値とパラメータを推定する点において、潜在変数を用いたモデルの柔軟性の高さが感じられる。

なぜ順序プロビットを使うのか

マーケティング業務において扱うデータにおいて、NPSやアンケートなど順序尺度の質的変数が多いので、それらのデータを二値データに落とし込んだり、そのまま基数データとして扱うのではなく、適切に扱いたいというモチベーションがあります。加えて、順序尺度の質的変数をもとに予測する際は普通のOLSだと、今回のケースで3を下回ったり、8を超えたりする可能性があり、予測結果として使いにくいです。
アンケートの点数をそのまま被説明変数として回帰しているケースは、データアナリティクスにこだわりの無いメンバーとかであればままあることなので、順序プロビットの民主化というか、布教していきたいと思います。

今回扱うデータ

勝手ながら大人のirisだと思っているワインデータです。今回は赤ワインに絞って、品質に関する順序変数を被説明変数として、各変数との相関を見ていきます。
まずはGGallyパッケージのggpairs()関数を適用して傾向を掴みます。見にくいので是非コードを回して確かめてください。

データに関する説明はワインの味(美味しさのグレード)は予測できるか?(1)で丁寧になされていますので、ご確認ください。

モデル

データセットに含まれる全部を含めて順序プロビットで回帰してみようと思います。
つまり、「酒石酸濃度、酢酸濃度、クエン酸濃度
残留糖分濃度、塩化ナトリウム濃度、遊離亜硫酸濃度
総亜硫酸濃度、密度、pH、硫酸カリウム濃度、アルコール度数」
の全てを使って赤ワインの質への影響を見ていきます。

Stanコード

最初に、Stanのユーザーガイド2.17の138ページにあるOrdered Probitのサンプルコードを使ってみたのですが、
収束しなかったので、初期値を設定するか弱情報事前分布を導入するかの判断が必要となります。
そこで、jabranhamさんが係数が平均0で分散10の正規分布に従うとするサンプルコードを書かれていたので、そちらを使って推定します。
書き換えているところはデータの制約くらいです。

StanをキックするためのRコード

推定結果の可視化を行うためのcommon.Rは松浦さんのGitHubにあるものになります。

結果

まず、MCMCが収束したかどうかの判断ですが、ShinyStanに従うものとします。

ShinyStanによる収束診断をクリアできています。

続いて、推定したパラメータです。

係数の符号がはっきりと分かれている、赤ワインの品質に影響を与えそうな変数としては、volatile acidity(酢酸濃度)、chlorides(塩化ナトリウム濃度)、total.sulfur.dioxide(総亜硫酸濃度)、sulphates(硫酸カリウム濃度)、alcohol(アルコール度数)のようです。

最後に、推定した閾値です。

1と2の閾値が近く、2と3の開きが大きく、あとは比較的均等のようです。

比較

質的変数をそのまま重回帰した際の結果ですが、符号やその大小はあまり変わらないです。やはり予測の際にどの順序尺度の値に対応するかがわかるのが使う利点だと思います。

考察

マーケティングにおいて、順序尺度の質的変数を扱う際に順序プロビットを積極的に使っていきたいと思いますが、アンケート分析を行う際は、ユーザーごとの評点の癖が点数に影響を与えている可能性があります。
そのため、点数の付き方がユーザーごとに違うとする階層モデルへの拡張を今後行っていくのが面白いと思いますし、実際に研究されている論文があります。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
stan-examples/limited-dv/oprobit.stan
stan-dev/stan users-guide-2.18.0.pdf
Stanによる順序ロジット回帰
第9章 順序選択モデル: 年金投資選択問題
Wine Quality Data Set
ワインの味(美味しさのグレード)は予測できるか?(1)
世界一簡単な収束[シナイ]Stanコード
RStanとShinyStanによるベイズ統計モデリング入門

おまけ

数式をブログに載せる際は、こちら
Online LaTeX Equation Editor – create, integrate and download
でインタラクティブに数式を作成し、その結果を
QuickLaTex Publish Math on the Web without compromising quality
に貼り付けて画像を出力しています。

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

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

データのおさらい

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

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

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

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

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

モデル

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

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

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

今回のモデルのDAG

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

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

stanコード

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

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

実行結果

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

K=1のケース

K=3のケース

K=5のケース

K=10のケース

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

おわりに

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

参考情報

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

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

はじめに

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

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

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

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

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

目的

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

データ

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

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

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

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

モデル

  • 家計ごとに異なる、マーガリン価格に対する反応を想定。各マーガリンのブランドの価格に対するパラメータの数は家計の数だけある。
  • 価格に対する反応は家計の属性によっても決まる。
    という前提に立ち、以下のセッティングで推論していきます。

  • 6つのブランド選択に関する多項ロジスティックモデル(カテゴリカル分布とsoftmax関数の適用)

  • 1階層目はブランドごとの価格を説明変数とし、価格に対する反応係数をかけ合わせたものを多項ロジスティックモデルの入力とする。
  • 2階層目はブランドの価格に対する反応係数が家計ごとの定数項と属性データに属性ごとの係数をかけ合わせたものからなる。
  • 家計ごとの定数項は平均0、分散V_betaの正規分布に従う。
  • 属性ごとの係数は平均vec(delta_bar)、分散V_betaクロネッカーのデルタA^(-1)の正規分布に従う。
  • 分散V_betaは平均υ、分散Vの逆ウィシャート分布に従う。
  • A = 0.01、υ = 6 + 3 = 9、V = υI(Iは単位行列)

コード

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

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

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

実行結果

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番目のブランドの係数が他のブランドに比べて小さいことがわかります。

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

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

多峰性などはなく、正規分布に従っているようです。他のブランドと比較して、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ライブラリを用いて行います。

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

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

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

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

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

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

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

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

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

Rでもやってみる

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

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

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

トピックの吐き出し

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

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

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

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

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

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

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

とりあえず、種族値合計(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版

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に変換)に、日時を整数に変更しています。

FMの実行

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

結果

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

こちらは、この中で性能の良かった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について調べてみた

References   [ + ]

1, 2, 3, 4, 5. predict(model,data.test)-target.test)^2

非負値行列因子分解(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年の新定番!ユーザーエクスプローラーをもっと活用するための簡単な方法

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

データ確認

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

行列因子分解(Matrix Factorization)

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

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

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

レコメンド結果の確認

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

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

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

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