アンケートや調査時に回答カテゴリごとの回答数（回答頻度）を集計したい

こんにちは。

お久しぶりです。

いきなりですが、アンケートや質問調査などを集計するときに便利な関数を備忘録的にまとめます。

base内のtabulate関数についてです。

やりたいこと

複数の変数（回答者・項目）の回答カテゴリの回答数（回答頻度）を集計したい

→ table関数で解決です

→ ただ、全変数の回答カテゴリ数が揃っていないと結構面倒（だと思っていました）

こんなときにどうすればよいのか？

はい、やっていきます。

簡単な場合

今回はltmパッケージのScienceデータを使います

回答者 x 項目 の回答データ(今回は4件法)になっています.

library(ltm)
data("Science")

dat <- 
  Science |>  
  lapply(as.numeric) |> 
  as.data.frame()

## Work項目に関する回答を集計したい場合
table(dat$Work)

 #  1   2   3   4 
 # 33  98 206  55 

複数列の場合 1

全ての変数（項目）で全ての回答カテゴリが使われている場合

dat |> 
  apply(2,table) |> 
  as.data.frame() ## stringsAsFactors = TRUE となっているので、注意が必要な場合があります

#   Comfort Environment Work Future Technology Industry Benefit
# 1       5          29   33     14         18       10      21
# 2      32          90   98     72         91       47     100
# 3     266         145  206    210        157      173     193
# 4      89         128   55     96        126      162      78

このように、項目ごとの回答カテゴリ（1~4）への回答頻度を集計したデータフレームを返してくれます。

複数列の場合 2

一部の変数（項目）で一部の回答カテゴリが使われていない場合

まず、ゆがんだデータを作成します.

dat2 <- dat

## わざとデータをゆがませる(Work列で1の回答を全て2に上書き)
dat2$Work <- ifelse(dat2$Work==1,2,dat2$Work) 

table(dat$Work)
 #  1   2   3   4 
 # 33  98 206  55 

table(dat2$Work)
#   2   3   4 
# 131 206  55 

ゆがんだデータdat2ができました。

ここで先ほどの集計方法を使ってみます。

dat2 |> 
  apply(2,table) |> 
  as.data.frame()

 # (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE,  でエラー: 
 #   引数に異なる列数のデータフレームが含まれています: 4, 3 

このように変数（項目）間のカテゴリ数があわなくなるので、データフレームができないと言われます。

私はこの問題に結構前から悩まされていて、いろいろと面倒なことをやって解決していたのですが、

base内のtabulate関数で解決できました。

dat2 |> 
  apply(2,tabulate, nbins=4) ## 項目ごとに頻度を計算してくれます

#      Comfort Environment Work Future Technology Industry Benefit
# [1,]       5          29    0     14         18       10      21
# [2,]      32          90  131     72         91       47     100
# [3,]     266         145  206    210        157      173     193
# [4,]      89         128   55     96        126      162      78

## 回答者(行ごと)の場合も同様に可能(返り値は、カテゴリ数 x 回答者の行列)
dat2 |> 
  apply(1,tabulate, nbins=4)

nbinsは、最大カテゴリ数（今回の例では4件法なので4）を指定します。

この関数によって、回答頻度が0件もきちんと補ってくれて集計することができました。

注意すべき点

ヘルプに書いてあるとおりですが、table関数とは違って、負の数や0は集計しません(以下の3番目の例)。

また、nbinsと最大カテゴリ数が合っていないと（nbins<最大カテゴリ数の場合）、nbins以降は無視されます(以下の4番目の例)。

参照先：ヘルプの例の実行例

tabulate(c(2,3,5))
# > [1] 0 1 1 0 1

tabulate(c(2,3,3,5), nbins = 10)
# > [1] 0 1 2 0 1 0 0 0 0 0

tabulate(c(-2,0,2,3,3,5))  # -2 and 0 are ignored
# > [1] 0 1 2 0 1

tabulate(c(-2,0,2,3,3,5), nbins = 3)
# [1] 0 1 2

以上です。

このような集計を使う人がどれだけいるかわかりませんが。何かの役に立てば。

久々のEnjoy！

Stan Advent Boot Camp 2020 第10日目； おすしでインドスカルしてたら､推しメンが結婚した話｡

ういっす｡この記事は､Stan Advent Calendar 2020 10日目の記事です｡

多次元尺度構成法(Multidimensional scaling; MDS)についてStanで実行してみます｡

初めに書いておくと､良い/面白い結果が出ていないので､やってみたぐらいのノリです｡

qiita.com

~余談~

1ヶ月前のある日､突如某slackでスコココと通知音が鳴った｡

なにかと思えば､来月に控えたStan Advent Calendarの埋まり具合が悪いらしい｡それを危惧した某メンバーが最初の10日間で初心者のためのStan Boot Campをやることを提案した｡僕はその様子を静観していた｡しかし､それは許されなかった｡そして､ここから悲劇は始まる｡気付いたら､Stan Boot Campの担当者になっていたのだ｡これは道を歩いていたら虎に襲われるぐらいの衝撃である｡しかし､驚くべきはこれだけではなかった｡最初は自身の研究領域であるIRTについて執筆する予定だったのだが､｢専門の人が専門領域の記事を書いても面白くない｣などと言った謎の意見により､気付いたらMDSの記事を書くことになっていた｡これはまさに道を歩いていたら気が狂って目の前にいる虎を襲ってしまうぐらいの衝撃である｡そして､今日に至る｡｡｡

