項目反応理論のモデル群

正規累積モデル

\(i\)番目の被験者の\(j\)番目の項目の値\(y_{ij}\)が二値\(\{0, 1\}\)であるとする（例えば正解・不正解だったり、アンケートの「あてはまる」「あてはまらない」という2件法など）。

\(y_{ij}\)の背後には潜在的な能力の連続量\(\theta_i \in \mathbb{R}\)が存在し、\(\theta_i\)が閾値\(b_j\)を超えていたら1、超えていなければ0が観測されるとする。つまり\(y_{ij}\)が以下のように決まるとする。

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } \theta_i < b_j\\ 1 & \text{ if } \theta_i \geq b_j\\ \end{cases} \end{split}\]

1パラメータ正規累積モデル

しかし、実際には被験者\(i\)の体調や運（たまたま正解できた）などにより、常にこのようにきれいに正解・不正解が決まるわけではないと考えられる。こうした誤差を表すパラメータ\(\varepsilon_{ij} \sim N(0, \sigma^2_{\varepsilon})\)も追加して

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\theta_i - \varepsilon_{ij}) < b_j\\ 1 & \text{ if } (\theta_i - \varepsilon_{ij}) \geq b_j\\ \end{cases} \end{split}\]

とする。誤差が確率変数のため、\(y_{ij}\)のとる値も確率変数として考えることができるようになる。\(b_j\)を移項すると

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\theta_i - \varepsilon_{ij} - b_j) < 0\\ 1 & \text{ if } (\theta_i - \varepsilon_{ij} - b_j) \geq 0\\ \end{cases} \end{split}\]

となる。\(\theta_i - \varepsilon_{ij} - b_j \sim N(\theta_i - b_j, \sigma^2_{\varepsilon})\)である。 \(\varepsilon_{ij}\)を移項すれば

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\theta_i - b_j) < \varepsilon_{ij}\\ 1 & \text{ if } (\theta_i - b_j) \geq \varepsilon_{ij}\\ \end{cases} \end{split}\]

でもあるので「\(y_{ij}=1\)となるのは誤差\(\varepsilon_{ij}\)が\(\theta_i - b_j\)以下のとき」とわかる。

仮に\(\varepsilon_{ij}\)が標準正規分布（\(\sigma^2_{\varepsilon} = 1\)の正規分布）に従うならば、特性値\(\theta_i\)の人が項目\(j\)に当てはまると回答する確率は

\[\begin{split} \begin{aligned} P(y_{ij} = 1) &= P(\varepsilon_{ij} \leq \theta_i - b_j)\\ &= \int_{-\infty}^{\left(\theta_i-b_j\right)} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^2}{2}\right) d z \end{aligned} \end{split}\]

となる（最後のは、\(\varepsilon_{ij}\)が従う標準正規分布のうち \(-\infty\) から \(ｋ\theta_i-b_j\) までの範囲の面積が\(P(\varepsilon_{ij} \leq \theta_i - b_j)\)ということ）。

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
epsilon = np.linspace(-4, 4, 1000)
y = norm.pdf(epsilon)
fig, ax = plt.subplots(figsize=[4,2])
ax.plot(epsilon, y, 'k-')

c = 1

ax.fill_between(epsilon, y, where=(epsilon < c), color='blue', alpha=0.3, label=r'$P(y_{ij} = 1)$')
ax.fill_between(epsilon, y, where=(epsilon >= c), color='red', alpha=0.3, label=r'$P(y_{ij} = 0)$')
ax.axvline(x=c, color='grey', linestyle='--', label=r'$\theta_i - b_j$')
ax.set(
    title=r'$\varepsilon_{ij} \sim N(0, 1)$',
    xlabel=r"$\varepsilon_{ij}$",
    ylabel="density"
)
ax.legend()
ax.grid(True)
fig.show()

2パラメータ正規累積モデル

