【書評】PsychoPyでつくる心理学実験

本記事について

朝倉書店様と訳者の先生方より、2020年9月1日出版の、PsychoPyでつくる心理学実験をご恵投いただきました！

www.asakura.co.jp

PsychoPyの作者自身が筆者であるため、本当に素晴らしい教科書でした。まずPsychoPyについて簡単に紹介したのち、書評を記します。

※ただ書いてから気づきましたが、書評というより要約みたいになってしまいました。まあいいか。

そもそもPsychoPyとは？

University of NottinghamのJonathan Peirce氏が作成された、「心理学、神経科学、言語学のための多彩でダイナミックな実験を作成するための、オープンソース（無料）のソフトウェアパッケージ」です（本書より直接引用）。

端的にいえば、「心理学、神経科学、言語学でよく用いられる実験課題を無料で作成できるソフト」です。Pythonをベースとしているので、名前に「Py」が付いています。

PsychoPyのダウンロードはこちらから。

PsychoPyの作者は海外の研究者ですが、愛媛大学・十河宏行先生が日本語化や、日本語ドキュメントの執筆を行ってくださっています。

インストール方法や基本的な使用方法は、関西学院大学・小川洋和先生の解説が最も分かりやすいと思います。

これまでにも同様の目的を実現する商用ソフトはありました。しかしPsychoPyは、それらと遜色ない、いやむしろ今やそれら以上に充実した機能を備えたソフトであり、しかもタダ。作者自身が実験心理学者であるため、痒い所にちゃんと手が届く仕様。

特徴的なのは、以下の3通りの方法で実験課題を作成できること。

1. Builder（プログラミングのコードを一切書く必要がないGUIで、マウス操作が主体）
2. Coder（Pythonのコードを書く）
3. BuilderとCoderのハイブリッド（大半はBuilderで作成し、一部の機能をCoderで実現する）

このうち、Builderの操作が本当に分かりやすく、それでいて多機能なので、学部生向けの実験実習から、プロ研究者の実験まで、幅広く用いることができます。僕も学部生向けの実験実習で長年使用していますが、上記の関西学院大学・小川洋和先生の解説に沿って実演するだけで、30分程度で基本的な機能は習得してもらえます。

僕自身でも、長年にわたりBuilderで実験課題を作成してきましたが、最近はちょっと複雑な実験課題を組むことがあるため、柔軟性の高いCoderに移行しています。PsychoPyはPythonベースのソフトと書きましたが、上述のように「パッケージ」なので、心理実験等でよく組み込まれる処理（例：刺激を200msec呈示）を行うためのPythonコードが、関数として多数用意されています。

十河先生が執筆された「心理学実験プログラミング: Python/PsychoPyによる実験作成・データ処理」を丁寧に写経すれば、Python初心者でも十分にCoderで実験課題を作成できます。

www.asakura.co.jp

この記事も非常におすすめです。

qiita.com

書評

ではここから書評を記します。僕が本書の特長だと思った点は、以下の4点です。

1. Builder操作に関する解説が主体
2. いくつかの代表的な実験課題（例：視覚探索課題）を例に挙げながら、機能を解説
3. 実験心理学の教科書（←これが良い！）
4. 外部機器との連携

1つずつ説明していきます。

1. Builder操作に関する解説

上で、「Coderのほうが柔軟性が高いから、僕自身では、最近はCoderを使っている」と書きました。え、じゃあBuilderって、GUIで簡単に操作できるとはいえ、できることが限られるんじゃないの？　Coderの解説書の方が需要が高いんじゃないの？　と思うかもしれません。

その疑問に対する答えは、「Yes、でもそれでいいっ！　それがBEST！」です（リサリサ先生の声で）。

Builderの設計思想は、GUIで任意の実験課題を作れるようにすることではありません。Builderは、GUIで「標準的な実験課題」を作れるように用意されたもの。Builderには多くのメニューが用意されていますが、意図的に機能を網羅していません。なぜなら、多すぎるメニューは利便性を低下させるし、何より初心者を混乱させるからです。幅広いユーザーが、簡便に標準的な実験課題を作れるように調整されたインタフェースが、Builderなのです。

したがって、初心者は難なく使用方法を習得できるし、プロ研究者であっても標準的な実験課題を用いる場合にはBuilderを使って何の問題もありません。むしろ素早く実験を作成できるのだから、その方がいいじゃないか。デバグも楽になるし。

あ、熱く語りましたが、これは本書に記載されていることです。僕はここを読んだとき、ハートが震えてヒートが燃えつきました（いや燃えつきたらあかん）。PsychoPyの作者自身が執筆した書籍なので、こういう作者の思いが垣間見えるのがいいですね。

ちなみに、後述もしますが、PsychoPyはオンライン実験にも対応しています。PsychoPyで作成した実験課題を専用のサーバにアップロードすることで、インターネットブラウザを介して実験を受けられるようになるということです。本記事の執筆時点（2020/8/18）では、オンライン実験用のプログラムは、原則としてCoderを使用せず、Builderで作成しなければならないと思います。そういう意味でも、Builderで課題を作成する方法の解説には意義があります。

※なお本書では、オンライン実験のための手続きは掲載されていません

2. 代表的な実験課題（例：視覚探索課題）を例に挙げながら、機能を解説

PsychoPyが見据えている研究領域（心理学、神経科学、言語学）の特徴は、「こういうデータを測定したかったら、こういう実験がよく用いられるよね」という実験パラダイムが存在することかと思います。

例えば抑制能力（特定の情報処理を我慢して、別の情報処理を行う）を測定したかったら、ストループ課題が例に挙がると思います。ストループ課題とは、「あか」「あお」のように、単語の意味とインクの色が一致する刺激（あか）と、一致しない刺激（あお）を呈示し、インクの色を回答させるという実験です（「あか」でも「あお」でも、インクの色は赤）。

あるいは、特定の情報を探し出す速さを測定したかったら、視覚探索課題が用いられると思います。例えば画面上に多くの「L」が散りばめられており、その中に「T」が1つだけ紛れ込んでいるか否かを判断するという実験です。

本書では、このような実験課題を例に挙げながら、Builderの機能を広く解説しています。汎用的な機能も学びつつ、特定の実験課題を作り上げられるので、一石二鳥ですね。

しかも、冒頭で書いたような、「作者自身が実験心理学者であるため、痒い所にちゃんと手が届く仕様」の紹介も随所に存在しています。

実験では、刺激を呈示して、それに対する被験者の反応を測定することが基本ですが、それだけすれば十分というわけではありません。ときに、測定精度を高めるための”工夫”を交えることがあります。

例えば、通常は本番の実験の前に練習を行い、被験者に実験手続きを理解していただくことが必要になります。被験者が実験手続きを理解していなかったら、本番中に測定されるデータ（例：反応の速度）は、どのような情報処理を反映しているのか分かりにくくなるからです（例：反応時間が長くても、それは刺激が処理しにくかったのか、反応方法が分からなかったのか、区別できない）。

そこで、例えば練習試行の最後の5試行ぶんだけ抽出し、反応時間や反応の正確性を確認したいというモチベーションが生まれます。そこで極端に反応が遅かったり、正答率が低かったりしなければ、被験者が実験手続きを理解したと推測することができます（もちろん、これだけで十分という意味ではありませんが）。

じゃあPsychoPyでどう実装するのか。ちゃんと実例が紹介されているんですね、これが。

このように、基本的な機能の紹介にとどまらず、「実践的テクニック」も随所に盛り込まれているのが、素晴らしい点だと思います。

3. 実験心理学の教科書

さて、僕が何より感動したのが、この点です。

ここまで、「PsychoPyを使えば簡単に実験課題が作れるよ！　実例も載ってるよ！」と、”楽ができる”ことを推してきました。しかし一般的に、「実験の準備」という手続き自体は、楽ではありません。

実験を実施する前に、丁寧に時間をかけて確認しなければならないこと、調整しなければならないことが、たくさんあります。そのために知らなければならないことが、たくさんあります。

PsychoPyは、「作る」という作業においてユーザーに楽を与えてくれます。しかし、それは諸手続きを軽視してよいことを意味しません。

例えば、現在ではPCの液晶ディスプレイを用いて実験刺激の呈示を行うことが多いと思います。一般的には、画面は60Hzで更新されることが多いでしょう（= 1秒間に60回、画面が更新される）。いわゆる、リフレッシュレートというやつです。

画面に、0.2秒だけ刺激を呈示したいとします。これは可能です。1回の更新が1/60秒なのだから、12*1/60 = 0.2であり、12回画面を更新すればよいからです。しかし同様の理屈により、画面に0.22秒だけ正確に刺激を呈示することは困難です。13.2*1/60 = 0.22なので、画面を13.2回更新すればよいわけですが、そんな半端なことはできません。

PsychoPy上で、刺激を0.22秒呈示すると設定することはできます。しかし、それは実際に刺激が0.22秒呈示されていることを意味しません。

というわけで、環境（例：使用するディスプレイ）によって実験手続きに制約が実は存在しているわけです。え...そんなの知らんかった...言うといてや...。安心してください、書いてますよ。

そもそも、液晶ディスプレイって実験で使っていいのだろうか。あれ、そういえば昔受けた実験だと、ブラウン管（CRT）モニターを使用していた気がするぞ。そもそも液晶ディスプレイとCRTモニターの違いはどうなっているのだろう？　安心してください、書いてますよ。

他にも、気を払わなければならない（かもしれない）ポイントはたくさんあります。刺激の色や輝度を厳密にコントロールしたい場合もあるでしょう。特定の情報を探す速さを測定する視覚探索課題の場合、当たり前ですが、物理的に目立つ刺激は見つけられやすくなります。こういった事情を気にする場合、キャリブレーションやガンマ補正などのキーワードが頭に浮かびます。安心してください、書いてますよ。

刺激の呈示ルールはどうすればいいだろう。とりあえずランダムな順序で呈示されるようにしておけばいいかな。でもランダム呈示って本当に万能なのかなあ。安心してください、書いてますよ。

被験者が特定の情報を認識できる／できない境界線（閾値）を知りたい。テクニックがあるのかなあ。そういえば階段法とか聞いたことがあるな...。安心してください、書いてますよ。

と、このように、実は本書は実験心理学の教科書としても非常に有用です。PsychoPyの作者自身が実験心理学者である強みが、まさにここに表れています。巻末には、三角関数の簡単な紹介も載っていて（刺激を特定の位置に呈示するためには、sinやcosの計算が必要だから）、本当に...痒い所に手が届く...。

4. 外部機器との連携

僕自身がこの用途を体験したことがないので、ここはさらっとした紹介になってしまいますが、脳活動関係（fMRI、EEG）の測定をする場合や、眼球運動の計測をする場合における、外部機器との連携が求められる発展的な機能も解説されています。これらの測定やデータの解析と相性の良い商用ソフトもありますが、フリーのPsychoPyでも実現できるのは嬉しいところだと思います。

結び

ここまでレビューしてきたように、本書は、初心者から熟練者まで、幅広いユーザーが待ち望んでいた「実験心理学の教科書」なのではないかと思います。

PsychoPyの作者自身が執筆した書籍であることによる充実度はもちろんのこと、訳者も全員、実験心理学を専門とし精力的に活躍する研究者であるため、翻訳の精度も高いと思います（まだ隅々まで精読したわけではありませんが、ざっと読んだ限りでは、非常に読みやすいと思いました）。

皆様もぜひ！　Enjoy！

Windows上でzipファイルを解凍するときの、「ファイル名文字化け」への対処

本記事について

自分用メモです。

僕は論文ファイルを保存するとき、ファイル名を論文の概要にする習慣があります（この記事に感銘を受けました）。

たまった論文をフォルダごと移動させようと思って、ZIPにした後で解凍したところ、ことごとく日本語のファイル名が文字化けしました。どうやらWindowsだとこういうことが起きやすいみたいですね（理由）。

楽な解決方法が見つかったので、記録しておきます。

前提

（最近の）Windows10ではデフォルトで用意されている、Windows Subsystem for Linux（WSL）の機能を用いるなどして、Windowsマシン上でUbuntuを動かせる環境にあることを前提とします。この環境構築方法については、以下の記事が参考になります。

kscscr.com

Ubuntu上でzipファイルの解凍

あとはもうこっちのもんです。

Ubuntu上でunarコマンドを実行すれば、解凍時に、文字コードを自動検出して良きに計らってくれると。

qiita.com

まずはunarをインストールして、次にunarコマンドで解凍するだけ。

sudo apt-get install unar

unar ファイル名.zip

a-zumi.net

補足

WSLとWindows間でファイルをやり取りする場合は、この方法で。

一番楽なのは、解凍対象のzipファイルを

\\wsl$\Ubuntu-(バージョン)\home\(ユーザー名)\

のディレクトリに入れた後で、上記のunarコマンドを実行することかな。

qiita.com

