rstanでお手軽にGLMMができる関数，glmmstan（）を作りました。

 

今回はRの話です。

社会心理学会の方法論セミナーでもGLMMをとりあげましたが，階層ベイズの話も久保先生のトークの中にありました。

GLMMでは，変量効果の数が増えると最尤法だと推定が難しくなったりするので，ベイズ推定のほうが向いています。

しかし，GLMMを直接ベイズ推定してくれる商用ソフトもあまりない（あっても機能が部分的）ので，stanなどのフリーのソフトに頼らざるを得ません。しかし，stanは初心者にはなかなかとっつきにくいので，今回はだれでも簡単にGLMMがベイズ推定できる関数を作ってみました。

実は過去に，同様にGLMMを簡単にstanで走らせてくれるglmer2stanというパッケージを紹介したことがあります（こちら）。しかし，glmer2stanはあとで挙げるようにいろいろ使い勝手が悪いところもあり，自分用に使いやすいものを作ろうと思ったのがはじまりです。

追記：

SapporoRで発表したスライドをスライドシェアにアップしています。こちらの資料のほうが詳しく書いているので，こちらも合わせて参照してください。

rstanで簡単にGLMMができるglmmstan()を作ってみた

 

glmmstan（）とは

R上で動くstanの代筆関数で，rstanパッケージが入っていることが前提となります。

こちらのページをのコードを保存して自由に使ってください。あるいは，

 

を走らせてもOKです。

lme4の文法でモデルを書けば，あとは自動的に内部でstanコードを作成し，MCMCを走らせます。

またすべてのモデルについてWAICを出力するので，モデル比較を柔軟に行うことができます。（ただ，対数正規分布のモデルについてだけは出力しません）

この関数については，６月１３日のSapporoR#4で発表するので興味ある方はぜひ来てください。

 

glmer2stanとの違い

glmer2stanは便利ですが，コードの微調整がしにくかったり，個人的に出力が気に食わなかったり，謎のエラーがでたりと，使いにくいところもあります（それがこの関数の作成の動機です）。

以下にglmer2stanとの違いを書いておきます。

１．glmmstanのほうが1.5～2倍は計算速度が速い
２．stanコードが変数名や固定効果の数に依存しないのでコンパイルを避けられる
３．stanのコードを簡単にいじれる
４．使えるfamilyが違う（後述）
５．標準偏差パラメータの事前分布にcauch分布を使っている
６．ポアソンと負の二項分布にオフセット項を利用できる
７．交互作用項がある場合，単純効果をパラメータとして推定できる

などがあります。

 

glmmstan()の基本的な使いかた

glmmstan()は，基本的にはlme4と同じ文法でモデルを入力してMCMCで推定することができるものです。その点についてはglmer2stanと全く同じです。

使いかたは，下のように，

 

という感じで使います。

x1，x2は固定効果の推定に用いるものです。（）内のモデルは変量効果で，groupのところに個人や集団，項目といった変量効果の単位を識別する変数を入れます。1は切片を，z1などは傾きの変動を見たい変数を指定します。もちろんz1とx1が同じ変数でも構いません（だいたいは同じ事のほうが多い）。

なお，変量効果は指定しなくてもOKで，その場合はただの一般化線形モデルになります。

familyは目的変数の残差が従う分布を指定します。選べるのは正規分布（gaussian)，ベルヌーイ分布(bernoulli)，二項分布（binomial）,ポアソン分布（poisson），負の二項分布(nbinomial)，ガンマ分布(gamma)，対数正規分布(lognormal)，ベータ分布（beta），順序カテゴリカル（ordered）です。今後増やしていこうと思うので，これを追加してほしい！とかあれば連絡ください。

戻り値はstanfitクラスのオブジェクトで，stanが返すものと同じです。なのでそのままMCMCリストを取り出したりができます。

