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)

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