実験参加者内要因の条件間でパラメータを比較する（パラメータ間の相関を考慮した比較）

はじめに

本記事は、Stan Advent Calendar 2019　22日目の記事です。

本記事の内容はたのしいベイズモデリング　事例で拓く研究のフロンティア第15章の内容とほぼ同じですが、ちょっとだけ発展的な内容を含んでいます。

次のような手続きで実験を行ったとしましょう。PCのディスプレイ上に、「O」と「D」がランダムに呈示されます。実験参加者は、

• 「O」が出たと思ったら右のボタンを押す
• 「D」が出たと思ったら左のボタンを押す

という風に反応することを求められます。刺激が呈示されてから、実験参加者がボタンを押すまでの反応時間が記録されます。一見簡単に見えますが、これをずーーーっと続けていると、時折間違えた反応をしてしまうこともあるでしょう。また、反応が徐々に遅くなっていくと予想されます。すなわち反応の正確性や反応時間は、「情報処理が円滑に行われている程度」を反映する指標と考えられます。

モデリングの方針

このように、実験参加者に二通りの反応が許されているとき（今回は、左右いずれかのボタンを押す）、その結果として得られる反応時間のデータを、drift diffusion modelによってモデリングしてみましょう。そのようなモデリングを行う利点については、Stan Advent Calendar 2019　16日目の記事を参照してください。

本記事の目的

複数の条件を、異なる実験参加者グループが経験している場合（実験参加者「間」計画）と、同じ実験参加者が経験している場合（実験参加者「内」計画）で、モデリングの仕方が変わります。

実験参加者「内」計画の場合、同じ実験参加者がすべての条件を経験するので、条件同士で反応に相関関係が認められると予想されます。例えば反応時間を測定する実験の場合、慎重な人はすべての条件で反応が遅いかもしれないし、せっかちな人はすべての条件で反応が早いかもしれません。

本記事ではこのような実験参加者内計画における反応時間のデータを、条件間の相関を考慮したうえでモデリングする方法を紹介します。

実践

サンプルデータの作成

drift diffusion modelでは、反応時間がWiener分布に従うと考えます。そこで、Wiener分布に従う乱数を発生させましょう。brmsパッケージには、brms::rwiener()という関数が用意されています。これは内部ではRWienerというパッケージを利用しています。もし以下のコードを実行したとき「RWienerというパッケージは無い」というエラーが出た場合には、追加でRWienerパッケージをインストールしておいてください。

# load packages ----------------------------
library(tidyverse)
library(brms)
library(rstan)

# sample data ------------------------------
trial = 1000 # 一つの条件における試行数
N_subject = 5 # 実験参加者数
tau = 0.5 # Wiener分布の非決定時間パラメータ
beta = 0.5 # Wiener分布の開始点パラメータ

set.seed(123)

# 以下は順に、実験参加者1の条件1、条件2、実験参加者2の条件1、...、実験参加者5の条件2
dat_all = rbind(brms::rwiener(n = trial, alpha = 0.6, tau = tau, beta = beta, delta = 2),
                brms::rwiener(n = trial, alpha = 0.75, tau = tau, beta = beta, delta = 2.5),
                brms::rwiener(n = trial, alpha = 0.5, tau = tau, beta = beta, delta = 1.7),
                brms::rwiener(n = trial, alpha = 0.55, tau = tau, beta = beta, delta = 2.4),
                brms::rwiener(n = trial, alpha = 0.65, tau = tau, beta = beta, delta = 2.2),
                brms::rwiener(n = trial, alpha = 0.8, tau = tau, beta = beta, delta = 2.5),
                brms::rwiener(n = trial, alpha = 0.5, tau = tau, beta = beta, delta = 1.8),
                brms::rwiener(n = trial, alpha = 0.6, tau = tau, beta = beta, delta = 2.2),
                brms::rwiener(n = trial, alpha = 0.7, tau = tau, beta = beta, delta = 1.6),
                brms::rwiener(n = trial, alpha = 0.75, tau = tau, beta = beta, delta = 2.1)) %>% 
  dplyr::mutate(ID = rep(x = 1:N_subject, each = 2 * trial),
                condition = rep(1:2, each = trial, time = N_subject))

Wiener分布は4つのパラメータを持ちます。このうち、非決定時間tauと、情報蓄積の開始点betaについては、今回はすべて共通としています。

反応が起こるために必要な情報蓄積の閾値alphaと、情報蓄積の速さdeltaが、実験参加者や条件によって異なるサンプルデータです。このサンプルデータには、以下のように4つの変数があります。

• q
  • 反応時間（秒）
• resp
  • 右のボタンを押したか（1）、左のボタンを押したか（0）
• ID
  • 実験参加者のID。5人いるので、1~5。
• condition
  • 「O」条件（1）、「D」条件（2）

head(dat_all)

          q resp ID condition
1 0.5569206    1  1         1
2 0.6966415    0  1         1
3 0.5378940    1  1         1
4 0.5503652    1  1         1
5 0.6447536    1  1         1
6 0.7017261    0  1         1

tail(dat_all)
              q resp ID condition
9995  0.6605005    0  5         2
9996  0.6017658    1  5         2
9997  0.6058292    1  5         2
9998  0.5311541    1  5         2
9999  0.7406046    1  5         2
10000 0.5725912    0  5         2

Stanコード

Stanにはwiener_lpdf()という分布関数が用意されているので、基本的にはここに目的変数とパラメータを書き入れるだけです。なお16日目の記事とはやや違う書き方をしています。

パラメータ名は以下のように名付けなおしています。

• 閾値alpha　→　boundary parameter a
• 非決定時間tau　→　t_raw及びt_new
  • これら二つの違いは後述
• 情報蓄積の速さdelta　→　drift rate parameter d

data {
  int<lower=0> N;                         //全データ数
  int<lower=0> N_ID;                      //被験者数
  int<lower=0> N_cond;                    //被験者内要因の条件数
  
  int<lower=0> ID[N];                     //被験者を識別するIDが格納されたベクトル
  int<lower=0> cond[N];                   //被験者内条件を識別するラベルが格納されたベクトル
  
  vector<lower=0>[N] RT;                  //反応時間ベクトル
  real<lower=0> minRT;                    //全被験者のなかでの最小RT（推定されたtがこれを超えるとエラーになるため）
  
  int Resp[N];                            //反応のベクトル
  
  real<lower=0, upper = 1> beta;         //starting point (0.5に固定)
}

parameters {
  vector[N_cond] a[N_ID];                 //各条件のboundary parameterを計算する材料。被験者内条件数の長さで、そのベクトルが被験者数ある。
  vector[N_cond] mean_mu_a;               //各条件のboundary parameterが従う、多変量正規分布の平均ベクトル
  corr_matrix[N_cond] omega_mu_a;         //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列を作るための、相関行列
  vector<lower=0>[N_cond] tau_mu_a;       //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列を作るための、標準偏差ベクトル
  
  real t_raw;                             //全被験者に共通するnon-decision time parameterを求めるための材料
  
  vector[N_cond] d[N_ID];                 //各条件のdrift rate parameterを計算する材料。被験者内条件数の長さで、そのベクトルが被験者数ある。
  vector[N_cond] mean_mu_d;               //各条件のdrift rate parameterが従う、多変量正規分布の平均ベクトル
  corr_matrix[N_cond] omega_mu_d;         //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列を作るための、相関行列
  vector<lower=0>[N_cond] tau_mu_d;       //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列を作るための、標準偏差ベクトル
}          

transformed parameters{
  real<lower=0, upper=1> t_new;           //全被験者に共通するnon-decision time parameter
  cov_matrix[N_cond] mu_a_covmatrix;      //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列
  cov_matrix[N_cond] mu_d_covmatrix;      //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列
  
  t_new = inv_logit(t_raw) * minRT;       //ロジスティック変換で0～1に変換し、それに最小RTを掛けているため、
                                          //Wiener分布から推定するt_newは最小RTを超えないように工夫
  
  mu_a_covmatrix = quad_form_diag(omega_mu_a, tau_mu_a); //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列
  mu_d_covmatrix = quad_form_diag(omega_mu_d, tau_mu_d); //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列
}

model {
  // model ----------------------------------------------------
  for(n in 1:N){
    if(Resp[n] == 1){
      target += wiener_lpdf(RT[n] | exp(a[ID[n], cond[n]]), t_new, beta, exp(d[ID[n], cond[n]])); //exp(a_raw), exp(v_raw)により、両parameterが多変量「対数」正規分布から生成されたことを表現
    }
    else{
      target += wiener_lpdf(RT[n] | exp(a[ID[n], cond[n]]), t_new, beta, -exp(d[ID[n], cond[n]]));
    }
  }
  
  // boundary parameter
  target += multi_normal_lpdf(a | mean_mu_a, mu_a_covmatrix);
  
  // drift rate parameter
  target += multi_normal_lpdf(d | mean_mu_d, mu_d_covmatrix);
  
  
  // priors -----------------------------------------------------
  target += normal_lpdf(mean_mu_a | 0, 5);
  target += lkj_corr_lpdf(omega_mu_a | 2);
  target += student_t_lpdf(tau_mu_a | 3, 0, 3) - student_t_lccdf(0 | 3, 0, 3);
  
  target += normal_lpdf(t_raw | 0, 5);
  
  target += normal_lpdf(mean_mu_d | 0, 5);
  target += lkj_corr_lpdf(omega_mu_d | 2);
  target += student_t_lpdf(tau_mu_d | 3, 0, 3) - student_t_lccdf(0 | 3, 0, 3);
}

ポイントは以下の点です。

1. wiener_lpdf()のなかで、aとdのパラメータに指数変換を施している
   • これらのパラメータは正の値を取ります。ということは、両パラメータにlower=0と下限を指定するだけでもいいような気がしますが、そうすると収束が不安定になりがちでした。
   • aやd自体には下限を指定せず、exp()によって正の値に変換しています。
2. 非決定時間パラメータ
   • 反応時間の最小値を超えてしまうと推定できなくなるので、それを抑える工夫です。詳しくはコード内のコメントを見てください。
3. 各条件のaとdのパラメータが、多変量正規分布に従うと仮定している
   • 今回は2条件しかないので、実質的に2変量正規分布。
   • 多変量正規分布は、分散共分散行列をパラメータに持つので、ここで条件間の相関を表現できる
   • なお今回は、aとdを最終的に指数変換しているため、結果的にこれらのパラメータは多変量対数正規分布に従うというモデルになっている

1.と2.は推定効率を高めるための工夫ですが、3.が本記事の目的である、相関のあるパラメータを推定するために必要な部分です。

実験参加者「内」条件間でパラメータを比較したいとき、この多変量（対数）正規分布の平均ベクトルを比較することが出来ると思います。今回のモデルではmean_mu_aやmean_mu_dです。

もっとも、そういったやり方ではなく、条件間でパラメータが等しいというモデルと、パラメータが異なることを仮定するモデルを、例えばベイズファクター等を用いてモデル比較するという方法もあると思います。

推定

Rコードは以下の通りです。自己相関が高かったので、thin = 10としています。もうちょっとモデルを工夫できるのかもしれません。

# model -------------------
sm = rstan::stan_model("ddm.stan")

# data --------------------
sdata = list(N = nrow(dat_all),
             N_ID = dat_all %>% dplyr::distinct(ID) %>% nrow(),
             N_cond = dat_all %>% dplyr::distinct(condition) %>% nrow(),
             ID = dat_all$ID,
             cond = dat_all$condition,
             RT = dat_all$q,
             Resp = dat_all$resp,
             minRT = min(dat_all$q),
             z = 0.5)

# sampling --------------------
fit = rstan::sampling(object = sm,
                      data = sdata,
                      seed = 1234,
                      iter = 8000,
                      warmup = 3000,
                      thin = 10)

結果

パラメータが多いので、ここでは結果の要約は載せず、パラメータリカバリできたかどうかに注目します。推定されたaやdを指数変換したものが、Wiener分布のパラメータであることに注意してください。

直線上に点が載っていたら、パラメータの真値と推定結果（パラメータの事後期待値）が一致していることを意味します。閾値aについてはうまくリカバリできたようです。

しかし情報蓄積の速さdについては、真値に近い値は得られているものの、ちょっと誤差があります。このパラメータの推定はもうちょっと工夫が要るのかもしれません。

おわりに

今回は2肢選択事態における反応時間データを、drift diffusion modelによってモデリングしました。さらにその際に、実験参加者内要因、つまり同じ実験参加者が全ての実験条件を経験したことにより、条件間でパラメータに相関が仮定される状況を取り上げました。

実験心理学ではこのような実験計画を組むことが多いので、参考になれば幸いです。

Enjoy!

R初心者の館（RとRStudioのインストール、初期設定、基本的な記法など）

本記事について

R Advent Calendar 2019　2日目の記事です。

本記事執筆のモチベーション

