廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

正規累積項目反応曲線のギブスサンプリングによる推定. Rcpp による実装例 (2)

今日の川柳

Albert (1992): 正規累積項目反応曲線のギブスサンプリングによる推定. Rcpp による実装例. - 廿TT のつづきです.


Albert, J. H. (1992). Bayesian estimation of normal ogive item response curves using Gibbs sampling. Journal of educational statistics, 17(3), 251-269. https://www.jstor.org/stable/1165149 を参考にしていますがちょっといじっています.

モデルと分析対象

正解, 不正解の2値に明確に分類できるような試験問題を考えます.

被験者 ij 問目への回答の正解, 不正解を変数 y_{i,j} で表します(1が正解, 0が不正解).

6人の被験者が6問に答えた試験の場合、下表みたいなデータがとれます(行が被験者, 列が問題).

0 0 0 1 1 0
0 0 1 1 1 0
1 0 1 1 1 1
0 0 0 1 1 0
0 0 1 1 1 0
1 1 0 1 1 0

2パラメータ正規累積モデル(2PL normal ogive model)では, 被験者 i (i=1,\ldots,n) の j (j=1,\ldots,m) 問目に正解する確率を,

 \displaystyle p_{i,j}=\Phi(\alpha_j(\theta_i-\gamma_j);0,1)

とモデル化します.

(Albert (1992) では

 \displaystyle p_{i,j}=\Phi(\alpha_j\theta_i-\gamma_j;0,1)

としていますが, ギブスサンプリングに使う条件付き分布の導出がちょっと面倒なので変えました.)

\Phi(x;\mu,\sigma^2) は平均 \mu, 分散 \sigma^2正規分布の分布関数を表す記号で, \Phi(x;0,1) は標準正規分布の累積分布関数です.

 \theta_i は能力パラメータと呼ばれ, 被験者の能力を表し, \gamma_i は困難度パラメータと呼ばれ, 問題の難しさを表すと解釈されます.

\alpha_j は識別パラメータとよばれます.

そこで  \theta_i の事前分布に, 平均 0, 分散 1 の正規分布, \gamma_j の事前分布にも平均 0, 分散 1 の正規分布を, \alpha_j の事前分布にパラメータ 1, 1 のガンマ分布を仮定して, ベイズ法でパラメータを推定します.

潜在変数を用いた解釈

\Phi(x) は単に p_{i,j} を 0 から 1 の範囲に押さえ込むためのリンク関数と捉えることもできますが, Albert (1992) はちょっとおもしろい解釈を推奨しています.

テストの背後には被験者の発揮するパフォーマンスみたいな潜在変数 Z_{i,j} があり, それが平均 \theta_i-\gamma_j正規分布に従うとします.

パフォーマンスが閾値0を超えた場合に正解にたどり着き, 0以下のとき誤答すると考えると,

正解のとき: \Phi(0;\theta_i-\gamma_j,1/\alpha_j)^{y_{i,j}}
不正解のとき: \{1-\Phi(0;\theta_i-\gamma_j,1/\alpha_j)\}^{y_{i,j}}

となり, 正規累積モデルの尤度と一致します.

離散的なデータの背後に連続的な潜在変数を考えるのはアルバート先生の得意技らしく, 順序データの分析の本でも似た議論が出てきます.

Ordinal Data Modeling (Statistics for Social and Behavioral Sciences)

Ordinal Data Modeling (Statistics for Social and Behavioral Sciences)

識別パラメータ \alpha_j は潜在変数の分散をコントロールする集中度パラメータです. \alpha_j が大きいほどパフォーマンスを精度よく発揮できる問題ということです.

これが項目反応理論のおもしろいところだと思います. 被験者を評価すると同時に問題も評価されているのです.

ギブスサンプリング

ギブスサンプリングそのものについては、

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

などを参照してください.

潜在変数 Z_{i,j} の実現値 z_{i,j} が得られているとします.

\alpha_j, \gamma_j が得られたときの \theta_i の事後分布は以下に比例します. パラメータ\theta_i に依存しない因子は無視するのがコツです.


 \displaystyle  \exp\left(-\frac{\sum_{j=1}^m(\alpha_j ( z_{i,j}-(\theta_i-\gamma_j) ))^2}{2}\right) \\
 \displaystyle \quad \times \exp\left(-\frac{\theta_i^2}{2}\right)\\
 \displaystyle =  \exp\left(-\frac{\sum_{j=1}^m\alpha_j\{(\theta_i-\gamma_j)-z_{i,j})\}^2}{2}\right)\\
