[R]ボージョレ・ヌーボーのコメントに対してLDATSパッケージを使って時系列トピックモデルを扱う

はじめに

先日、某勉強会でLTをしました。その際に10秒だけ紹介したRのパッケージについて記事を書いてみようと思います。

LDATSパッケージについて

時系列でのトピックモデルを推定することができるパッケージです。
やっていることとしてはLDAでトピックを推定して次元を減らし、そのトピックの多変量時系列に関してベイズ手法による変化点検知のためのパラメータ推定を行っているようです。GitHubの該当しそうなソースコードに多変量のデータに対するsoftmax関数での回帰をやっているとの記述がある。(multinomial Bayesian Time Series analysis

元となっている論文を見る限り、BoW(Bag of Words)を想定して作っておらず、20~30程度のグループからなるデータに対して適用するのがちょうど良いです。アクセスログのページカテゴリや、マーケティングの顧客セグメントであればそんなに数は多くないので扱いやすいと思います。

データ

Webサイトから集めてきたボージョレ・ヌーボーのキャッチコピー14年分を今回は扱います。実は販売店側のキャッチコピーとワイン委員会が決めた評価が存在します。私の知っている世界は販売店側のキャッチコピーだけでした。

試してみた

今回はとにかく動くことだけを考えて、汚いコードとなっております。やっていることとしては、キャッチコピーを販売側とワイン委員会側のものを一つにつないで、数字を正規表現で「数字」に変換し、RMeCabで形態素解析をし、LDATS向けの形式のデータを作成していきます。
途中で、日本語の文字化け問題を回避するためにGoogle翻訳を使って単語名を置き換えています。
1時系列につき1文書となるようにデータを作っていく必要があるのですが、今回はボージョレ・ヌーボーのキャッチコピーなので最初から1時系列につき1文書となっているため都合が良いです。
データとソースコードはこちら

こちらは論文の図と同じものだとドキュメントの説明にあったので、論文の説明を見る限り、表すものとしては以下のようです。

  • 一番上の積み上げグラフはトピックごとの単語の割合を表しています。
  • 二番目の折れ線グラフはLDAによって推定されたトピックの時系列推移です。
  • 三番目のヒストグラムは二番目の時系列における変化点を集計したものです。
  • 四番目の折れ線グラフはモデルが推定したトピック割合の変化点の前後での推移です。

今回の図では文字が潰れていて見にくいですが、

  • トピック1はボキャブラリーが比較的リッチなコメント(「フルーティー」「フレグランス」「複雑」)
  • トピック2は数字を用いたコメント(「何年に一度の!」みたいな)
  • トピック3はボキャブラリーが貧相なコメント(「すごい!」みたいな)

のようです。
二番目の折れ線グラフを見る限り、周期的に数字を用いたコメントが現れているように思われます。四番目の折れ線グラフの変化点を見る限り、近年は数字を用いたコメントが相対的に減ってきて、リッチなボキャブラリーになってきているようです。

おわりに

時系列トピックモデルをカジュアルに試せる面白そうなパッケージだなと思い、LDATSパッケージを触ってみましたが、そもそもBoWなどを想定して作られているパッケージではないので、単語数が多いような分析ではそもそも可視化ができず使いにくいだろうなと思いました。マーケティングなどでユーザーのセグメントの推移を分析したい場合などにちょうど良いのだろうと思われます。

参考情報

[1] Long‐term community change through multiple rapid transitions in a desert rodent community
[2] Latent Dirichlet Allocation coupled with Bayesian Time Series analyses
[3] Package ‘LDATS’

2019年に読んだデータ分析系の本の振り返り(21+1冊)

はじめに

2020年、あけましておめでとうございます。年末に自分自身を振り返ろうと思ったのですが、結局データ分析と勉強しかしていないわけで、書籍を振り返ろうと思うに至りました。私の知り合いのデータサイエンティストはだいたい全冊持っているであろうと思われますが、良い本だと思うので思い出していただければ幸いです。

1.『ベイズモデリングの世界』(岩波書店)

基本的に階層ベイズモデルを使って、個体ごとの異質性を考慮した分析手法が提案されています。前半はオムニバス形式で様々な先生がモデルの適用について執筆されており、後半では伊庭先生による階層ベイズモデルの講義になっています。途中でスタイン統計量による縮小推定の話があげられ、柔軟なモデリングのためには「階層化した方が少なくとも望ましい推定量が得られる」という数学的証明を捨てることもやむを得ないと書かれています。

2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』(コロナ社)

この本はトピックモデルの教科書というよりも、ベイズ推定の教科書という側面が強い印象があります。途中で出てくる数式は流し読みするのは難しく、最低2冊以上のノートが別途必要になると思います。一度でもLDAのパラメータを導出してみたいという方には良い教科書だと思います。疑似コードが提供されているので、それをもとにRやPythonでコーディングしていけば、一番シンプルなLDAが非常に短い行で実行できてしまうことに驚かれるかもしれません。人間が手を動かして推定アルゴリズムを導出しているからこそ、短いコードで済むということを実感できるはずです。

3.『構造的因果モデルの基礎』(共立出版)

グラフィカルなアプローチで因果推論を扱っている書籍です。Judea Pearl流の因果推論アプローチについて記すことを目的に書かれています。基礎と書かれていますが決して簡単ではありません。ただ、扱われる数学のレベルとしては確率と線形代数がわかれば大丈夫だと思われます。余談ではありますが、1章の相関関係と因果関係の事例紹介で「おむつとビールの話」が都市伝説ではなくきちんと記事としてWall Street Journalという雑誌に掲載されていたことが明らかにされています。

4.『現場で使える!PyTorch開発入門 深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』(翔泳社)

PyTorchを触ったことがないが、深層学習の手法について知っている層を対象とした本です。6章まではGoogleのColabで動かせるのでGoogleに課金することなく深層学習による回帰、CNN、GAN、RNN、Encoder-Decoderモデル、ニューラル行列因子分解をPyTorchで試すことができます。写経したものはこちら。転移学習や高解像度化や画像生成、文章のクラス分類、文書生成、機械翻訳などもできるので、PyTorchでこれくらいの量をコーディングしたらこれくらいのことができるのかという学びや、他の人の書いたPyTorchコードを読みやすくなるなどの便益は十分にあると思いました。

5.『作ってわかる! アンサンブル学習アルゴリズム入門』(シーアンドアール研究所)

会社で行っているPythonもくもく会用に買った本で、scikit-learnを使わずに機械学習のアルゴリズム(アンサンブル系)をコーディングするための本です。pythonのコードについて逐次、細かい解説が行われているわけではないので、1行1行自分でコメントを加えながら写経をしていけば力が付くという本かなと思われます。sklearnはそれはそれで素晴らしいですが、こういう本でフルスクラッチで修行できるのはいいですね。

6.『数理統計学―基礎から学ぶデータ解析』(内田老鶴圃)

統計検定1級を合格された方のブログで紹介されていた教科書です。理系の大学生レベルの数学知識があれば、数理統計学の基礎を学べると思います。中心極限定理の証明や、様々な分布の期待値や分散、様々な分布の性質について数式を用いてしっかり理解することができます。数式もほどよく端折られているので、無論ですがノートが数冊必要になります。各章毎にある練習問題も解くことで力が付くと思います。日本の大学の授業の教科書がこれだったらジェノサイド(再履修者の大量発生)が起きるんだろうなと思ってしまった。

7.『44の例題で学ぶ統計的検定と推定の解き方』(オーム社)

統計の検定に関してだけ扱った珍しい本です。第3部までは統計学の普通の教科書ですが、それ以降であらゆる検定の例題が44件も載せられています。パラメトリックな検定から、ノンパラメトリックな検定まで幅広く扱われています。一番気にいっているのは仮説検定法の分類の表です。これさえあれば、どのデータに対してどの検定を行えばいいかが一目瞭然です。

8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』(技術評論社)

機械学習の原理を手早く数式を交えて学べる本です。かゆいところに手が届いていると言うか、既出の教科書では捨象されがちな、条件付き確率における2変数以上の条件づけでの表現に紙面を割いていたりしてくれるのが嬉しいです。ある程度数学の話はわかるが、だいぶ忘れているビジネスパーソンには大変にありがたいコンテンツと言えると思います。ベイズ線形回帰に関しても行列を用いた、わかりやすい導出方法が紹介されています。またコラムで紹介されている、測度論にどう向き合えばいいかの著者の見解は参考になります。

9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

R言語とstanを用いてベイズ統計学を入門レベルから学べる本です。各トピックごとにそれなりの紙面が割かれています。例題も豊富にあるので、線形回帰・MCMC・情報量基準・階層ベイズモデルまで、ベイズ統計学を基礎から応用までしっかりと学べると思います。youtubeで著者の講義も配信されているので、留学気分を味わえます。

10.『scikit-learnとTensorFlowによる実践機械学習』(オライリージャパン)

2019年に日本で開かれたML SummitでTFの開発者がおすすめしていた教科書です。前半部分で機械学習の入門から応用までをわかりやすい説明で学ぶことができます。数式は少ないですが、図とソースコード(Python)がちりばめられており、手を動かして理解を進めることができます。後半部分はTensorFlowを用いた深層学習の基礎を同様に手を動かして学ぶことができます。ただ、TFのバージョンも変わってきているので前半の説明をアテにして読むのも良いと思います。

11.『AIアルゴリズムマーケティング 自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)

