こんにちは。前回・前々回と項目反応理論（Item Response Theory, IRT）のスクラッチ実装に関する記事でしたが、今回もその続きというか発展編です。これまではベイズ推定を用いるのが通例でしたが、やはりベイズは気長に待たなければならないという悲しい点があります。そこで、今回はより高速なEMアルゴリズムを利用して項目反応理論を速く推定してみましょう～！という内容になります。

ここでは、Polson et al. (2013)*1 らによるPolya-Gamma分布を利用したロジットモデル推定のスキームを応用していきます。彼らの論文では、Polya-Gamma分布を用いてロジットモデルでのギブスサンプリングを効率的に実現するということが提案されていました。それを応用することで、ロジットモデルでも効率的にEMアルゴリズムによる推定を行うことができるようになります。ゆっくりと解説していきます。

Polya-Gamma分布による定式化については、株式会社Nospareさんのqiita記事がより詳しいのでそちらも参照してください。 qiita.com

また、今回使用するコードはすべてGithub上にアップしてますので、そちらも併せてご覧ください。

github.com

セットアップ
Polya-Gamma data augmentationを利用したIRTモデルの推定
- EMアルゴリズムへの応用
- パラメトリックブートストラップによる信頼区間の導出
応用例：第106回連邦議会における米国上院議員の投票行動分析
おわりに

セットアップ

まずは前回と同様にIRTモデルのセットアップから行います。被験者 $i=1,...,I$と問題$j=1,...,J$があるとします。この時、$i$が$j$で正解する（$y_{ij}=1$, otherwise $0$）確率を以下で定義します：

\begin{equation} \text{Pr}(y_{ij}=1) = \frac{\exp(\psi_{ij})}{1 + \exp(\psi_{ij})}, \ \ \ \psi_{ij} = \alpha_j + \beta_j \theta_i \end{equation}

ここで、$\alpha_j, \beta_j$はそれぞれ問題$j$の解答困難度・識別度であり、$\theta_i$は被験者$i$の能力です。上記のモデルに対してベイズ推定などを用いることで潜在的なパラメータである$\alpha_j, \beta_j, \theta_i$を推定します。

仮に上記のモデルをベイズ推定する場合、事後分布が上手く定まらないのでMetropolis Hastingsアルゴリズムなどに頼る必要があります。しかし、これはギブスサンプリングに比べると非効率であり、そのためギブスサンプリングが簡単に行える別の二項選択モデルであるプロビットモデルが好まれていました*2。しかし、Polson et al. (2013)が提案したPolya-Gamma data augmentationのスキームを用いれば、プロビットモデルと同様に簡単に効率的なギブスサンプラーを構築できるのです。

Polya-Gamma data augmentationを利用したIRTモデルの推定

ここからは、Goplerud (2019) *3 によるPolya-Gamma分布を用いたIRTの定式化に従っていきます。まず、最初にPolya-Gamma分布の定義を見ていきます。以下の場合、確率変数$X$は$b>0, c\in \mathop{\mathbb{R}}$を持つPolya-Gamma分布に従うとします ($X\sim PG(b, c))$：

\[X = \frac{1}{2\pi^2}\sum\limits_{k=1}^{\infty}\frac{g_k}{(k - 1/2)^2 + c^2/(4\pi^2)}\]

また、$g_k \sim Ga(b, 1)$です。ここでIRT式について、$\psi_{ij} \in \mathop{\mathbb{R}}$で$\omega_{ij} \sim PG(1, 0)$の場合：

\[ \frac{\exp(\psi_{ij})^{y_{ij}}}{1 + \exp(\psi_{ij})}= \frac{1}{2}\exp (k_{ij}\psi_{ij}) \int\limits_{0}^{\infty} \exp(-\omega_{ij}\psi_{ij}^2/2) f(\omega_{ij}) d\omega_{ij}, \ \ \ k_{ij} = y_{ij} - 1/2 \]

として書くことができます。この場合、Polya-Gammaの確率変数$\omega_{ij}$によるdata augmentationを行った後の完全データ尤度関数$f(y_{ij}, \omega_{ij}\mid \alpha_j, \beta_j, \theta_i)$は