~余談おわり~

MDSとは､｢距離から地図を作る方法である(小杉, 2018)｣と説明がされています｡

www.amazon.co.jp

距離と言うと､実際に東京から大阪といった物理的な距離をイメージされることが多いと思いますが､ここでは必ずしも物理的な距離のみを指すわけではありません｡

• たとえば､心理学にはパーソナルスペースという言葉があります｡これは､他人に入ってきて欲しくない(入られると不快になる)距離を指します｡電車に乗っていて､そんなに混んでいないのにもかかわらずこの人めっちゃ近いなみたいな場面を想像して貰えば良いと思います｡このような心理的距離(測定すれば物理距離として定義もできますが)も､距離として扱うことができます｡
• また､類似度は距離として解釈されることがあります｡松屋･すき家･𠮷野家の値段の類似度を出したときに､松屋とすき家間は同価格帯で類似しており距離が近いとすると､松屋と𠮷野家は前者間よりも価格に開きがあり(類似していなく)距離が遠いと定義することができます｡これも物理的な距離ではありませんが距離として扱うことが可能です｡
• そして､対象への選好(好み)も距離として捉えることができます｡りんごとれもんのどちらがどれだけ好きかを考え､れもんが好きだとしましょう｡このとき､好きな者ほど､距離が近いと定義すれば､ぼくとれもんは､ぼくとりんごよりも近いと解釈することができ､これは様々なフルーツで行うことで､僕の好きなフルーツを地図で表すことができます｡このように選好も距離と考えることができます｡
• また､これらの距離は必ずしも対象とは限りません｡例えば､仮に僕がある女優Tさんが好きだったととしましょう｡あくまでも仮の話です｡このとき僕の思うこの気持ち a.k.a 想いは､ 距離です｡しかし､残念なことに女優Tさんは僕の存在すら認識していませんので､女優Tさんから私に対する気持ち a.k.a 想いは距離ですが､未観測なのでNULLです｡つまり､この場合､この私の気持ち a.k.a 恋は､片思いであり､すなわち非対称距離になるのですね｡MDSでは､このような非対称距離を扱うこともできるそうです(この例の場合女優から私への距離はNULLなので厳密に非対称距離とはなっていませんかもしれませんが)｡

以上のように､単純に距離と言っても多種多様です｡

では､今回は､対象への選好を距離と捉えて､その関係性を表す地図を作って行きたいと思います｡

今回用いるデータは､寿司の好みについて評定したデータです｡

こちらをお借りしました｡

www.kamishima.net

これは全国の日本人参加者(N=5000)が各寿司ネタ(全100ネタ)の好みを評定した大変貴重なデータです｡

今回はここから､｢えび｣｢まぐろ｣｢とろ｣｢かっぱ巻き｣｢はまち｣の5ネタを評定した回答者(N=5)に絞り､分析を実施しました｡

※本当は､10ネタぐらいでやりたかったのですが､10ネタすべてを評定した回答者がいなかったので(いずれかのネタが欠損になってしまう)､その兼ね合いで5ネタとしました｡正直この少なさでMDSを使うのは向いていない気もします｡また欠損値を含むMDSができれば良かったのですが､MDSの初心者なのでよくわからなかったです｡検索したら僕のボスの博論が出てきました｡

モデルや分析コードは以下の小杉先生の章を参考にしました｡

www.kitaohji.com

詳細なモデルの数理等については関連書籍を確認してください（まとめる余裕なかった）。そんなに難しくないです。真のネタ間の距離が、個人ごとの選好の違いと掛け合わされて、誤差とともに出てきたものが今回の評定データであるというモデルになっています。

このようなモデルは、INdivisual Differences multidimensional SCALing(INDSCAL)と呼ぶそうです。なお、INDSCALはインドスカルではなく、インスカルと読むらしい。

また､モデルは5ネタでN=5という通常あり得ないサンプルサイズでの解析なので本当にやってみた程度の結果です｡

これらをもとに作成した地図がこれです(点はEAP,エラーバーは50%信用区間)。

うーん、やはりネタ数が少なすぎて解釈できないですね。

いろいろ寿司ネタを変えたり､試したんですが､無理でした｡このデータだと厳しかったかもしれませんね｡

別のデータでリベンジしたいですね｡

締まらない感じ&オーバータイムになってしまいましたが､とりあえず｡

コード等は以下

https://github.com/dastatis/dastatis.github.io/tree/master/pdf/StanAdvent2020

日本視覚学会@2021年冬季大会にてOpen Science Frameworkの解説を行います｡

こんにちは｡

この記事は宣伝です｡

以下のAdvent カレンダーの2日目の記事です｡

adventar.org

オープンサイエンスやっていますか?

やらねばいけないとわかっていながら､結局手を出していないひと｡

やりたいけれども､何から始めれば良いかわからないひと｡

