PsyStat.

Just for fun.

マニアックな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)

f:id:dastatis:20191212015749p:plain
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)

f:id:dastatis:20191212020032p:plain
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)

f:id:dastatis:20191212020111p:plain
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)

f:id:dastatis:20191212020322p:plain
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)

f:id:dastatis:20191212020454p:plain
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

f:id:dastatis:20191212180552p:plain
5PLa

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

f:id:dastatis:20191212181052p:plain
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を騒がせていた、

 

 

 

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

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

 

結果から書くと、

 

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)です。

もちろん、通常?のlinuxAWS,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!!

 

 

 





事後最頻値をとる関数(get_posterior_map)を作りました。

どうも。

この記事は、StanAdventも何も関係ありませんw

 

Stanで事後平均値(EAP)をとってくるときにget_posterior_mean()があると思いますが、事後分布の最頻値(MAP)をとってきたい場合もたまにあるかと思います。

 

僕はその都度R codeを書いていたのですが、めんどくさくなったので、MAP推定値をとってくる関数を書きました。
それ以上でもそれ以下でもない記事ですw



依存パッケージは、dplyrだけです(%>%が使いたいだけ)。

構造:get_posterior_map(stanfit, pars)
stanfit : Stan結果オブジェクト
pars : 取りたいパラメータ

注意:抽出したいパラメータの次元が4次元まで対応しています。

使用したい場合、以下のコードをコピーして実行してください。

get_posterior_map <- function(stanfit,pars) {
  if(!require("dplyr",character.only = T)) install.packages("dplyr", dependencies = T); require("dplyr",character.only = T)
  maxpost <- function(par){ density(par)$x[which.max(density(par)$y)]}
  r <- rstan::extract(stanfit,pars=pars)
  Dim <- r[[pars]] %>% dim %>% length - 1
  if (Dim == 1 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2]))
    for (l in 1:dim(r[[pars]])[2])
      map_est[l] <- r[[pars]][,l] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 2 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        map_est[l,m] <- r[[pars]][,l,m] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 3 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3],dim(r[[pars]])[4]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        for (n in 1:dim(r[[pars]])[4]) 
          map_est[l,m,n] <- r[[pars]][,l,m,n] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 4){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3],dim(r[[pars]])[4],dim(r[[pars]])[5]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        for (n in 1:dim(r[[pars]])[4]) 
          for (o in 1:dim(r[[pars]])[5]) 
            map_est[l,m,n,o] <- r[[pars]][,l,m,n,o] %>% maxpost(.)
    return(map_est)
  }else{
    print("Check Parameter Dimensions")
  }
}


って感じです。
もっとうまく書く方法はあると思いますが、僕はこれでうまく取ってこれるので問題無いです:)

Enjoy!

OSF (Open Science Framework) で、オープンな心理学研究を目指して

はじめに 

 

こんにちは

 

現在、The 51st Annual Meeting of the Society for Mathematical Psychology

(通称: MathPsych/ICCM 2018 )に参加しています。

 

この学会では、Bayesian Cognitive Modelingに関する研究発表が盛んに行われています。

一方で、この学会に参加されている心理学者の多くは、ベイズ統計モデリングはもちろん、心理学再現性の問題に対しても関心が高く、方法論ないしはそれよりも広い視点から心理学現象を理解するために研究を行っています。

 

その中でも、OSFを用いた研究の"可視化"はマストになりつつあります。

現に、このMathPsychで発表されるスライドやポスターを1カ所で共有し、オープンにする取り組みも行われています。 

 

ここでは、論文執筆(もしくは研究発表)を行うときにOSFをどのように用いるかを簡単にまとめてみます。

(そもそも、この記事を書いた理由は、最近OSFをいじっていた際に査読用の匿名機能があることを知ったので、自分への備忘録を兼ねて書きます。)

 

まず、OSFをご存じない方は、こちらのスライドをご覧ください。

www.slideshare.net

 

論文をsubmitするまで

研究の一連の流れとして、たとえば以下が考えられます。

  1. 実験(調査)を計画し、
  2. 分析を行い、
  3. 論文を執筆する。

我々が、上記の過程で研究の質を高めるために行うべきこととして、プレレジ*1などが挙げられますが、ここでは話題に挙げません。

 

今回は特に、2.の分析で出た、結果の補足もしくは分析コード、生データをOSFで公開するための方法を紹介いたします。

 

まず、OSFにアカウントの登録します(これは説明を省きます)。

 

そして、ログインし、研究プロジェクト(もしくは論文ごと)のページを作成します。

"Create new project" を押します。

f:id:dastatis:20180724181439p:plain

 

ここでは、試しに「論文A」というプロジェクトを作成します。

f:id:dastatis:20180724181606p:plain

 

できました!

f:id:dastatis:20180724181916p:plain

あとは、Filesに追加したい分析コードや生データを追加するだけです。

必要な記述・ライセンス・タグなどは必要に応じて記述してください。

 

このプロジェクトページを公開する場合、右上の"Make Public"を押すことで公開することができます。

一度公開したものを非公開にすることも可能ですが、もちろん非公開にするまでにダウンロードされてしまうことはありますので注意が必要です。

 

 論文に記載すべきもの