ゼミや講義でRを使いたいことがあります。しかし、インストールや初期設定、基本的な記法についての説明で時間を使ってしまうのはもったいないと思い、「これを事前に読んできて」と言えば済むような資料を用意したいと思いました（もちろんすでに、ネット上には有用な記事がたくさんあります）。もし同様の要望をお持ちの方がいらっしゃったら、本記事をご活用いただければ幸いです。  

そういうわけで、本記事では、Rをまったく触ったことがない初心者を読者に想定しています。また、筆者の環境がWindowsであるため、同環境を事例として説明しています。

目次

• RとRStudioのインストール
• RStudioの初期設定
• RStudioの機能
• Rの基本的な記法・使用方法
  • 四則演算
  • 代入演算子
  • その他の演算子
  • データの型・構造
  • ベクトルの作成
  • matrixの作成
  • data.frameの作成
  • パッケージのインストールと読み込み
  • 外部ファイルの読み込み

RとRStudioのインストール

Rとは

公式ページの説明の通り、Rはフリーの、オープンソースのプログラミング言語で、様々な環境上で動きます（Linux系、Windows、MacOS）。統計解析や作図に強い点が特徴です。

Rのデフォルトの機能だけでも相当のことができますが、Rはパッケージ（package）を活用することで、どんどんその機能を拡張することが出来ます。

ピカチュウはでんきタイプのポケモンなので、レベルアップに伴い、基本的にはでんき系の技を中心に覚えていきます。「でんじは」とか「かみなり」とか。 しかし、わざマシンを使うと、ピカチュウがその技を使用できるようになります。例えば「ねむる」のわざマシンを使うと、自力では覚えないはずの「ねむる」を使用可能になります。

パッケージはこのように、本来備わっていない機能を、後から追加するための道具と考えてください。詳しい説明は後述します。

Rのインストール

Rにまつわるエトセトラは、CRAN（しーらん：Comprehensive R Archive Network）からダウンロードします。日本では、以下のいずれかのミラーサイトからダウンロードするとよいでしょう。

• 統計数理研究所のCRANミラーサイト
• 山形大学のCRANミラーサイト  

自分の環境（Linux、MaxOS、Windows）に適したリンクを選択し...

「install R for the first time」をクリックしましょう（下図の反転している部分）。最新バージョンのRをダウンロードするための画面に移動します。「Download～」をクリックするとインストーラーがダウンロードされます。

インストーラーを起動したら、「OK」と「次へ」ボタンだけを押しましょう。他の設定は一切変更しないでください。

※もし変更する場合は、インストール場所に注意してください。過去に、「パスに日本語が含まれる場所」にインストールしたり（例：C:/Users/hogehoge/Documents/分析）、同期しているクラウドストレージ（onedrive等）上にインストールしたりしたりしたため、トラブった例がありました。

これでRがインストールできました。このままRを立ち上げて操作してもいいのですが、RStudio経由でRを操作したほうが圧倒的に便利です。そこで次節ではRStudioのインストール方法を説明します。  

RStudioとは

RStudioは、Rをより使いやすくするための統合開発環境（IDE）です。よく、「RStudioがあれば、Rは不要ですか？」と聞かれることがありますが、そうではありません。Rを直接操作するのではなく、RStudio経由でRを操作するわけなので、RもRStudioもそれぞれインストールする必要があります。

RStudioのインストール

ローカル環境にインストールして使用するRStudio Desktopと、サーバ上でRStudioを操作するRStudio Serverがありますが、ここではRStudio Desktopの利用方法を説明します。

RStudioの公式サイトから、最新版のRStudioインストーラーをダウンロードしてください。Installers for Supported Platformsの中から、各自の環境に合ったインストーラーを選んでください。

インストーラーを起動したら、「次へ」と「インストール」ボタンだけを押しましょう。他の設定は一切変更しないでください。

※もし変更する場合は、インストール場所に注意してください。過去に、「パスに日本語が含まれる場所」にインストールしたり（例：C:/Users/hogehoge/Documents/分析）、同期しているクラウドストレージ（onedrive等）上にインストールしたりしたりしたため、トラブった例がありました。

これでRStudioのインストールも完了です。RStudio経由でRを操作できるようになったので、以後はRStudioだけを起動させればよいです。Windowsの場合はスタートメニューから探してください。

RStudioの初期設定

次は、RStudioを本格的に使い始める前の初期設定を整えます。以後の初期設定はすべて、Tools　→　Global Options...　の中で行います。

デフォルトの作業ディレクトリ（Working Directory）の設定

通常、ファイルを参照する際に、「Cドライブの中の、R_practiceというフォルダの中の、practice_1.csvという名前のファイル（C:/R_practice/practice_1.csv）」というようにファイルが置かれた場所を絶対パスで正確に指定する必要があります。いわば、郵便を送るときに、「東京都港区～～丁目の・・・さん」のように住所を書くようなものです。これは結構めんどくさいです。

しかし、作業ディレクトリを設定しておけば、その作業ディレクトリを起点とした相対パスを書けば済みます。例えば作業ディレクトリが「Cドライブの中の、R_practiceというフォルダ」であれば、いきなり「practice_1.csv」というファイル名を指定するだけで、そのファイルを特定することができます。

RStudioには「プロジェクト」という機能があり、目的ごとに作業ディレクトリを設定することができます。しかしプロジェクトが指定されていない時は、デフォルトの作業ディレクトリが適用されます。プロジェクトについては今回は説明を省略しますが、本格的にRを使うようになると必須の機能なので、その時に改めて調べてみてください。

Browse...ボタンを押して、任意のフォルダを作業ディレクトリに設定しましょう。ここではCドライブ内に「R_practice」というフォルダを新しく作り、そこをデフォルトの作業ディレクトリに指定しました。

作業ディレクトリは任意の場所で構いませんが、なるべくパスに日本語が入らないようにしたほうが安心です。次に説明するように、文字化けしたら困るので。

   

文字コードの設定

プログラム内に、日本語でコメントを書き入れることがあると思います。しかし日本語のようにマルチバイト文字を使用する場合、文字コードによっては文字化けが生じる恐れがあります。そこでまず文字コードを設定してしまいましょう。

左側のカラムから「Code」、上部のタブから「Saving」を選び、以下の画面を開きます。Change...ボタンを押すことで、任意の文字コードを指定できます。ここではUTF-8を選びましょう。

”見栄え”の設定

次は好みの問題ですが、好きな見栄えを選びましょう。例えばフォントサイズを変更することが出来ます。

また、RStudio上でRコードを書いた場合は、予約語やコメントに色が付いてハイライトされます。その色のパターンを、Editor themeの中から選ぶことが出来ます。

パッケージをダウンロードするCRANリポジトリの設定

先述のように、Rは次々にパッケージをインストールしていくことによって、その機能を拡張していくことが出来ます。パッケージは、世界中に存在するCRANリポジトリからインターネット経由でダウンロード＆インストールすることが出来ます。その都度リポジトリを選ぶことも出来ますが、ここで設定したほうが楽です。

Change...ボタンを押すとリポジトリを変更可能です。日本国内では、以下の機関がCRANリポジトリを管理しているので、いずれかを選ぶと良いでしょう。

• 統計数理研究所（Japan (Tokyo) [https] - The Institute of Statistical Mathematics, Tokyo）
• 山形大学（Japan (Yonezawa) [https] - Yamagata University）

他にも細かい設定を行うことが可能ですが、最低限、以上の設定が終わればよいでしょう。

RStudioの機能

一度RStudioを閉じて、もう一度立ち上げなおしてください。下図の赤枠で囲ったところに、先ほど設定した作業ディレクトリが表示されていることがわかります。

さて、RStudioは上図のように、4つの領域に役割分担されています。それぞれの領域を、ペイン（Pane）と呼びます。各ペインの配置場所は変更可能ですが、ここではデフォルトの配置を前提として説明していきます。  

Sourceペイン（左上）

いわゆる、スクリプトを書くところです。実はRStudioは、RもPythonもStanも...様々なプログラミング言語のコードエディタとして機能します。どのような言語であろうと、このSourceペインでスクリプトを書きます。演劇でいうところの、台本だと思ってください。

Sourceペインでは、スクリプトをただ書くだけではなく、任意のコードを即座に実行することが出来ます。実行結果は、次に説明するConsoleペインに表示されます。

コードを実行するためには、以下の操作を実行してください。MacではCtrlではなくCmdキーを使用します。

• 任意の行のコードを実行したい
  • その行にカーソルを移動させて、Ctrl + Enter
• 複数行を同時に実行したい
  • 複数行を選択したうえでCtrl + Enter
• あるファイル内の全ての行を実行したい
  • Ctrl + Shift + Enter

Consoleペイン（左下）

上図を見るとわかるように、Sourceペインと同じコードと、その実行結果が表示されていますね。演劇でいうところの、舞台だと思ってください。

実はConsole上に直接コードを入力し、実行することも出来ます。いわば、台本にない台詞をアドリブで舞台上で披露するようなものです。しかし、分析結果の再現・再生可能性を高めるために、後で参照する可能性があるコードは全て、Sourceペインに記入したほうがよいでしょう。

Environmentペイン、Historyペイン（右上）

これらはタブでペインを切り替えます。Environmentペインには、作成したオブジェクトの一覧が表示されます。

例えば上図では、

• a <- 5
• b <- 3

により、aというオブジェクトに5を、bというオブジェクトに3を代入しています（後述しますが、Rでは代入演算子として、矢印を示す<-を使用します）。どのオブジェクトが作成されたか、それにはどのようなデータが入っているかを、Environmentペインで確認できるというわけです。

Historyペインでは、過去に実行したコマンドの履歴が確認できます。

Filesペイン、Plotsペイン、Packagesペイン、Helpペイン、Viewerペイン（右下）

やはりこれらもタブで切り替えます。

• Filesペイン
  • 作業ディレクトリ内のファイル一覧が表示されます。ここから任意のファイルを開くことも可能です。
• Plotsペイン
  • 作図結果がここに表示されます。
• Packagesペイン
  • 新しいパッケージを追加したい場合に、ここから選びます。また、すでにインストールされているパッケージを確認することもできます。
• Helpペイン
  • 関数やパッケージ等についてのサポート情報を確認することが出来ます。
• Viewerペイン
  • ちょっと特殊な出力（例えばJavaScriptベースのグラフ）が表示されます。  

Rの基本的な記法・使用方法

ではいよいよ、コードを書いていきましょう。左上のアイコンを押して、新規ファイルを立ち上げます。ここではRのコードを書いていくので、一番上の「R Script」を選びましょう。このファイルの拡張子は「.R」になります。

四則演算

以下のの演算子を用いることで、四則演算が可能です。割り算の余りを返したい場合の演算子は%%です。

※注意：以後、Consoleペイン上に表示される結果を表示しています。実際のコードでは冒頭に「>」を入力する必要はありません。
※Rでは、「#」以後に書かれたものはコメントアウトされます。つまり、「#」以後に書いたものは、プログラムのコードとは認識されないため、様々なメモを残すことが可能ということです。

> 5 + 3   #足し算
[1] 8

> 5 - 3   #引き算
[1] 2

> 5 * 3   #掛け算
[1] 15

> 5 / 3   #割り算
[1] 1.666667

> 5 %% 3  #割り算の余り
[1] 2

代入演算子

先述のように、Rにおける代入演算子は、矢印を表す<-です（半角の「<」の後に、半角の「-」を続けて打つ）。この代入演算子の面白いところは、->のように左右反転させても機能することです。この場合、->の左側の情報が右側のオブジェクトの中に代入されます（もっとも、この用法はあまりおすすめはしません）。
代入演算子<-は、Alt + =でショートカット入力させることもできます（+は入力する必要ありません）。

多くのプログラミング言語では、代入演算子は=ですが、Rでも使用可能です。現在のRではどちらの演算子を用いてもほぼ問題ありません。ごく稀に、<-でしか機能しない状況が存在するため、そういう意味では常に<-を使用したほうが汎用的ではあります。

> a <- 5    #aに5を代入
> a  #aの中身を表示
[1] 5
 
> b = 3     #bに3を代入
> b  #bの中身を表示
[1] 3
 
> 7 -> c    #cに7を代入
> c  #cの中身を表示
[1] 7

なお筆者はここ数年、代入演算子は=を用いています。Rを教える際に、他の言語（C、Python...）をすでに学んでいる／これから学ぶ学生が混乱しないように配慮していたのですが、いつの間にか筆者自身が=に慣れてしまいました。しかし特に困ったことはありません。

その他の演算子

Rには非常に多くの演算子があり、すべてをいきなり理解するのは難しいので、徐々に覚えていけばいいと思います。より高度な演算子は、以下のサイトを参考にしてください。

Rの演算子特集  

データの型・構造

とてもとてもとてもとても大事なところです。ただし、ここをちゃんと説明しようと思ったらかなり大変なので、参考となるサイトを紹介するにとどめます。ぜひ以下のリンク先を熟読してください。英語のサイトですが、非常に端的にまとまっています。ここを読むのがおすすめです。