1/20にOSFの使い方について､実際のハンズオンも含め､日本視覚学会冬季大会2021にてトークします｡

visionsociety.wixsite.com

オンラインですので､ぜひとも聞きに来てくださいませ｡

僕の講演の要旨は以下です(後ほど､HP上で公開されると思います)｡

近年､心理学を初めとする社会科学の諸領域において､研究の再現性が問題になっている｡研究の再現性を高めること､これを示すことは決して受動的にできることではない｡これからの時代は､自らの研究が信頼されうる研究であることを積極的に示すことが必要になる｡本セッションでは､再現性の問題について簡単に振り返りつつ､自身の研究がReproducible(再現可能)な研究であること､すなわち信頼されうる研究であることを示す方法について解説し､実践する｡具体的には､Open Science FrameworkというWebサイトを使用し､研究の目的･方法の再現性と結果の再現性を高める方法について解説し､これらを実際に動かし､その簡便性を体感してもらう｡最終的に本セッションの参加者が､この日から再現性を高める研究を実践できることを目指す｡

また､企画セッションの他の講演では､様々なコンテンツを学ぶことができそうです｡

｢オンライン実験･調査をする際の諸々｣や｢計算論的アプローチの最前線研究｣､どちらもとても面白そうですね｡

早期申込は12/7までのようですよ｡

ということで､ 簡単な宣伝で申し訳ないですが､Enjoy!

余力があれば､データ匿名化の話をこのアドベントカレンダーの何日目かに書きます｡

お楽しみに｡

【Tips】cmdstanr使うときの諸々

こんにちは。

自分用にcmdstanrを使うときの諸々のTipsを随時更新していきます。

環境

WindowsからWSLでRstudio Serverからcmdstanrを動かすことを想定。

普段はrstanを使用。

cmdstanrのiterationの指定がrstanと異なる

rstanは、warmupを含む反復回数をiterとして指定

cmdstanrは、warmupを含まない（after warmup）をnum_sumplesとして指定(JAGS等と同一)

cmdstanrの推定結果(stan_csv)の保存先を変えたい

cmdstanrのsample()のオプションでoutput_dirを指定する

WindowsとWSL間でファイルのやりとりがしたい

cmdstanrそのものと関係ありませんが、次の記事を参考にしました。

https://qiita.com/quzq/items/1

qiita.com

096c638c0d86795be13

(Windows → WSL) のやりとりをしたい

想定例：WSLに保存した分析結果にアクセスしたい

以下のパスからアクセス可能。

//wsl$

なので、WSLのRstudio Server内のファイルへのアクセスは例えば以下から。

//wsl$/Ubuntu-18.04/home/rstudio 

(WSL → Windows) のやりとりをしたい

想定例：WSLのRstudio Serverでの分析結果を直接Windowsに保存したい

以下のパス指定で可能(dドライブへのパスの場合)

/mnt/d/

例えば、デスクトップにggplotの図を保存したい。

ggsave(filename = "/mnt/c/Users/dhojo/Desktop/demo.png")

※指定の際のポイントは、ドライブの文字は小文字じゃないとうまくいかない

※先述のoutput_dirと組み合わせれば、推定結果をWindowsに保存可能だと。

デフォルトでwarmup情報は保持されない。

save_warmup = FALSEがデフォルトなので、必要に応じてこのオプションをTRUEとする。

vbの動作方法

コンパイルしたモデルのオブジェクトにvariationalをつける

例

model <- cmdstanr::cmdstan_model("demo.stan")

fit_vb <- model$variational

== 2020/04/29 追加 ==

MCMCサンプルを取り出したい（rstan::extractをしたい）

$draws()で取り出せる。

fit$draws

> fit$draws()
# A draws_array: 250 iterations, 4 chains, and 641 variables
, , variable = lp__

         chain
iteration     1     2     3     4
        1 -3790 -3797 -3837 -3816
        2 -3786 -3786 -3813 -3784
        3 -3794 -3788 -3809 -3781
        4 -3822 -3789 -3803 -3781
        5 -3828 -3811 -3798 -3805

, , variable = theta[1,1]

         chain
iteration      1      2     3    4
        1  0.725  0.448 0.849 1.06
        2  0.261 -0.176 0.625 0.12
...



fit.vb$draws

> fit_vb$draws()
# A draws_matrix: 1000 draws, and 642 variables
    variable
draw  lp__ lp_approx__ theta[1,1] theta[2,1] theta[3,1] theta[4,1]
  1  -3936        -306     -0.101      0.657     -0.095      -1.08      
  2  -3861        -240      0.218     -0.866      0.602       0.55      
  3  -3890        -260      0.390      0.189      0.865      -0.46      
  4  -3896        -268     -0.462     -1.276     -0.112      -0.23      
  5  -3910        -274     -0.179     -1.564      0.194       0.72      
  6  -3906        -281      0.167     -0.087     -0.710      -0.85      
  7  -3909        -271      0.453     -0.735     -0.389      -1.21      

...

==2020/05/01追記== 

bridgesamplingパッケージを使いたい

現状、cmdstanrのオブジェクトそのままでは使えない？と思うので、

以下の要領で使えるようにする。