ただ，いくつかglmmstan用の出力をくっつけているので，それらを参照したい場合は以下の出力に関する関数を利用してください。

 

 出力に関する関数

stanfitオブジェクトをある程度要約した結果を知りたい場合は，output_result()を使います。

 

戻り値はリストで，betaは固定効果の結果を，tau_sdは変量効果の標準偏差，tau_corrは変量効果の相関行列，tauは共分散行列を返します。

それ以外に，WAICも出力します。2N*WAICはAICやDICなどに近くなるような定義で計算されたものです。

結果の一部を知りたい場合は次のように書くといいでしょう。たとえば固定効果だけを知りたい場合は，

 

という感じです。$WAICとすればWAICだけをとりだせます。

大局パラメータのみのrhatなどを確認したい場合は，output_stan()関数を使うと便利です。

 

これを使うと，固定効果，残差，変量効果の標準偏差，相関，共分散の要約統計量やrhatを確認できます。

あと，stanコードを確認したい場合は，output_code()を使います。

 

これで走らせたstanのコードが見れます。

あと，output_ggmcmc()を使うと，固定効果についてはパラメータ名を変数に対応するように変更してggmcmc用のオブジェクトを返します。

 オプション

stan（）のオプションはすべてglmmstan（）でも入力できます。とくによく使うものを紹介しておくと，iter=でサンプリング回数，warmup=でバーンイン回数，chains=でMCMCの数，thin=で何回ごとに採択するかを変更できます。glmmstanのデフォルトはiter=2000, chains=2,thin=1です。warmupは指定しないとiterの半分になります。

glmmstanはstanの結果を返すだけでなく，いろいろ出力できます。

オプションで，dataonly=Tとすればモデルに従って作成されたstanに入力する用のデータをリスト型で返します。

codeonly=Tとすればモデルに従って作成されたstanコードを返します。

modelonly=Tとすれば，コンパイルだけしたstanmodelクラスのオブジェクトを返します。あとで書くように，これを作っておくと固定効果を変えるだけなら毎回コンパイルの必要がありません。

 

 モデルについてのオプション

glmmstan()では，心理学などの社会科学者向けに交互作用効果を検討した場合，単純効果もパラメータとして推定できます。また，center=Tとすれば説明変数を全体平均で中心化します。

 

と書けば，x1とx2の交互作用効果を中心化した状態で作成し，slice=で指定した変数で群分けした場合の単純効果の事後分布が計算できます。

さらに，offset =とすればオフセット項を追加できます。いまのところオフセット項の設定ができるのはポアソンと負の二項分布のみです。

 

 コードを少しだけいじって使う

stanに使い慣れている人は，自分なりのカスタマイズがしたいと思います。glmmstanではstancode=のあとにstanコードが入ったオブジェクトをいれると入力されたコードを使ってstanを回します。

 

というように書けば，内部でコードを作らず，入力されたもので分析します。

同様に，standata =を使えばstanに入力する用のデータを優先して入力できます。また，すでにコンパイルしたstanmodelがあれば，

 

と書くと，コンパイルなしにMCMCサンプリングが始まります。またstanfitクラスのオブジェクトをstanfit =でいれてやれば，同様にコンパイルなしで分析できます。

なお，glmmstanではコードは固定効果については変数の数，名前に依存しないものを使っていますので，固定効果を増やしたり違う変数を使ってもstanfit=などを使えばコンパイルなしで分析できます。

 

変量効果については，いくつか制限があって，次のようなときにコードが変わります。

１．変量効果の単位を増やす
２．変量効果が１つの場合と，２つ以上の場合（つまり（）内で指定した変数が１だけの場合と，＋z１・・を書いた場合でコードが違う）

ただし，これらも変数名には依存しません。

今回は以上です。今後いろいろ改良しようと思いますので，他にも改良してほしい点などがあれば，連絡ください。

あと，まだまだバグがあるかと思うので，変だな，と思ったら同じく連絡いただければと思います。