項目反応理論やってみよう

さて,Stan advent calender 2020 7日目です。本日は項目反応理論ですね。

1 項目反応理論て?

項目反応理論自体の説明は,加藤・山田・川端 (2014)が非常に丁寧にしてくれています。書籍なので購入する必要がありますが,お薦めです。また,項目反応理論やitem response theoryでググるとタダで読めるいろんな資料もヒットしますのでそれらを読むとよいでしょう。

特にStanで項目反応理論を実施する様々な例を清水さんがご紹介下さっているので、そちらも参照ください。

項目反応理論をStanで実行するときのあれこれ

展開型項目反応理論をStanで推定する

2 いろいろな項目反応理論

この記事では,Leo & jiao (2017) に沿って,項目反応理論の中でもちょっと複雑な,かゆい所に手が届くかもしれない以下のモデルについて紹介します。

  • 3パラメタ(以下,PL)モデル

  • 段階反応モデル(graded response model: GRM; Samejima, 1969)

    • IRTは通常2値の回答になるものに対して実施します。段階反応モデルではリッカート尺度のような多値の回答にIRTを適用したモデルです。段階反応モデルは,あるカテゴリーからあるカテゴリーに回答する差分を推定していくことになります。
  • 名義反応モデル(nominal response model: NRM; Bock, 1972)

    • 名義反応モデルも段階反応モデル同様に多値の回答にIRTを適用したモデルです。各カテゴリーにおける困難度などを推定するモデルです。
  • ラッシュテストレットモデル(Rasch testlet model, Wang & Wilson, 2005)

    • これはラッシュモデル(識別力は1に,当て推量固定を0に固定し,各項目の困難度を推定)を多次元に拡張したものになります。例えば,学力を測定する問題が50問あったとしたときに,その50問を10問ずつの5セット(こうしたものをテストレットとよぶ。訳語わかんない)に分けて別の人たちに実施したような場合に適用可能です。
  • マルチレベル3PLモデル

    • 3PLモデルのマルチレベル版です。例えば,ある県の児童の学力を調査した場合には,その児童の学力は学級によっても違うかもしれないし,学校によっても違うかもしれません。そうしたときに児童は,学級にネストされ,さらに学校にもネストされていると言ったりします。こうした状態をそのまま分析する方法がマルチレベル分析で,IRTにマルチレベル分析を適用したものです。

    • 行列式読めるよって人は,マルチレベルIRTについてはFox & Glas (2001) を読むとよいかもしれない。

よくあるIRTの説明では,1PLモデルやラッシュモデルから入ることが多いと思います。しかし,この記事では,いきなり3PLモデルから説明していきます。なぜあえて複雑な3PLモデルから説明するのかというと,論文の構成がそのようになっているからラッシュモデル,1PLモデルや2PLモデルは,3PLモデルに制約を加えたモデルとして理解できるからです。

この記事では、主に3PLモデルを中心に、モデルの解説とStanコードと結果の抽出方法について紹介します。

2.0.1 IRTの3PLモデルとその他のモデルの関係性

あとでそれぞれが何なのかはもう少しだけ詳しく書きますが,3PLモデルでは

  • 困難度
  • 識別力
  • 当て推量

の3つのパラメタを推定します。このうち,ラッシュモデルと1PLモデルでは困難度を,2PLモデルでは困難度と識別力を推定します。3PLモデルを使って2PLモデルに一致させるには,2PLモデルで推定されない当て推量について,そのパラメタの値を0にしてしまえば,2PLモデルと一致することになります。

また,3PLモデルで,1PLモデルと一致させるには,当て推量の値を0にし,識別力を項目間で等値制約を置けばOKです。明示的に1に固定すればラッシュモデルとなるとLeo & jiao (2017)の中では書かれています12

3 3PLモデル

3PLモデルは,以下のように表現できる。

\[ p_{ij} (u_{ij} = 1| \theta_i, a_j, b_j, c_j) = c_j + \frac{1-c_j}{1 + {\rm exp}(-a_j(\theta_i - b_j))} \]

  • \(u_{ij}\): \(項目_j\)に対する\(回答者_i\)の反応

  • \(\theta_i\): \(回答者_i\)の能力(潜在変数)

  • \(b_j\): \(項目_j\)の困難度PL。この値が変わると能力値\(\theta\)によって以下の図(このような曲線を項目特性曲線(item characteristic curve: ICC)と呼ぶ)のような感じで正答率が平均的に変わっていきます。同じ能力値\(\theta\)でも変わってきます。下の図は,3PLモデルの式に,\(a\)を1に,\(c\)の値を0に固定して,\(b\)の値のみを変えています。

