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コードと結果の抽出方法について紹介します。

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の値が同じだったとしても，それらの項目がθを同程度に識別できるということを意味するわけではない。

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

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

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

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 );
    }
  }
}

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);
  }
}

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))

library(knitr)
kable(datastan$mu[1:20,])

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

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

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))

#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))

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

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モデルを適用できるように一般化したものであるといえます³。GRMでは多値の反応をモデリングするために，あるカテゴリ以下の反応とそれ以上の反応の差分をモデリングします。以下のような式で表すことができます⁴。

\[ 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モデルと一緒です。

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]);
   }}
  }

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ページくらいからをご参照ください。

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

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]);
 }}
 }

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モデルのときと同じです。