では、ここからは、論文に記載する際にどうすれば良いのかについて説明します。

 

結論から書きます。

先ほどのプロジェクトページのURLを論文の本文もしくは脚注に記載します。

今回で「論文A」プロジェクトでいえば、https://osf.io/bz79v/ です。

 

OSFでは、自動的に短縮URLとなるため、論文に記載しても場所をあまり取りません。 以上で、完了です。驚くほど、簡単です!

公開する際は、"Make Public"にするのを忘れないでくださいね!

 

ただ、これで一件落着かと思ったら、そうではありません。

論文を執筆し、submitしたら査読されます。

このとき、論文は著者情報を匿名化し審査されることが一般的です。

 

しかしながら、現段階では、プロジェクトの Contributorsに自分の名前が記載されたままで名前がばれてしまいます。

このような場合どうすれば良いでしょう?

 

プロジェクトの匿名化のために

A. 査読の際は、DropBox、GoogleDriveといったファイル共有サービスに一時的にアップし、匿名リンクを論文に貼ればよい。

→これももちろん間違いではありません。現に私はこのような方法を用いていました。しかし、アクセプト後にOSFにファイルを移して、色々やるのは手間ですし、何よりここでミスが起きてしまうと元も子もありません。

 

しかし、OSFはそれすらも解決してくれます。

なので、推奨する方法はこれです。

A. OSFの査読用機能で匿名化する

 

まず、先ほどの「論文A」プロジェクトページ上部の"Settings"に行きます。

f:id:dastatis:20180724185246p:plain

 

そして、"View-only Links"というタブの"+Add"を選択します。

f:id:dastatis:20180724185356p:plain

 

 このような画面がでますので、「Anonymize...」のチェックボックスにチェックを入れて、Createします。必要があればLink nameも設定できます。

f:id:dastatis:20180724185536p:plain

 

Createすると先ほどの画面にこのようなリンクが作成されます。

f:id:dastatis:20180724185805p:plain

これのShared project linkに書いてあるリンクをコピーして論文に貼り付ければ、匿名化されたプロジェクトのページを査読者が確認することができます。

ここで作られる匿名用のリンクは元々のリンクであったhttps://osf.io/bz79v/ が少し長くなるだけなので、アクセプトされた後も直すのか簡単です。

また、査読の際には、事前に一言、「元々のリンクは著者情報を含むため、一時的にこちらのリンクにコード(付録)を置いた」的なことを書くとスムーズかと思います。

 

実際に、匿名化されたプロジェクトページはこんな感じになります。

f:id:dastatis:20180724190054p:plain

先ほどまでContributorsに僕の名前が入っていましたが、これがAnonymous Contributorsに変わります(かっこいい)。

 

さて、これで自分の研究がだいぶオープンになりました!

OSFは本当に便利ですね!論文執筆の際の参考に!

 

 

またどこかで!Enjoy!

 

 

*1:pre-registrationの略。事前に実験(調査)計画を登録しておくことで、データ収集後の仮説の変更や不都合な結果がでた場合に公表されないパブリケーションバイアスを防ぐとされている

Tokyo.R #70 「Bayesian Sushistical Modelling」を発表

タイトルのイベントで話しました。

他の方の資料は以下よりご覧になることができます。

第70回R勉強会@東京(#TokyoR) : ATND

 

 私は、「Bayesian Sushistical Modelling」というタイトルで話しました。

 

www.slideshare.net

 

ベイズモデリングの話をお願いします。」と運営の方に頼まれたわけですが、

内容を見てもらってわかるようにほぼ記述統計とIRT(項目反応理論)の基礎のお話です。

 

実は裏テーマがあって、それは「寿司データで学ぶ項目反応モデル」です。

マーケティングとかでは、よく調査が行われると思うのですが、

実態は単純集計のみの公表で終わることが多々あると思います(偏見です)。

(集計して、可視化することは非常に重要なので、それ自体が悪いことでまったくありません)

 

ただ、このスライドに出てくるように、5件法の調査項目が存在するときに、

僕の「4.すきです」とあなたの「4.すきです」は同じように解釈して良いでしょうか?

 

極端な例ですが、このように素点(そのままの点数)で結果を解釈することは

僕の「4.すきです」と宇宙人の「4.すきです」を同じように解釈していることになります。

また、猫ちゃんがたまたまボタンを押した「4.すきです」も同じように解釈することになります。(極論ですけどね。。。そして、猫ちゃんはかわいいからゆるす。)

 

昨日は、そこから一歩踏み出して、Let's ベイズ項目反応モデル!という感じでした。

もう1歩や2歩飛び出すとそこにはベイズモデリングの世界が広がっているといった印象です(明確な定義はありませんが)。 

 

私自身もまだまだ勉強中ですが、

IRTモデル+ベイズモデリングのエッセンスで

社会の様々な問題が解決(真正な人材採用や人材評価が)できるかもしれません。

 

なので、これを機に、項目反応モデルに興味を持って頂ける方が一人でも増えて、そこにベイズモデリングの要素を組み込んだ、"より良い社会にむけたベイズモデリング"が広がると良いなと思っております。

 

最後に、

昨日、参加された方、鑑賞された方、発表された方、運営された方、

みなさまありがとうございました🍣

 

またどこかで! Enjoy!