curve((1 + 0)/(1 + exp(-1.0*(x-0.0))), -3, 3, col = "black", lty = 1, ylab = "probability") # a = 1.0, b = 0.0, c = 0.0
curve((1 + 0)/(1 + exp(-1.0*(x-1.0))), -3, 3, col = "red", lty = 2, add = TRUE) # a = 1.0, b = 1.0, c = 0.0
curve((1 + 0)/(1 + exp(-1.0*(x-(-1.0)))), -3, 3, col = "blue", lty = 3, add = TRUE) # a = 1.0, b = -1.0, c = 0.0
legend("topleft", 
       legend = c("a = 1.0, b = 0.0, c = 0.0", "a = 1.0, b = 1.0, c = 0.0", "a = 1.0, b = -1.0, c = 0.0"), 
       col = c("black", "red", "blue"),
       lty = 1:4, cex = 1.0)

  • \(a_j\): \(項目_j\)の識別力。この値が変わると以下の図のような感じで,S字の曲線の傾きが変わります。つまり,識別力が変わると(テストであれば)正答率が同じ能力値\(\theta\)でも変わってくることを意味します(ただし正答率50%となる\(\theta\)の値は変わらない)。下の図は,3PLモデルの式に,\(b\)\(c\)の値を0に固定して,\(a\)の値のみを変えています。
curve((1 + 0)/(1 + exp(-1.0*(x-0.0))), -3, 3, col = "black", lty = 1, ylab = "probability") # a = 1.0, b = 0, c = 0
curve((1 + 0)/(1 + exp(-0.5*(x-0.0))), -3, 3, col = "red", lty = 2, add = TRUE) # a = 0.5, b = 0, c = 0
curve((1 + 0)/(1 + exp(-2.0*(x-0.0))), -3, 3, col = "blue", lty = 3, add = TRUE) # a = 2.0, b = 0, c = 0
legend("topleft", 
       legend = c("a = 1.0, b = 0.0, c = 0.0", "a = 0.5, b = 0.0, c = 0.0", "a = 2.0, b = 0.0, c = 0.0"), 
       col = c("black", "red", "blue"),
       lty = 1:4, cex = 1.0)

  • \(c_j\): \(項目_j\)の当て推量。この値が変わると以下の図のような感じで,S字の曲線の下限値が変わってきます。つまり,当て推量の値の大きさは,たまたま正答してしまう可能性や誰もがある反応をしてしまう可能性をモデルに組み込んだものといえます。下の図は,3PLモデルの式に,\(a\)を1,\(b\)の値を0に固定して,\(c\)の値のみを変えています。
curve(0.0 + ((1 - 0.0)/(1 + exp(-1.0*(x-0.0)))), -3, 3, col = "black", lty = 1, ylab = "probability") # a = 1.0, b = 0, c = 0
curve(0.5 + ((1 - 0.5)/(1 + exp(-1.0*(x-0.0)))), -3, 3, col = "red", lty = 2, add = TRUE) # a = 1.0, b = 0, c = 0.5
curve(0.1 + ((1 - 0.1)/(1 + exp(-1.0*(x-0.0)))), -3, 3, col = "blue", lty = 3, add = TRUE) # a = 1.0, b = 0, c = 0.1
legend("bottomright", 
       legend = c("a = 1.0, b = 0.0, c = 0.0", "a = 1.0, b = 0.0, c = 0.5", "a = 1.0, b = 0.0, c = 0.1"), 
       col = c("black", "red", "blue"),
       lty = 1:4, cex = 1.0)

なお,上述した加藤他 (2014, p.78-79) では,この\(c\)が推定されることによって,3PLモデルの解釈は2PLモデルとは同じ解釈ができなくなることを指摘しています。具体的には,

3PLモデルでは,

  1. 3PLモデルと2PLモデルでのaの値が同じだったとしても,それらの項目がθを同程度に識別できるということを意味するわけではない。
  1. 3PLモデルと2PLモデルでのbの値が同じだったとしても,それらの項目が同程度に難しいということを意味するわけではない。

ということが生じます。これは上記の2つの事柄が,どちらもcの影響を受けるためです。項目の正答確率からcの分を取り去った残りである1-c に対して,2PLモデルのICCを当てはめている状況です。