マーケティングへのデータサイエンスの適用に関する珍しい書籍です。ソースコードはついていないですが、業務で使う際のアイデアが手に入ることもあります。一般的な回帰、生存時間分析、オークション、アトリビューション分析、アップリフトモデリング以外にも、情報検索やレコメンデーションやトピックモデルなどマーケティングながら学際的なトピックも扱われています。レコメンドなどで使われる、ランク学習に関して詳しく書かれた書籍をあまり知らないので、この本はその点においてもありがたい本でもあります。

12.『入門 統計的因果推論』(朝倉書店)

ほぼ全ての章でグラフィカルなアプローチで因果推論を扱っています。例題も豊富なので、一つ一つ丁寧にやれば理解が捗ります。おそらく、例題の多さを含め一番丁寧にd分離性、do演算子、バックドア基準、フロントドア基準に関する説明をしてくれている本なのかなと思いました。グラフでの因果推論に関して初めての人でも、確率さえ知っていれば読み進めることができるはずです。また、途中で操作変数法の紹介もされ、経済学出身者としては読みやすい。ただ、傾向スコアのくだりや、DIDなどのくだりはあまり出てきません。あと、やってないですが章末の練習問題に対するSolution Manualが提供されているようです。

13.『実践 ベイズモデリング -解析技法と認知モデル-』(朝倉書店)

