2018の5日目の記事です。今日は書く予定ではなかったんですが、カレンダーが空いてたので埋めてみよう、という記事です。そんなに準備はできてないです。竹林さんの3日目の記事、tidybayesは便利そうですね。tidyなみなさんはぜひこちらを使ってほしいんですが、どうも僕は研究ではそもそもデータがきれいなものを扱うことが多いので、そんなにtydyverseにお世話になっていないです。いうこともあり、dplyrやggplot2などの文法に慣れていないので、Stanを使うとこもapplyとplotを使うことが多いです。いうわけで、今回は、僕的にStanの結果を非tidyverseで見るための工夫をメモ代わりに書いておこうと思います。知ってることもあると思いますし、tidybayesパッケージとかのほうが便利なこともあるかと思いますが、つらつら書いておきます。基本的に「俺流」なので、あまり参考にならないかもですが、こういう記事もいいかなと思って。◆今回の例で使うモデルこの記事の読者はむしろ心理統計にあまり関心がない人が多いかもしれませんが、今回の例では項目反応理論を挙げておきます。GPCM)と呼ばれる、順序尺度データに対応した項目反応理論を紹介します。GPCMはある人iが項目jのカテゴリc(0,1,...K)に反応する確率を次のような関数で表せるとします。ここで、alphaは識別力パラメータといって、項目がテスト全体の傾向を反映してるかどうかの指標です。betaは閾値パラメータといって、各反応カテゴリの閾値を表しています。このGPCMをStanでモデルどおりに書くと次のようになります。}}走らせるRコードは以下です。GPCMに従って乱数を生成して、それをStanで推定している感じです。sampling(model.irt5,data=datastan,iter=2000,warmup=1000,chains=4,cores=4)◆僕が使っているパッケージ僕がStanを使うときに使うパッケージは以下です。僕は実はあまりパッケージを使いたくないタイプでして(仕様が変わって使い勝手が変わったり、再現できないのが嫌なタイプ)、極力素のRでコードを書く人間です。library(rstan)library(psych)library(magrittr)library(bridgesampling)rstanを使えばggplot2は勝手に使えるので書いてません。dplyrは使わないのですが、stanfitオブジェクトがアホみたいに容量がでかいので、あまり複製をしたくないため、パイプ演算子は必須です。magrittrだけ使うことが多いです。psychパッケージはlogistic関数とか、基本的なのがはいってるので常に入れてます。bridgesamplingはベイズファクターの計算で使います。◆前もって作っておくと便利な関数MCMCの結果を見るときに、Rにははじめから入ってなかったり、あるけど重くて使いにくいのがいくつかあるので、ここで紹介しておきます。MAP推定値をかえすmap()という関数を作ります。これはStanユーザーでは有名な処理です。}カーネル密度推定を利用して、密度の最大値を推定するものです。続いて、WAICを返す関数です。looパッケージにwaic関数があるんですが、これ、よくRstudioが爆発するんです(やっと推定終わったーさぁWAICでも見るかーボカーンという地獄)。いうわけで自分で定義したものを僕は使ってます。return(list(waic=waic,lppd=lppd,p_waic=p_waic))}◆mcmcの結果を見るではさっそくMCMCの結果を見てみましょう。僕はパイプ演算子が好き、かつ、一つの処理を1行で書きたいという謎の癖を持っています。Rhat=1).僕はprintが好きです(謎)。digits=で簡単に桁数も変えられます。rhatだけを見たい場合は、stan_rhatが便利です。stan_rhatうまくいってそうです。自己相関の程度を見るのはessが便利。stan_essStan2.18からESSの計算方法が変わって、1を超えることがあるようですね。WAICは次のようにやります。さっき作ったwaic2関数を使います。waic2◆パラメータを取り出す識別力母数alphaは項目ごとにあるので、applyを使います。別に何も特別なことではないですけども。ここでは先程作ったmap関数を使って、MAP推定値を取り出してみます。data.frame(alpha=.)縦に表示したいがために、最後にdata.frameを使ってます。別にここはいらないです。mean)とかになります。1.92最後にprintで桁数を変えるのもポイント。識別力母数を項目ごとに事後分布を見る場合は次のようにやってます。plot出た!シンプルな図!ploter魂が震えますね。シンプルすぎるので、たまにはこんなのもします。rstan内蔵の関数を使います。T)収束してるのもよくわかってよい。traceplotは普通にrstan内蔵のものをよく使ってます。stan_trace(pars="alpha[1]")潜在変数を取り出します。apply(2,mean)hist(theta)味わい深いですね。HAD11.00をアップしました。HAD12.0をアップしました。Framework.。