\[f(y_{ij}, \omega_{ij}\mid \alpha_j, \beta_j, \theta_i) \propto \exp\left(k_{ij}\psi_{ij} - \frac{\omega_{ij}(\psi_{ij})^2}{2}\right) f(\omega_{ij}\mid 1, 0)\]

となります。パラメータに対して正規分布を事前分布として定めることで、事後分布も正規分布として得ることができるためギブスサンプリングが可能となります。めっちゃ便利ですね。ちなみに、Polya-Gamma data augmentationを使ってIRTのギブスサンプリングをしている例がJiang & Templin (2019) の研究です*4。しかし、今回はEMアルゴリズム実装に興味があるので、ギブスサンプリングの詳細は省きます。

EMアルゴリズムへの応用

次に、Polya-Gamma分布を用いた拡大法のEMアルゴリズムへの適用の仕方を見ていきます。まず、パラメータの事前分布を

\[\alpha_j \sim N(\alpha_0, A_0), \ \ \ \beta_j \sim N(\beta_0, B_0), \ \ \ \theta_i \sim N(0, 1)\]

と設定します。$\xi, \xi^*$をそれぞれ現在のイテレーション・一つ前のイテレーションのパラメータの集合とします。この場合、E-stepにおいてQ関数は以下のように求めます：

E-step

\begin{align} Q(\xi, \xi^*) &\propto \sum\limits_{i=1}^I\sum\limits_{j=1}^J k_{ij}\psi_{ij} - \frac{\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij}, \xi^*]\psi_{ij}^2}{2}\\\nonumber &- \sum\limits_{i=1}^I\frac{\theta_i^2}{2} - \sum\limits_{j=1}^J \frac{(\alpha_j-\alpha_0)^2}{2A_0} - \sum\limits_{j=1}^J\frac{(\beta_j-\beta_0)^2}{2B_0} \end{align}

また、$k_{ij} = y_{ij} - 1/2; \ \ \ \psi_{ij} = \alpha_j + \beta_j \theta_i$であり、$\omega_{ij}$の一つ前のイテレーションにおけるパラメータと観測される$y_{ij}$を条件とした場合の条件付き期待値は：

\[\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij}, \xi^*] = \frac{\tan(\psi_{ij}^*/2)}{2\psi_{ij}^*}\]

と表せます。式中の$tan$はハイパボリックタンジェントのtanhを打ちたかったのですが、はてなリンクが自動に付与されてしまう関係で数式が表示されなくなってしまうので仕方なく$tan$にしてます。タイポではないです、ご了承ください。本題に戻りますが、M-stepでは単純にQ関数を任意のパラメータについて最大化することでパラメータをアップデートするだけです。パラメータは$\theta \rightarrow \alpha \rightarrow \beta$の順にアップデートしていきます：

M-step

その1. $\theta_i$のアップデート

\begin{align} \frac{\partial Q}{\partial \theta_i} &= \sum\limits_{j=1}^J \beta_j^* k_{ij} - \alpha_j^*\beta_j^* \mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] - (\beta_j^*)^2\theta_i\mathop{\mathbb{E}}[\omega_ij\mid y_{ij},\xi^*] - \theta_i = 0\\ \therefore \theta_i &= \left(\sum\limits_{j=1}^J(\beta_j^*)^2\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] + 1\right)^{-1}\left(\sum\limits_{j=1}^J \beta_j^*k_{ij} - \alpha_j^*\beta_j^*\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*]\right) \end{align}

その2. $\alpha_j$のアップデート

\begin{align} \frac{\partial Q}{\partial \alpha_j} &= \sum\limits_{i=1}^Ik_{ij}- \alpha_j\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] - \beta_j^*\theta_i\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] - \frac{\alpha_j - \alpha_0}{A_0}= 0\\ \therefore \alpha_j &= \left( A_0^{-1} + \sum\limits_{i=1}^I\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*]\right)^{-1}\left(\alpha_0A_0^{-1} + \sum\limits_{i=1}^I k_{ij} - \beta_j^*\theta_i\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] \right) \end{align}

その3. $\beta_j$のアップデート

\begin{align} \frac{\partial Q}{\partial \beta_j} &= \sum\limits_{i=1}^I k_{ij}\theta_i - \alpha_j\theta_i\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] - \beta_j \theta_i^2 \mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] - \frac{\beta_j-\beta_0}{B_0}=0\\ \therefore \beta_j &= \left(B_0^{-1} + \sum\limits_{i=1}^I\theta_i^2\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] \right)^{-1}\left(\beta_0B_0^{-1}+\sum\limits_{i=1}^I k_{ij}\theta_i - \alpha_j\theta_i\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*]\right) \end{align}