と述べられています。こうしたことにより,3PLモデルを具体的に解釈することは,非常に困難であり,実際の研究でも2PLモデルくらいまでで,それよりも複雑なモデルが用いられていることは少ないように思います。

3.1 3PLモデルをStanでやってみる

3.1.1 Stan コード

data{
  int N; // 回答者数
  int n; // 項目数
  int<lower=0, upper=1> mu[N,n]; 
  // 各回答者の各項目への回答
}

parameters{
  real<lower=0, upper=5> a[n]; //識別度パラメータ
  real<lower=-4, upper=4> b[n]; // 困難度パラメータ
  real<lower=0, upper=1> c[n]; //  当て推量パラメータ
  real theta[N]; // 能力値
}

model{
  for(i in 1:N){
    for(j in 1:n){
      p= inv_logit( a[j] * (theta[i] - b[j]) );
      mu[i,j] ~ bernoulli( c[j] + (1-c[j]) * p );
      }
    }

  /// 事前分布
  a ~ lognormal(0,1); 
  b ~ normal(0,1);  
  c ~ beta(5,23);  
  theta ~ normal(0,1);
}


generated quantities {
  vector[n] log_lik[N]; ///// 対数尤度
  int<lower=0, upper=1> mu_rep[N,n]; // 予測分布

  for(i in 1:N){
    for(j in 1:n){
      p = inv_logit( a[j] * (theta[i] - b[j])) 
      log_lik[i,j] = bernoulli_lpmf(mu[i,j] | c[j] + (1-c[j])*p );
      mu_rep[i,j] = bernoulli_rng( c[j] + (1-c[j]) *p );
    }
  }
}

3.1.2 3PLモデル用のデータ生成Stan/Rコード

data{
  int P;
  vector[P] a;
  vector[P] b;
  vector[P] c;
}

parameters{
  real temp;
}

model{
  temp ~ uniform(0,1);
}

generated quantities{
 int y[P];
 real theta;
 theta = normal_rng(0,1);
  for(p in 1:P){
    real z;
    z=inv_logit( a[p] * (theta - b[p]));
    y[p] = bernoulli_rng( c[p] + (1-c[p]) *z);
  }
}

3.1.3 3PLモデル用のデータを生成するRコード

library(rstan)

# セッティング
P = 6  #項目数
a_true = rep(1.7,P)  #識別力パラメータ
b_true = seq(3,-3,length.out = P) #困難度パラメータ
c_true = seq(0.2,0.5,length.out = P) #あて推量母数

# 3PLモデルのstanコードを読み込み
model.random_3pl <- stan_model("stan_code/3PL_random_normal.stan")

# データの生成

N <- 300 #サンプルサイズ
set.seed(123)
data.random=list(P=P,a=a_true,b=b_true,c=c_true)
rdata <- sampling(model.random_3pl,data=data.random,iter=N,warmup=0,chains=1)
dat <- as.data.frame((rstan::extract(rdata)$y))
theta_true <- as.data.frame((rstan::extract(rdata)$theta))
N <- nrow(dat)
P <- ncol(dat)
datastan <- list(N=N,n=P,mu=as.matrix(dat))

3.1.4 できたデータセット

こんなデータができあがります。 上から20人分のデータを示しています。

library(knitr)
kable(datastan$mu[1:20,])
V1 V2 V3 V4 V5 V6
0 0 1 1 1 1
0 0 0 1 1 1
1 1 0 0 1 1
0 0 1 1 1 1
0 1 1 1 1 1
0 0 1 0 1 1
0 0 1 1 1 1
0 1 0 0 1 1
0 0 1 1 1 1
0 0 0 1 1 0
0 0 0 0 1 1
1 1 1 1 1 1
0 1 1 1 0 1
1 0 0 0 1 1
1 0 1 1 1 1
0 0 0 1 1 1
0 0 0 1 1 1
0 0 0 1 1 1
0 0 0 1 0 1
0 0 1 1 1 1

3.1.5 推定の実行

# stanコードのコンパイル
model.3PL <- stan_model("stan_code/ThreePL.stan")
fit.mcmc <- sampling(model.3PL,data=datastan)

3.1.6 推定結果の確認

tidybayesを使うことでmcmcの結果を短いコードでサクッと整然化できます。 使い方の詳細は下記をご参照ください。