ベイズモデリングを様々な事例に適用する方法がオムニバス形式で記された本です。ワイブル分布や異質性を考慮した二項分布、無制限複数選択形式のアンケートデータに対する手法、トピックモデル、項目反応理論などが扱われています。マーケティングの実務で使える事例が多いように感じました。こちらはサポートサイトでRコードとstanコードが提供されています。あと、appendixにあるプレート表現の見方も参考になります。

14.『機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習などで用いるベイズ推論を扱った教科書です。入門とありますが、入門者は書かれた数式をそのまま見ていても頭に入らないのではないでしょうか。手を動かしてなんぼの本だと思います。ノート2冊は絶対に必要です。たぶん、数式の展開を丁寧に記すと倍以上の厚みの本になると思います。各々のモデルに関してグラフィカルモデルが記されているのや、サンプルコードとしてGitHubにJuliaで書かれたソースコードが提供されているのも良いです。

15.『その問題、数理モデルが解決します』(ベレ出版)

物語形式で、様々な問題に対して数理モデリングのアプローチが紹介されています。途中でマッチング理論やゲーム理論やオークションなども登場することから、経済学出身者も喜ぶ内容かもしれません。社会人になってからナッシュ均衡という言葉が書かれた本は中々出会って来なかった。

16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』(CCCメディアハウス)

2013年と結構古い本ですが、データ分析を様々な事象に対して適用した事例紹介本です。アップリフトモデリングへの言及もあり、こういったものに関して日本は何年も遅れてブームが来るんだなという実感を与えてくれた本でもありました。appendixに分析事例が147個ほどあげられているのも参考になります。

17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』(北大路書房)

主にstanを用いたベイズモデリングによる分析事例が1と2で38本もオムニバス形式で載っています。ほとんどの事例で階層ベイズモデルが扱われています。2では若干マーケティングに近い内容の題材も扱われ、データサイエンティストの人にも嬉しい内容かもしれません。もちろんデータとstanとRのコードがサポートサイトで提供されています。

18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』(共立出版)

状態空間モデルで時系列予測を行うための手法が記されている本です。RのKFASパッケージが全面に渡って扱われています。トレンドを考慮したり、カレンダー効果を追加したり、共変量を追加したりなど様々なアプローチが紹介されコードも伴っているわけですから、業務でも抜群に役に立ちました。

19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』(SBクリエイティブ)

自分のいる会社で最低限の数学がわかると思われる若いメンバーに買ってもらうように言っている本です。微積分・線形代数だけでなく、カルシュ・キューン・タッカー条件(最適化数学)に関しても扱ってくれているので、ここで出てくる数学がわかれば大体の論文に立ち向かえると思います。さらに、Pythonの基礎もこれで学ぶことができるので一石二鳥な素敵な本ですね。また、最後の方でスクラッチでアルゴリズムを書くパートがあり、こちらも勉強になります。

20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

機械学習における前処理の指針を与えてくれる本です。Pythonのコードが提供されています。例えばですが、「テストデータにだけある、新しい単語は取り除いてしまえばいい」などの細かいアドバイスが何気に嬉しいです。「Effectコーディング」「特徴量ハッシング」「ビンカウンティング」「バックオフ」「leakage-proof統計量」などは読むまで知らないところだったので勉強になりました。