1．warmup情報を保存しておく

上記に書いたとおり、warmup保存オプションを書いておく。

> save_warmup = TRUE

2．rstan::stan()でモデルを空回しする

1．の準備を行っても、「モデルが妥当ではない」という旨のエラーをはかれるので、

## 事前に

## fit_cmdstan <- rstan::read_stan_csv(fit$output_files())

## をしておく。

fit_rstan <- rstan::stan("hoge.stan",data,iter=0)

bridgesampling::bridge_sampler(fit_cmdstan,fit_rstan)

これで一応回る。

ただ、これだと結局rstanでモデルのコンパイルをすることになります笑

まぁ、モデル評価の時ぐらいはコンパイル時間を待ちましょうかね！

また、随時追加します。

参考：

https://qiita.com/quzq/items/1096c638c0d86795be13

https://mc-stan.org/cmdstanr/reference/cmdstanr-package.html

マニアックなIRTモデルのパラメータリカバリ性能をひたすらチェックするよ(泣)

こんにちは｡

この記事は､Stan Advent Calendarの12日目の記事になります｡

本当は､今年書く予定はなかったのですが､某先生に名指しで｢12日､空いてるよ｣と強い圧をかけられたので書くことになりました笑

最初は何言ってんのかわからなくて､デートにでも誘われているのかと思いました｡

そんな時期です｡

今回の記事は書くことは一昨日決まったので､割とやっつけですがすみません｡

最近､心理学等の社会科学の分野において

(Stanで)ベイズ統計モデリングを行うことが注目されつつあるかと思います｡

データをとれば自由な統計モデルで推定することができるようになりました｡

適当に複数モデル作って､こっちの方がよさそうみたいなことも言えると思います(例えば､情報量規準やBFをもちいて)｡

ただ､これらの統計モデリングをやる上で

まず､大事なことはそのモデルの挙動を知ることかなと思います｡

挙動を知るとは､

モデルのあるパラメータの値を動かすと(変化させると)､

どのようなデータが生成されたり､他のパラメータにどのような影響を与えるかを検討することをここでは指しています(他にもあるかもしれませんが)｡

そして､モデルの挙動を知るために､

統計学の多くのモデル研究でまず検討されることは､

パラメータリカバリ性能の確認かなと思います｡

心理学で(ベイズ)統計モデリングが今後さらに流行っていく上で､

絶対に必要な研究手続きになると思います｡

なので､今回は､IRTモデル(説明変数が潜在変数である特殊なロジスティック回帰モデル)を用いて考えていきます｡

改めて､

パラメータリカバリの性能の確認とは何でしょうか?

我々が一般的な(応用的な)心理学研究を行う際は､

データを測定し､統計モデルにおけるパラメータ推定を行います｡

なので､端的に表すと

実データ → 統計モデル → パラメータ
という流れになると思います｡

ただ､ここで統計モデルがあまりよくない場合に､

統計モデル → パラメータ

がうまくいかないこともあります｡

パラメータリカバリの性能の確認とは､これを確認するシミュレーションになります｡

基本的な流れとしては､

パラメータの真値設定 → 統計モデル → 疑似データ → 統計モデル → 推定されたパラメータ

を行い､最初に設定したパラメータ(真値)と推定されたパラメータの差異等を確認することでそのモデル(や推定法)の推定精度をみることです｡

認知モデルなんかでは､これを行った後に､その真値パラメータを動かしたときのデータの生成挙動や他のパラメータの挙動を確認すべきと思います｡

こういうのも所謂､感度分析に近いものだと思います｡

では､本編に入ります｡
IRTの2値のロジスティックモデルにおいて､
最もよく使われるのが1パラメータもしくは2パラメータかなと思います｡
当て推量を考慮した3パラメータロジスティックモデルも知っている人はいるかと思いますが､
実はもう少しマニアック?なロジスティックモデルとして､4パラメータ､5パラメータモデルも存在します｡
今回は､この二つのモデルのパラメータリカバリを行ってみたいと思います｡

ここらへんのモデルは荘島先生のサイトが端的にまとまっています｡
http://www.rd.dnc.ac.jp/~shojima/exmk/index.htm

一応､整理すると以下のような感じです｡

Parameter	解釈	1PL	2PL	3PL	4PL	5PL
$$\theta_i$$	回答者iの特性 / 能力 ( trait / ability )	〇	〇	〇	〇	〇
$$\beta_j$$	項目jにおける困難度 ( difficulty )	〇	〇	〇	〇	〇
$$\alpha_j$$	項目jにおける識別力 ( discriminaion )		〇	〇	〇	〇
$$\gamma_j$$	項目jにおける当て推量 / 下限			〇	〇	〇
$$\upsilon_j$$	項目jにおける当て推量 / 上限				〇	〇
$$\zeta_j$$	項目jにおける非対称性					〇

モデル

1 Parameter Logistic Model ( Rasch Model )

\begin{align}
P(Y_{ij} | \theta_i, \beta_j)
&= logit^{-1}( \theta_i - \beta_j) \\
&= \frac{1}{1 + \exp( \theta_i - \beta_j)}
\end{align}

項目特性曲線(1PL)

1PL