\displaystyle \quad \times \exp\left(-\frac{\theta_i^2}{2}\right)\\
 \displaystyle \propto  \exp\left(-\frac{\theta_i^2 (\sum_{j=1}^m\alpha_j+1)-2\sum_{j=1}^{m}\alpha_j(\gamma_j+z_{i,j})\theta_i}{2}\right)\\
 \displaystyle =  \exp\left(-\frac{\theta_i^2 -2\theta_i\{\sum_{j=1}^{m}\alpha_j(\gamma_j+z_{i,j})\}/(\sum_{j=1}^m\alpha_j+1)}{2/(\sum_{j=1}^m\alpha_j+1)}\right)

これは正規分布の密度関数です.

\alpha_j, \theta_i が得られたときの \gamma_j の事後分布は以下に比例します.


 \displaystyle \exp\left(-\frac{\sum_{j=1}^n\alpha_j(z_{i,j}-(\theta_i-\gamma_j))^2}{2}\right) \\
 \displaystyle \quad \times \exp\left(-\frac{\gamma_j^2}{2}\right)\\
\displaystyle \propto \exp\left(-\frac{(\gamma_j^2-2\gamma_j \alpha_j \sum_i^n(\theta_i-z_{i,j})/(\alpha_j (n+1) )}{2/(\alpha_j (n+1))}\right)

正規分布の密度関数です.

\theta_i, \gamma_j が得られたときの \alpha_j の事後分布は以下に比例します.


 \displaystyle  \alpha_j ^{-n/2} \exp\left(-\alpha_j \frac{\sum_{i=1}^n( z_{i,j}-(\theta_i-\gamma_j))^2}{2}\right)\\
\displaystyle \quad \times \alpha^{1-1}\exp\left(-1 \alpha)\right)\\
\displaystyle = \alpha_j ^{-n/2+1}\exp\left(-\alpha_j \frac{\sum_{i=1}^n( z_{i,j}-(\theta_i-\gamma_j))^2+2}{2}\right)

これはガンマ分布の密度関数です.

アルゴリズムを Rcpp で書くと以下のようになりました.

#include <Rcpp.h>
using namespace Rcpp;

double rnorm_u (double eta, double sigma){
  return R::qnorm(R::runif(R::pnorm(0,eta,sigma,true,false),1),eta,sigma,true,false);
}

double rnorm_l (double eta, double sigma){
  return R::qnorm(R::runif(0,R::pnorm(0,eta,sigma,true,false)),eta,sigma,true,false);
}

// [[Rcpp::export]]
List Gibbs1PLNO(int iter, IntegerMatrix y){
  int n;
  int m;
  n = y.nrow();
  m = y.ncol();
  NumericMatrix theta(iter,n);
  NumericMatrix gamma(iter,m);
  NumericMatrix alpha(iter,m);
  alpha(0,_) = rgamma(m,1,1);
  NumericMatrix Z(n,m);
  for(int k=1; k<iter; k++){
    for(int i=0; i<n; i++){
      for(int j=0; j<m; j++){
        if(y(i,j)==1){
          Z(i,j) = rnorm_u((theta(k-1,i)-gamma(k-1,j)),1/sqrt(alpha(k-1,j)));
        }else{
          Z(i,j) = rnorm_l((theta(k-1,i)-gamma(k-1,j)),1/sqrt(alpha(k-1,j)));
        }
      }
    }
    for(int j=0; j<m; j++){
      alpha(k,j) = R::rgamma(n/2.0+1,2.0/(sum(pow(Z(_,j)-(theta(k-1,_)-gamma(k-1,j)),2))+2));
      gamma(k,j) = R::rnorm(sum(alpha(k,j)*(theta(k-1,_)-Z(_,j)))/(n*alpha(k,j)+1),1/sqrt(n*alpha(k,j)+1));
    }
    for(int i=0; i<n; i++){
      theta(k,i) = R::rnorm(sum(alpha(k,_)*(Z(i,_)+gamma(k,_)))/(sum(alpha(k,_))+1),1/sqrt(sum(alpha(k,_))+1));
    }
  }
  return List::create(Named("ability") = theta, _["item"]=gamma, _["alpha"]=alpha);
}

シミュレーション

乱数でデータを作ってちゃんと推定できるか確かめてみます.

100人に20問の試験を行ったことにしました.