\(\sigma^2_{\varepsilon}\)が項目ごとに異なる場合を考える。\(\sigma^2_{\varepsilon}=1/a_j\)とすると、誤差の確率分布は

\[ \varepsilon_{ij} \sim N\left(0, \frac{1}{a_j}\right) \]

となる。両辺を\(a_j\)倍すると、\(a_j \varepsilon_{ij} \sim N(0, 1)\)と表すことができ、引き続き標準正規分布を使うことができる。そのためモデルは\(a_i\)が追加され

\[\begin{split} \begin{aligned} P(y_{ij} = 1) &= P(a_j \varepsilon_{ij} \leq \theta_i - b_j)\\ &= \int_{-\infty}^{ a_j (\theta_i-b_j)} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^2}{2}\right) d z \end{aligned} \end{split}\]

となる。

パラメータの意味

\(b_j\)が大きくなると\(\theta_i - b_j\)の値は小さくなり、\(P(y_{ij} = 1)\)の面積が小さくなる。\(y_{ij} = 1\)が正解を表しているとするなら、正答率が低くなる方向に作用する。そのため\(b_j\)は 項目困難度（item difficulty） と呼ばれる。

項目特性曲線（ICC）との関係でいうと、ICCが0.5になる点が\(b_j = \theta_i\)となる。つまり、困難度が高いほど正答率50%の点にくる\(\theta_i\)も高くなる

また\(a_j\)は値が大きくなると\(\varepsilon_{ij}\)の分散を下げて分布がより尖っていく。また横軸に\(\theta_i - b_j\)、縦軸に\(P(y_{ij} = 1)\)のグラフを書くとき、この曲線の傾きを急にして、\(\theta_i\)が低い人と高い人の間で\(P(y_{ij} = 1)\)の変化を大きくする。例えば、\(\theta_i\)のある点を境にして正解/不正解がはっきり分かれるデータならICCの傾きは急になり、\(a_j\)は高くなる。 そのため\(a_j\)は 項目識別力（item discrimination） と呼ばれる。

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
x = np.linspace(-4, 4, 1000)

fig, axes = plt.subplots(figsize=[8, 2], ncols=2)


a = 1
for beta in [-1, 0, 1]:
    axes[0].plot(x, norm.cdf(x, loc=beta, scale=1/a), label=r"$b$ = " + f"{beta}")
axes[0].set(xlabel=r"$\theta_i$", ylabel=r"$P(y_{ij} = 1)$", xticklabels=[], yticklabels=[])
axes[0].legend()
axes[0].grid(True)

beta = 0
for a in [0.5, 1, 2]:
    axes[1].plot(x, norm.cdf(x, loc=beta, scale=1/a), label=r"$a$ = " + f"{a}")
axes[1].set(xlabel=r"$\theta_i$", ylabel=r"$P(y_{ij} = 1)$", xticklabels=[], yticklabels=[])
axes[1].legend()
axes[1].grid(True)

fig.show()

なお横軸に\(\theta_i\)、縦軸に\(P(y_{ij}=1)\)をとったグラフは 項目特性曲線 （item characteristic curve: ICC） と呼ばれる。

因子分析と2P正規累積モデルは等価

因子分析と2パラメータ正規累積モデルは数学的に等価であると知られている。標準化していない（切片が0でない）1因子モデルは、切片\(\tau_j\)、因子負荷量\(a_j\)、因子スコア\(f_i\)で

\[ y_{ij} = \tau_j + a_j f_i - \varepsilon_{ij} \quad (\varepsilon_{ij}\sim N(0, \sigma_{\varepsilon})) \]

となる（IRTの説明に合わせて誤差の符号をマイナスにしている）

カテゴリカル因子分析では「離散的な観測変数\(y_{ij}\)はその背後にある潜在的な連続変数によって決まる」という考え方をするため、IRTの冒頭の説明と同じで

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\tau_j + a_j f_i) < \varepsilon_{ij}\\ 1 & \text{ if } (\tau_j + a_j f_i) \geq \varepsilon_{ij}\\ \end{cases} \end{split}\]