• Understanding Basic Data Types and Data Structures in R

日本語で書かれたサイトでは、以下の記事が参考になります。

• Ｒのデータ型、データ構造

ベクトルの作成

上記のサイトにも書かれている内容のうち、重要な内容を一部抜粋して説明します。まずはベクトルの作成についてです。

ここではベクトルとは、「複数の情報が並んだもの」くらいの理解で大丈夫です。例えばトランプから無作為に5枚のカードを抜き取ったところ、数字は順に「3, 1, 2, 5, 8」だったとしましょう。これもベクトルです。

Rでは関数c()によってベクトルを作成します。要素をカンマで区切って並べます。文字列を要素とする場合は、引用符（""）で囲う必要がある点に注意です。

> v1 <- c(3, 1, 2, 5, 8)
> v1
[1] 3 1 2 5 8

> v2 <- c("abc", "あいう", "123", "春", "-----")
> v2
[1] "abc"  "あいう"  "123"  "春"  "-----" 

ベクトル内の要素にアクセス

あるベクトルの中の、n番目の要素にアクセスするときは、以下のように書きます。

> #1番目（最初）の要素を表示
> v1[1]
[1] 3

> # length()でベクトルの要素数を取得
> # 取得した要素数をインデックスにすることで、最後の要素を表示
> v1[length(v1)]
[1] 8

ベクトル名の直後に[]を書いて、その中に数字を入力することで、対応する要素にアクセスすることが出来ます。多くのプログラミング言語では、インデックスは0から始まると思いますが（最初の要素はv1[0]、2番目の要素はv1[1]）、Rでは1から始まることに注意が必要です。

matrixの作成

ベクトルが縦、あるいは横にずらっと並んだものが行列（matrix）のイメージです。例えば上で作成した2つのベクトル、v1とv2を縦方向に並べてみると...

#rbind()はベクトルを「縦方向」に連結する関数
#各ベクトルは左から右に向かってデータが並ぶ、すなわち「行（row）」になるので、頭文字r
> m1 <- rbind(v1, v2)
> m1
   [,1]  [,2]     [,3]  [,4] [,5]   
v1 "3"   "1"      "2"   "5"  "8"    
v2 "abc" "あいう" "123" "春" "-----"

横方向に並べてみると...

#cbind()はベクトルを「横方向」に連結する関数
#各ベクトルは上から下に向かってデータが並ぶ、すなわち「列（column）」になるので、頭文字c
> m2 <- cbind(v1, v2)  
> m2
     v1  v2      
[1,] "3" "abc"   
[2,] "1" "あいう"
[3,] "2" "123"   
[4,] "5" "春"    
[5,] "8" "-----" 

いずれの場合も、個々のベクトルの行数、または列数が等しくないと連結できないことに注意してください。

> m3 <- cbind(c(1, 2, 3),
              c(4, 5, 6, 7, 8)
              )
 警告メッセージ: 
 cbind(c(1, 2, 3), c(4, 5, 6, 7, 8)) で: 
  number of rows of result is not a multiple of vector length (arg 1)

matrixの要素にアクセスする

matrixには行と列が存在するため、ある要素にアクセスしたい場合には、行番号と列番号をそれぞれ指定する必要があります。例えば上で作成した、m2というmatrixから「春」というデータを参照したい場合には、以下のように書きます。[行番号, 列番号]という順番であることに注意してください。

> m2[4, 2]
  v2 
"春"

また、特定の行番号だけ、あるいは特定の列番号だけを指定すれば、該当するベクトルを参照することも出来ます。

> # 2行目だけ取得。カンマの右側が空白であることに注意
> m2[2,]
      v1       v2 
     "1" "あいう" 

> # 1列目だけ取得。カンマの左側が空白であることに注意
> m2[,1]
[1] "3" "1" "2" "5" "8"

さらに、c()によって複数の行番号や列番号を指定することもできます。

> # 3行目と5行目だけを取得 
> m2[c(3, 5),]
     v1  v2     
[1,] "2" "123"  
[2,] "8" "-----"

data.frameの作成

matrixは、複数のベクトルが縦方向または横方向にずらっと並び、長方形の形になったものですが、Rでは同じ形をしているデータフレーム（data.frame）というクラスもあります。Rでは、データをdata.frameクラスにしたうえで扱うことが多いです。

data.frame()関数に複数のベクトルを与えると、data.frameを作成することが出来ます。

> df1 <- data.frame(v1,
                    v2 = c("a", "b", "c", "d", "e"),
                    v3 = 10:14) #c(10, 11, 12, 13, 14)に等しい

> df1
  v1 v2 v3
1  3  a 10
2  1  b 11
3  2  c 12
4  5  d 13
5  8  e 14

ここで、ベクトルv1はすでに作成していたものを利用しています。v2とv3は新たに作成したベクトルです。作成されたdata.frameでは、個々のベクトルに名前が与えられ（列名）、上から下にデータが並んでいます。

あたかも、cbind()でベクトルを結合した場合と同じように見えますが、クラスは異なります。

> class(m2)
[1] "matrix"

> class(df1)
[1] "data.frame"

data.frameでもmatrixと同様に、行番号や列番号を指定することで、要素にアクセスすることが出来ます。さらにdata.frameでは列名が存在するので（もしdata.frame作成時に列名を明示しなかった場合には、適当な列名がつけられます）、列名を明示して特定の要素にアクセスすることもできます。

例えば上で作成したdf1というdata.frameで、v1列をごっそり抽出したかったら、data.frame名$列名と書きます。

> df1$v1
[1] 3 1 2 5 8

これはベクトルとして抽出されているので、さらにベクトル内の位置を指定することで、該当する要素にアクセスできます。

> df1$v1[3]
[1] 2

パッケージのインストールと読み込み

インストール

PackagesペインのInstallボタンを押すと...

下図（左）のウィンドウが現れます。インストールしたいパッケージ名を入力すると、候補が現れるので、適切なパッケージを選びましょう。ここでは、データの要約や因子分析や媒介分析や...とにかく様々なことが可能になる便利なパッケージ、「psychパッケージ」をインストールしてみます。Installボタンを押すとインストールが開始されます。

あるいは以下のようにinstall.packages()関数を用いて、コードによりインストールするパッケージを指定することも出来ます。

install.packages("psych")

一度パッケージをインストールすれば、Rのバージョンが変わらない限り（つまり別のバージョンのRを再インストールしない限り）、そのPCでは同じパッケージをインストールする必要はありません。

ただし、パッケージによっては開発が続けられており、バージョンアップが行われることがあります。もし最新版のパッケージがCRAN上で配布されていれば、再び上記の手続きをとることにより、最新版のパッケージをインストールすることが出来ます。

「もし最新版のパッケージがCRAN上で配布されていれば」というところがポイントで、パッケージによっては、最新版が現時点ではCRANにまだ登録されていなかったり（いずれは登録されるようになるでしょう）、そもそもパッケージ自体がCRANに登録されていなかったりすることがあります。よくあるのは、「GitHub上ではパッケージが公開されているけれど、CRANには（まだ）登録されていない」というケースです。

その場合、

1. devtoolsというパッケージをインストールして、GitHubからパッケージをインストールする機能をRに拡張させる（すでにdevtoolsパッケージをインストール済みならこの手順は省略可）
2. devtoolsパッケージの関数を用いて、GitHubから、本来インストールしたいパッケージをインストールする

という手順になります。

例えばggplot2というパッケージは、CRANに登録されていますが、開発中のバージョンはGitHubから直接インストールすることが出来ます（install_github()の中は、当然パッケージごとに異なるパスを書く必要があります）。

# install.packages("devtools")
devtools::install_github("tidyverse/ggplot2")

ほとんどの場合は上記のようなコードがパッケージの公式ページに記載されているはずなので、コピペすればOKです。

パッケージの読み込み

さて、パッケージはポケモンで言うところの「わざマシン」のように、後付けでRの機能を拡張するための道具だと説明しました。ただし「わざマシン」と違うところがあります。それは、パッケージは使うたびに読み込みが必要というところです。

「わざマシン」であれば、ひとたびピカチュウに「ねむる」のわざを覚えさせると、ポケモンセンターに寄ろうが、パーティーメンバーから外そうが、バトルで倒れようが、覚えたわざを使い続けることが出来ます。

しかしパッケージは、Rのセッションを立ち上げるたびに（Rを起動するたびに）、パッケージを読み込みなおす必要があります。インストールではありません。パッケージをそのセッション内で使用する場合は、読み込みが必要ということです。

例えば、psychパッケージには、ある量的変数に関する様々な要約統計量を計算してくれる、describe()という関数があります。この関数を、mtcars$mpgという変数に適用してみましょう。

> describe(mtcars$mpg)

 describe(mtcars$mpg) でエラー: 
   関数 "describe" を見つけることができませんでした 

できひんのかーい。

それもそのはず。このRセッションでは、まだpsychパッケージを読み込んでいなかったからです。R側からすると、「describe()...？いえ、知らない関数ですね」ということになります。

ではパッケージを読み込みましょう。それには、library()という関数を用います。

> library(psych)
> describe(mtcars$mpg)

   vars  n  mean   sd median trimmed  mad  min  max range skew kurtosis   se
X1    1 32 20.09 6.03   19.2    19.7 5.41 10.4 33.9  23.5 0.61    -0.37 1.07

パッケージを読み込みさえすれば、ちゃんとそのパッケージ内の関数を使えるようになりますね。

さて、ここで注意があります（次に述べることは筆者個人の意見です）。それは、Rにデフォルトで備わっている関数ではなく、パッケージ内の関数を用いる場合は、その関数がどのパッケージの関数かを常に明示すべきということです。

なぜかというと、プログラムのコードは、自分一人が見るものではなく、多くの場合は「今の自分以外」と共有するからです。「今の自分以外」には、他者も、「数年後の自分」も含まれます。

上記の例では、「describe()という関数が、psychパッケージの関数であることを知らない人」がいたら、混乱することでしょう。そこで、筆者はどのような関数であれ、パッケージ内の関数を使用する場合は、以下のように書くことを自分ルールとしています。

> library(psych)
> psych::describe(mtcars$mpg)

   vars  n  mean   sd median trimmed  mad  min  max range skew kurtosis   se
X1    1 32 20.09 6.03   19.2    19.7 5.41 10.4 33.9  23.5 0.61    -0.37 1.07

違いがわかるでしょうか。psych::describe(mtcars$mpg)と、関数名の前にパッケージ名を明記した点だけが違います。Rではこのように、パッケージ名を明記する場合は、「パッケージ名::パッケージ内の関数」と書きます（半角のコロン2つ）。

これには別の利点もあります。Rでは、複数のパッケージが、同じ名前の関数を持っていることがあります。例えば、psychパッケージは、psych::mediate()という関数を持っていますが、別のmediationパッケージも、mediation::mediate()という関数を持っています。
関数名しか書かなかったとしたら、どちらのパッケージのmediate()関数なのか、区別がつきません。

外部ファイルの読み込み

csvファイルの読み込み

Rで解析をする際に、実験データなどがまとめられたファイルをRに読み込むことが多いでしょう。外部ファイルを読み込むための関数は、そのファイルがどのような形式でデータをまとめているかによって変わります。

典型的には、csvファイルとしてまとめられていることが多いと思います。csvとはComma Separated Valueの略で、下図のようにカンマでデータが区切られた、テキスト形式のファイルです。

Excelでデータをまとめている場合にも、「名前を付けて保存」する際に、「ファイルの種類」をcsvに変更すれば、csvファイルとして保存されます（拡張子が.csvになっていることを確認してください）。

csvファイルを読み込む方法は、主に2つです。

一つ目は、EnvironmentペインのImport Datasetボタンから「From Text」を選択し、GUIにより読み込む方法です。「From Text」には以下の2つがあります。

• From Text (base)
  • デフォルトで使用可能なread.csv()関数によってデータを読み込む
• From Text (readr)
  • readrパッケージのread_csv()関数によってデータを読み込む

どちらを使うのが良いかは話すと長くなるので、ここでは割愛します。もし数値だけが格納されたデータセットを読み込むのならば、どちらでも構わないと思います。

二つ目の方法は、直接read.csv()やreadr::read_csv()という関数をエディタやコンソール上で実行することです。実は上記のGUIを用いた場合にも、結局これらの関数が実行されることになります。

tmp <- read.csv("C:/R_practice/sample_data.csv")

注意してほしいことは、原則的に、読み込むデータセットは、長方形でなければならないということです。つまり、列によって行数が違うような凸凹したデータセットは読み込めないということです。

原則的に、と書いたのは、凸凹したデータセットでも読み込む方法があるからです。

例えば、データがこんな風にまとめられていたとしましょう。

このデータセットは、以下の悩ましい問題を含んでいます。