上記のE-stepとM-stepを収束するまで繰り返します。ここで、収束判断の基準については現在のパラメータと一つ前のイテレーションでのパラメータとの相関係数を用い、$\min(1 - cor(\xi, \xi^*)) < 1e-6$となる場合に収束したと判断します。

パラメトリックブートストラップによる信頼区間の導出

EMアルゴリズムは確かに速いのですが、点推定値しか求めてくれません。そこで、Lewis & Poole (2004) *5 が提案したパラメトリックブートストラップを用いて信頼区間を求めていきます。手順は以下の通りです：

元モデルでのパラメータ推定値によって予測確率を求め、ベルヌーイ分布から$\hat{y}_{ij}$をサンプリングする
$\hat{y}_{ij}$を用いて再度IRTを行い、パラメータを推定し保存。

上記を任意の回数繰り返します。これで得られた推定値を利用して信頼区間を求めます。また、信頼区間を求める際にバイアス補正を行い、具体的には任意のパラメータ$\delta$についての区間推定値を求める場合：

\[ \left[\delta_{\text{boot}}^{\text{lower}} + \hat{\delta} - \bar{\delta}_\text{boot}, \ \ \delta_{\text{boot}}^{\text{upper}} + \hat{\delta} - \bar{\delta}_\text{boot}\right] \]

とします。ここで、$\hat{\delta}$は元モデルで得られた推定値、$\bar{\delta}_\text{boot}$はブートストラップによって得られた推定値の平均、$\delta_{\text{boot}}^{\text{lower}}$, $\delta_{\text{boot}}^{\text{upper}}$はその下側・上側パーセント点を表します。この操作をすることによって、元モデルからの推定値に対応した区間を求めることができます。

応用例：第106回連邦議会における米国上院議員の投票行動分析

今回も、{MCMCpack}パッケージに収録されているアメリカ連邦議会上院の第106議会での点呼式投票データ*6を使って各議員の理想点推定をしていこうと思います。政治学の文脈では、空間投票モデルとIRTの整合性が非常に高いので、例えば政治アクターの投票データ等を用いて彼らの潜在的な選好・イデオロギーを推定する際にIRTは頻繁に用いられています。アメリカ議会から最高裁、国際連合総会など様々な投票が起こる政治アリーナの分析では重宝されています。

アメリカ議会の投票行動を分析すると「共和党 vs. 民主党」のような議員の党派性がくっきりと表れた投票行動が為されています。それらをIRTを使ってリカバーしていこうというのがここでのモチベーションです。

実際のところ、私自身が政治学畑の人間なのでこういった投票行動の分析に関する応用が多いとは思いますが、元を辿ればIRTはテストデータ等の分析に使われるので、当然そちらの方にも適用できます。ただ、政治学の文脈だと、IRTは人々のイデオロギーを測る際に有用であり、目に見えない人の選好を見ようとする試みが気持ち悪くて結構好きです。ともあれ、早速やっていきましょう。

まずは、データを読み込みましょう。

data(Senate, package = "MCMCpack")

データを確認すると、点呼式投票 (rc~列) 以外の議員に関するデータが含まれていることがわかります。

Senate[1:10, 1:10]

              id statecode   state party     member rc1 rc2 rc3 rc4 rc5
SESSIONS   49700        41 ALABAMA     1   SESSIONS   1   0   0   0   1
SHELBY     94659        41 ALABAMA     1     SHELBY   1   0   0   0   1
MURKOWSKI  14907        81  ALASKA     1  MURKOWSKI   1   0   0   0   1
STEVENS    12109        81  ALASKA     1    STEVENS   1   0   0   0   1
KYL        15429        61 ARIZONA     1        KYL   1   0   0   0   1
MCCAIN     15039        61 ARIZONA     1     MCCAIN   1   0   0   0   1
HUTCHINSON 29306        42 ARKANSA     1 HUTCHINSON   1   0   0   0   1
LINCOLN    29305        42 ARKANSA     0    LINCOLN   1   0   0   1   0
BOXER      15011        71 CALIFOR     0      BOXER   1   1   1   1   0
FEINSTEIN  49300        71 CALIFOR     0  FEINSTEIN   1   1   1   1   0

