こんにちは。
この記事は、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)
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)
曲線の傾きが項目毎に変わるのが特徴です。
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)
曲線の下限が項目毎に変わるのが特徴です。
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)
曲線の上限が項目毎に変わるのが特徴です。
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)
わかりにくいですが、
曲線の対称性が項目毎に変わるのが特徴です。
(反応確率が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
ちなみに、同真値で1PLモデルのβを見てみると下のようになります。
美しい!!!
こんなに綺麗に相関していれば、ちゃんと実践でも使えそうですね。
最後、ほんまにやっつけですみません。。。
忙しくてw
詳細コードは以下に置きました。
Enjoy!!