1. 1行目に、1セルだけ、分析とは関係ない情報（記録日時）が掲載されている
2. 2行目が空白
3. Scoreの列に空白のセルがある

でも大丈夫。外部ファイルを読み込むための関数には、データを柔軟に読み込むための様々な引数が存在するのです。

tmp <- read.csv("C:/R_practice/sample_data.csv", skip = 2, na.strings = "", header = TRUE)

> tmp
  ID Age Score
1  1  25    80
2  2  30    NA
3  3  22    68

うまく読み込むことができました。それぞれの引数を説明します。

• skip = 2
  • 上から何行を無視するかを決めることができます。今回は、2行目までを無視しました。
• na.strings = ""
  • データが存在しないセルを、RではNAで表します。どのようなデータをNAとみなすかを指定できます。
  • 今回は空白のセルをNAとしたかったので、引用符の中に何も記述しませんでした。
  • もしデータが欠測していたセルを.で表していたなら、na.strings = "."とすればよいです。
• header = TRUE
  • 読み込む対象の範囲のうち、一番上の行をヘッダ（列名）とみなすか、データの一部とみなすか、を指定できます。
  • もし一番上の行がいきなりデータから始まっていたら、header = FALSEとしましょう。

他にも色々と引数があります。

Excelファイルの読み込み

readxlパッケージのread_excel()という関数によって、Excelファイル（拡張子が.xlsxなど）を読み込むことも出来るようになりました。直接この関数をエディタやコンソール上で打ち込んでもいいし、やはりGUIも用意されています。

readxl::read_excel()の引数は、read.csv()やreadr::read_csv()とはちょっと異なっています。詳しくはコンソール上で?readxl::read_excelと入力して、ヘルプを確認してください。

あ、そうそう。このように、?の後ろに関数名を書くと、RStudioの右下にあるHelpペインに解説が表示されます。積極的に活用してください。

おまけ：データ構造の確認

自分でdata.frameなどを作成する場合には、どのようなデータを扱っているか（特に、数値か文字列か）について把握しやすいですが、外部からデータを読み込んだ場合は、本当に自分が想定するデータの構造になっているか、確認することが必須です。というか常に確認してください。

例えば、lme4というパッケージを読み込むと参照できるようになる、sleepstudyというデータセットがあります。これは、睡眠時間が削られることで、刺激に対する反応時間がどう変化していくかを記録したデータセットです。以下の通り3列が存在し、いずれも数値が格納されているように見えます。

> #install.packages("lme4")
> library(lme4)

> data(lme4::sleepstudy)
> head(lme4::sleepstudy) #data.frameの冒頭6行だけ表示する

  Reaction Days Subject
1 249.5600    0     308
2 258.7047    1     308
3 250.8006    2     308
4 321.4398    3     308
5 356.8519    4     308
6 414.6901    5     308

ここでstr()という関数を用いて、データ構造（structure）を確認してみましょう。

> str(lme4::sleepstudy)

'data.frame':  180 obs. of  3 variables:
 $ Reaction: num  250 259 251 321 357 ...
 $ Days    : num  0 1 2 3 4 5 6 7 8 9 ...
 $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...

反応時間を表すReactionや、実験経過日数を示すDaysには、numと書かれています。本記事ではデータの型や構造について説明を省略したので、何のことかわからないかもしれませんが、要は「数値」と見なされているということです。numはnumericの略です。

一方、実験参加者のID番号を示す、Subjectには、Factorと書かれています。よく見ると、"308", "309", "310",...と、各数字が引用符で囲われていますね。すなわち、見た目は数値であっても、これらは文字列のように見なされているということです。

このようなデータ構造をしっかりと確認しないと、分析時などに思わぬミスにつながります。特に外部ファイルを読み込んだ場合や、他者が作ったデータセットを使用する場合は、必ずstr()でデータ構造を確認するようにしてください。

おわりに

他にも書きたいことは山ほどありますが、キリがないのでこのあたりにします。あくまで今回の記事の目的は、「まったくRを触ったことが無い人に、事前に読んできてもらう資料」を作成することでした。これらを踏まえたうえで、他の基本的な機能や、より発展的な機能を勉強していただければ幸いです。

Enjoy!

「RとStanではじめる ベイズ統計モデリングによるデータ分析入門」書評

この記事について

著者の馬場真哉様より、2019年7月10日に講談社より発売の、「RとStanではじめる　ベイズ統計モデリングによるデータ分析入門」をご恵投いただきました。ありがとうございます！！

www.kspub.co.jp

事前に献本をいただけるということを伺っていたので、その時から「ご恵投いただきました！」とTwitterで報告するだけでなく、簡単にでも読んでみた感想を書こうと決めていました。

まだざっと読んだ段階で、コードを実際に走らせてもいないのですが（もちろん後でじっくり読みながら実行します）、感想や関連書籍との比較をしていきたいと思います。

本記事の方針

本書の「はじめに」の部分やサポートページには、以下のような方を対象読者としていると書かれています。本記事も、そのことを念頭に書いていこうと思います。

統計学の基礎やベイズの定理などの基本事項を学んでみたものの、その有効性がピンとこない、という方を対象としています。
例えば、ベイズの定理の導出では終わらない、より実用的なベイズ統計学について学んでみたいと思った方、あるいは統計モデリングに興味が出てきた方は、本書の対象読者であるといえます。理系文系は問わず、データを分析する必要性がある人はもちろん、データサイエンスが専門でないエンジニアの方でも読めるような内容を目指して執筆しました。

逆に、R やStanを使ったベイズ統計モデリングにすでに明るいという方は、読者として想定していません。ベイズ統計学の数理的な側面もある程度省略しました。基本的な事項を、実装を何度も繰り返すことを通じて学ぶ構成になっています。
ベイズ統計モデリングの上級者をさらなる高みへと持っていく本ではないことに、留意してください。

また、RとStanでベイズ統計モデリングを行う入門書というと、真っ先に頭に浮かぶのが、松浦健太郎著「StanとRでベイズ統計モデリング（通称、アヒル本）」だと思います。僕自身も、恐らく人生で最も読み返した回数が多い本ではないかと思うくらい、お世話になっています。ページをめくったり書き込みしたりでボロボロになってきたので、最近2冊目を買いました。

www.kyoritsu-pub.co.jp

それぞれの本の強みを紹介するため、要所要所で、アヒル本と本書を比較していこうと思います。そのため本記事の想定読者も、アヒル本を一度読んだことがある人、としています。

なお結論を先取りすると、それぞれの本で強調されている点が異なるので、両方読むといいと思います！ということになります。

第1部：【理論編】ベイズ統計モデリングの基本

本書第1部では、ベイズ統計モデリングとは何かや、それを行う上で必要となる数学的基礎知識について紹介しています。特に、ベイズ統計モデリングが基礎としている「確率」について、紙面が割かれています。

これらの基礎的なパートについては、本書とアヒル本で「力を入れている」箇所が異なるように思いました。

アヒル本

アヒル本では、確率、確率密度関数、確率質量関数、尤度などについては、ある程度背景知識を持っていることを前提として説明が進むように思います。もちろん一通り説明はあるのですが、（初心者にとっては）難しい言葉もしばしば登場します。

その代わり、モデル、確率モデル、統計モデルとは何か。そもそもなぜこれらを考える意義があるのか、といったところは紙面が割かれていて、統計モデリングに対する著者の熱い思いを感じます。「いいからMCMCだ！」ではなく、その前にすべきこと（例えばメカニズムの想像）について、慎重に指南してくれる本です。

また、アヒル本では1章まるまる使って、様々な確率分布を丁寧に紹介しているパートがあります。正規分布、ポアソン分布、二項分布などの「メジャーな」確率分布はもちろんのこと、カテゴリカル分布やディリクレ分布、多変量正規分布など、発展的なモデルを考える際にお世話になる分布についても、網羅しています。

アヒル本では、基本的な部分を順を追って丁寧に説明することを重視しているものの、基礎を身に着けた読者が自身の発想でモデリングできるよう、その先の道しるべも用意してくれているような印象を持っています。

本書

本書では、アヒル本よりもより初心者を想定しているように思います。そのため確率、確率密度関数、確率質量関数、尤度などについて、この本を読みながら知識を得ていくことが可能です。馬場さんが過去に書かれた「Pythonで学ぶあたらしい統計学の教科書」でもこのあたりは言及されています。

www.shoeisha.co.jp

紹介されている確率分布も、正規分布、ベルヌーイ分布、二項分布、ポアソン分布、一様分布といった「メジャーな」ものに限定されており、より初心者にフォーカスした説明を意識されているのかなと思いました。

第2部：RとStanによるデータ分析

第2部は本書の特徴的なパートだと思います。「RとStan」でベイズ統計モデリングを行う際に必要となるRの基礎や、使用するパッケージの記法について、紹介しています。

アヒル本

アヒル本では、Rの利用方法については特に説明はなく、外部パッケージも極力使わない方針で解説されています。可視化において大活躍するggplot2パッケージは、rstanパッケージ自体が依存しているのでアヒル本でも活用されていますが、その記法は説明されていません。その代わりStanやrstanの説明に注力している感じです。

説明の順序についても、書籍間で違いがあります。アヒル本では、ページをめくるたびに、一つずつ技術を身に着けていくような説明の仕方になっている気がします。

例えばモデルの評価方法の一つである事後予測チェックは、アヒル本では簡単なモデルのパラメータ推定ができるようになった後で言及されます。

まずmodelブロックまで一通り書けるようになることを目指し、その次にgenerated quantitiesブロックの説明をする中で、事後予測チェックの方法を紹介しています（generated quantitiesブロック内で乱数を発生させる、という技術）。

本書 

RからStanにデータを渡すためには、それなりにRのテクニックが必要になります。例えばlistであったり、データの型であったり。本書では最初にRの主要な機能を解説しているため、初めてRを触る人であっても問題ありません。

さらに本書は積極的に拡張パッケージを活用する方針のため、簡単ではありますがggplot2の記法を解説しています。またそれだけでなく、MCMCの結果を可視化する際にbayesplotパッケージも活用し、その使用方法も説明しています。例えば診断の類はbayesplot等の拡張パッケージに任せてしまうため、まだgenerated quantitiesブロックの話が出てきていない段階で、事後予測チェックの話が出てきます（もちろんその後でgenerated quantitiesブロックの説明も出てきます）。

※2019/7/9追記：すいません見落としていましたが、generated quantitiesブロックの中でモデルに従う乱数を生成する方法も載っていました（136ページ目）

ちなみに手前味噌になりますが、bayesplotパッケージについては、以前本ブログで記事を書いたのでご参照ください。

das-kino.hatenablog.com

個人的には、bayesplotについて本書で紹介しているのは、かなり重要なポイントだと思います。診断は非常に大事な手続きだと思うのですが、診断に必要な手続きが複雑であれば、「めんどくさい」と思ってスキップしてしまうかもしれません。

実際のところ、診断の観点は様々あるのでこれ自体容易ではないのですが、少なくとも「視覚的に」診断するうえで、bayesplotはユーザの負荷を大きく下げてくれます。

第3部：一般化線形モデル

アヒル本

アヒル本では、データやパラメータをint型やreal型で宣言し、modelブロック内でforループを用いて書くシンプルな記法から始め、けっこう長い間この記法を使い続けます。ベクトルや行列を用いた書き方が登場するのは、全12章中、9章目と遅めです。

恐らくこれは、繰り返すように、「少しずつ技術を身に着けていく」ことを狙っているからではないかと思います。

本書

一方本書では、第2部でStanの基本的な記法を簡潔に紹介したのち、第3部の実践編ではいきなりベクトル・行列を活用した記法に切り替わります。説明が丁寧なので問題ないのですが、不慣れな人はここでちょっと段差が高くなった印象を受けるかもしれません。

brmsパッケージ

本書の最大の特徴と言っても過言ではないのが、brmsパッケージの活用です。brmsパッケージを用いると、ユーザ自身がStanコードを書く必要なく、シンプルな記法で線形モデルのパラメータ推定を実現することができます。

なお本ブログでも以前に記事を書いたので、ご参照ください。

das-kino.hatenablog.com

このbrmsを活用したパートは、本当に現時点では稀有だと思います。僕自身が↑のブログ記事を書いたモチベーションの一つも、日本語のチュートリアル記事がほとんどないということでした。少なくとも現時点では、brmsパッケージの使用方法を事例とともにここまで丁寧に説明した和書は、他にないのではないかと思います（僕が知らないだけかもしれませんが...）。

本書第3部では、分散分析モデルやGLMの推定をbrmsパッケージに思い切って任せてしまうことで、Stanコード自体の説明はアヒル本より少なくなるけれど、そのぶん結果の可視化も含めて様々なケースを紹介しています。

特に、交互作用について独立した節を作って事例を紹介しているのは、まさに痒い所に手が届くという感じです。