行名は議員の名前になっており、属する州や政党のインデックス (1 = Republican, 0 = Democrat)が入ってます。とはいえ、IRTをするうえでは不要なので全部切り捨てて、またSenateデータはデータフレームなので行列に変換します。

Y <- as.matrix(Senate[, 6:length(Senate)])

EMアルゴリズムをプログラムする

次にEMアルゴリズムを回す関数やパラメータをアップデートする関数を定義します。ここで、推定をより高速にするためにC++を利用していきます。$\texttt{R}には$RcppやRcppArmadilloといった効率的にC++を利用できる環境があるので思う存分使っていきましょう。

$\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij}, \xi^*]$を求める関数

E-stepの箇所で説明した通り、$\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij}, \xi^*]$は一つ前のイテレーションでのパラメータを基にして計算します。式とコードは以下の通りです（sourceCpp(draw_E_Omega.cpp)で読み込み）：

\[\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij}, \xi^*] = \frac{\tan(\psi_{ij}^*/2)}{2\psi_{ij}}\]

//[[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

//[[Rcpp::export]]
arma::mat draw_E_Omega(arma::vec alpha, arma::vec beta, arma::vec theta) {
  // Generate I length row vector filled with 1
  arma::vec intercept(theta.n_rows);
  intercept = intercept.ones();
  // 2 by J item parameter matrix 
  arma::mat beta_tilde = join_horiz(alpha, beta).t();
  // I by 2 individual parameter matrix
  // 1st column is intercept
  arma::mat Theta =  join_horiz(intercept, theta);
  // Calculate psi
  arma::mat psi = Theta * beta_tilde;
  // Calculate E(omega)
  arma::mat Omega = arma::tanh(psi / 2) / (2 * psi);
  return(Omega);
}

$\theta_i$をアップデートする関数

M-stepの箇所で紹介した通りにコードを書いていきます。式とコードは（sourceCpp(draw_theta.cpp) で読み込み）：

\[\theta_i = \left(\sum\limits_{j=1}^J(\beta_j^*)^2\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] + 1\right)^{-1}\left(\sum\limits_{j=1}^J \beta_j^*k_{ij} - \alpha_j^*\beta_j^*\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*]\right)\]

#include <Rcpp.h>
using namespace Rcpp;

//[[Rcpp::export]]
NumericVector draw_theta(NumericMatrix Y, NumericVector alpha, NumericVector beta, 
                         NumericMatrix Omega, int constraint) {
  int I = Y.nrow();
  int J = Y.ncol();
  
NumericVector theta(I);
  for (int i = 0; i < I; i++) {
    double left_part = 0;
    double right_part = 0;
    for (int j = 0; j < J; j++) {
      if (NumericVector::is_na(Y(i, j))) {
        left_part += 0;
        right_part += 0;
      } else {
        double omega_ij = Omega(i, j);
        double k_ij = Y(i, j) - 0.5;
        left_part += std::pow(beta[j], 2.0) * omega_ij;
        right_part += beta[j] * k_ij - alpha[j] * beta[j] * omega_ij;
      }
    }
    left_part = left_part + 1;
    theta[i] = (1 / left_part) * right_part;
  }
  
  if (theta[constraint - 1] < 0) {
    theta = -theta;
  }
  
  return(theta);
}

ここでconstraintという変数を用いることで、自分が指定したユニットの$\theta_i$が常に正の値を取るように設定してます。IRTのようなモデルでは解の符号が反転することがあるのでそれらを防ぐのに効果的かと思われます。

$\alpha_j$をアップデートする関数

$\alpha_j$も同様にしてコードを書いていきます（sourceCpp(draw_alpha.cpp) で読み込み）：

\[ \alpha_j = \left( A_0^{-1} + \sum\limits_{i=1}^I\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*]\right)^{-1}\left(\alpha_0A_0^{-1} + \sum\limits_{i=1}^I k_{ij} - \beta_j^*\theta_i\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] \right)\]

#include <Rcpp.h>
using namespace Rcpp;