つまり標準正規分布に従う誤差\(\varepsilon_{ij}\)が\(\tau_j + a_j f_i\)より小さいときに\(y_{ij}=1\)となると考えるため、その確率\(P(y_{ij}=1)\)は

\[\begin{split} \begin{aligned} P(y_{ij} = 1) &= P(\varepsilon_{ij} \leq \tau_j + a_j f_i)\\ &= \int_{-\infty}^{ (\tau_j + a_j f_i)} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^2}{2}\right) d z \end{aligned} \end{split}\]

となる。2つのモデルのパラメータを\((f_i, a_j, \tau_j) = (\theta_i, a_j, a_j b_j)\)と対応させると

\[\begin{split} \begin{aligned} \tau_j + a_j f_i &= a_j \left( f_i + \frac{\tau_j}{a_j} \right) \\ &= a_j (\theta_i - b_i) \end{aligned} \end{split}\]

となる。

ロジスティックモデル

正規累積モデルはプロビット回帰と同様のことをするので、コンピュータで積分計算をするのがやや難しいという問題がある。そこでロジスティック分布に置き換えたものが使われる。

ロジスティック分布の確率密度関数と累積分布関数は

\[ f(x)=\frac{\exp (-x)}{[1+\exp (-x)]^2}, \quad F(x)=\frac{1}{1+\exp (-x)} \]

となる。とくに\(x\)を約1.7倍したロジスティック分布は累積分布関数が正規分布と非常に近くなることが知られている。

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, logistic
x = np.linspace(-10, 10, 1000)

fig, axes = plt.subplots(figsize=[12, 3], dpi=90, ncols=2)
ax = axes[0]
ax.plot(x, norm.pdf(x), 'k-', label='Normal Distribution')
ax.plot(x, logistic.pdf(x),  'orange', alpha=0.9, label='Logistic Distribution')

ax.set_title('Normal Distribution vs Logistic Distribution')
ax.set_xlabel('x')
ax.set_ylabel('Probability Density')
ax.legend()
ax.grid(True)

ax = axes[1]
ax.plot(x, norm.cdf(x), 'k-', label=r'CDF of $Normal(x)$')
ax.plot(x, logistic.cdf(x), 'orange', alpha=0.9, label=r'CDF of $Logistic(x)$')
ax.plot(x, logistic.cdf(x * 1.704), 'red', alpha=0.9, label=r'CDF of $Logistic(1.7 x)$')
ax.set_title('Normal Distribution vs Logistic Distribution (CDF)')
ax.set_xlabel('x')
ax.set_ylabel('Probability')
ax.legend()
ax.grid(True)
plt.show()

2PLモデル

正規分布の代わりにロジスティック分布を使った 2パラメータロジスティック（2PL）モデル は以下のように表される。

2PLモデル

\[ P(y_{ij} = 1) = \frac{1}{1+ \exp(-D a_j \big( \theta_i - b_j) \big)} \]

• \(a_j\)：項目識別力

• \(b_j\)：項目困難度

なお\(D\)はロジスティック・シグモイド関数を正規累積モデルの関数に近づけるための定数（通常は\(D=1.7\)か\(D=1\)にする）なので、正規累積モデルと比較する必要がなければ不要（\(D=1\)でいい）。

3PLモデル

例えば4択問題では、正解がわからなくて適当に選んだとしても1/4は当たることになる。こうした影響を「当て推量」パラメータ\(c_j\)として取り入れたモデル。

3PLモデル

\[ P(y_{ij} = 1) = c_j + \frac{1 - c_j}{1 - \exp(-a_j \big( \theta_i - b_j) \big)} \]

• \(a_j\)：項目識別力

• \(b_j\)：項目困難度

• \(c_j\)：当て推量