tidybayesパッケージで推定結果の整然化

3.1.6.1 事後分布の要約統計量の確認

パラメータa, bの中央値と95%確信区間の上限値と下限値の推定結果は下記のように確認できます。

library(tidyverse)
library(tidybayes)


fit.mcmc %>%
  spread_draws(a[n],b[n]) %>%
      median_qi(.width = c(.95)) ->ab
ab
## # A tibble: 6 x 10
##       n     a a.lower a.upper      b b.lower b.upper .width .point .interval
##   <int> <dbl>   <dbl>   <dbl>  <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
## 1     1 1.95    0.824    4.56  1.66    1.12    2.59    0.95 median qi       
## 2     2 1.30    0.505    4.19  1.74    1.09    2.82    0.95 median qi       
## 3     3 1.07    0.361    3.07  0.111  -0.344   0.733   0.95 median qi       
## 4     4 0.971   0.424    2.56 -0.932  -1.91   -0.346   0.95 median qi       
## 5     5 1.93    1.00     4.45 -1.92   -2.92   -1.38    0.95 median qi       
## 6     6 2.80    1.43     4.84 -2.55   -3.48   -2.02    0.95 median qi

パラメータa, bの中央値と95%,80%,50%確信区間をプロットするためには下記のようにするとできます。

fit.mcmc %>%
  spread_draws(a[n],b[n]) %>%
      median_qi(.width = c(.95, .8,.5)) %>% 
        ggplot(aes(y = n, x = a, 
                   xmin = a.lower, 
                   xmax = a.upper, size = -.width))+
        geom_pointinterval(interval_size_range = c(0.5, 2))

3.1.6.2 項目特性曲線の描画

#dataframe for ggplot2
theta <- seq(-5,5,length.out=201)
P<-datastan$n
ggdf <- data.frame(matrix(ncol=,nrow=201))


for(i in 1:P){
  logistic_f<-function(x) 1/(1+exp(-x))
  ggdf[,i] <- logistic_f(ab$a[i]*(theta-ab$b[i]))
}

ggdf$theta <- theta

#item characteristic curve
ggdf %>% 
  tidyr::gather(key=var,value,-theta,factor_key=TRUE) %>%
  ggplot(aes(x=theta)) + geom_line(aes(y=value,color=var))

3.1.6.3 項目情報曲線, テスト情報曲線

ggdf2 <- data.frame(matrix(ncol=P,nrow=201))

for(i in 1:P){
  p_theta= 1/(1+ exp(-ab$a[i]*(theta - ab$b[i])))
  q_theta=1-p_theta
  ggdf2[,i]<-(ab$a[i]^2)*p_theta*q_theta
}

ggdf2$theta <- theta

# item info
ggdf2 %>% 
  tidyr::gather(key=var,value,-theta,factor_key=TRUE) %>%
  ggplot(aes(x=theta))+ geom_line(aes(y=value,col=var))+
  ylab("Item information") ->Iinfo


#testinfo
ggdf2$testinfo<-apply(ggdf2[,1:6],1,sum)

ggdf2 %>% 
  ggplot(aes(x=theta))+ geom_line(aes(y=testinfo)) +
  ylab("Test information")->Tinfo

library(patchwork)
Iinfo|Tinfo

3.1.6.4 情報量基準

モデル選択のための指標としてWAICやWBICがありますが、generated quntitiesで対数尤度を抽出するようにコードを書いたので、looパッケージなどでそれらを求めることができます。

library(loo)
log_lik <- extract_log_lik(fit.mcmc)
waic(log_lik)
## 
## Computed from 4000 by 1800 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -780.4 21.3
## p_waic       131.9  7.0
## waic        1560.9 42.6
## 
## 29 (1.6%) p_waic estimates greater than 0.4. We recommend trying loo instead.

あとはさらっと紹介します。

4 段階反応モデル (graded response model: GRM)

GRMは,項目に対する多値の反応(例えば質問紙調査にみられるようなリッカート尺度に対する反応)でも,2PL IRTモデルを適用できるように一般化したものであるといえます3。GRMでは多値の反応をモデリングするために,あるカテゴリ以下の反応とそれ以上の反応の差分をモデリングします。以下のような式で表すことができます4