//[[Rcpp::export]]
NumericVector draw_alpha(NumericMatrix Y, NumericVector beta, NumericVector theta, 
                         NumericMatrix Omega, double a0, double A0) {
  int I = Y.nrow();
  int J = Y.ncol();
  
  NumericVector alpha(J);
  for (int j = 0; j < J; j++) {
    double left_part = 0;
    double right_part_1 = 0;
    double right_part_2 = 0;
    for (int i = 0; i < I; i++) {
      if (NumericVector::is_na(Y(i, j))) {
        left_part += 0;
        right_part_1 += 0;
        right_part_2 += 0;
      } else {
        double omega_ij = Omega(i, j);
        double k_ij = Y(i, j) - 0.5;
        left_part += omega_ij;
        right_part_1 += k_ij;
        right_part_2 += theta[i] * omega_ij;
      }
    }
    left_part = left_part + 1/A0;
    double right_part = a0/A0 + right_part_1 - beta[j] * right_part_2;
    alpha[j] = (1 / left_part) * right_part;
  }
  
  return(alpha);
}

$\beta_j$をアップデートする関数

最後に$\beta_j$の関数を書きます（sourceCpp(draw_beta.cpp) で読み込み）：

\[\beta_j = \left(B_0^{-1} + \sum\limits_{i=1}^I\theta_i^2\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*] \right)^{-1}\left(\beta_0B_0^{-1}+\sum\limits_{i=1}^I k_{ij}\theta_i - \alpha_j\theta_i\mathop{\mathbb{E}}[\omega_{ij}\mid y_{ij},\xi^*]\right)\]

#include <Rcpp.h>
using namespace Rcpp;

//[[Rcpp::export]]
NumericVector draw_beta(NumericMatrix Y, NumericVector alpha, NumericVector theta, 
                        NumericMatrix Omega, double b0, double B0) {
  int I = Y.nrow();
  int J = Y.ncol();
  
  NumericVector beta(J);
  for (int j = 0; j < J; j++) {
    double left_part = 0;
    double right_part_1 = 0;
    double right_part_2 = 0;
    for (int i = 0; i < I; i++) {
      if (NumericVector::is_na(Y(i, j))) {
        left_part += 0;
        right_part_1 += 0;
        right_part_2 += 0;
      } else {
        double omega_ij = Omega(i, j);
        double k_ij = Y(i, j) - 0.5;
        left_part += std::pow(theta[i], 2.0) * omega_ij;
        right_part_1 += k_ij * theta[i];
        right_part_2 += theta[i] * omega_ij;
      }
    }
    left_part = left_part + 1/B0;
    double right_part = b0/B0 + right_part_1 - alpha[j] * right_part_2;
    beta[j] = (1 / left_part) * right_part;
  }
  
  return(beta);
}

EMアルゴリズムを回す関数

アルゴリズムを回す関数は$\texttt{R}$で定義します。なぜかというと、単純に私の力不足でC++でのコーディングができなかっただけです。コードは以下の通りです：