21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』(オライリージャパン)

データ分析の仕事をする上で最低限必要な知識を幅広く抑えることができる本です。数式は少なく、ところどころ出てくるコードはR言語です。参考文献などがブログだったりするため厳密さがめちゃあるわけではないですが、業務で使う分には問題ないと思います。分類問題において、AUCなどの評価指標だけでなく、予測値自体の探索的分析のすすめなどが書かれており、参考になりました。また、特徴量エンジンとしてのk-NN法の話も面白いと思いました。

[+α]『プログラマのためのGoogle Cloud Platform入門 サービスの全体像からクラウドネイティブアプリケーション構築まで』(翔泳社)

Google Cloud Platformを初めて触るデータ分析者にはちょうど良い本です。説明もわかりやすいので、いきなりアカウントを作ってドキュメントを解読するよりかは戸惑いは減るはずです。この本を土台に、GCS・GCEを駆使してML系のAPIを呼び出して使うなどの最低限の操作は私でもできるようになりました。GCPの画面や機能もどんどん変わっていくので書籍を買ってもアレなんですが、歴史的な背景も若干記述されているので、それはそれで勉強になります。ただ、エンジニアにこの本を買うべきか聞いた際にネガティブな意見があったのですが、たぶん現役プログラマからすると簡単過ぎるからなんだろうなと思います。

終わりに

2019年もぼちぼち勉強できましたが、2020年もこれまで同様にノートとペンを大事にする勉強を続けていき、コーディングも分析ももっともっと数をこなして会社や社会に求められるようなデータ分析官を目指していこうと思います。あぁ、英会話などの勉強をする時間を作るのが難しい。

[Stan]ロジスティック回帰の階層ベイズモデルとk-foldsクロスバリデーション

はじめに

stanは意思決定のための分析などでのパラメータ推定に使うことが多く、機械学習のために扱うことはありませんでした。ただ、過去にリク面などでお話したデータサイエンティストの方はstanで機械学習していて、クロスバリデーションもしているとの発言をされていました。
先日、記事を漁っていたらstanでクロスバリデーションを行うためのコードを書いている方を見つけたので、その方のコードをもとに大人のirisであるwineデータを用いて、質の高いワインかどうかを分類するために階層ベイズでのロジスティック回帰モデルを回してみたいと思います。

データについて

UCI Machine Learning Repositoryにある、赤ワインの評価と成分のデータです。データに関する説明はワインの味(美味しさのグレード)は予測できるか?(1)で丁寧になされていますので、ご確認ください。今回は6点以上であれば1を、そうでなければ0を取るものを教師データとしています。

分析方針

今回は階層ベイズモデルを扱うことから、グループごとにロジスティック回帰のパラメータが異なるという仮定を置きます。そのために、citric.acidというデータを3つのカテゴリデータに変換して、それをグループとして扱います。モデルでは、今回のデータセットの変数を全て回帰項として使います。もちろん、回帰用の式からはcitric.acidは除外します。
全体の80%を訓練データに、20%をテストデータとして、10foldsクロスバリデーションでのstanによる予測結果の平均AUCを評価指標とします。最後に、テストデータを用いた予測のAUCを確かめます。また、階層ベイズモデルではないモデルでの10foldsクロスバリデーションでのAUCとも比較します

分析概要

・データ整形
・訓練データとテストデータの分割
・クロスバリデーション用のデータの作成
・stanの実行
・クロスバリデーション結果の出力
・テストデータでの予測
・非階層モデルとの比較

全体のコード以下のリンクにあります。
kick_logistic_regression_allowing_k_hold_cross_validation_hierachical.R
logistic_regression_allowing_k_fold_cross_validation_hierachical.stan

データ整形

階層ベイズで扱うグループをcitric.acidから作っています。

訓練データとテストデータの分割

ワインの質に関するバイナリーデータをこちらで作成し、80%を訓練データに、20%をテストデータに分割しています。

クロスバリデーション用のデータの作成

こちらのコードでは任意の数でクロスバリデーション用のデータを作成し、stanで扱う訓練用データのlistに追加しています。
また、参考にしているブログより転用したstan_kfoldという関数を定義しています。k分割した際のstanの推定結果をリストに格納するための関数です。

stanの実行

こちらのstanのコードでは、M個のグループごとにパラメータが異なるというモデルを書いています。modelブロックの途中でholdoutを入れることで一部のデータを推定に使わないようにしています。

こちらはstanをキックするためのコードです。いつもと違い、先程定義したstan_kfoldを用いています。