\(c_j\)は項目特性曲線の下限となる。\(\theta_i\)がどんなに低い人でも必ず\(c_j\)以上の\(P(y_{ij} = 1)\)になるということ。

4PLモデル

項目特性曲線の上限を表すパラメータ\(d_j\)を追加したもの。\(\theta_i\)がどんなに高い人でも100%の正答率にはできない高難度な状況（運ゲー）を想定したモデル。

4PLモデル

\[ P(y_{ij} = 1) = c_j + \frac{d_j - c_j}{1 - \exp(-a_j \big( \theta_i - b_j) \big)} \]

• \(a_j\)：項目識別力

• \(b_j\)：項目困難度

• \(c_j\)：当て推量。項目特性曲線の下限

• \(d_j\)：項目特性曲線の上限

4PLMになるとかなりモデルが複雑になりパラメータの推定も不安定になるので、1~3PLMほど一般的ではない。対応していないライブラリも多い。

5PLモデル

「非対称性」のパラメータ\(e_j\)を追加したもの。4PLまでは項目特性曲線の動き方が0.5を中心に対称になっている。5PLでは「最初は\(\theta_i\)があがるほど急激に\(P(y_{ij}=1)\)が上がるが、徐々に上がりにくくなる」などの状況を表すことができる。

5PLモデル

\[ P(y_{ij} = 1) = c_j + \frac{d_j - c_j}{\left[ 1 - \exp(-a_j \big( \theta_i - b_j) \big) \right]^{e_j}} \]

• \(a_j\)：項目識別力

• \(b_j\)：項目困難度

• \(c_j\)：当て推量。項目特性曲線の下限

• \(d_j\)：項目特性曲線の上限

• \(e_j\)：非対称性

多値型モデル

項目の反応\(y_{ij}\)が多値になった場合のモデルも存在する。

段階反応モデル

段階反応モデル（graded response model: GRM） は複数の二値IRTモデルを組み合わせて多値反応を表現する。

回答者\(i\)の項目\(j\)に対する回答\(y_{ij} = k \quad (k=1,2,\dots,K)\)について、「\(k\)以上のカテゴリを選ぶ確率」を考えると、これはまだ「\(k\)未満 or \(k\)以上」の二値なので2PLなどで表せる。例えば以下のようになる。

\[ P(y_{ij} \geq k) = \frac{1}{1+ \exp \big(-a_j ( \theta_i - b_{jk}) \big)} \]

なお、困難度は項目\(j\)のカテゴリ\(k\)ごとに用意されるため\(b_{jk}\)に変更している。

このモデルを組み合わせると、「ちょうど\(k\)番目のカテゴリを選ぶ確率」は

\[ P(y_{ij} = k) = P(y_{ij} \geq k) - P(y_{ij} \geq k + 1) \]

と表すことができる。ただし端のカテゴリは\(P(y_{ij} \geq 1) = 1, P(y_{ij} \geq K + 1) = 0\)とする。また確率100%の困難度は低くて当然なので\(b_{j1} = -\infty\)とする。

段階反応モデル

\[ P(y_{ij} = k) = \frac{1}{1+ \exp(-a_j \big( \theta_i - b_{jk}) \big)} - \frac{1}{1+ \exp(-a_j \big( \theta_i - b_{jk+1}) \big)} \]

名義反応モデル

順序のないカテゴリカル変数を名義尺度（nominal scale）という。名義尺度の反応を扱えるのが 名義反応モデル（nominal response model: NRM） 。softmax関数のような形で多値化する。

\[ P(y_{ij} = k) = \frac{ \exp (a_{jk} \theta_i + \gamma_{jk}) }{ \sum_{\kappa=1}^K \exp (- a_{j\kappa} \theta_i + \gamma_{j\kappa}) } \]

連続反応モデル

連続反応モデル（Continuous Response Model, CRM）は連続的な反応（例：時間、強度、割合）を扱うことができる。