# Loading wrapper function ====================
pg_emirt <- function(Y, init = NULL, prior = NULL, constraint = NULL, maxit = 500, tol = 1e-6, verbose = 10) {
  
  cat("===========================================================================\n")
  cat("2PL Item Response Model with EM Algorithm and Polya-Gamma Data Augmentation\n")
  cat("===========================================================================\n")
  
  # Measuring implementation times
  stime <- proc.time()[3]
  
  if (is.null(init)) {
    cat("Initial values is not supplied, so computing from the data.....")
    theta_init <- as.numeric(scale(rowMeans(Y, na.rm = TRUE)))
    alpha_init <- beta_init <- rep()
    for (j in 1:ncol(Y)) {
      c <- suppressWarnings(coef(glm(Y[, j] ~ theta_init, family = "binomial")))
      alpha_init[j] <- ifelse(abs(c[1]) > 10, 0.1, c[1])
      beta_init[j] <- ifelse(abs(c[2]) > 10, 0.1, c[2])
    }
    init <- list(theta = theta_init, alpha = alpha_init, beta = beta_init)
    cat("DONE!\n")
  }
  if (is.null(prior)) {
    cat("Prior is not supplied, so automatically setting:\n")
    cat("alpha ~ N(0, 100), beta ~ N(0, 100)\n")
    prior <- list(a0 = 0, A0 = 100, b0 = 0, B0 = 100)
  }
  
  # Get initial values
  theta <- init$theta
  alpha <- init$alpha
  beta <- init$beta
  
  # If theta constraint is not supplied
  if (is.null(constraint)) constraint <- which.max(theta)
  if (theta[constraint] < 0) theta <- -theta
  
  # Get priors
  a0 <- prior$a0
  A0 <- prior$A0
  b0 <- prior$b0
  B0 <- prior$B0
  
  # Convergence check statistics store (correlation between current and previous parameters)
  all_cor <- matrix(NA, maxit, 3)
  colnames(all_cor) <- c("alpha", "beta", "theta")
  
  # EM
  cat("Start EM:\n") 
  for (g in 1:maxit) {
    ## Previous parameters
    theta_old <- theta
    alpha_old <- alpha
    beta_old <- beta
    
    ## E-step (calculate E[omega])
    Omega <- draw_E_Omega(alpha_old, beta_old, theta_old)
    
    ## M-step (theta -> alpha -> beta)
    theta <- draw_theta(Y, alpha_old, beta_old, Omega, constraint)
    alpha <- draw_alpha(Y, beta_old, theta, Omega, a0, A0)
    beta <- draw_beta(Y, alpha, theta, Omega, b0, B0)
    
    ## Convergence stat
    all_cor[g, 1] <- cor(alpha_old, alpha)
    all_cor[g, 2] <- cor(beta_old, beta)
    all_cor[g, 3] <- cor(theta_old, theta)
    
    ## If converged
    if (min(1 - all_cor[g, ]) < tol) {
      ### Conv stat
      all_cor <- na.omit(all_cor)
      attr(all_cor, "na.action") <- NULL
      
      ### Parameters
      names(theta) <- rownames(Y)
      parameters <- list(alpha = alpha, beta = beta, theta = theta)
      
      ### Organize
      L <- list(parameters = parameters,
                converge = TRUE,
                conv_check = all_cor,
                iter = g,
                data = list(Y = Y,
                            init = init,
                            prior = prior,
                            constraint = constraint,
                            maxit = maxit,
                            tol = tol))
      
      ### Record time
      etime <- proc.time()[3]
      cat("Model converged at iteration", g, ": total time", round(etime - stime, 3), "sec\n")
      return(L)
    }
    
    ## If failed to converge
    if (g == maxit) {
      ### Parameters
      names(theta) <- rownames(Y)
      parameters <- list(alpha = alpha, beta = beta, theta = theta)
      
      ### Organize
      L <- list(parameters = parameters,
                converge = FALSE,
                conv_check = all_cor,
                iter = g,
                data = list(Y = Y,
                            init = init,
                            prior = prior,
                            constraint = constraint,
                            maxit = maxit,
                            tol = tol))
      
      ### Record time
      etime <- proc.time()[3]
      cat("Model failed to converge", ": total time", round(etime - stime, 3), "sec\n")
      warning("The result may not be reliable; return the result as it is.")
      return(L)
    }
    
    ### Print the status
    if (g %% verbose == 0) {
      mtime <- proc.time()[3]
      cat("Iteration", g, "eval =", min(1 - all_cor[g, ]), ": time", round(mtime - stime, 3), "sec\n")
    }
  }
}

引数は以下の通りです：

Y: 解答（投票）行列
init: 初期値
prior: 事前分布
constraint: 常に$\theta_i$が正の値を取るユニットの番号
maxit: 最大イテレーション数（デフォルトは500）
tol: 収束基準（デフォルトは1e-6）
verbose: 何イテレーション毎にメッセージを出力するか（デフォルトは10）

パラメトリックブートストラップを行う関数

ブートストラップの関数も$\texttt{R}$で書きました：