2 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j)
&= logit^{-1}(a_j \theta_i - \beta_j) \\
&= \frac{1}{1 + \exp(a_j \theta_i - \beta_j)}
\end{aligned}
$$

項目特性曲線(2PL)

2PL

曲線の傾きが項目毎に変わるのが特徴です｡

3 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j, \gamma_j)
&= \gamma_j + ( 1 - \gamma_j ) \times logit^{-1}(a_j \theta_i - \beta_j) \\
&= \gamma_j + ( 1 - \gamma_j ) \times \frac{1}{1 + \exp(a_j \theta_i - \beta_j)} \\
&= \gamma_j + \frac{( 1 - \gamma_j )}{1 + \exp(a_j \theta_i - \beta_j)}
\end{aligned}
$$

項目特性曲線(3PL)

3PL

曲線の下限が項目毎に変わるのが特徴です｡

4 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j, \gamma_j, \upsilon_j)
&= \gamma_j + ( \upsilon_j - \gamma_j ) \times logit^{-1}(a_j \theta_i - \beta_j) \\
&= \gamma_j + ( \upsilon_j - \gamma_j ) \times \frac{1}{1 + \exp(a_j \theta_i - \beta_j)} \\
&= \gamma_j + \frac{( \upsilon_j - \gamma_j )}{1 + \exp(a_j \theta_i - \beta_j)}
\end{aligned}
$$

項目特性曲線(4PL)

4PL

曲線の上限が項目毎に変わるのが特徴です｡

5 Parameter Logistic Model

$$
\begin{aligned}
P(Y_{ij} | \theta_i,\alpha_j,\beta_j, \gamma_j, \upsilon_j, \zeta_j)
&= \gamma_j + \left[( \upsilon_j - \gamma_j ) \times logit^{-1}(a_j \theta_i - \beta_j)\right]^{\zeta_j} \\
&= \gamma_j + \left[( \upsilon_j - \gamma_j )\times \frac{1}{1 + \exp(a_j \theta_i - \beta_j)}\right]^{\zeta_j} \\
&= \gamma_j + \left[\frac{( \upsilon_j - \gamma_j )}{1 + \exp(a_j \theta_i - \beta_j)}\right]^{\zeta_j}
\end{aligned}
$$

項目特性曲線(5PL)

5PL

わかりにくいですが､
曲線の対称性が項目毎に変わるのが特徴です｡
(反応確率が0.5を境に非対称になっている)

意外とありましたね､IRT2値ロジスティックモデル｡
非対称\zeta_jを混ぜた2パラメータロジスティックモデルとか使ったら楽しそう(ここではやらない)｡

では､ここから､5PL(パラメータロジスティックモデル)でパラメータリカバリの性能を確認します｡
上の各モデルの項目特性曲線含めて､今回は以下のパラメータの真値からデータを生成しました｡
詳しいコードは末尾につけるgithubリンクを見てください｡

パラメータの真値

I <- 1000 ##  回答者数
J <- 20 ## 項目数

alpha <- rlnorm(J, meanlog = 0, sdlog = 1) ## discrimination
beta <- stats::rnorm(J, mean = 0, sd = 2) ## difficulty
gamma <- stats::runif(J, min = 0, max = 0.2) ## guessing 
upsilon <- stats::runif(J, min = 0.8, max = 1) ## 
zeta <- extraDistr::rhnorm(J, sigma = 1) ## asymmetric
theta <- rnorm(I, mean = 0, sd = 1) ## trait or ability

統計モデル(Stan)

data{
  int<lower = 0> I;
  int<lower = 0> J;
  int<lower = 0, upper = 1> y[I,J];
  int<lower = 0, upper = 1> mode;
}

transformed data {
  real<lower = 0, upper = 1> watanabe_beta;
  if ( mode == 0 ) {
    watanabe_beta = 1.0;
  }
    else { // mode == 1
    watanabe_beta = 1.0/log(I);
  }
}

parameters {
  vector<lower = 0>[J] alpha;
  vector[J] beta;
  vector<lower = 0.0, upper = 0.2>[J] gamma;
  vector<lower = 0.8, upper = 1.0>[J] upsilon;
  vector<lower = 0>[J] zeta;
  vector[I] theta;
}

transformed parameters {
  real<lower = 0, upper = 1> prob[I,J];
  for (i in 1:I) {
    for (j in 1:J) {
      prob[i,j] = gamma[j] + (upsilon[j]-gamma[j]) * inv_logit(alpha[j] * (theta[i] - beta[j]))^zeta[j];
    }
  }
}

model{
  target += lognormal_lpdf(alpha | 0 , 1) - lognormal_lccdf(0 | 0, 1);
  target += normal_lpdf(beta | 0 , 5);
  target += uniform_lpdf(gamma | 0.0 , 0.2) - uniform_lccdf(0 | 0.0 , 0.2) -  uniform_lcdf(0.2 | 0.0 , 0.2);
  target += uniform_lpdf(upsilon | 0.8 , 1.0) - uniform_lccdf(0.8 | 0.8 , 1.0) -  uniform_lcdf(1.0 | 0.8 , 1.0);
  target += normal_lpdf(zeta | 0 , 1) - normal_lccdf(0 | 0, 1);
  target += normal_lpdf(theta | 0 , 1);
  for (i in 1:I) {
    for (j in 1:J) {
      target += watanabe_beta * bernoulli_lpmf(y[i,j] | prob[i,j]);
    }
  }
}

