廿TT

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

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

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

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

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

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

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

Albert (1992) のモデルではここにさらに識別パラメータというのが入るんだけど, 今回は省略します.

このモデルはこのままでは未知パラメータどうしの引き算により, 識別不能になります.

たとえば,  \theta_i が 2 で \gamma_j が 1 だとすれば引き算すれば 1 です。平行移動して \theta_i が 3 で \gamma_j が 2 だとしても引き算すれば 1 で, 尤度はまったく一緒です. 尤度がまったく一緒になるパラメータの組み合わせがいくつも存在するので最尤推定量は一致性を持ちません.

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

潜在変数を用いた解釈

\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)^{y_{i,j}}
不正解のとき: \{1-\Phi(0;\theta_i-\gamma_j,1)\}^{y_{i,j}}

となり, 1パラメータ正規累積モデルの尤度と一致します.

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

Ordinal Data Modeling (Statistics for Social and Behavioral Sciences)

Ordinal Data Modeling (Statistics for Social and Behavioral Sciences)

ギブスサンプリング

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

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

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

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

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

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


 \displaystyle  \exp\left(-\frac{\sum_{j=1}^m(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\{(\theta_i-\gamma_j)-z_{i,j}\}^2}{2}\right)\\
\displaystyle \quad \times \exp\left(-\frac{\theta_i^2}{2}\right)\\
 \displaystyle =  \exp\left(-\frac{\sum_{j=1}^m\{(\theta_i-(\gamma_j+z_{i,j})\}^2}{2}\right)\\
\displaystyle \quad \times \exp\left(-\frac{\theta_i^2}{2}\right)\\
 \displaystyle = \exp\left(-\frac{(m\theta_i-\sum_{j=1}^m(\gamma_j+z_{i,j}))^2}{2}\right)\\
 \displaystyle \quad \times \exp\left(-\frac{\theta_i^2}{2}\right)\\
 \displaystyle =\exp\left(-\frac{\{(m+1)\theta_i-\sum_j(\gamma_j+z_{i,j})\}^2}{2}\right)\\
 \displaystyle =  \exp\left(-\frac{(\theta_i-\sum_j(\gamma_j+z_{i,j})/(1+m))^2}{2(1+m)^2}\right)

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

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


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

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

Albert (1992) のアルゴリズムをまとめると以下のようになります.

  1. 初期値を適当に決める(今回は 0 とした)
  2. \theta_i, \gamma_j が得られたとき潜在変数 Z_{i,j}を切断正規分布からサンプリング
  3. Z_{i,j}, \gamma_j が得られたときの \theta_i正規分布からサンプリング
  4. Z_{i,j}, \theta_i が得られたときの \gamma_j正規分布からサンプリング
  5. 2, 3, 4 を繰り返す

Rcpp で書くと以下のようになりました.

#include <Rcpp.h>
using namespace Rcpp;

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

double rnorm_l (double eta){
  return R::qnorm(R::runif(0,R::pnorm(0,eta,1,true,false)),eta,1,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 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));
        }else{
          Z(i,j) = rnorm_l(theta(k-1,i)-gamma(k-1,j));
        }
      }
    }
    for(int i=0; i<n; i++){
      theta(k,i) = R::rnorm(sum(Z(i,_)+gamma(k-1,_))/(m+1),1/sqrt(m+1));
    }
    for(int j=0; j<m; j++){
      gamma(k,j) = R::rnorm(sum(theta(k,_)-Z(_,j))/n,1/sqrt(n));
    }
  }
  return List::create(Named("ability") = theta , _["item"]=gamma);
}

シミュレーション

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

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

library(Rcpp)
library(parallel)
#データを生成
set.seed(2)
n=100;m=20
theta <- rnorm(n,0,1)
gamma <- rnorm(m,0,1)
y <- matrix(NA,n,m)
for(i in 1:n){
  for(j in 1:m){
    y[i,j] <- rbinom(1,1,pnorm(theta[i]-gamma[j]))
  }
}

MCMC が 1 系列しかないとちゃんと同じところに収束するか不安なので, 並列して4つの chain を回すことにします.

#推定
sourceCpp("albert1992.cpp")
set.seed(2)
fitlist <-mclapply(1:4,function(i){
  Gibbs1PLNO(2000,y)
},mc.cores = detectCores())
a_samples <-lapply(fitlist, function(x)x$ability)
i_samples <-lapply(fitlist, function(x)x$item)

収束をトレースプロットで確かめてみます.

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)

f:id:abrahamcow:20180210230118p:plain

f:id:abrahamcow:20180210230127p:plain

大丈夫そうです.

全部のパラメータのトレースプロットを見るのは大変なので,  \hat R 指標も利用します.

 \hat R 指標についてはたしか,

Rによる計算機統計学

Rによる計算機統計学

に解説があります.

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

最初の1000個のサンプルはバーンインとして捨てます.

ex_a_samples_l <-lapply(a_samples, function(x){x[-c(1:1000),]})
ex_i_samples_l <-lapply(i_samples, function(x){x[-c(1:1000),]})
> all(Rhat(ex_a_samples_l)<1.1)
[1] TRUE
> all(Rhat(ex_i_samples_l)<1.1)
[1] TRUE

シミュレーションで設定したパラメータの真値と推定値(事後期待値)をプロットしてみます.

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

hattheta <-apply(ex_a_samples,2,mean)
hatgamma <-apply(ex_i_samples,2,mean)
ran1 <-range(hattheta,theta)
ran2 <-range(hatgamma,gamma)
plot(theta,hattheta,xlim=ran1,ylim=ran1)
abline(0,1,lty=2)
plot(gamma,hatgamma,xlim=ran2,ylim=ran2)
abline(0,1,lty=2)

f:id:abrahamcow:20180211134132p:plain

f:id:abrahamcow:20180211134148p:plain

おおむね近い値が求まっていそうです.

項目反応理論の人はあまり区間推定には関心がないようですが, 95%信用区間が真のパラメータを含んだ回数を数えると, 名目上の数字と一致することがわかります.

> 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] 0.95

付録

同じモデルを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;
}
model{
  theta ~ normal(0,1);
  for(i in 1:n){
    for(j in 1:m){
     if(y[i,j]==1){
       target += normal_lccdf(0|theta[i]-gamma[j],1);
     }else{
       target += normal_lcdf(0|theta[i]-gamma[j],1);
     }
    }
  }
}