# Loading parametric bootstrap function ===================================
pg_em_boot <- function(model, boot = 100, verbose = 10) {
  
  cat("===========================================\n")
  cat("Parametric bootstrap for Polya-Gamma EM IRT\n")
  cat("===========================================\n")
  
  stime <- proc.time()[3]
  
  # For suppressing the messages
  quiet <- function(x) { 
    sink(tempfile()) 
    on.exit(sink()) 
    invisible(force(x)) 
  } 
  
  # Extract the estimated parameters
  alpha <- model$parameters$alpha
  beta <- model$parameters$beta
  theta <- model$parameters$theta
  
  # Storage
  theta_boot <- matrix(NA, length(theta), boot)
  rownames(theta_boot) <- names(theta)
  alpha_boot <- beta_boot <- matrix(NA, length(alpha), boot)
  
  # Bootstrap
  for (b in 1:boot) {
    # Generating random binomial variates (prediction of Y)
    set.seed(b)
    p <- plogis(cbind(1, theta) %*% rbind(alpha, beta))
    y_pred <- matrix(rbinom(length(theta) * length(alpha), 1, p),
                     length(theta), length(alpha))
    
    # Estimation
    fit_boot <- quiet(pg_emirt(Y = y_pred,
                               init = model$data$init,
                               prior = model$data$prior,
                               constraint = model$data$constraint,
                               maxit = model$data$maxit,
                               tol = model$data$tol,
                               verbose = model$data$maxit + 1)) 
    
    # Saving
    theta_boot[, b] <- fit_boot$parameters$theta
    alpha_boot[, b] <- fit_boot$parameters$alpha
    beta_boot[, b] <- fit$parameters$beta
    
    if (b %% verbose == 0) {
      mtime <- proc.time()[3]
      cat("Bootstrap", b, "done : time", round(mtime - stime, 3), "sec\n")
    }
  }
  
  # Return the results
  etime <- proc.time()[3]
  L <- list(alpha = alpha_boot, beta = beta_boot, theta = theta_boot)
  cat("Bootstrap finished : total time", round(etime - stime, 3), "sec\n")
  return(L)
}

初期値・事前分布の設定

実装するにあたり、初期値と事前分布を設定します。まずは、初期値の設定手順は以下の通りです：

投票行列に対して行ごとの平均値を計算し標準化することで$\theta_i$の初期値を得る
$\theta_i$の初期値を投票行列の各項目に対する投票（0 or 1）に回帰させる（ロジット回帰をJ回行う）
切片の推定値を$\alpha_j$の初期値、また$\theta_i$初期値の係数値を$\beta_j$の初期値とする

コードは以下の通りです：

## Generate initial values, constraint and priors ===========
# Calculate the standardized value of the probability of yea response 
# as initial value of theta
theta_init <- as.numeric(scale(rowMeans(Y, na.rm = TRUE)))

# With initial value of theta, get initial values of alpha and beta
# by regressing Y on theta
alpha_init <- beta_init <- rep()
for (j in 1:ncol(Y)) {
  c <- coef(glm(Y[, j] ~ theta_init, family = "binomial"))
  alpha_init[j] <- ifelse(abs(c[1]) > 10, 0.1, c[1])
  beta_init[j] <- ifelse(abs(c[2]) > 10, 0.1, c[2])
}

# Aggregate
init <- list(theta = theta_init,
             alpha = alpha_init,
             beta = beta_init)

また、constraintはMCMCpack::MCMCirt1dの関数説明で紹介されている通りに従い、共和党のHELMSに定め、仮にHELMSの$\theta_i$の初期値が負の値を取っている場合には初期値全体を反転させます：

# Set constraint such that theta of HELMS always takes a positive value
constraint <- which(rownames(Y) == "HELMS")
if (theta_init[constraint] < 0) theta_init <- -theta_init

次に、事前分布を以下の通りに設定します：

\[\alpha_j \sim N(0, 25), \ \ \ \beta_j \sim N(0, 25), \ \ \ \theta_i \sim N(0, 1)\]

この場合、$\alpha_0=0, A_0=25, \beta_0 = 0, B_0 = 25$となります。よって：

# Priors
prior <- list(a0 = 0, A0 = 25, b0 = 0, B0 = 25)

これで準備が整いました！

実装

それでは実装していきます：

## Run PG-EM IRT =============================================
model <- pg_emirt(Y = Y, init = init, prior = prior, constraint = constraint,
                  maxit = 500, tol = 1e-6, verbose = 10)

===========================================================================
2PL Item Response Model with EM Algorithm and Polya-Gamma Data Augmentation
===========================================================================
Start EM:
Iteration 10 eval = 5.770342e-05 : time 0.093 sec
Iteration 20 eval = 9.77759e-06 : time 0.169 sec
Iteration 30 eval = 1.749291e-06 : time 0.199 sec
Model converged at iteration 35 : total time 0.215 sec

