【R/Rcpp】項目反応理論 (IRT) をスクラッチで実装する
タイトル通りです。IRTをスクラッチでMCMC実装しました。検索をかけてもあまり実例が出てこなかったため、自分のメモ帳替わりということも含めて記事にします。
本記事では、Wim J. van der Linden (2016) の``Handbook of Item Response Theory Volume 2: Statistical Tools" 15章を参考にしていますので、より詳しく知りたい方はそちらをご覧ください。
また、当方も勉強中のため間違っているところが多々あると思うので、その都度お知らせいただけると幸いです。
準備
ざっくり書いていきます。個人 $i = 1,...,I$ が 問題 $j = 1,...,J$ にて正解する ($y_{ij} = 1$; otherwise 0) 確率を以下で与えます:
\begin{align*}
p_{ij} = \text{Pr}(y_{ij}=1\ |\ \alpha_j, \beta_j, \theta_i) = \text{logit}^{-1}(\beta_j \theta_i - \alpha_j).
\end{align*}
ここで、$\alpha_j$, $\beta_j$はそれぞれ問題$j$の難易度 (difficulty) と識別度 (discrimination)、$\theta_i$は個人$i$の能力 (ability) で、$\text{logit}^{-1}$は逆ロジット関数です。Likelihood $f(\cdot\ |\ \cdot)$は
\begin{align}
f(y\ |\ \alpha, \beta, \theta) &= \prod\limits_{i=1}^{I}\prod\limits_{j=1}^J \left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right]\\
\end{align}
です。次に、各パラメータの事後分布を考えていきます。事前分布として、以下を与えます:
\begin{align*}
\alpha_j &\sim \mathcal{N}(a_0, A_0),\\
\beta_j &\sim \mathcal{N}(b_0, B_0), \\
\theta_i &\sim \mathcal{N}(0, 1).
\end{align*}
ベイズの定理より、その他のパラメータを所与としたうえでの各パラメータのconditionalな事後分布は、それぞれ以下で与えられます:
\begin{align*}
\pi(\alpha_j\ |\ y, \beta_j, \theta_i) &\propto f(y | \alpha, \beta, \theta)\pi(\alpha)\\
&\propto \prod\limits_{i=1}^{I} \left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right] \cdot
N(\alpha_j \ |\ a_0, A_0) \forall j = 1,\dots J,\\
\pi(\beta_j\ |\ y, \alpha_j, \theta_i) &\propto f(y | \alpha, \beta, \theta)\pi(\beta)\\
&\propto \prod\limits_{i=1}^{I}\left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right] \cdot
N(\beta_j \ |\ b_0, B_0) \forall j=1,...J,\\
\pi(\theta_i\ |\ y, \alpha_j, \beta_j) &\propto f(y | \alpha, \beta, \theta)\pi(\theta)\\
&\propto \prod\limits_{j=1}^J \left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right] \cdot
N(\theta_i \ |\ 0, 1) \forall i=1,...,I.
\end{align*}
超ざっくり要点だけさらっていきましたが、こんな感じです。あとはサンプラーを作るだけです。
MCMCする (Metropolis-Hastings Algorithm)
通常のギブスサンプリングだと上記の目標分布からのサンプリングは困難ということで、Metropolis-Hastings Algorithmを用います。
ここで、ノテーションの簡略化のために任意のパラメータを$\delta_k$と定めます ($\delta_{1} = \alpha$, $\delta_2 = \beta$, $\delta_3 = \theta$)。まず、各パラメータの候補$\delta_k^*$を以下のようなランダムウォークな分布からサンプリングします。
\[\delta_k^* \sim N(\delta_k^{(t-1)}, \kappa_\delta)\]
$\delta_k^{(t-1)}$は1イテレーション前の$\delta_k$の値であり、$\kappa_\delta$はMH法でいうところのチューニングパラメータです。また、採択確率 (acceptance probability) は以下で計算します:
\[
ap = \min\left\{\frac{f\left(\delta_k^{*}\ | \ y, \delta_{-k}^{(t-1)}\right) g\left(\delta_k^{(t-1)}\ |\ \delta_k^{*} \right)}{f\left(\delta_k^{(t-1)}\ | \ y, \delta_{-k}^{(t-1)}\right) g\left(\delta_k^*\ |\ \delta_k^{(t-1)} \right)}, 1\right\}
\]
$g(\cdot |\ \cdot)$は提案分布です。また、とあるrandom number $u$を$U(0, 1)$からサンプリングして、$u < ap$ならばパラメータのサンプル候補である$\delta_k^{*}$を採択します。そうでない場合は、$\delta_k^{(t-1)}$をサンプルとして採択します。
しかし、コンピューター的な問題で$ap$等の確率の値がPCが認識できないレベルまで小さくなる場合があります (例えば、1000個の観測とそれぞれの確率が0.001の場合、つまり0.0011000の場合を想像してみてください)。そこで、上記の採択確率を対数に変換します。すなわち
\begin{align*}
\log ap = \min&\{\log f\left(\delta_k^{*}\ | \ y, \delta_{-k}^{(t-1)}\right) + \log g\left(\delta_k^{(t-1)}\ |\ \delta_k^{*} \right) \\
&- \left[\log f\left(\delta_k^{(t-1)}\ | \ y, \delta_{-k}^{(t-1)}\right) + \log\left(\delta_k^*\ |\ \delta_k^{(t-1)} \right)\right], 0\}
\end{align*}
であり、$\log(u) < \log(ap)$であれば$\delta_k^*$を採択します。そうでない場合は、$\delta_k^{(t-1)}$をサンプルとして採択します。これらのステップを任意のイテレーション数まで繰り返します。
サンプラーを作る
以下でサンプリングのための関数を定義していきます。また、高速化のためにC++で一部関数を書きました。
対数尤度を計算する関数
この関数では対数尤度を計算します。各パラメータのサンプリングの際に対数尤度を用いるので補助的に定義します。また、forループを使用するのでC++で書きました。ファイル名はlog_likelihood.cppです。
//log_likelihood.cpp #include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] NumericMatrix log_likelihood(NumericMatrix Y, NumericVector alpha, NumericVector beta, NumericVector theta) { int I = theta.length(); int J = alpha.length(); NumericMatrix prob(I, J); NumericMatrix log_p(I, J); for (int j = 0; j < J; j++) { for (int i = 0; i < I; i++) { prob(i, j) = std::exp(beta[j] * theta[i] - alpha[j]) / (1 + std::exp(beta[j] * theta[i] - alpha[j])); log_p(i, j) = Y(i, j) * std::log(prob(i, j)) + (1 - Y(i, j)) * std::log((1 - prob(i, j))); } } return(log_p); }
また、この関数を利用するために{Rcpp}パッケージを読み込んで、コンパイルしましょう。そうすればlog_likelihood関数として利用できます。
library(Rcpp) sourceCpp("log_likelihood.cpp")
$\alpha$をサンプリングする関数
alpha_samplingという関数を定義します。
alpha_sampling <- function(Y, old_data, a0, A0) { # sum function with na.rm=TRUE sum_ <- function(x) sum(x, na.rm = TRUE) # Get old parameter alpha_old <- old_data$alpha MH_alpha <- old_data$MH$alpha J <- length(alpha_old) # Generate candidate parameter: alpha_can ~ N(alpha_old, tuning_param_alpha) alpha_star <- suppressWarnings(rnorm(J, alpha_old, MH_alpha)) # Calculate acceptance probability ## Log posterior density ### Candidate: log(f(Y | theta_old, alpha_star, beta_old)) + log(pi(N(alpha_star | a0, A0))) log_cc_star <- apply(log_likelihood(Y, alpha = alpha_star, beta = old_data$beta, theta = old_data$theta), 2, sum_) + log(dnorm(alpha_star, a0, A0)) ### Old: log(f(Y | theta_old, alpha_old, beta_old)) + log(pi(N(alpha_old | a0, A0))) log_cc_old <- apply(log_likelihood(Y, alpha_old, old_data$beta, old_data$theta), 2, sum_) + log(dnorm(alpha_old, a0, A0)) ## Proposal density ### Candidate: g_m(alpha_star | alpha_old) log_prop_star <- log(dnorm(alpha_star, alpha_old, MH_alpha)) ### Old: # g_m(alpha_old | alpha_star) log_prop_old <- log(dnorm(alpha_old, alpha_star, MH_alpha)) ## Acceptance probability log_ac <- pmin(log_cc_star + log_prop_old - log_cc_old - log_prop_star, 0) ## Acceptance indicator - if log(u) < log(acceptance prob), then return 1 (otherwise 0) acc_new <- ifelse(log(runif(J)) < log_ac, 1, 0) # Save current_data <- old_data ## If acceptance indicator == 1, accept alpha_star current_data$alpha <- ifelse(acc_new == 1, alpha_star, alpha_old) ## Calculate acceptance rate current_data$acc$alpha_acc <- current_data$acc$alpha_acc + acc_new return(current_data) }
$\beta$をサンプリングする関数
同様に、beta_samplingという関数を定義します。
beta_sampling <- function(Y, old_data, b0, B0) { sum_ <- function(x) sum(x, na.rm = TRUE) # Get old parameter beta_old <- old_data$beta MH_beta <- old_data$MH$beta J <- length(beta_old) # Generate candidate parameter: beta_can ~ N(ln beta_old, tuning_param_beta) beta_star <- suppressWarnings(rnorm(J, beta_old, MH_beta)) # Calculate acceptance probability ## Log posterior density ### Candidate: log(f(Y | theta_old, alpha_old, beta_star)) + log(pi(N(beta_star | b0, B0))) log_cc_star <- apply(log_likelihood(Y, alpha = old_data$alpha, beta = beta_star, theta = old_data$theta), 2, sum_) + log(dnorm(beta_star, b0, B0)) ### Old: log(f(Y | theta_old, alpha_old, beta_old)) + log(pi(N(beta_old | b0, B0))) log_cc_old <- apply(log_likelihood(Y, alpha = old_data$alpha, beta = beta_old, theta = old_data$theta), 2, sum_) + log(dnorm(beta_old, b0, B0)) ## Proposal density ### Candidate: g_m(beta_star | beta_old) log_prop_star <- log(dnorm(beta_star, beta_old, MH_beta)) ### Old: # g_m(beta_old | beta_star) log_prop_old <- log(dnorm(beta_old, beta_star, MH_beta)) ## Acceptance probability log_ac <- pmin(log_cc_star + log_prop_old - log_cc_old - log_prop_star, 0) ## Acceptance indicator - if log(u) < log(acceptance prob), then return 1 (otherwise 0) acc_new <- ifelse(log(runif(J)) < log_ac, 1, 0) # Save current_data <- old_data ## If acceptance indicator == 1, accept beta_star current_data$beta <- ifelse(acc_new == 1, beta_star, beta_old) ## Calculate acceptance rate current_data$acc$beta_acc <- old_data$acc$beta_acc + acc_new return(current_data) }
$\theta$をサンプリングする関数
theta_samplingという関数を定義します。パラメータのサンプリング関数はそれぞれ事前分布の部分を書き換えただけで、そこまで難しくはないと思います。
theta_sampling <- function(Y, old_data) { sum_ <- function(x) sum(x, na.rm = TRUE) # Get old parameter theta_old <- old_data$theta MH_theta <- old_data$MH$theta I <- length(theta_old) # Generate candidate parameter: theta_can ~ N(theta_old, tuning_param_theta) theta_star <- suppressWarnings(rnorm(I, theta_old, MH_theta)) # Calculate acceptance probability ## Log posterior density ### Candidate: log(f(Y | theta_star, alpha_old, beta_old)) + log(pi(N(theta_star | 0, 1))) log_cc_star <- apply(log_likelihood(Y, alpha = old_data$alpha, beta = old_data$beta, theta = theta_star), 1, sum_) + log(dnorm(theta_star, 0, 1)) ### Old: log(f(Y | theta_old, alpha_old, beta_old)) + log(pi(N(theta_old | 0, 1))) log_cc_old <- apply(log_likelihood(Y, alpha = old_data$alpha, beta = old_data$beta, theta = theta_old), 1, sum_) + log(dnorm(theta_old, 0, 1)) ## Proposal density ### Candidate: g_m(theta_star | theta_old) log_prop_star <- log(dnorm(theta_star, mean = theta_old, sd = MH_theta)) ### Old: # g_m(theta_old | theta_star) log_prop_old <- log(dnorm(theta_old, mean = theta_star, sd = MH_theta)) ## Acceptance probability log_ac <- pmin(log_cc_star + log_prop_old - log_cc_old - log_prop_star, 0) ## Acceptance indicator - if log(u) < log(acceptance prob), then return 1 (otherwise 0) acc_new <- ifelse(log(runif(I)) < log_ac, 1, 0) # Save current_data <- old_data ## If acceptance indicator == 1, accept theta_star current_data$theta <- ifelse(acc_new == 1, theta_star, theta_old) ## Calculate acceptance rate current_data$acc$theta_acc <- old_data$acc$theta_acc + acc_new return(current_data) }
MHサンプラー
最後にMH法によるサンプラーMH_irtを定義します。R映えというかサンプリング状況を見やすくするために、イテレーション数の出力など色々加えてます。また、サンプリングの順番は$\theta\rightarrow\alpha\rightarrow\beta$です。
MH_irt <- function(datamatrix, alpha_init = NULL, beta_init = NULL, theta_init = NULL, MH_alpha = 0.4, MH_beta = 0.4, MH_theta = 0.4, iter = 2000, warmup = 1000, thin = 1, refresh = 100, seed = NULL, a0 = 1, A0 = 10, b0 = 0, B0 = 10) { # Rename matrix as Y for convenience Y <- datamatrix # Extract the # of individuals and items I <- nrow(Y) J <- ncol(Y) cat("\n================================================================\n") cat("Run Metropolis-Hasting Sampler for 2PL Item response model...\n\n") cat(" Observations:", I * J,"\n") cat(" Number of individuals:", I, "\n") cat(" Number of items:", J, "\n") cat(" Total correct response:", sum(as.numeric(Y), na.rm = TRUE), "/", I * J, "[", round(sum(as.numeric(Y), na.rm = TRUE) / (I * J), 2) * 100, "%]","\n\n") cat(" Priors: \n") cat(" alpha ~", paste0("N(", a0, ", ", A0, "),"), "beta ~", paste0("N(", b0, ", ", B0, "),"), "theta ~ N(0, 1).\n") cat("================================================================\n\n") # Preparation ## Measure starting time stime <- proc.time()[3] ## Set seed if (is.null(seed)) { seed <- as.integer(runif(1, 1, 50000)) } set.seed(seed) ## Define the number of iterations total_iter <- iter + warmup ## Generate storage alpha_store <- matrix(NA, J, iter/thin) beta_store <- matrix(NA, J, iter/thin) theta_store <- matrix(NA, I, iter/thin) ## Define initial number if (is.null(alpha_init)) { alpha_init <- rep(0, J) } if (is.null(beta_init)) { beta_init <- rep(1, J) } if (is.null(theta_init)) { theta_init <- rep(0.1, I) } ## Generate data current_data <- list(alpha = alpha_init, beta = beta_init, theta = theta_init, MH = list(alpha = MH_alpha, beta = MH_beta, theta = MH_theta), acc = list(alpha_acc = 0, beta_acc = 0, theta_acc = 0)) ## Generate priors priors <- list(a0 = a0, A0 = A0, b0 = b0, B0 = B0) # SAMPLING PHASE for (ii in 1:total_iter) { # Sampling theta parameter current_data <- theta_sampling(Y, current_data) if (all(is.na(current_data$theta))) { cat("Error at", ii, ": All values are NAs.\n") break } # Sampling alpha parameter current_data <- alpha_sampling(Y, current_data, priors$a0, priors$A0) if (all(is.na(current_data$alpha))) { cat("Error at", ii, ": All values are NAs.\n") break } # Sampling beta parameter current_data <- beta_sampling(Y, current_data, priors$b0, priors$B0) if (all(is.na(current_data$beta))) { cat("Error at", ii, ": All values are NAs.\n") break } if (ii <= warmup & ii %% refresh == 0) { cat("Warmup:", ii, "/", total_iter, "[", round(ii/total_iter, 2) * 100, "%]", round(proc.time()[3] - stime, 1), "sec\n") } if (ii > warmup) { if ((ii - warmup) %% thin == 0) { iii <- (ii - warmup) / thin alpha_store[, iii] <- current_data$alpha beta_store[, iii] <- current_data$beta theta_store[, iii] <- current_data$theta } if (ii %% refresh == 0) { cat("Sampling:", ii, "/", total_iter, "[", round(ii/total_iter, 2) * 100, "%]", round(proc.time()[3] - stime, 1), "sec\n") } } } # Calculate acceptance rate alpha_acc <- current_data$acc$alpha_acc / total_iter beta_acc <- current_data$acc$beta_acc / total_iter theta_acc <- current_data$acc$theta_acc / total_iter # Generate variable labels label_iter <- paste0("iter_", 1:(iter/thin)) alpha_lab <- paste0("alpha_", 1:J) beta_lab <- paste0("beta_", 1:J) theta_lab <- paste0("theta_", 1:I) colnames(alpha_store) <- colnames(beta_store) <- colnames(theta_store) <- label_iter rownames(alpha_store) <- alpha_lab rownames(beta_store) <- beta_lab rownames(theta_store) <- theta_lab # Redefine quantile function lwr <- function(x) quantile(x, probs = 0.025, na.rm = TRUE) upr <- function(x) quantile(x, probs = 0.975, na.rm = TRUE) mean_ <- function(x) mean(x, na.rm = TRUE) median_ <- function(x) median(x, na.rm = TRUE) # Calculate statistics alpha_post <- data.frame(variable = alpha_lab, mean = apply(alpha_store, 1, mean_), median = apply(alpha_store, 1, median_), lwr = apply(alpha_store, 1, lwr), upr = apply(alpha_store, 1, upr)) beta_post <- data.frame(variable = beta_lab, mean = apply(beta_store, 1, mean_), median = apply(beta_store, 1, median_), lwr = apply(beta_store, 1, lwr), upr = apply(beta_store, 1, upr)) theta_post <- data.frame(variable = theta_lab, mean = apply(theta_store, 1, mean_), median = apply(theta_store, 1, median_), lwr = apply(theta_store, 1, lwr), upr = apply(theta_store, 1, upr)) # Aggregate result <- list(summary = list(alpha = alpha_post, beta = beta_post, theta = theta_post), sample = list(alpha = alpha_store, beta = beta_store, theta = theta_store), ac_rate = list(alpha = alpha_acc, beta = beta_acc, theta = theta_acc)) etime <- proc.time()[3] cat(crayon::yellow("Done: Total time", round(etime - stime, 1), "sec\n")) return(result) }
さあ、準備は整いました。実際に動かしてみましょう。
応用：アメリカ最高裁判事の理想点推定
政治学では、議会における点呼式投票や裁判所における審判記録などを用いて各アクターの理想点 (ideal point) を推定します。理想点とは簡単に言えば「イデオロギー」のようなものであり、例えばリベラル・保守のようなものです。IRTは理想点推定に頻繁に用いられるテクニックであり、Clinton, Jackman and Rivers (2004)*1で詳しく説明されています。
一次元の政策空間$R$上に政治アクター$i=1,...,I$がいるとして、点呼式投票 or 判決 $j=1,...,J$について賛成票 (Yea) を投じる場合は$y_{ij}=1$となり、反対票 (Nay) を投じる場合は$y_{ij}=0$となります。ここで、棄権や欠席は欠損値とします*2。この時、賛成票を投じる確率は先ほどのIRTの式で表されます*3
\begin{align*}
\text{Pr}(y_{ij}=1\ |\ \alpha_j, \beta_j, \theta_i) = \text{logit}^{-1}(\beta_j \theta_i - \alpha_j).
\end{align*}
ここで、$\alpha_j, \beta_j$はそれぞれ政策争点のcutpoint (cutpointを超える$\theta$を持つ$i$は投票する)、分極度 (polarity, saliency) として捉えられます。また、$\theta_i$が$i$の理想点です。
先ほども述べましたが、アメリカ最高裁判事の理想点が推定されている研究は多々あります。どの判事がリベラルで保守なのか等様々に議論されています。また、皆さん大好き{MCMCpack}の作者であるMartinさんとQuinnさんもその一例です。(Martin & Quinn 2002)*4。
嬉しいことに、{MCMCpack}には2000年期の全会一致の判決以外の判決43個と判事9人を記録したデータが提供されています。今回はそれを利用します。
data(SupremeCourt, package = "MCMCpack") Y <- t(SupremeCourt) Y[1:9, 1:5]
1 2 3 4 5 Rehnquist 0 0 0 0 1 Stevens 1 1 1 0 1 O Connor 1 0 0 0 0 Scalia 0 0 0 0 0 Kennedy 1 0 0 0 1 Souter 1 1 1 0 0 Thomas 0 0 0 0 0 Ginsburg 1 1 1 0 0 Breyer 1 1 1 1 0
上記のように、点呼式投票データや判決データは教育学や心理学でいうテストデータと同じ構造になっていることがわかります。それでは早速MCMCしましょう。今回は$\alpha_j$の事前分布に$N(0, 10)$、また$\beta_j$の事前分布に$N(1, 0.2)$を与えます。
# Initial values alpha_init <- rep(0.1, ncol(Y)) beta_init <- rep(0.1, ncol(Y)) theta_init <- rep(0.1, nrow(Y)) # Set anchors theta_init[which(rownames(Y) == "Ginsburg")] <- -2 # more liberal theta_init[which(rownames(Y) == "Scalia")] <- 2 # more conservative # RUN fit <- MH_irt(datamatrix = Y, # roll-call matrix alpha_init = alpha_init, # init for alpha beta_init = beta_init, # init for beta theta_init = theta_init, # init for theta MH_alpha = 2, # tuning param for alpha MH_beta = 0.6, # tuning param for beta MH_theta = 0.9, # tuning param for theta a0 = 0, # prior mean for alpha A0 = 10, # prior sd for alpha b0 = 1, # prior mean for beta B0 = 0.2, # prior sd for beta iter = 10000, # num of sampling warmup = 5000, # num of warmup thin = 10, # thining refresh = 1000, # Output the status of sampling every 1000 steps seed = 1) # seed
================================================================ Run Metropolis-Hasting Sampler for 2PL Item response model... Observations: 387 Number of individuals: 9 Number of items: 43 Total correct response: 190 / 387 [ 49 %] Priors: alpha ~ N(0, 10), beta ~ N(1, 0.2), theta ~ N(0, 1). ================================================================ Warmup: 1000 / 15000 [ 7 %] 2.4 sec Warmup: 2000 / 15000 [ 13 %] 4.9 sec Warmup: 3000 / 15000 [ 20 %] 7.5 sec Warmup: 4000 / 15000 [ 27 %] 10.3 sec Warmup: 5000 / 15000 [ 33 %] 12.5 sec Sampling: 6000 / 15000 [ 40 %] 14.6 sec Sampling: 7000 / 15000 [ 47 %] 16.7 sec Sampling: 8000 / 15000 [ 53 %] 18.7 sec Sampling: 9000 / 15000 [ 60 %] 20.8 sec Sampling: 10000 / 15000 [ 67 %] 22.3 sec Sampling: 11000 / 15000 [ 73 %] 23.5 sec Sampling: 12000 / 15000 [ 80 %] 24.8 sec Sampling: 13000 / 15000 [ 87 %] 26 sec Sampling: 14000 / 15000 [ 93 %] 27.2 sec Sampling: 15000 / 15000 [ 100 %] 28.5 sec Done: Total time 28.5 sec
終わりました。収束診断しますが、パラメータが多いので$\theta_i$だけ見てみます。
# Extract estimated thetas theta_iter <- fit$sample$theta %>% t() %>% as_tibble() # Convert into longer data theta_iter <- theta_iter %>% gather() %>% mutate(iter = rep(1:1000, times = nrow(Y))) # Plot theta_iter %>% ggplot(aes(iter, value)) + geom_line() + theme_bw() + xlab("Iteration") + ylab("") + facet_wrap(~ key, scales = "free_y")
大体イイ感じかと。それでは係数プロットを使って判事の理想点を可視化します。また、任命された政党のデータも追加してみました。
# Extract estimated theta theta <- fit$summary$theta # Plot theta %>% mutate(name = rownames(Y), party = c("Republican", "Republican", "Republican", "Republican", "Republican", "Republican", "Republican", "Democrat", "Democrat") %>% factor(levels = c("Republican", "Democrat"))) %>% ggplot(aes(y = reorder(name, mean), x = mean, color = party)) + geom_pointrange(aes(xmin = lwr, xmax = upr), size = 1) + theme_light() + xlab("Estimated Ideal Point") + ylab("") + theme(legend.title = element_blank(), legend.position = "bottom", legend.direction = "horizontal")
符号が反転してる感がありますが（通常、保守は右側でリベラルは左側）、それなりに妥当な理想点が推定できていると思います。
おわりに
余力があれば順序モデルや多項モデルに拡張していきたいですね。また、先ほどのMartin & Quinn (2002) で紹介されていたダイナミックな時変的な推定も実装してみたいですね。
実はCppで対数尤度関数を書く前に、全部Rで書いていたのですが結構遅くて、Cppに切り替えたら10倍以上速くなって感動しました。以上です。何か間違っていることがあればお知らせください。
*1:Clinton, J., Jackman, S., & Rivers, D. (2004). The statistical analysis of roll call data. American Political Science Review, 98(2), 355-370.
*2:Rosas and Shomer (2015) がこれらの取り扱いについて詳しく述べています。
*3:実際には、空間投票モデルやランダム効用理論を用いてIRTとの類似性を議論する必要がありますがここでは割愛します。
*4:Martin, A. D., & Quinn, K. M. (2002). Dynamic ideal point estimation via Markov chain Monte Carlo for the US Supreme Court, 1953–1999. Political analysis, 10(2), 134-153.