クロスバリデーション結果の出力

以下は、k個ずつ手に入ったクロスバリデーションでの推定結果から、各パラメータの平均値を計算し、ロジスティック回帰モデルで2値の予測を行い、平均AUCを計算するコードです。

平均AUCは0.675となりました。すごくいいわけではないですが、手抜きモデルとしてはまずまずと言ったところでしょうか。

テストデータでの予測

以下のコードで最初に分けていたテストデータでの予測結果を返します。

実行の結果、AUCは0.665と、クロスバリデーションでの平均AUCと比べてあまり下がりませんでした。

非階層モデルとの比較

非階層モデルでも同様に10foldsクロスバリデーションの平均AUCを計算しました。非階層モデルよりもAUCが1%ポイントくらいは高いようです。

おわりに

現時点において、stanでの柔軟なモデリングを機械学習に活かす作法について紹介されている文献はあまりなく、選手人口はどれくらいいるのか気になるところですが、発見したブログのやり方でクロスバリデーションをカジュアルに行えるので、より多くの方がstanでの機械学習にチャレンジしうるものだなと思いました。ただ、このレベルの階層ベイズだとrstanarmで簡単にできてしまうので、より深く分析してモデルをカスタムしていきたいですね。

参考情報

[1]Lionel Hertzog (2018), “K-fold cross-validation in Stan,datascienceplus.com”
[2]Alex Pavlakis (2018), “Making Predictions from Stan models in R”, Medium
[3]Richard McElreath (2016), “Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)”, Chapman and Hall/CRC
[4]松浦 健太郎 (2016), 『StanとRでベイズ統計モデリング (Wonderful R)』, 共立出版
[5]馬場真哉 (2019), 『実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門』, 講談社

ベイジアン線形回帰モデルの式変形とRでのギブスサンプリングの適用

今回は特に目新しい手法というわけでもなく、線形回帰モデルのギブスサンプリングについて忘備録として残しておきたいと思います。
ベイジアン線形回帰モデルはプログラミング言語で言う、Hello World!的なものなので、あえてブログで取り上げる必要があると考えていないのですが、導出をしては忘れの繰り返しが嫌なので自分のために残しておこうと考えました。加えて、Stanのありがたみを感じられ、Stanへのコミットメントが増すのではないかとも期待しています。

・モデル
・数式の展開
・Rのコードの紹介
・おわりに
・参考情報

モデル

東北大学の照井教授の『ベイズモデリングによるマーケティング分析』に載せられている表記に従い、以下のように記します。

説明変数の数がk個の正規線形モデル

を考える。その場合、尤度関数は

となる。

係数パラメータの事前分布や条件付きの誤差分散の事前分布は以下のように設定する。(βは正規分布に従い、σ2|βは逆ガンマ分布に従う。)

数式の展開

私が大学院生だった時に、数式の展開をどう進めるかを手っ取り早く知る方法としては、「ネットに上がっている海外の大学院の講義資料を漁る」という作戦を取っていました。こうすることで数学のセンスがそれほど高くなくても、理解し進めることができました。今回に関してもおそらく、わかりやすく解説している海外の研究者がいるはずだと思い漁ってみたところ、コロンビア大学の機械学習の講義資料を見つけることができました。

資料はこちらのPDF(Course Notes for Bayesian Models for Machine Learning)で126ページにもなっていますが、導出のステップなどが非常に丁寧に書かれています。

それでは、今回の講義ノートを参考にしながら、線形モデルにおいて、ギブスサンプリングを行うところまでの導出を行いたいと思います。

まず、同時事後分布を以下の左辺のように置き、ベイズの定理を用いて右辺のように表記する。

次に、条件付き確率の定義と先程の尤度関数から以下のようになる。

yが与えられたもとでのp(y)は一定のため、比例している分子だけを残すと以下のようになる。

同時事後分布に事前分布の関数を代入していくと、

となる。両辺について対数を取ると、

となる。ここでβやσ2についての事前分布の形状から、同時事後分布におけるβやσ2について整理するための目標となる形状を確かめる。
まず、βはp(β)の定義より、対数を取りβについて整理すると、

となる。つまり、1/B0や1/B0・β0に該当する表現を先程の対数を取った同時事後分布から得ることを目標とする。
他方、σ2についても同様に、p(β|σ2)の定義より対数を取りσ2について整理すると、

となる。つまり、ν0やδ0に該当する表現を、同じく対数を取った同時事後分布から得ることを目標とする。

以上のパラメータごとの目標とする形状になるように各々のパラメータについて、対数を取った同時事後分布を整理する。

まずはβについてまとめ、関係のない項をconst.にする。