鬼高速ですねｗ推定に0.2秒ほどしかかかりませんでした。次にブートストラップを100回ほどしていきます：

## Boostrap ==================================================
boot <- pg_em_boot(model)

===========================================
Parametric bootstrap for Polya-Gamma EM IRT
===========================================
Bootstrap 10 done : time 1.36 sec
Bootstrap 20 done : time 2.455 sec
Bootstrap 30 done : time 3.541 sec
Bootstrap 40 done : time 4.685 sec
Bootstrap 50 done : time 5.778 sec
Bootstrap 60 done : time 6.804 sec
Bootstrap 70 done : time 7.851 sec
Bootstrap 80 done : time 8.937 sec
Bootstrap 90 done : time 9.979 sec
Bootstrap 100 done : time 11.138 sec
Bootstrap finished : total time 11.139 sec

こちらも11秒ほどで終了しました。計11.3秒とかなり速いです。前回の記事でMCMCpack::MCMCirt1dを用いて同様のデータで分析した際にburninが5000、mcmcが10000で設定したら大体55秒くらいかかってましたが、それと比べると中々の快挙ですね。とはいえ、EMは収束したらアルゴリズムを勝手に止めるという感じなので、そりゃ速いやろって感じですが。

結果

しっかり推定できているか確認します。信頼区間を導出する際には、バイアス補正をした形で行っていきます：

## Result ====================================================
# Bias correction
## Calculate bootstrapm mean
boot_mean <- rowMeans(boot$theta)

## Bias (estimated mean - boot mean)
bias <- model$parameters$theta - boot_mean

## 95% CI
lwr <- apply(boot$theta, 1, quantile, probs = .025) + bias
upr <- apply(boot$theta, 1, quantile, probs = .975) + bias

## Plot
tibble(name = Senate$member, 
       party = Senate$party, 
       theta = model$parameters$theta,
       lwr = lwr,
       upr = upr) %>% 
  mutate(party = if_else(party == 1, "Republican", "Democrat") %>% 
           factor(levels = c("Republican", "Democrat"))) %>% 
  ggplot(aes(x = theta, y = reorder(name, theta), 
             color = party)) +
  geom_pointrange(aes(xmin = lwr, xmax = upr), size = 0.4) +
  xlab("Ideal Point") + 
  theme_light() +
  theme(legend.title = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_text(size = 3),
        legend.position = "bottom",
        legend.direction = "horizontal")

赤は共和党で青は民主党を示します。党派性がくっきりとでた図になっていると思います。今までのベイズを使った推定値とほぼ同じ結果になっていると思います。前述の通り、試験問題解答データを使っても上手く推定できるかと思いますので、ぜひ活用してください。

おわりに

今回はPolya-Gamma data augmentationのスキームを用いたEMアルゴリズムによる項目反応理論モデルの推定を紹介しました。IRTモデルにはベイズ推定が用いられることが殆どですが、このようにEMアルゴリズムを用いた方が圧倒的に速いので割とオススメします。紹介した論文により詳細に手続きが記されているので興味がある人は見てみてください。私自身まだまだ勉強中ですので間違い等があれば指摘してください。ここまでご覧いただきありがとうございました。

再度にはなりますが、今回使用したコードはGithub上にまとめていますのでそちらもチェックしてみてください。

github.com

*1:Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Pólya–Gamma latent variables. Journal of the American statistical Association, 108(504), 1339-1349.

*2:項目反応理論を用いた政治学研究でもプロビットモデルが好まれています。

*3:Goplerud, M. (2019). A Multinomial Framework for Ideal Point Estimation. Political Analysis, 27(1), 69-89.

*4:Jiang, Z., & Templin, J. (2019). Gibbs samplers for logistic item response models via the Pólya–Gamma distribution: A computationally efficient data-augmentation strategy. psychometrika, 84(2), 358-374.

*5:Lewis, J. B., & Poole, K. T. (2004). Measuring bias and uncertainty in ideal point estimates via the parametric bootstrap. Political Analysis, 12(2), 105-127.

*6:詳しいことはこちら https://cran.r-project.org/web/packages/MCMCpack/MCMCpack.pdf

砂になった人

砂になった人の趣味や勉強についての備忘録

【R/Rcpp】EMアルゴリズムで項目反応理論を実装するー Polya-Gamma Data Augmentationの応用