library(Rcpp)
library(parallel)
sourceCpp("~/Documents/albert1992.cpp")
set.seed(1)
n=100;m=20
theta <- rnorm(n,0,1)
gamma <- rnorm(m,0,1)
alpha <- rgamma(m,1,1)
y <- matrix(NA,n,m)
for(i in 1:n){
  for(j in 1:m){
    y[i,j] <- rbinom(1,1,pnorm(alpha[j]*(theta[i]-gamma[j])))
  }
}

###
fitlist <-mclapply(1:4,function(i){
  Gibbs1PLNO(4000,y)
},mc.cores = detectCores())
a_samples <-lapply(fitlist, function(x)x$ability)
i_samples <-lapply(fitlist, function(x)x$item)
alpha_samples <-lapply(fitlist, function(x)x$alpha)

col4<-c("royalblue","salmon","forestgreen","orange2")
plot(a_samples[[1]][,1],type="l",col = col4[1], lty=1,ylab="theta")
for(i in 2:4)lines(a_samples[[i]][,1],type="l",col = col4[i], lty=1)

plot(i_samples[[1]][,1],type="l",col = col4[1], lty=1,ylab="gamma")
for(i in 2:4)lines(i_samples[[i]][,1],type="l",col = col4[i], lty=1)

plot(alpha_samples[[1]][,1],type="l",col = col4[1], lty=1,ylab="alpha")
for(i in 2:4)lines(alpha_samples[[i]][,1],type="l",col = col4[i], lty=1)

ex_a_samples_l <-lapply(a_samples, function(x){x[-c(1:2000),]})
ex_i_samples_l <-lapply(i_samples, function(x){x[-c(1:2000),]})
ex_alpha_samples_l <-lapply(alpha_samples, function(x){x[-c(1:2000),]})

Rhat <- function(samples_l){
  samples <- simplify2array(samples_l)
  term2 <-apply(apply(samples, 2:3, mean),1,var)
  varW <-apply(apply(samples, 2:3, var),1,mean)
  n.r <-nrow(samples)
  var1 <-((n.r-1)/n.r)*varW+term2
  sqrt(var1/varW)
}

all(Rhat(ex_a_samples_l)<1.1)
all(Rhat(ex_i_samples_l)<1.1)
all(Rhat(ex_alpha_samples_l)<1.1)

ex_a_samples <-do.call("rbind",ex_a_samples_l)
ex_i_samples <-do.call("rbind",ex_i_samples_l)
ex_alpha_samples <-do.call("rbind",ex_alpha_samples_l)

hattheta <-apply(ex_a_samples,2,mean)
hatgamma <-apply(ex_i_samples,2,mean)
hatalpha <-apply(ex_alpha_samples,2,mean)

ran1 <-range(hattheta,theta)
ran2 <-range(hatgamma,gamma)
ran3 <-range(hatalpha,alpha)
plot(theta,hattheta,xlim=ran1,ylim=ran1)
abline(0,1,lty=2)
plot(gamma,hatgamma,xlim=ran2,ylim=ran2)
abline(0,1,lty=2)
plot(alpha,hatalpha,xlim=ran3,ylim=ran3)
abline(0,1,lty=2)
a_CI <-t(apply(ex_a_samples,2,quantile,probs=c(0.025,0.975)))
mean(a_CI[,1]<=theta & theta <=a_CI[,2])
#[1] 0.95
i_CI <-t(apply(ex_i_samples,2,quantile,probs=c(0.025,0.975)))
mean(i_CI[,1]<=gamma & gamma <=i_CI[,2])
#[1] 1
alpha_CI <-t(apply(ex_alpha_samples,2,quantile,probs=c(0.025,0.975)))
mean(alpha_CI[,1]<=alpha & alpha <=alpha_CI[,2])
# [1] 0.8

なんとなく近い値が求まっていそうなことがわかります.

f:id:abrahamcow:20190807030417p:plain

f:id:abrahamcow:20190807030443p:plain

f:id:abrahamcow:20190807030500p:plain

付録

同じモデルをStanで書くとこうなります.

data{
  int<lower=1> n;
  int<lower=1> m;
  int<lower=0,upper=1> y[n,m];
}
parameters{
  vector[n] theta;
  vector[m] gamma;
  vector<lower=0>[m] alpha;
}
model{
  theta ~ normal(0,1);
  gamma ~ normal(0,1);
  alpha ~ gamma(1,1);
  for(i in 1:n){
    for(j in 1:m){
     if(y[i,j]==1){
       target += normal_lccdf(0|alpha[j]*(theta[i]-gamma[j]),1);
     }else{
       target += normal_lcdf(0|alpha[j]*(theta[i]-gamma[j]),1);
     }
    }
  }
}