先程もとめた目標の形状を当てはめると以下のようになる。

よって、βの事後分布は以下のようになる。

他方、σ2についても同様に、関係のない項をconst.にし、目標の形状にまとめると以下のようになる。

目標の形状と比較すると以下のようになる。

よって、σ2の事後分布は以下のようになる。

Rのコードの紹介

条件付き事後分布からβやσ2の従う分布の形状がわかったので、それらを使ってRでギブスサンプリングを行います。先日、たまたま見つけた線形回帰モデルのギブスサンプリングのRのソースコードを拝借しようと思います。

ギブスサンプリングでは、先程導出した条件付き分布からβ→σ2と交互にサンプリングしていきます。それを記述したRコードは以下の通りです。

先程導出したβの事後分布である正規分布からのサンプリングの後(15~17行目)、そのサンプリングしたβを用いて、同じく先程導出したσ2の事後分布である逆ガンマ分布からサンプリングし(19~21行目)、それを指定した回数だけ繰り返し、所定の数まではバーンインとして除外します。(25~27行目)こうして導出した数式と、ギブスサンプリングのコードを見比べると理解が捗ると思いました。

実際に、先程のGitHubのソースコードを回してみると、以下のようにギブスサンプリングのイタレーションのプロットや、パラメータの事後分布を確認できます。

全体のコードはこちらです。

おわりに

シンプルなモデルですらこれだけ導出に手間がかかるということからも、Stanなどの確率的プログラミング言語のありがたみは非常に大きいなと思いました。こうして残すことで今後忘れたとしてもすぐに思い出せる気がします。
しかしながら、Stanでは事前分布と尤度を指定してしまえば、事後分布を計算し、知りたいパラメータについて解いた条件付き分布からサンプリングしてくれるわけですから、研究者の寿命を伸ばしたと言っても過言ではないかもしれません。

参考情報

[1]John Paisley (2016), “Course Notes for Bayesian Models for Machine Learning”, Columbia University
[2]照井伸彦 (2008), 『ベイズモデリングによるマーケティング分析』, 東京電機大学出版局
[3]須山敦志 (2017), 『機械学習スタートアップシリーズ ベイズ推論による機械学習入門』, 講談社
[4]stablemarkets,BayesianTutorials/MultipleLinearReg/multiplelinearreg.r

[Stan]生存時間分析のコードと便利なデータセットについて

はじめに

仕事で生存時間分析を使うことは結構あるのですが、マーケティングの良いデータセットがない印象でブログにしにくいと感じていました。また、Stanでの生存時間分析の事例もあまり把握していません。そこで使えそうなデータセットやStanのコードを探して、そのデータに対して生存時間分析を適用してみたいと思います。

目次
・生存時間分析とは
・生存時間分析で使えるデータ
・生存時間分析をマーケティングで使う際の用途
・先行研究
・生存時間分析で使えるデータセット
・Stanでの実行例
・おわりに
・参考文献

生存時間分析とは

生存時間分析は、ある時点から興味のあるイベント(マーケティングだと解約など)が発生するまでの期間を分析対象としています。データを手に入れた時点で、すでに解約して、真の累積の契約期間が判明している人と、解約しておらず今後いつ解約するかわからない中での累積の契約期間が残されている人のようなデータを扱うことが多いです。ここでの後者をcensoring(打ち切り)されたデータと呼びます。

生存時間分析をマーケティングで使う際の用途

ブログなどを読み漁る限り、以下の用途で生存時間分析を活用できるようです。

  • 顧客のサービス離脱率の予測、離脱原因の特定
  • 顧客がマーケティングキャンペーンに反応するまでの期間の長さ
  • 故障率の予測、故障原因の特定

先行研究

Stanを用いた分析事例は、調べた限りですが以下のモデルがあるようです。

  • 指数分布のモデル
  • Weibull(ワイブル)分布による比例ハザードモデル
  • ハザードの対数値についてのランダムウォークモデル
  • 2階差分のマルコフ場モデル(生存時間の確率分布は正規分布)
  • 1階差分のランダムウォークモデル(生存時間の確率分布は正規分布)

生存時間分析で使えるデータセット

事例を調べる過程で見つけた、生存時間分析に適したデータセットは以下の通りです。

どうやら、マーケティング、HR、離婚、再犯と幅広いデータがオープンソースで手に入るようです。

Stanでの実行例

今回は、「Using Customer Behavior Data to Improve Customer Retention」のデータセットを用いて、先行研究のソースコードにより生存時間分析をしてみようと思います。データは電話会社の顧客の解約に関するもので、様々な顧客の履歴データなどが用意されています。
先行研究のソースコードはWeibull分布を想定した比例ハザードモデルです。今回は決済の電子化の有無と離脱の関係を確かめてみます。なお、今回の打ち切りデータは契約期間となります。