CRMは段階反応モデル（GRM）をさらに多段階に拡張していって連続値を扱う。

受験者\(i\)が項目\(j\)において\(x_j\)点以上をとる確率は、次のようにロジスティック関数で表すことができる

\[ p_i^*(x_j)=\frac{1}{ 1 + \exp \big(- a_j (\theta_i - b_{x_j}) \big)} \]

これを使うと、受験者\(i\)が項目\(j\)において\(x_j\)点をとる確率 \(p_i(x_j)\) を次のように定義できる

\[ p_i(x_j) = \lim _{\Delta x_j \to 0} \frac{p_i^*(x_j)-p_i^*(x_j+\Delta x_j)}{\Delta x_j} \]

\(b_{x_j}\)は項目\(j\)において\(x_j\)以上の点を取る困難度で、ロジスティック関数の逆関数を用いて以下のように定義される

\[ b_{x_j} = \beta_j + \frac{1}{\alpha_j} \log \frac{x_j}{K_j-x_j} \]

ここで\(\alpha_j, \beta_j\)はロジスティック関数の逆関数における識別力と困難度を表現するパラメータ、\(K_j\)は\(x_j\)のとる最大値（\(x_j \in [0, K_j]\)）

参考

• 初出：Samejima, F. (1973). Homogeneous case of the continuous response model. Psychometrika, 38(2), 203-219.

• 応用研究の例：秋山隆, 尾崎幸謙, & 豊田秀樹. (2012). 項目反応理論における連続反応モデルを用いた顧客の価格感度測定—ID-POS データへの適用例を通じて—. 行動計量学, 39(2), 93-101.

多次元モデル

多次元因子分析モデルでは項目反応\(y_{ij}\)を

\[ y_{i j} = \tau_j + \sum_{t=1}^T f_{i t} b_{t j} - \varepsilon_{ij} \]

のように表す。ここから考えると、多次元IRTモデル（の2PLM）の項目反応関数は

\[ P(y_{ij}=1) = \frac{1}{ 1 + \exp \left[ -\sum_{t=1}^T \alpha_{j t} \theta_{i t} - \tau_j \right]} \]

のように考えられる。\(\tau_j\)は切片のパラメータで、これを困難度に変換するには、次元ごとの識別力を一つにまとめた 多次元識別力 である

\[ \operatorname{MDISC}_j=\sqrt{\sum_{t=1}^T \alpha_{t j}^2} \]

を使って

\[ \beta_j=\frac{-\tau_j}{\mathrm{MDISC}_j} \]

と変換する（豊田 2013 『項目反応理論［中級編]』） 。

補償型と非補償型

ある次元の特性値が低くても、別の次元の特性値がとても高い場合は反応確率が高くなるようなモデルは 補償型 (compensatory) モデル と呼ばれる。例えば先述の

\[ P(y_{ij}=1) = \frac{1}{ 1 + \exp \left[ -\sum_{t=1}^T \alpha_{j t} \theta_{i t} - \tau_j \right]} \]

は補償型モデルである。

一方で複数の次元の特性値がそろわないと高くならない多次元IRTモデルは 非補償 (noncompensatory) モデル あるいは 部分補償(partially compensatory)モデル と呼ばれる（例：「総合学力」に対して「数学力スコア」と「英語力スコア」はそれぞれ独立していると考える）。最もシンプルな非補償モデルの項目反応関数は

\[ P(y_{ij}=1) = \prod_{t=1}^T \frac{1}{1+\exp \left[-\alpha_{j t}\left(\theta_{i t}-\beta_{j t}\right)\right]} \]

である。それぞれの次元の特性値\(\theta_{i t}\)で説明される確率の積にして同時確率のようになっている。

推定の難しさなどの観点から補償型のほうがよく使われるらしい

拡張されたIRTモデル

Latent Regression Model