generated quantities {
  int<lower = 0, upper = 1> y_pred[I,J];
  vector[I] log_lik;
  vector[J] log_lik_temp;
  for (i in 1:I) {
    for (j in 1:J) {
      log_lik_temp[j] = bernoulli_lpmf(y[i,j] | prob[i,j]);
      y_pred[i,j] = bernoulli_rng(prob[i,j]);
    }
    log_lik[i] = log_sum_exp(log_lik_temp);
  }
}

細かいところを見てくださったかたは､お気づきかと思いますが､
事前分布の設定と真値の設定を微妙に変えています(項目特性曲線のときとも少し値変わってます)｡
本来は､同様の分布が望ましいと思いますが､今回はデモなのでこうしています｡

はい､Stanでベイズ推定しました(雑っ!!!)｡

5PLのalphaの真値と推定値の相関

はい､あんまりうまく推定できてないね(推定出来ているはずだったw)｡
二値データからパラメータが5つあるモデルはやはり難しいね｡
ほら､パラメータリカバリの確認は大事でしょ?w

5PLa

ちなみに､同真値で1PLモデルのβを見てみると下のようになります｡

1PLb

美しい!!!
こんなに綺麗に相関していれば､ちゃんと実践でも使えそうですね｡

最後､ほんまにやっつけですみません｡｡｡
忙しくてw

詳細コードは以下に置きました｡

github.com

Enjoy!!

【ポエム】どうして僕がOpen Scienceをやるのか｡

こんにちは｡

この記事は､｢Open and Reproducible Science Advent Calendar 2019｣の3日目の記事です｡そして､ポエムです｡

気分を害したくないなら､読まない方が身のためです｡

adventar.org

じゃあ､早速ポエムを始めましょうw

僕がOpenScienceを行う理由｡

それはとても簡単です｡僕は統計ができないからです｡

だから､私は常にOpenScienceに大きな関心と注意を注いでいます｡

現在､僕は心理統計学を研究しています｡学部3年のゼミから"専攻"していますので､だいたい､学習歴は6年ぐらい?ですかね｡

6年といえば､小学校に在籍する期間と一緒です｡

ただ､それだけ学習しても(心理)統計はさっぱりわかりません｡

僕の人一倍の要領の悪さといい加減な学習方法が､この原因だと思っています｡

ただ､さすがに6年も専門的にやっているので､

心理学研究を行っている母集団の中では､ほんとにほんとにほんの少しだけ統計がわかるかもしれません｡

だけど､やはり統計学を専攻している､所謂､理系の皆様よりはるかに統計ができません｡

唯一の救いはわからないことを知っていることです｡

だから､僕はOpenScienceを行うことを心がけています｡

その理由､ちゃんと説明していきます｡

ここからの話は､特に心理学研究を想定しています｡

（量的/定量的/データを用いた）研究を行う上で､統計学及び､それに関連する知識は必須です｡絶対にこれから逃れることが出来ません｡

これは､理論･応用研究を問わずにそうです｡

上でいう､統計学に関連する知識とは､

たとえば､

• 機械学習の知識
• 因果推論の知識

などが該当するのかなと思っています｡

場面にもよりますが､

僕は､(心理学に問わず)社会科学を研究する上では､これらの知識は絶対に必須だと考えており､完璧に漏らすことなく知っている必要があると考えています｡

あたりまえですが､そんなの無理です｡全知全能の神でも無い限り｡

少し大げさに書きました｡*1

僕には､これが不可能だからOpenScienceを行います｡

OpenScienceは､統計学のできない僕みたいなのこそ､

真っ先に身につけ､実践すべき研究法です｡

ここで改めて､OpenScienceとは何かを考えます｡

研究全体もしくはその過程を可視化することがOpenScienceです(正確な定義ではないかもしれませんが)｡

僕がOpenScienceに関して､真っ先に思いつくのが､研究で収集､分析したデータの公開です｡みなさんも同じものを想像したかもしれません｡

ただ､これだけがOpenScienceではありません｡

• 調査･実験に至るまで
• 方法(データ収集､手順の詳細)
• 解析手法(論文に載せていないものも含む)
• 等々

上記を正確に記述し､公開することもOpenScienceの一つだと思います｡

これに実施することによって､科学は発展します(と考えています)｡

少なくとも､科学の流れを止めません｡止めるのがどうかみたいな話はしません｡

たとえば､

僕は､OpenScienceの作法に則って､データを集め､解析Aを行い､論文を書きました｡

しかし､統計学マスターの謙介さん(仮名)*2が､論文を読んだところ､今回の実験に関して､(解析Aでも間違っていないが)解析Bがより適していることに気付きました｡

※この状況では､実際､解析Bを行うことがより正しい､わかることが多いとします｡

ここで､僕がOpenScienceに取り組んでいなければ､ 科学の発展は終わりです｡少なくとも停滞します｡