まずはStanのコードはこちらです。Xobs_bgに説明変数が来るようにデータを用意しておきます。

続いて、このStanコードをキックするためのRのソースコードです。 元のデータが7043件と多すぎるのでランダムサンプリングしています。サンプリング数を8000、チェイン数を4にして実行します。(なお、可視化のソースコードもあるので結構長くなっていますので。頑張ってスクロールしてください。)

Rhatは全て1.05以下になっています。

traceplotを見る限り、重なり合っているので問題なさそうです。

各パラメータごとの自己相関係数です。こちらも問題なさそうです。

推定したパラメータの分布です。

横軸は推定した継続期間です。決済の電子化をしていない消費者は、契約期間の短い際の確率密度が低い傾向があるようです。

どうやら離脱率に関して決済の電子化をしていない消費者は、そうでない消費者よりも低いようです。

こちらは95%で取りうる範囲をそれぞれプロットしたものです。

おわりに

Stanで生存時間分析を行うという事例はそんなに多くはないものの、業界の長たちが良いコードを作成してくれていました。また、面白そうなデータセットも見つけることができました。このようなデータがもっと広まっていけば、マーケティングにおける生存時間分析がより活発に行われるのかもしれません。

参考文献

[1] 豊田秀樹 (2017) 『実践 ベイズモデリング -解析技法と認知モデル-』 朝倉書店
[2]生存時間解析入門
[3]比例ハザードモデルはとってもtricky!
[4]Stanで生存時間解析(Weibull 回帰)
[5]生存時間分析をStanで実行してみた
[6]階層ベイズ生存解析を用いたwebサイトの訪問者分析に関するStanでの実装
[7]生存時間分析 – ハザード関数に時間相関の制約を入れる
[8]Bayesian Survival Analysis 1: Weibull Model with Stan
[9]Bayesian Inference With Stan ~062~
[10]生存時間解析について – 概要編
[11]Survival Analysis for Employee Attrition ※kaggleで提供されているHR系のデータをサバイバル分析に用いている。
[12]Survival Analysis with R※Random Forests Modelによる生存時間の推定がなされている。
[13]Survival Analysis with R and Aster ※服役後の犯罪に関する分析や、離婚の分析などをしている。
[14]Survival Analysis of Mobile Prepaid Customers Using the Weibull Distribution(ダウンロード注意)

[Stan]項目反応理論(IRT)の段階反応モデルでbaysemのアンケートデータの分析をしてみる

はじめに

stanのユーザーガイドを見ていて、項目反応理論(IRT)についての章があり気になりました。勉強会のLTなどで手法の名前をちらっと聞いたことはあったのですが、使い道について調べていませんでした。ビジネスにおける実活用もしやすそうだと思ったので、カジュアルに分析して備忘録として残したいと思います。

目次
・項目反応理論(Item Response Theory:IRT)とは
・ビジネスでの適用可能性について
・データ
・モデルの推定
・結果の解釈
・おわりに

項目反応理論(Item Response Theory:IRT)とは

関西学院大学の教授のブログによると、

項目反応理論とは、テストについての計量モデルで、問題に対する正解・不正解のデータから、問題の特性や、回答者の学力を推定するためのモデルです。

とあります。また、Wikipediaによると、TOEFLの問題の評価のために使われているそうです。

主に、バイナリーと順序変数のモデルがあるようで、以下の母数がモデルに想定されています。どちらもほぼ同じです。

回答が2値変数のモデル

  • 2母数のロジスティックモデル
    • 特性値(例えば、広告配信の満足度とか)
    • 識別度母数(項目特性曲線の傾き)
    • 困難度母数(項目特性曲線の切片)
    • 定数

回答が順序変数のモデル(まずい < まぁまぁ < おいしい)

  • 段階反応モデル
    • 特性値(例えば、広告配信の満足度とか)
    • 識別度母数(項目特性曲線の傾き)
    • 困難度母数(項目特性曲線の切片)

※項目特性曲線は横軸に特性値、縦軸に質問の正答率を取ったものです。

ビジネスでの適用可能性について

  • 顧客のアンケート結果の解釈
    • 異質な集団間の得点を比較可能
    • 異なる尺度間の得点を比較可能(昔のアンケートだと5段階、今のアンケートは7段階などの状況はビジネスデータでありうる。)
  • 人事評価のバイアスの統制
    • 採用面接時の個人特性の正当な評価
  • アンケート項目の項目削減によるアンケートコストの低減
    • 各アンケート項目が理解されたかどうかを分析し、一つ一つのアンケート項目の精度を高める