被験者の潜在能力 \(\theta\) を観測共変量 \(x\) の関数としてモデル化する拡張。能力分布を固定（例：\(\mathcal N(0,1)\)）とせず、説明変数で条件付ける。

\[ \theta_i=\beta_0+\boldsymbol{\beta}^{\top} \mathbf{x}_i+\varepsilon_i, \quad \varepsilon_i \sim \mathcal{N}(0, \sigma_\theta^2) \]

• \(\mathbf{x}_i\)：年齢・性別・学歴・SES 等の共変量

• \(\boldsymbol{\beta}\)：能力に対する回帰係数

• \(\sigma_\theta^2\)：共変量で説明しきれない能力分散

実装方法

例えばRの mirt パッケージで推定できる

library(mirt)
mod <- mirt(data, 1,
            covdata = cov_df,
            formula = ~ x1 + x2)

Multilevel IRT

Multilevel IRT（階層IRT／マルチレベルIRT） は、被験者や項目が階層構造（クラス・学校・施設・テストフォーム等）にネストされている状況を明示的にモデル化するIRTの拡張。

種類

1. 能力パラメータ\(\theta\)の階層化

受験者\(i\)の能力がグループ\(g\)（例：学校）ごとに異なる

\[ \theta_{i g}=\mu_g+\eta_{i g}, \quad \mu_g \sim \mathcal{N}\left(\mu_0, \tau^2\right), \eta_{i g} \sim \mathcal{N}\left(0, \sigma^2\right) \]

2. 項目パラメータ\(a,b\)の階層化

項目困難度・識別力がグループごとに異なるとするもの（例えばオンラインのテストで、同じ設問でもPCとスマホで見るのとで違うとか？？）

\[ a_{j g} \sim \mathcal{N}\left(\bar{a}_j, \omega_a^2\right), \quad b_{j g} \sim \mathcal{N}\left(\bar{b}_j, \omega_b^2\right) \]

3. 両方

被験者×項目×群の三層（例：学生×問題×学校）

実装方法

• ベイズ（MCMC）が主流（階層が深いほど有利）。

• Rのmirtパッケージでの最尤推定： multipleGroup() 関数

適合度の評価

局所独立性の確認

• \(\chi^2\)統計量

個人適合度

「\(\theta_i\)が高い人なのに困難度が低い項目で間違えるのはおかしい」といった考え方から、特性\(\theta_i\)と困難度\(b_j\)の関係性を見る。

• \(z_h\)（\(l_z\)）統計量 ：ある回答者の反応パターン\(\boldsymbol{y}_i = (y_{i1}, y_{i2}, \dots, y_{iJ})\) と「考えられる全反応パターン」での尤度を比較する

IRTの前提条件・仮定

局所独立性（local independence）の仮定

\(\theta_i\)で条件づけたとき、項目\(j\)と項目\(l\) （\(j \neq l\)） への回答は完全に独立であるという仮定。最尤推定の計算の簡易化のために必要になる。

仮定を満たすか確認する方法：

• 例えばデータから\(\theta_i\)が特定の値の人を集めて（\(\theta_i\)で条件づけたデータを用意して）、項目\(j,l\)の回答のクロス表を作ったとき、関連性がない場合は局所独立性が満たされていることがわかる。

一次元性

もともとのIRTは単一のテストに対する単一の能力を測定することが主たる目的のため、推定する因子数は1である（因子数を拡張したモデルも存在するが、基本は1）。そのため、すべての項目が同じ能力・特性のみを反映していることが必要となる。

仮定を満たすか確認する方法：

• 例えば、1因子の確認的因子分析を行い、適合度指標（RMSEAやCFA、AGFIなど）が許容範囲内かどうかで確認する。

• DIMTEST（Stout et al., 1996）や DETECT（Zhang & Stout, 1999）といった手法もある

参考

• Chapter 8 項目反応理論 1

• 項目反応理論やってみよう