\[ p_{ij} (u_{ij} = k| \theta_i, a_j, b_{jk}) = \frac{1}{1 + {\rm exp}(-a_j(\theta_i - b_{jk}))} - \frac{1}{1 + {\rm exp}(-a_j(\theta_i - b_{jk + 1}))} \]

  • \(k\): これはカテゴリ数を表す。例えば5段階のリッカート式の尺度への回答であれば,\(k = 5\)となる。

  • \(b_{jk}\): \(項目_j\)のあるカテゴリに対する困難度。

あとのパラメタは3PLモデルと一緒です。

4.1 段階反応モデルのStanコード

論文に紹介されているStanコードです。

data{
   int<lower=2, upper=4> K; //反応カテゴリ数
   int <lower=0> N; //回答者数
   int <lower=0> n; //項目数
   int<lower=1,upper=K> mu[N,n];
 }

 parameters {
   vector[N] theta;
   real<lower=0> alpha [n];
   ordered[K-1] kappa[n]; //カテゴリーの困難度
   real mu_kappa; //カテゴリーの困難度の事前分布の平均
   real<lower=0> sigma_kappa; //カテゴリーの困難度の事前分布の標準偏差
 }

 model{
   alpha ~ cauchy(0,5);
   theta ~ normal(0,1);
 for (i in 1: n){
   for (k in 1:(K-1)){
     kappa[i,k] ~ normal(mu_kappa,sigma_kappa);
   }}
 mu_kappa ~ normal(0,5);
 sigma_kappa ~ cauchy(0,5);
 for (i in 1:N){
  for (j in 1:n){
    mu[i,j] ~ ordered_logistic(theta[i]*alpha[j],kappa[j]);
 }}
 }
 
 generated quantities {
 vector[n] log_lik[N];
   for (i in 1: N){
     for (j in 1: n){
     log_lik[i, j] = ordered_logistic_lpmf (mu[i, j]|theta[i]*alpha[j],kappa[j]);
   }}
  }

4.1.1 GRMの実行

下記のようにモデルをコンパイルして、ltmパケージに含まれているScienceデータから4項目選んで、GRMを適用すると下記のようになります。

model.GRM <- stan_model("stan_code/GRM.stan")

library(ltm)
data<-sapply(Science[c(1,3,4,7)],as.numeric)

N<-nrow(data)
n<-ncol(data)
K<-length(table(data[,4]))

datastan<-list(N=N,n=n,K=K,mu=as.matrix(data))
fit.mcmc_grm <- sampling(model.GRM,data=datastan)

5 名義反応モデル (nominal response model: NRM)

NRMは,GRMと同様に多値型の反応に対してIRTを適用できるようにしたモデルで,ある問題に対して,正答と複数の誤答があり,誤答には特に難易度が設定されていないような場合に適用できるモデルとなっているらしいです。NRMは,以下のように表現できます。

\[ p_{ij} (u_{ij} = k| \theta_i, a_{jk}, \gamma_{jk}) = \frac{{\rm exp}(a_{jk}\theta_i + \gamma_{jk})}{\sum_{h = 1} ^{c} {\rm exp} (-a_{jh}\theta_i + \gamma_{jh}))} \]

ソフトマックス関数のようになっていますね。

  • \(a_{jk}, \gamma_{jk}\): \(カテゴリ_k\)におけるICCのような曲線の傾きと切片に影響するパラメタ。

よくわからなかったので,\(\gamma_{jk}\)を一定の値にして,\(a_{jk}\)をいろいろと変更してみたら,曲線の傾きが主に影響を受ける様子なので,\(カテゴリ_k\)における識別力みたいなものを意味していそうです。

curve(exp(1 * x + 1)/(exp(1.0*x + 1) + exp(0.2*x + 1) + exp(0.1*x + 1)), -3, 3, col = "black", lty = 1, ylab = "probability") # a_11 = 1.0, a_12 = 1.0, gamma_11 = 1.0, gamma_12 = 1.0
curve(exp(0.4 * x + 1)/(exp(0.4*x + 1) + exp(1.0*x + 1) + exp(1.0*x + 1)), -3, 3, col = "red", lty = 2, ylab = "probability", add = TRUE) # a_11 = 1.0, a_12 = 0.4, a_13 = 0.1, gamma_11 = 1.0, gamma_12 = 1.0, gamma_13 = 1.0,
curve(exp(-0.5 * x + 1)/(exp(-0.5*x + 1) + exp(1*x + 1) + exp(0.5*x + 1)), -3, 3, col = "blue", lty = 2, ylab = "probability", add = TRUE) # a_11 = 1.0, a_12 = 0.4, a_13 = 0.1, gamma_11 = 1.0, gamma_12 = 1.0, gamma_13 = 1.0,
curve(exp(1.0 * x + 1)/(exp(1.0*x + 1) + exp(0.3*x + 1) + exp(0.3*x + 1)), -3, 3, col = "green", lty = 2, ylab = "probability", add = TRUE) # a_11 = 1.0, a_12 = 0.4, a_13 = 0.1, gamma_11 = 1.0, gamma_12 = 1.0, gamma_13 = 1.0,