データ

今回扱うデータはbaysemパッケージに入っているデータセットです。Yellow Pagesの広告プロダクトにおける満足度サーベイの回答データで、全ての回答は1から10のスケールで点数が付けられています(1がPoorで10がExcellent)。質問数は10個で、回答数は1811件です。

各質問の内容(baysemパッケージのドキュメントに載っていました。)
q1:全体の満足度

価格について
q2:競争的な価格設定
q3:昨年と同じ広告の最小値に対しての価格の引き上げ
q4:消費者の数に対しての適切な価格設定

効果について
q5:広告の購入の潜在的な影響
q6:広告を通じて自身のビジネスへの集客ができたか
q7:多くの消費者にリーチしたかどうか
q8:年間を通じて消費者に対する長期での露出があったか
q9:多くの家計やビジネスを必要としている人に届いたかどうか
q10:ビジネスを必要としている地理上のエリアに届いたかどうか

今回のIRT適用における特性値は、「広告プロダクトに関する満足度の傾向」としてみたいと思います。

モデルの推定

今回は教科書にならって以下の段階反応モデルを用います。

ここでaは識別力(広告の満足度が高まりやすいかどうか)、bは境界パラメータ(回答カテゴリ間の境界値)、θは特性(回答者がどれだけ広告に満足しているか)を表しています。Dは定数項で、以下では1とおいています。cはアンケートの回答のカテゴリ番号です。今回の例では10段階の評価が入ることになります。最後に、uは反応を、jは質問の番号を表しています。

実践 ベイズモデリング -解析技法と認知モデル-

こちらの本のサポートサイトからダウンロードできるzipファイルにstanのコードやRコードがありますので、そちらを利用しています。

モデルですが、以下のような設定となっています。

こちらをキックするためのRコードです。

処理時間としては、2014年末モデルのMacbook Proのcorei5、メモリ8GBで数時間程度でした。(正確な時間はわかりませんが、寝て起きたら計算が終わっていました。)


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


Rhatも1.1未満におさまっています。

結果の解釈

まず、推定した特性値の値のユーザーごとの平均値を求めて、ヒストグラムを描いてみます。どうやら、上限周辺にやたらと高い評価をしてそうなユーザーがいるようです。

最後に、項目特性曲線を質問ごとに、そして回答ごとに描いてみようと思います。

質問1~10に関して、10段階の回答ごとの項目反応曲線を以下に描いています。上まで戻るのが面倒なので、質問内容を再掲します。

q1:全体の満足度
q2:競争的な価格設定
q3:昨年と同じ広告の最小値に対しての価格の引き上げ
q4:消費者の数に対しての適切な価格設定
q5:広告の購入の潜在的な影響
q6:広告を通じて自身のビジネスへの集客ができたか
q7:多くの消費者にリーチしたかどうか
q8:年間を通じて消費者に対する長期での露出があったか
q9:多くの家計やビジネスを必要としている人に届いたかどうか
q10:ビジネスを必要としている地理上のエリアに届いたかどうか

これらの傾向から、9〜10点を獲得するにはある程度は特性値が高まる必要がある質問としては、q1〜q6のように見えます。価格や購買など自身のビジネスに直結しそうな質問が多い印象です。逆にふわっとした質問であるq7~q10は特性値が低くても9〜10点を取れる可能性が高い傾向があります。

おわりに

Stanのユーザーガイドを読むことで、普段自分が業務で扱っているアプローチなどが如何に限定的であることが実感できました。今回はIRTのアンケートデータへの適用事例を知れ、そこから様々な文献や便利なコードに至ることができました。社内のアンケートデータへの適用は面白そうだと思いますので業務で使ってみようと思います。

参考情報

[1] 豊田秀樹 (2017) 『実践 ベイズモデリング -解析技法と認知モデル-』 朝倉書店
[2] Yoshitake Takebayashi (2015) 「項目反応理論による尺度運用」 SlideShare
[3] 持主弓子・今城志保 (2011) 「IRTの組織サーベイへの応用」
[4] 清水裕士 (2017) 「項目反応理論をStanで実行するときのあれこれ」 Sunny side up!
[5] 清水裕士 (2016) 「Stanで多次元項目反応理論」 Sunny side up!
[6] 小杉考司 (2013) 「項目反応理論について」
[7] Daniel C. Furr et al. (2016) “Two-Parameter Logistic Item Response Model”
[8] daiki hojo (2018) “Bayesian Sushistical Modeling” Tokyo.R#70
[9] abrahamcow (2017) 「[RStan]項目反応理論の応用でフリースタイルダンジョン登場ラッパーの強さをランキングしてみた」

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コードです。