第4部：一般化線形混合モデル

アヒル本

アヒル本では、階層モデルについて、「階層という発想を導入する必要性」やStanコードの書き方も含め、詳しく説明されています。Stanでは、階層モデルはよりシンプルなモデルの自然な拡張として書けるうえ、実用的機会が多いため習得することが重要だと思います。

僕自身、アヒル本を読んで初めて、「そういうことだったのか！」と目の前の霧が晴れたように理解できた部分が多かったです。

本書

一方本書では、意外なことに紙面が少ないです。それは、brmsを用いることで、階層モデル（ここでは一般化線形混合モデル）も極めて簡単に推定できてしまえるからです（lmerやlmerTestとほぼ同じ記法で書けます）。もちろんランダム切片モデルとかランダム係数モデルについての説明もあるのですが、アヒル本と対比すると説明はシンプルです。

第5部：状態空間モデル

出ました！という感じです。なにせ馬場さんといえば、時系列分析に特化した書籍、「時系列分析と状態空間モデルの基礎：RとStanで学ぶ理論と実装」の著者なのですから。

www.amazon.co.jp

アヒル本でも状態空間モデルについて非常に分かりやすく解説されているのですが、本書では第5部がまるまる時系列データの分析に割かれています。

その代わりアヒル本では、混合分布モデルや、空間構造データを用いたモデリングなど、異なる観点で詳しくモデリング事例を紹介しています。

まとめ

冒頭にも書いた通り、アヒル本と本書は、狙いが恐らく異なっており、それぞれ強調されているポイントが違うので、両方読めば様々なケースに対応できるようになるのではないかと思いました。

あくまで僕の勝手な推測ですが、

アヒル本は、ユーザが自由なモデリングをすることを支援しているように思いました。そのために、Stanの記法について丁寧に解説し、収束しなかった場合の対処やトラブルシューティング、打ち切りや外れ値への対処、階層モデルという発想、離散パラメータを扱うテクニックなど、上級者へ至るためにいつかは通らなければならないポイントを、あらかじめ道しるべとして残していてくれているように思います。

世の中にはこういうモデルがあるんだよ、ユーザは必要なモデルを選べばいいよ、ではなく、君は君らしくモデリングする自由があるんだ（以下略）という印象です。

一方本書でも、もちろんユーザのモデリングを支援しているのですが、恐らく本書の想定読者は、アヒル本よりも初心者寄りであることから、様々なモデリングの事例をメニューとして用意してあげて、それらを導入するハードルを下げてやろう、そのためにbrmsやbayesplotなどの便利なパッケージもどんどん紹介していこう、という狙いなのかなあと思いました。

どちらか1冊読めば全方位に対応できる、ということではなく、それぞれの本が異なる強みを持っているので、これから初めて勉強する人も、すでにある程度触れている人も、両者を行き来することでさらにもう一段高みへ上れるのではないかと思います。

駆け足で読んだ段階なので、まだ見落としている部分があるかもしれません。もしかしたら今後追記するかもしれませんが、ひとまず感想まで。

ggplot2でヒストグラムを描くときに外枠だけに線を引きたい

タイトルの通り、「ggplot2でヒストグラムを描くときに外枠だけ線を引きたい」という質問を受けました。図にすると以下のような状況です。

枠線を付けたいときにパッと思いつくのは、引数colourを指定することですが、そうすると下図・右のようになってしまいます。

library(ggplot2)

bw <- diff(range(faithful$waiting))/20

# 左図 -------------------------------------------------
ggplot(data = faithful, mapping = aes(x = waiting)) +
  geom_histogram(binwidth = bw, fill = "skyblue") +
  labs(title = "通常のヒストグラム")

# 右図 -------------------------------------------------
ggplot(data = faithful, mapping = aes(x = waiting)) +
  geom_histogram(binwidth = bw, fill = "skyblue", colour = "darkblue") +
  labs(title = "全ての棒に枠線を付けたヒストグラム")

で、どうやるんだろうと調べてみました。安直に「ggplot2　histogram　only　outline」というググり方をしましたが、幸いにもすぐにヒットしました。ggplot2の調べ物は英語でググるのが吉です。
参考：can I draw a histogram that has an outline around the entire shape, but not the individual bars?

ヒストグラムを描くときに、Y軸にマッピングする変数を省略するのではなく、mapping = aes(y = ..density..)としてkernel density estimateをマッピングすることを明示したうえで、geom_step()を用いてヒストグラムの輪郭だけを描くんですね。