今度は,上の図の\(a_{jk}\)の値のまま,\(\gamma_{jk}\)をいろいろと変更してみたら,曲線の傾きもちょっと変わっている気もするが主に切片が影響を受けている様子。うーん,たぶん困難度のようなパラメタを意味しているのだろうか。。よくわからない。

curve(exp(1 * x + 1)/(exp(1.0*x + 1) + exp(0.2*x + 1) + exp(0.1*x + 1)), -3, 3, col = "black", lty = 1, ylab = "probability") # a_11 = 1.0, a_12 = 1.0, gamma_11 = 1.0, gamma_12 = 1.0
curve(exp(0.4 * x + 0.1)/(exp(0.4*x + 0.1) + exp(1.0*x + 1) + exp(1.0*x + 1)), -3, 3, col = "red", lty = 2, ylab = "probability", add = TRUE) # a_11 = 1.0, a_12 = 0.4, a_13 = 0.1, gamma_11 = 1.0, gamma_12 = 1.0, gamma_13 = 1.0,
curve(exp(-0.5 * x - 1)/(exp(-0.5*x - 1) + exp(1*x + 1) + exp(0.5*x + 1)), -3, 3, col = "blue", lty = 2, ylab = "probability", add = TRUE) # a_11 = 1.0, a_12 = 0.4, a_13 = 0.1, gamma_11 = 1.0, gamma_12 = 1.0, gamma_13 = 1.0,
curve(exp(1.0 * x + 0.5)/(exp(1.0*x + 0.5) + exp(0.3*x + 1) + exp(0.3*x + 1)), -3, 3, col = "green", lty = 2, ylab = "probability", add = TRUE) # a_11 = 1.0, a_12 = 0.4, a_13 = 0.1, gamma_11 = 1.0, gamma_12 = 1.0, gamma_13 = 1.0,

NRMについてはちゃんとBock (1972) を確認してみないといけないですね。

NRMは3PLモデルやGRMと比べると知らない人が多いかもしれませんが、使い方次第でとても面白いリサーチクエスチョンの検討を可能にしてくれます。例えば、回答の選択肢間に順序性を仮定していないという仮定を生かして、順序性の仮定があるモデル(一般化部分採点モデル: GPCM)と比較することで、回答選択肢の中間反応の適切さを検討する際に用いられています。その例の紹介は下記のスライド43ページくらいからをご参照ください。

しなやかな項目反応モデル

5.1 NRMのStanコード

data{
   int<lower=2, upper=4> K;
   int <lower=0> N;
   int <lower=0> n;
   int<lower=1,upper=K> mu[N,n];
 }
 
 parameters {
   vector[K] zeta[n]; //freely estimated intercept
   vector[K] lambda[n]; //freely estimated slope
   vector[N] theta;
 }
 
 transformed parameters {
   vector[K] zetan[n]; //intercept with constraints
   vector[K] lambdan[n]; //slope with constraints
   for (k in 1:n) {
     for (l in 1:K) {
       zetan[k,l]=zeta[k,l]-mean(zeta[k]); //constrain intercept sum for each item to 0
       lambdan[k,l]=lambda[k,l]-mean(lambda[k]); //constrain slope sum for each item to 0
     }}
 }
 
 model{
   theta ~ normal(0,1);
   for (i in 1: n){
   zeta[i] ~ normal(0,2);
   lambda[i] ~ normal(0,2);
 }
 for (i in 1:N)
   for (j in 1:n)
     mu[i,j] ~ categorical_logit(zetan[j] + lambdan[j]*theta[i]);
 }
 
 generated quantities {
   vector[n] log_lik[N];
   for (i in 1: N){
     for (j in 1: n){
       log_lik[i, j] = categorical_logit_lpmf (mu[i, j]| zetan[j] + lambdan[j]*theta[i]);
 }}
 }

