すでに関数自体は作っていたglmmstanですが，それをパッケージ化しました。
関数名と同じ，glmmstanパッケージです。
インストール
glmmstanパッケージは，CRANに登録していないので，install.packages()ではインストール出来ません。
そこで，GitHubにおいてあるパッケージをインストールするためのdevtoolsパッケージをまずインストールする必要があります。
|
1 |
install.packages("devtools") |
続いて，次のコードを実行します。
|
1 2 |
library(devtools) install_github("norimune/glmmstan") |
すると，glmmstanパッケージがインストールされます。
使いかた
パッケージをインストールしたら，glmmstan()を使ってGLMMをMCMCで推定できるようになります。
詳しい使いかたについては，SapporoRで発表したスライドをスライドシェアにアップしています。こちらの資料に詳しく書いているので，参照してください。
rstanで簡単にGLMMができるglmmstan()を作ってみた
一応，この記事でもglmmstan関数の使いかたを書いておきます(前に書いた記事をコピペしただけですが）。
glmmstan()の使い方
glmmstan()は，基本的にはlme4と同じ文法でモデルを入力してMCMCで推定することができるものです。その点についてはglmer2stanと全く同じです。
使いかたは，下のように，
|
1 |
fit <- glmmstan(y~x1+x2+(1+z1|group),dat,family="gaussian") |
という感じで使います。
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()を使います。
|
1 |
output_result(fit) |
戻り値はリストで，betaは固定効果の結果を，tau_sdは変量効果の標準偏差，tau_corrは変量効果の相関行列，tauは共分散行列を返します。
それ以外に，WAICも出力します。2N*WAICはAICやDICなどに近くなるような定義で計算されたものです。
結果の一部を知りたい場合は次のように書くといいでしょう。たとえば固定効果だけを知りたい場合は，
|
1 |
output_result(fit)$beta |
という感じです。$WAICとすればWAICだけをとりだせます。
大局パラメータのみのrhatなどを確認したい場合は，output_stan()関数を使うと便利です。
|
1 |
ouput_stan(fit) |
これを使うと，固定効果，残差，変量効果の標準偏差，相関，共分散の要約統計量やrhatを確認できます。
あと，stanコードを確認したい場合は，output_code()を使います。
|
1 |
output_code(fit) |
これで走らせたstanのコードが見れます。
オプション
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とすれば説明変数を全体平均で中心化します。
|
1 |
fit <- glmmstan(y~x1+x2+x1:x2+(1|group),dat,family="gaussian",center=T, slice="x2") |
と書けば，x1とx2の交互作用効果を中心化した状態で作成し，slice=で指定した変数で群分けした場合の単純効果の事後分布が計算できます。
さらに，offset =とすればオフセット項を追加できます。いまのところオフセット項の設定ができるのはポアソンと負の二項分布のみです。
コードを少しだけいじって使う
stanに使い慣れている人は，自分なりのカスタマイズがしたいと思います。glmmstanではstancode=のあとにstanコードが入ったオブジェクトをいれると入力されたコードを使ってstanを回します。
|
1 |
fit <- glmmstan(y~x1+x2+(1+z1|group),dat,family="gaussian",stancode=code1) |
というように書けば，内部でコードを作らず，入力されたもので分析します。
同様に，standata =を使えばstanに入力する用のデータを優先して入力できます。また，すでにコンパイルしたstanmodelがあれば，
|
1 |
fit <- glmmstan(y~x1+x2+(1+z1|group),dat,family="gaussian",stanmodel=model1) |
と書くと，コンパイルなしにMCMCサンプリングが始まります。またstanfitクラスのオブジェクトをstanfit =でいれてやれば，同様にコンパイルなしで分析できます。
なお，glmmstanではコードは固定効果については変数の数，名前に依存しないものを使っていますので，固定効果を増やしたり違う変数を使ってもstanfit=などを使えばコンパイルなしで分析できます。
|
1 |
fit <- glmmstan(y~x1+x2+(1+z1|group),dat,family="gaussian", stanfit=fit1) |
変量効果については，いくつか制限があって，次のようなときにコードが変わります。
１．変量効果の単位を増やす
２．変量効果が１つの場合と，２つ以上の場合（つまり（）内で指定した変数が１だけの場合と，＋z１・・を書いた場合でコードが違う）
ただし，これらも変数名には依存しません。
今回は以上です。今後いろいろ改良しようと思いますので，他にも改良してほしい点などがあれば，連絡ください。
あと，まだまだバグがあるかと思うので，変だな，と思ったら同じく連絡いただければと思います。