• 僕がOpenScienceに取り組んでいることによって､

1. 謙介さん(仮名)*3が､同様の手順でデータを取り直すことができ､解析Bを行うことができる
2. もしデータが公開されているならば､データを取り直すことなく解析Bをその場で実行することができる

というように､僕よりも遙かに優秀な人が効率的に科学を発展させる(社会をより良くする)ことにつながるのです｡

もちろん､ここで僕もOpenScienceの作法に則り､貴重な研究基礎を提供したので､このようなコミュニティでは評価されることにつながると思います｡

一方で､

• 僕がOpenScienceに取り組んでいなかった場合､

1. 謙介さん(仮名)*4が､僕の論文同様の手順でデータを取ることが出来ないかもしれない*5｡これによって､解析Bが出来たとしても､その結果と元論文の分析Aの結果を比較可能なのか?という問題も出てくる可能性がある｡

というように､僕がOpenScienceを取り組んでいないがために､先ほどよりも研究の発展が難しくなるかもしれません｡ここで､私は謙介さん(仮名)の先行研究になり得ていますが､謙介さん(仮名)がOpenScienceの作法に則って研究を行い､解析Aと解析Bを行った場合､そちらの研究の方がより研究貢献をしているような気がします*6｡

従って､

僕がOpenScienceの作法に則り､研究を行えば､今後の誰かの研究の大きな一歩につながるかもしれないということです｡そして､自分の研究がその研究に大きく貢献する可能性があるということです｡

少し話しが逸れますが､

心理学の分野によってはOpenScienceが行いにくい分野ももちろんあります｡

これには様々な理由が考えられますが､

1. プライバシーの観点からデータ･手順を公開しにくい
2. 古くからの慣習
3. 大御所のおっさんが認めてくれない(多分あとで怒られる)

1.は臨床の分野だと多いかもしれませんね｡そして､それを理由にOpenScienceを避けているのをたまに見たような気がしますが､あれは夢だったかもしれません｡

まだ､心理学領域でほぼ語られていませんが､データを匿名化する方法はたくさん存在し､これらの方法を用いれば､たとえ臨床心理学で使われるような秘匿性の高いデータでも公開できるかもしれません｡

また､医療系では既に匿名化を行うことでカルテ等のデータを共有することを考えているという話も聞いています｡

おそらく､2-3年の内に(臨床)心理学でも匿名化の話が出てくると思います｡*7

ここまで書いてきたことを踏まえたうえで､ 

僕は､統計はできないけど(から)､OpenScienceは行います｡

*1:千里の道も一歩から､科学の発展は一日にしてならずという教えをボスから教わり､僕もそれを正しいと思っているので､始めから完璧な研究をやる必要は無いと思いますが､あくまで大げさに言っています｡

*2:深い意味は全くありませんw

*3:何度も書きますが､深い意味はありません

*4:何度でも書きますが､深い意味はありません

*5:なんか出来てしまいそうな気がしますがw

*6:あくまで1例で､一つの考え方です｡これによって､元論文の価値がなくなるとは思いませんが｡

*7:今が狙い目です｡

CmdStanRのコンパイルが速すぎたお話

こんにちは。

今朝、TwitterのStan好きTLを騒がせていた、

Holy shitballs CmdStanR is so fucking fast! https://t.co/stM9JnqoE0

— Dan Simpson (@dan_p_simpson) October 31, 2019

Help to test development version of CmdStanR, which is a lightweight interface to Stan for R usershttps://t.co/G3SOrxKF4O. CmdStanR misses some useful features of RStan, but has the benefit of getting more recent features faster. One example is the much faster model compilation. https://t.co/gOWVWFlm2Q

— Stan (@mcmc_stan) October 31, 2019

こちらについて簡単に検証してみました。

まだ、開発が始まったばかり程度なので、今後どうなるかわかりませんが、とても期待できそうであることがわかりました。

結果から書くと、

mc-stan.org

に記載のモデル(bernoulli.stan)では、

CmdStanRによるコンパイル時間
Time difference of 1.574471 secs
rstanによるコンパイル時間
Time difference of 44.28613 secs

という結果でした。本当にso fucking fast! で笑いました。

ということで、もう少し詳しく見ていこうと思います。

基本的には、上のGetting startedを見て、導入、実践してもらえればと思います。

まず、私の環境ですが、OSが Windows 10 Proです｡

そして､このおかげ?で､Windowsにうまくcmdstanを導入することが出来ませんでした｡

というのも､

devtools::install_github("stan-dev/cmdstanr")

library(cmdstanr)

 ここまではうまくいくのですが､

 install_cmdstan()

 これがうまくいきませんでした｡ 