5.1.0.1 NRMを実行するRコード、先ほどのGRMを適用したデータを使って実行してみましょう。

model.NRM <- stan_model("stan_code/NRM.stan")
fit.mcmc_nrm <- sampling(model.NRM,data=datastan)

6 ラッシュテストレットモデル

ラッシュテストレットモデルは多次元に拡張されたラッシュモデルです。ラッシュテストレットモデルは,以下のように表現できる。

\[ p_{i}(\theta_j) = \frac{1}{1 + {\rm exp}(-(\theta_i - b_j + \gamma_{jd(i)}))} \]

ラッシュモデルは識別力を1に固定し,困難度のみ推定するモデルであることを思い出すと,上の式はラッシュモデルに加えて新たに\(\gamma_{jd(i)}\)を推定しているモデルだということがわかります。この\(\gamma_{jd(i)}\)は回答者とテストレットの交互作用を表現したものになります。なお,\(\gamma_{jd(i)}\)の分散である\(\sigma_{\gamma_{jd(i)}}^2\)はテストレット効果の大きさを意味します。つまり,この値が大きければ,分割したテストによって結果がかなり異なってきていると思ってよいでしょう。

7 マルチレベル3PLモデル

最後は,マルチレベル3PLモデルです。マルチレベル3PLモデルは,以下のように表現できます。

\[ p_{i}(\theta_{gj}) = c_i + \frac{1 - c_i}{1 + {\rm exp}(-a_i(\theta_{gj} + \theta_g - b_j))} \]

これまでのモデルで出てこなかった新たなパラメタが何を意味しているかを説明していきます。

  • \(\theta_{gj}\): \(グループg\)にいる回答者\(j\)の能力値

  • \(\theta_{g}\): \(グループg\)の平均的な能力値

あとは3PLモデルのときと同じです。

7.1 最後に

ラッシュテストレットモデル、マルチレベル3PLモデルのStanコードもLeo & jiao (2017)に紹介されているので、関心のある方は是非論文をみて実行してみてください。

その他、項目反応理論は複数の因子(能力値)を仮定したモデルにも拡張可能です。それらの例は下記を参照ください。

Stanで多次元項目反応理論


  1. ただ,加藤他 (2014)では,二つのモデルは数学的には同一のモデルとして扱われています。加藤他(2014)では,1PLモデルを\(p_{j} (\theta) = \frac{1}{1 + {\rm exp}(-Da(\theta_i - b_j))}\)と表現しており,識別力(\(a\))と尺度因子(\(D\))を分けて考え,\(D\)を定数(\(D = 1.7\))としています。\(Da\)をまとめて1と置くこともできる述べられており,そのように考えた場合は1PLモデルとラッシュモデルの数学的な違いはないものとして考えられるのだと思われます。なお,加藤他 (2014)では,ラッシュモデルと1PLモデルは,モデルの理念や考え方に大きな違いがあると別の観点から説明されています。↩︎

  2. ただ,加藤他 (2014)では,二つのモデルは数学的には同一のモデルとして扱われています。加藤他(2014)では,1PLモデルを\(p_{j} (\theta) = \frac{1}{1 + {\rm exp}(-Da(\theta_i - b_j))}\)と表現しており,識別力(\(a\))と尺度因子(\(D\))を分けて考え,\(D\)を定数(\(D = 1.7\))としています。\(Da\)をまとめて1と置くこともできる述べられており,そのように考えた場合は1PLモデルとラッシュモデルの数学的な違いはないものとして考えられるのだと思われます。なお,加藤他 (2014)では,ラッシュモデルと1PLモデルは,モデルの理念や考え方に大きな違いがあると別の観点から説明されています。↩︎

  3. ただし,識別力は全てのカテゴリで共通↩︎

  4. 加藤他 (2014, p.221-226)を読むと,GRMはあるカテゴリへの反応とそれ以上の反応の累積確率曲線を用いると書いてあります。累積確率曲線とは,例えば,ある問題に対する得点が0,1,2,3の4つの得点があり得た場合に,2点以上となる確率を合算したもの。さらに,隣接するカテゴリの累積確率曲線の差を用いてカテゴリ確率曲線という2値反応に対するIRTのICCに相当するものを求めることができると書いてあります。今回参考にしているLeo & jiao (2017)では,Stanコード(Listing 2)見ても累積確率は使っていないように見えて,その辺りどうなっているのかよくわかっていない。↩︎