あとはヒストグラムのレイヤを重ねてやれば、見た目上は外枠だけに線が描かれたように見えるはずです。
ただし注意が必要なのは、ヒストグラムはY軸が度数（count）である点です。そこで以下のページを参考に、ひと工夫します。ポイントは、geom_step(mapping = eval(bquote(aes(y=..count..)))というところです。
参考：データのhistogramと確率密度関数を同時に描きたい

勝った、第3部完！

dens <- density(faithful$waiting)
bw <- diff(range(faithful$waiting))/20

ggplot(data = faithful, mapping = aes(x = waiting)) +
  geom_histogram(binwidth = bw, fill = "skyblue") +
  geom_step(mapping = eval(bquote(aes(y=..count..))),
            stat = "bin",
            binwidth = bw,
            colour = "darkblue") +
  xlim(range(dens$x))

ずれとるやんけ！

考えてみたらそりゃそうで、geom_step()が描く線の始点は、各階級幅の中央になるので、階級幅の半分だけずれることになります。仕方がないので、引数positionをいじって、X軸のずれを補正します。

dens <- density(faithful$waiting)
bw <- diff(range(faithful$waiting))/20

ggplot(data = faithful, mapping = aes(x = waiting)) +
  geom_histogram(binwidth = bw, fill = "skyblue") +
  geom_step(mapping = eval(bquote(aes(y=..count..))),
            stat = "bin",
            binwidth = bw,
            position = position_nudge(x = -bw/2), #←階級幅の半分だけ左にずらしている
            colour = "darkblue") +
  xlim(range(dens$x))

完成。
ただWarning messages出てるんですよね...。うーん、もうちょっと改善の余地がありそうです。もしもっと適切な方法があれば教えていただけたら嬉しいです。

Warning messages:
1: Removed 2 rows containing missing values (geom_bar). 
2: Removed 1 rows containing missing values (geom_path). 

---

追記：
atusyさん（@Atushi776）から、「stat_bin()を使うとよい　＆　引数pad = TRUEを書くとよい」とご提案いただきました。ありがとうございます！　確かにこちらのほうがシンプルで良いですね。そしてWarning messagesも出ない！

bw <- diff(range(faithful$waiting))/20

ggplot(data = faithful, mapping = aes(x = waiting)) +
  geom_histogram(binwidth = bw, fill = "skyblue") +
  stat_bin(binwidth = bw,
           colour = "darkblue",
           geom = "step",
           position = position_nudge(x = -bw / 2),
           pad = TRUE)

Enjoy !

brmsパッケージを用いたベイズモデリング入門

はじめに

本記事は、Stan Advent Calendar 2018 16日目の記事です。brmsパッケージを用いたモデリングの方法や、便利機能、ちょっとしたtipsを紹介します。

brmsパッケージのコンセプトは、
ベイズで（Bayesian）、回帰モデルを（Regression Model）、扱おう、　Stanを使って（using Stan）
です。 読み方はそのまま「ビー・アール・エム・エス」です（僕は勝手にブラムスと呼んでいますが）。

その名の通り、Stanのラッパーパッケージですが、ユーザは自身でStanコードを書く必要はなく、モデルを指定すると自動的に内部でStanコードが生成され、サンプリングが実行されます。

最近は、brmsを用いて分析した論文も見かけるようになってきました。例えば難波修史さんのこの論文とか

• Dynamic Displays Enhance the Ability to Discriminate Genuine and Posed Facial Expressions of Emotion

僕もたまに使っているのですが、あまりにも多機能なので、しっかりと把握できていない部分も多いです。もし本記事に間違いがあればご指摘いただければ幸いです。

brmsパッケージの良いところ

その1

回帰モデルと一言でいっても、実はその守備範囲はかなり広いです。線形モデルも非線形モデルも、その気になればSEMもできるそうです（それならblavaanでいいと思いますが）。

1. 一般線形モデル
2. 一般化線形モデル
3. 階層線形モデル
4. その他の自由なモデリング

嬉しいのは、上記の1~3は、多くのRユーザが慣れ親しんでいるであろう記法と同じ記法でモデルを指定できることです。 例えば

• 回帰モデルではlm(y ~ x, data = dat)
• ポアソン回帰モデルではglm(y ~ x, family = "poisson", data = dat)
• 階層線形モデルではlme4パッケージやlmerTestパッケ―ジを用いて、
  • lme4::lmer(y ~ x + (1 | group), data = dat)
  • lme4::glmer(y ~ x + (1 | group), family = "poisson", data = dat)

のように書きますよね。brmsパッケージでも、記法はほぼ同じになります。

その2

デフォルトで用意されている確率分布が、ものっっっっっっっすごく多いです。コンソールで

?brms::brmsfamily

と入力してみてください。用意されている確率分布の一覧が確認できます。 驚くのは、「ベルヌーイ分布とポアソン分布の混合分布である、ゼロ過剰ポアソン分布」等、異なる確率分布を組み合わせた、いわゆる混合分布も多く用意されていることです。通常、このような混合分布を用いたモデリングをStanで書くにはひと手間要るのですが、brmsパッケージでは確率分布を指定するだけで済みます。

また、様々な確率分布からの乱数を生成するための関数や...

確率（密度）を返す関数も用意されています。

その3

「内部で何をやっているのか」を知ることができます。

ベイズを用いる際に最も気になることの一つが、事前分布の設定でしょう。brmsでは、brms::get_prior()で指定したモデルの事前分布を知ることができます。brmsがデフォルトで採用している事前分布も分かるので、報告時に困ることはありません。

また、brms::make_stancode()で指定したモデルのStanコードを出力させることもできます。「この確率分布を使ったモデルをStanで書きたいとき、どうすればいいんだろう」と悩んだとき、brmsで同じ確率分布を使った適当なモデルでStanコードを出力させて、「なるほど、こう書けばいいのか」とヒントを得ることも可能です。

僕が驚いたのが、brms流のガンマ回帰の書き方です。parametersブロックと、modelブロックの一部を以下に抜粋します（事前分布のコードは省略しています）。

parameters { 
  vector[Kc] b;  // population-level effects 
  real temp_Intercept;  // temporary intercept 
  real<lower=0> shape;  // shape parameter 
} 

model { 
  vector[N] mu = temp_Intercept + Xc * b;
  for (n in 1:N) { 
    mu[n] = shape * exp(-(mu[n])); 
  } 
  target += gamma_lpdf(Y | shape, mu); 
} 

modelブロックでは、まず切片と説明変数の線形結合をmuに代入しています。 ガンマ分布にはshapeとrate、2つのパラメータがあり、平均はshape/rateになります。brmsではガンマ回帰におけるデフォルトのリンク関数は対数リンク関数（つまり逆リンク関数は指数関数）なので、exp(mu) = shape/rateの関係があります。 これを変形すると、rate = shape * exp(-mu)になります。

ところがparametersブロックではrateに相当するパラメータはそもそも宣言しておらず、一度別の用途で使ったmuという変数を再利用して上書きしています。あたかも、「もともとのmuには興味ないでしょ」と言わんばかりに。

さらには、そのmuはmodelブロック内で宣言されているため、最終的にrate（に相当する）パラメータを取り出すことができません。なかなか思い切ったことをするなあと思ったのですが、こんな風にbrmsはサンプリング効率を高めるようなStanコードの書き方をしているため、推定がかなり速いです。

その4

他のパッケージとの連携が充実しています。

例えばbrmsはbayesplotパッケージを利用して、brms::pp_check()という関数により、容易に事後予測チェックを行うことができます（詳しくはこちら）。

その他にも、bridgesamplingパッケージを利用してbayes factorを容易に計算できたり、欠測値がある場合にmiceパッケージで多重代入法による補完を行ったあとで、補完したそれぞれのデータセットからサンプリングを行えたりします（詳しくはこちらから）。

---

brmsパッケージの使用例

それではさっそく、使用例を紹介していきましょう。

インストール & ロード

brmsはCRANに登録されているので、以下でOKです。なお本稿執筆時点における最新バージョンは2.6.0です。

install.packages("brms")
library(brms)

上述のbayesplotパッケージやbridgesamplingパッケージ、それからStanをRから扱うためのrstanパッケージを含む、種々の依存パッケージが一緒にインストールされます。便利なものも多いので、初めてStanを使う人は、brmsをCRANからインストールするのが一番楽な気がします （Rtoolsは別途インストールする必要があります）。

2018/12/19追記：
記事を書いた直後にバージョン2.7.0にアップデートされました。"extending support for Gaussian processes"ですって...

一般線形モデル（回帰分析）

まずはベーシックな回帰分析をやってみましょう。サンプルデータは、プリセットのmtcarsです。32台の車の性能が収められたデータセットになります。

head(mtcars)

                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

モデルの指定

燃費（mpg）を、車体重量（wt）とマニュアル車かオートマ車か（am）で予測する重回帰分析を試みます。mtcars$amには0か1が格納されていますが、今回はこれらをfactor型に変換したうえで、交互作用項を含むモデルを指定することにします。

mtcars$am <- as.factor(mtcars$am)

stats::lm()を用いると以下のように書きますが...

fit <- lm(mpg ~ wt * am, data = mtcars)
summary(fit)

brmsでは、どのようなモデルであっても、統一的にbrm()という関数のなかでモデルを指定します。もちろん、rstan::sampling()でサンプリングする際に指定するような諸々の引数も、brm()は対応しています。例えば、サンプリング回数（iter）、ウォームアップ期間（warmup）、乱数のシード（seed）、チェイン数（chain）など。

fit <- brm(mpg ~ wt * am,
           data = mtcars,
           iter = 2000,
           warmup = 1000,
           seed = 1234,
           chain = 4)

summary(fit)

Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: mpg ~ wt * am 
   Data: mtcars (Number of observations: 32) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    31.42      3.17    24.90    37.59       2139 1.00
wt           -3.79      0.82    -5.37    -2.11       2150 1.00
am1          14.96      4.39     6.66    23.68       1682 1.00
wt:am1       -5.34      1.48    -8.37    -2.50       1828 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     2.72      0.41     2.09     3.61       1774 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

また、サンプリング結果を格納したオブジェクトをsummary()に入れると要約された結果が返されますが、brms::waic()やbrms::loo()に入れると、WAICやlooを計算することもできます（内部ではlooパッケージを使用して計算しています）。

brms::waic(fit)

Computed from 4000 by 32 log-likelihood matrix

          Estimate  SE
elpd_waic    -78.9 3.9
p_waic         4.1 1.0
waic         157.8 7.9
Warning message:
4 (12.5%) p_waic estimates greater than 0.4. We recommend trying loo instead. 

brms::loo(fit)

Computed from 4000 by 32 log-likelihood matrix

         Estimate  SE
elpd_loo    -79.0 4.0
p_loo         4.3 1.1
looic       158.1 7.9
------
Monte Carlo SE of elpd_loo is 0.1.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     31    96.9%   1187      
 (0.5, 0.7]   (ok)        1     3.1%   855       
   (0.7, 1]   (bad)       0     0.0%   <NA>      
   (1, Inf)   (very bad)  0     0.0%   <NA>      

All Pareto k estimates are ok (k < 0.7).
See help('pareto-k-diagnostic') for details.

事前分布の表示

brmsは、いくつかのパラメータについてデフォルトで事前分布を与えていることがあります。今回指定したモデルで、各パラメータにどのような事前分布が与えられていたのかを知るには、brms::get_prior()を使用します。この中にモデルを指定すると、以下のように一覧が表示されます。

brms::get_prior(mpg ~ wt * am, data = mtcars)

                 prior     class  coef group resp dpar nlpar bound
1                              b                                  
2                              b    am                            
3                              b    wt                            
4                              b wt:am                            
5 student_t(3, 19, 10) Intercept                                  
6  student_t(3, 0, 10)     sigma                 

Interceptの事前分布student_t(3, 19, 10)、19ってどういう根拠で設定されたんでしょうね...。brmsはこんな風に、いくつかのパラメータにデフォルトで弱情報の事前分布を与えていることがあるのですが、どういう根拠に基づいているのかは、まだよく分かっていません。

またbrmsでは基本的に（おそらく常に？）、回帰係数bには事前分布は明示されません。つまりデフォルトでは、無情報事前分布が適用されるということになります。これは後述する「bayes factorの計算」においては問題となるため、そのような場合にはユーザが事前分布を指定する必要がでてきます。

事前分布の指定

ユーザが任意の事前分布を指定したければ、brms::brm()のなかでprior =という引数にパラメータと事前分布をセットで指定します。指定方法には複数の記法があります。詳しくは、コンソール上で?brms::set_priorと入力して、ヘルプをご参照ください。

fit2 <- brm(mpg ~ wt * am,
            data = mtcars,
            iter = 2000,
            warmup = 1000,
            seed = 1234,
            chain = 4,
            prior = c(prior_string("normal(0, 10)", class = "b"),
                      prior(student_t(3, 19, 10), class = Intercept),
                      prior_(~student_t(3, 0, 10), class = ~sigma)
                      )
            )

Stanコードの表示

brms::make_stancode()のなかに、事前分布を含み任意のモデルを指定すると、そのStanコードをコンソールに表示してくれます。例えば先ほど事前分布を指定した、以下のモデルのStanコードを表示してみましょう。長いので、ここではmodelブロックだけ掲載します。

ちゃんと、ユーザが指定した事前分布も反映されていますね（target += normal_lpdf(b | 0, 10)という部分）。

brms::make_stancode(mpg ~ wt * am,
                    data = mtcars,
                    prior = c(prior_string("normal(0, 10)", class = "b"),
                              prior(student_t(3, 19, 10), class = Intercept),
                              prior_(~student_t(3, 0, 10), class = ~sigma)
                              )
                    )

# 以下、出力 --------------------------------------------
model { 
  vector[N] mu = temp_Intercept + Xc * b;
  // priors including all constants 
  target += normal_lpdf(b | 0, 10); 
  target += student_t_lpdf(temp_Intercept | 3, 19, 10); 
  target += student_t_lpdf(sigma | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10); 
  // likelihood including all constants 
  if (!prior_only) { 
    target += normal_lpdf(Y | mu, sigma);
  } 
} 

そうそう、brmsのStanコードは全て、target +=記法を用いて書かれています。

結果の可視化

plot()

brmsは、サンプリング結果の視覚化にも長けています。rstan::stan_trace()やrstan::stan_dens()のように、描きたいグラフに応じて個別の関数が用意されていることが多いと思いますが、brmsはbayesplotパッケージのbayesplot::mcmc_combo()を用いて結果を可視化します。bayesplotパッケージについては、Stan Advent Calendar 2018の2日目の記事や、11日目の記事をご覧ください。

bayesplot::mcmc_combo()は、各パラメータについて、事後分布とトレースプロットを並べて表示してくれます。

plot(fit)

繰り返しますがこのグラフはbayesplotパッケージを用いて出力されているので、同パッケージの関数を用いて見栄え（Aesthetic）を調整することができます。また、bayesplot::mcmc_combo()と同様に、bayesplotパッケージのmcmc_ シリーズの関数のうち、どの種類の関数を組み合わせて表示するかも選択することができます（mcmc_シリーズについては、Stan Advent Calendar 2018 2日目の記事を参照してください）。ただしmcmc_combo()で表示できるグラフの種類は限られているようです。

例えばこんな風に、色調をピンクに変更し、トレースプロット（bayesplot::mcmc_trace()）とチェインごとの事後分布（bayesplot::mcmc_dens_overlay()）を並べることもできます。

bayesplot::color_scheme_set("pink")
plot(fit, combo = c("trace", "dens_overlay"))

marginal_effects()

※注意：brms 2.11.0より、brms::marginal_effect()がbrms::conditional_effects()に名称変更されています。詳しくはこちら

さらにbrms::marginal_effects()を用いると、「主効果」や「交互作用」を可視化することもできます。今回は交互作用項を含むモデルを指定したので、いわゆる「単純傾斜」も自動的に可視化されます。下のグラフは、左から順に「wtの主効果」「amの主効果」「wtとamの交互作用」を示します。

引数effectsに何も指定しなければ、すべての項について可視化が行われますが、特定の項だけを指定することも可能です。また、交互作用項を指定する場合には、その順番を変えることで、どちらの変数についての単純効果を可視化するかも決めることができます。

marginal_effects(fit, effects = "wt:am") #左のグラフ
marginal_effects(fit, effects = "am:wt") #右のグラフ

離散変数の場合は（今回はam）、各水準における他方の変数の単純効果が可視化されます。連続変数の場合は（今回はwt）、平均±1SDのポイントにおける、他方の変数の単純効果が示されます。

mean(mtcars$wt) + sd(mtcars$wt) #平均+1SD
[1] 4.195707

mean(mtcars$wt) #平均
[1] 3.21725

mean(mtcars$wt) - sd(mtcars$wt) #平均-1SD
[1] 2.238793

事後予測チェック

事後予測チェックについては、Stan Advent Calendar 2018 11日目の記事も参照してください。Stanでモデルを書いた場合、事後予測チェックを行うためには、generated quantitiesブロックでモデルから乱数を生成させる手続きが必要になります。その乱数を抽出して、実際に観測されたデータとの類似性を把握するわけですが、なんとbrmsでは、サンプリング結果が格納されたオブジェクトをbrms::pp_check()に入れるだけで、事後予測チェックができます！

brms::pp_check(fit)

Stan Advent Calendar 2018 11日目の記事を読んだ方はお察しの通り、これは内部でbayesplot::ppc_dens_overlay()を使用しています。

他の事後予測チェックにも対応していますよ。type =という引数に、描画形式を指定してみます。

brms::pp_check(fit, type = "error_hist")

これはbayesplot::ppc_error_hist()を使用していますね。

一般化線形モデル

次は、正規分布以外の指数分布族を扱えるように拡張した、一般化線形モデルを例に挙げます。今度は目的変数を、オートマ車かマニュアル車かamにし、説明変数に燃費mpgと重量wtを投入してみます。

stats::glm()と同様に、目的変数が従う確率分布を、引数family =に指定しましょう。今回は、オートマ車かマニュアル車かamは2値の変数なので、family = "bernoulli"とベルヌーイ分布を指定するのがよいように思われます（余談ですが、stats::glm()では2値の目的変数に対しては二項分布（family = "binomial"）を指定しますが、brmsでは2値の目的変数に対して二項分布を指定すると、Only 2 levels detected so that family 'bernoulli' might be a more efficient choice.とサジェストされます）。

ただしbrmsでは、family = "binomial"を指定した際の記法がstats::glm()と異なるので、ここでは紹介のためあえて二項分布を適用した書き方をしてみます。

stats::glm()ではこのように書きます。

#amをfactor型に変換していたのを、いったんリセット
data(mtcars)

fit <- glm(cbind(am, 1 - am) ~ mpg + wt,
           family = "binomial",
           data = mtcars)
#目的変数が2値の場合は、
#fit <- glm(am ~ mpg + wt, family = "binomial", data = mtcars) でOK

summary(fit)

cbind()のなかで何をしているのかというと、総試行数中に占める、成功試行数（am）と失敗試行数（1 - am）の列を指定しています。今回は1台の車がオートマ車かマニュアル車かを予測するので、試行数は1であるため、失敗試行の列は1 - amとしています。

さて、brmsでは、目的変数が二項分布に従うと仮定するとき、以下のように書きます。|の左側には成功試行数を、右側にはtrials()の中に試行数を入れます。

#amをfactor型に変換していたのを、いったんリセット
data(mtcars)

fit <- brm(am | trials(1) ~ mpg + wt,
           family = "binomial",
           data = mtcars)
#目的変数が2値の場合は、
#fit <- brm(am ~ mpg + wt, family = "bernoulli", data = mtcars) でOK

summary(fit)

brms::marginal_effects()で主効果を可視化してみると、ロジットリンク関数によって、直線ではなく曲線になっていることがわかりますね。

豊富なbrms family

ゼロ過剰ポアソン分布

さて、ここからがbrmsの真骨頂です。stats::glm()ではfamilyに指定できる確率分布は以下の通りです。

• gaussian
• binomial
• Gamma
• inverse.gaussian
• poisson
• quasi

が、brmsでは豊富なfamilyが用意されているということを、冒頭で述べました。試しにちょっと珍しい分布を使ってみましょう。

Mr.Unadon氏がTwitter上で「あなたは、片方しかない靴下を何足もっていますか？」とアンケートをとった結果を拝借します。ネタ元はこちらの片方しかない靴下の数をゼロ過剰ポアソン分布でモデリングしてみたです。

dat <- data.frame(Y = c(rep(0, 746),
                        rep(1, 142),
                        rep(2, 142),
                        rep(3, 154),
                        rep(4, 23),
                        rep(5, 22),
                        rep(6, 9),
                        rep(7, 66)
                        )
                  )

set.seed(1234)
dat$X <- 50 + 2.5*dat$Y + rnorm(n = nrow(dat), mean = 0, sd = 10)

YがUnadon氏が集めた、片方しかない靴下の所持数のデータです。今回のサンプルの中では、ほとんどの人は靴下を片方だけなくすということはないようですね。僕が初めてこのアンケートを見た時、「靴下片方だけどうやって失くすねん」と思ったのですが、最近我が家で片方しかない靴下が発見されました。

さて、靴下の所持数はカウントデータなので、ポアソン分布と相性がよさそうに思えますが、上図のとおりゼロが圧倒的に多く、通常のポアソン分布の仮定を満たしそうにありません。そこで、ゼロが通常のポアソン分布よりも過剰であるという特徴を捉えた、ゼロ過剰ポアソン分布（zero-inflated poisson）という確率分布を考えます。

ゼロ過剰ポアソン分布は、ベルヌーイ分布とポアソン分布の混合分布です（以下の図はUnadon氏の記事から拝借しました）。靴下を失くすか否かがベルヌーイ分布に従うとと考え、もし失くさないなら当然、片方しかない靴下の所持数の回答は0。もし失くすなら、所持数はポアソン分布に従う、と考えます。なお、ポアソン分布に従う確率変数は0以上の整数をとりうることに注意してください。

出典：片方しかない靴下の数をゼロ過剰ポアソン分布でモデリングしてみた

一方Xは適当に作った変数ですが、仮に「日常の忙しさ指標（0-100点）」だとしておきましょう。仮定としては、忙しい人ほど靴下を片方なくしやすい、といった具合です。

さて、片方しかない靴下の所持数を、忙しさで予測する、シンプルなモデルを考えてみましょう。ここで、Stanでモデルを書くとしたら、Unadon氏の記事にあるように、target +=という記法を用いて、少し工夫した書き方が必要になります。

しかしbrmsでは、family = "zero_inflated_poisson"を指定するだけです。

fit <- brm(Y ~ X,
           data = dat,
           family = "zero_inflated_poisson",
           seed = 1234)

summary(fit)

Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: Y ~ X 
   Data: dat (Number of observations: 1304) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    -1.64      0.16    -1.98    -1.33       2039 1.00
X             0.04      0.00     0.04     0.05       2283 1.00

Family Specific Parameters: 
   Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
zi     0.47      0.02     0.43     0.50       1868 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

事後予測チェックしてみましょう。

brms::pp_check(fit)

いい感じですね。

階層線形モデル

本記事冒頭でも言及したように、lme4パッケージやlmerTestパッケ―ジを用いると、以下のように書きますよね。

• lme4::lmer(y ~ x + (1 | group), data = dat)
• lme4::glmer(y ~ x + (1 | group), family = "poisson", data = dat)

brmsでは、関数をbrm()に変えるだけなので、本記事では説明を省略します。

モデル比較

brmsパッケージを用いてサンプリングした結果を利用して、モデル比較を行ってみます。モデル比較には様々な観点がありますが、ここではbayes factorを指標とすることにします。bayes factorについては、北條大樹さんのこちらの記事や、清水裕士さんのスライドをご参照ください。

brmsを用いてモデル比較するために、kidneyというデータセットを用いて、試してみましょう。survivalパッケージにも同じ名前のデータセットがあるので、区別するためにbrms::kidneyと指定したほうが無難です。

話を簡単にするために、ここでは目的変数に「病気の再発までの時間time」、説明変数に「患者の年齢age」を投入した、単回帰分析を考えます。また、本当は右側打ち切り（right-censored）が起こっているため、それを考慮したモデリングをするべきですが、打ち切りが起こっていないデータだけを抽出しています。
なおbrmsでは、brm(y | cens(censored) ~ x, data = dat)のように書くことで、打ち切りデータのモデリングも容易に実現できます（cens()のなかに、打ち切りの有無を判別する変数を投入します）。

#データセットの作成 ---------------------------------------------
#1回目の再発までの時間 & 非打ち切りデータのみを使用
dat <- brms::kidney %>%
  dplyr::filter(recur == "1" & censored == 0) 

#データセットの行列数 -------------------------------------------
dim(dat) #38人分のデータが抽出された
[1] 38  7 

#冒頭6行分を表示 ------------------------------------------------
head(dat)

  time censored patient recur age    sex disease
1    8        0       1     1  28   male   other
2   23        0       2     1  48 female      GN
3   22        0       3     1  32   male   other
4  447        0       4     1  31 female   other
5   30        0       5     1  10   male   other
6   24        0       6     1  16 female   other

目的変数timeの分布はこのように左に歪んでいます。

正規分布を仮定したモデル

まずはベーシックな単回帰分析を実行してみます。目的変数が正規分布に従うと想定することになるので、family = "normal"と指定します（familyを指定しなければ自動的に正規分布を仮定することになるので、省略も可能です）。

さらに以下3点を対処する必要があります。

• 全パラメータに、無情報でない事前分布を指定
• iterを多めに
  • 正確には、有効なMCMCサンプルの数を多めに（2109/12/6追記）
• save_all_pars = TRUEと引数を指定

モデルはこのようになります。

fit_n <- brm(time ~ age,
             family = "normal",
             prior = prior(normal(0, 5), class = b),
             seed = 1234,
             iter = 100000,
             warmup = 5000,
             data = dat,
             save_all_pars = TRUE)

事後予測チェックを行ってみましょう。

brms::pp_check(fit_n)

実際に観測されたデータとズレがありますね。上のヒストグラムを見ればわかる通り、今回はtimeが正規分布に従うと仮定するのは、適切ではないかもしれません。

最後にbrms::bridge_sampler()を用いて、bridge samplingによる対数周辺尤度を計算してみましょう。これは内部ではbridgesampling::bridge_sampler()を利用しています。

brms::bridge_sampler(fit_n)

Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Bridge sampling estimate of the log marginal likelihood: -249.9636
Estimate obtained in 4 iteration(s) via method "normal".

ワイブル分布を仮定したモデル

brms::kidneyは生存時間解析に適したデータセットなので、比較対象として、パラメトリックな生存時間解析でよく用いられる、ワイブル分布を適用したモデルを考えてみます。

fit_w <- brm(time ~ age,
             family = weibull(),
             prior = prior(normal(0, 5), class = b),
             seed = 1234,
             iter = 100000,
             warmup = 5000,
             data = dat,
             save_all_pars = TRUE)

事後予測チェック。

brms::pp_check(fit_w)

先ほどよりマシになった気がしますね。

次はbridge samplingによる対数周辺尤度の計算。

brms::bridge_sampler(fit_w)

Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Bridge sampling estimate of the log marginal likelihood: -232.0824
Estimate obtained in 4 iteration(s) via method "normal".

正規分布 vs ワイブル分布

いよいよモデル比較です。各モデルの対数周辺尤度に、いったん-1を掛けます。つまり、自由エネルギーとして計算します。 「自由エネルギーの差」の指数をとるとbayes factorになります。

fe_n = -1 * brms::bridge_sampler(fit_n)$logml #正規モデルの自由エネルギー
fe_w = -1 * brms::bridge_sampler(fit_w)$logml #ワイブルモデルの自由エネルギー

exp(fe_n - fe_w) #bayes factor

[1] 58278089

今回のデータに対する当てはまりの観点では、正規分布を仮定したモデルに対して、ワイブル分布を仮定したモデルを支持する程度が、58278089倍...もう圧倒的に大きいと考えられます。

ただしbrmsは、bridgesampling::bayes_factor()を利用して直接bayes factorを計算する、brms::bayes_factor()という関数を用意しています。引数x1とx2に、各モデルのサンプリング結果を格納したオブジェクトを指定してください。

brms::bayes_factor(x1 = fit_w, x2 = fit_n)

Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Estimated Bayes factor in favor of bridge1 over bridge2: 58273598.05176

やはり今回のデータに対する当てはまりの観点では、正規分布を仮定したモデルよりも、ワイブル分布を仮定したモデルのほうが良いと考えられます。

こんな風に、brmsを用いてモデル比較をすることもできます。ただしモデル比較という手続きは非常に難しいので、実際にはもっと慎重に行う必要があると思います。 以下、北條さんの記事から引用します。

ベイズファクターの注意点
- ベイズファクターは事前分布の影響をかなり大きくうけます ｡
　→ 同一データ､同一モデルで事前分布のみを変化させた場合に､
　　どれだけベイズファクターが変化するかを検討するととてもよくわかります｡
- bridgesamplingパッケージでも計算がうまく行かないことはありますので､
　全てのモデル間のベイズファクターが計算できるわけではありません｡
- モデル比較を行うのはとても難しいです｡
　本来､我々が､比較すべきモデルは無限に存在します｡
- ベイズファクターは､2つのモデル間の比として考えるので､
　それ以外のモデルについては考えてはいません｡

2019/12/6追記
清水裕士さんによる、Stan Advent Calendar 2019 6日目の記事：brmsパッケージで安易にベイズファクターを使うと死ぬ話で、brmsでベイズファクターを求める際の注意点が書かれています。

その他のtips

Stanコードの保存

brms::brm()でサンプリングを行う際に、save_model = "hogehoge"と引数を追加してみましょう（もちろんhogehogeのところは任意の名前で構いません）。

fit_w <- brm(time ~ age, data = dat, save_model = "hogehoge")

するとカレントディレクトリに、hogehogeという名前のファイルができていると思います。R上でこのファイルを開くと、brms::make_stancode()で生成できるようなStanコードが記録されていることが分かります。

仮にbrmsのバージョンが変わったとしても、少なくとも今回のサンプリングを実行した時点でのStanコードは判明しているので、結果の再生可能性の向上に寄与するかもしれません。

サンプリング結果の保存

brms::brm()でサンプリングを行う際に、file = "fugafuga"と引数を追加してみましょう（もちろんfugafugaのところは任意の名前で構いません）。

fit_w <- brm(time ~ age, data = dat, file = "fugafuga")

するとカレントディレクトリに、fugafuga.rdsという名前でサンプリング結果が保存されています。 以下のようにreadRDS()という関数で、あとでロードすることができるので、「しまった、うっかりRStudioを閉じてしまった！」なんて心配も不要です。また、例えばサンプリング結果をRMarkdownでまとめたい場合などにも重宝すると思います。

fit_again <- readRDS("fugafuga.rds")

関数名の前に「brms::」と付けるクセをつける

あくまで個人的な意見ですが、何かしらのパッケージ内の関数を使用する場合には、パッケージ名を併記したほうがいいと思っています（例えばdplyr::select()とか）。これはいくつかの理由がありますが、一番の理由は、たまに複数のパッケージ間で同じ名前の関数が存在することがあるからです。例えばfilter()という関数は、少なくともstatsパッケージとdplyrパッケージに存在します。

また、他者とコードを共有するときに、何のパッケージの関数か分からないと混乱を招くと思うので、そういった意味でもパッケージ名を明示しておいた方がいいかなと思っています。

もちろん、文字数が増えるのでコードが見にくくなるという側面もありますが。

ただ、brmsを用いる際には、特に意識してbrms::と書いたほうがよいと思います。これは何故かというと、brmsは様々なパッケージに依存しているのですが、データセット名や関数名が衝突することがあるからです。例えば以下のようなデータセットや関数たち。

• brms::kidneyデータセットとsurvival::kidneyデータセット
• brms::bayes_factor()とbridgesampling::bayes_factor()
• brms::bridge_sampler()とbridgesampling::bridge_sampler()
• brms::bf()とbridgesampling::bf()

一番タチが悪いのが最後のbf()です。bridgesampling::bf()はbayes factorを計算する関数で、bridgesampling::bayes_factor()の略称なのですが、brms::bf()はbayes factorとは関係がありません。brms::bf()は、brms::brmsformula()の略称で、brmsを用いて高度なモデリングを可能にするための関数です。

分位点回帰を例に説明しましょう。brmsは、分位点回帰も容易に実行することができます。その場合、asymmetric Laplace分布を適用して、以下のように書きます。brm()のなかに、さらにbf()という関数を書いて、その中でモデルを指定します。

fit <- brm(bf(y ~ x, quantile = 0.2),
           data = dat,
           family = asym_laplace()
           )

これ、bridgesampling::bf()のほうが読み込まれてしまうと、当然エラーになります。したがって以下のように、brms::bf()であることを明示したほうが無難です。

fit <- brm(brms::bf(y ~ x, quantile = 0.2),
           data = dat,
           family = asym_laplace()
           )

あるいは、パッケージの優先順位をつけてもいいかもしれないですね。
参考記事：R でパッケージの優先順位を変えたい #rstatsj

おわりに

ここで紹介した機能は、brmsのほんの一部にすぎません。モデリングの自由度もかなり高いですし、ユーザも多いのでWeb上で色んな情報を得ることができます。 例えばこちらの記事では、brmsを用いて信号検出理論のパラメータを推定しています。

きっとこれからも、どんどん新しい機能が追加されていくと思うので、要チェックですね。

そういえば、この記事を書くためにbrmsについて調べていたら、brmsとggplot2とtidyverseを用いてStatistical Rethinkingを解説しているドキュメントを見つけました。
Statistical Rethinking with brms, ggplot2, and the tidyverse

ざっと見たところ、brmsの具体的な使用例についても紹介しているようです。参考にしてみてください。

Enjoy!!