Running bash "C:/Users/hogehoge/Documents/R/win-library/3.6/cmdstanr/make_cmdstan.sh" \ -j2 -w
Error in rethrow_call(c_processx_exec, command, c(command, args), stdin, :
Command not found @win/processx.c:977

 という感じのエラーがはかれます｡

恐らく､実行に際し､make_cmdstan.shというファイルを実行して､cmdstanをインストールするようなのですが､そこでこけてしまっています｡

.shファイルは､bashと呼ばれる､linuxとかで開く前提の拡張子らしいので?(あっている?)そもそも､windowsで開けないっぽい｡

windowsで.shファイルを開く方法をググったらできそうな感じはありましたが､それを実行するのが面倒だったので window subsystem for linux (ubuntu LTS 18.04) を使って環境を構築しました｡なので､環境は､仮想linux(ubuntu)です｡

もちろん､通常?のlinuxやAWS,AWS,GCP等のlinux環境でも普通に実行できると思いますが｡

また､Rからではなく､自力でcmdstanを導入することもできると思いますが､今回は行いませんでした(面倒だった)｡

github.com

まずは､仮想linuxにRやRstudio Serverを構築しました(ここはCmdStanR関係ないので省きます)｡ 

調べればたくさん出てきますので､そちらを参考にしてください｡

そしたら､Rで(僕は､Rstudio Serverでやりました｡)

devtools::install_github("stan-dev/cmdstanr")
library(cmdstanr)

install_cmdstan()

先ほど同様のコードを走らせます｡

すると､cmdstanがインストールされます｡

その後､cmdstanのpathを指定します(動作環境によって以下でうまくいかないようですが､その際は公式( https://mc-stan.org/cmdstanr/articles/cmdstanr.html )を見てください)｡

set_cmdstan_path()

続いて､.stanモデルを指定します｡

今回は､パッケージインストールと同時に導入される例題モデルをコンパイルしてみます｡

## .stanファイルのpathを指定する｡

file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")

start <- Sys.time()
mod <- cmdstan_model(file) ##モデルのコンパイル
end <- Sys.time()
print(end-start)

上記を実行させれば､モデルがコンパイルされます｡

僕の環境では､

Time difference of 1.574471 secs

でした｡

ちなみに､コンパイルされた例題モデルは

mod$print()

で確認することが出来ます｡

上記では､

data {
int<lower=0> N;
int<lower=0,upper=1> y[N];
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1,1); // uniform prior on interval 0,1
y ~ bernoulli(theta);
}

が､出力されると思います｡

あとは､

data_list <- list(N = 10, y =c(0,1,0,0,0,0,0,0,0,1))
fit <- mod$sample(data = data_list, seed = 123, num_chains = 2)

を実行すれば､MCMCのサンプリングを行ってくれます｡

VBも出来るそうですが､今回は実行しませんでした(内容の本質ではありませんので)｡

結果オブジェクトは､rstanの結果オブジェクトとは異なります｡

class(fit)

[1] "CmdStanMCMC" "R6" 

 という形式で保存されています｡

そのため､結果の取り出し方も少し違います｡

## 要約統計量

fit$summary()

## MCMCの診断
fit$diagnose()

## MCMC samples (rstanのextract関数)

fit$draws()

 ここらへん､少し不便かなと思ったのですが､以下のように実行してあげると､rstanの結果オブジェクトと同様に扱ってあげることができます｡

stanfit <- rstan::read_stan_csv(fit$output_files())
print(stanfit)

やっていることとしては､

cmdstanrの結果オブジェクトは､.csvファイルで保存されている(される?)ようなので､

これらをrstanに標準装備されているread_stan_csv()関数で読み込んであげれば､結果オブジェクトを復元することができます｡

あとは､いつものようにrstanの結果オブジェクトと同様に扱うことが可能です｡

ここらへんは､今後改善されるてきなことが公式にも記載されていたような気がします｡

それでは､少し話しを戻して､cmdstanrでコンパイル時間を計算できたので､比較してみましょう｡

先ほどのモデルコードを事前にbernoulli.stanのようにどこかに保存しておいてください｡

比較するので､念のため､僕はRを再起動しました｡

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

start <- Sys.time()
mod <- stan_model("Documents/R/model/bernoulli.stan")
end <- Sys.time()
print(end-start)

あとは､こんな感じで書いてあげれば､コンパイル時間を計測することが出来るかと思います｡

Time difference of 44.28613 secs

という結果でした｡

比較してみると､

cmdstanRによるコンパイル時間
Time difference of 1.574471 secs
rstanによるコンパイル時間
Time difference of 44.28613 secs

so fucking fast! でした(言いたいだけ)｡

経験上､モデルの違いでコンパイル時間は大きく変化しない気がするので(主にCPUのスペックだと思います)､他のモデルでも早いのではないでしょうか｡

ちなみに､なぜこんなに早いのかについても公式で言及されています｡

The RStan interface (rstan package) is an in-memory interface to Stan and relies on R packages like Rcpp and inline call C++ code from R. On the other hand, the CmdStanR interface does not directly call any C++ code from R, instead relying on CmdStan for compilation, running algorithms, and writing results to output files.

https://mc-stan.org/cmdstanr/articles/cmdstanr.html

Rcppを一切介さず諸々をやるので速いようです｡

実際に計算しているのも､CmdStanくんが頑張っているようです｡

速いっていいですね｡

今回は､サクッと比較したので､記事に誤字､脱字多いかもしれません｡すみません｡

あと､Windowsで簡単に出来る方法(CmdStanをインストールさえしてしまえばという感じですが)をご存じの方は教えていただけると嬉しいです｡

Enjoy!!