目次
- 目次
- ZEN Studyにおける学力定量化の試みと狙い
- 項目反応理論とベイズ推定のおさらい
- NumPyroの概要
- NumPyroによるIRTモデルの実装
- モデル定義と計算実行
- 事後サンプリングの可視化
- 事後予測モデルチェック
- 垂直尺度化のためのモデル拡張
- 落ち穂拾い【MCMCの実用上の注意点】
- おわりに
- We are hiring!
ZEN Studyにおける学力定量化の試みと狙い
こんにちは。ドワンゴ教育事業本部のデータサイエンティストの板宮です。主に学力測定基盤の開発を担当しています。
ドワンゴの教育事業本部では、ZEN Studyと呼ばれる、N高グループ・ZEN大学向け教育用プラットフォームを開発しています。ZEN Studyにおける学びの重要なテーマのひとつが、基礎学力の着実な習得と学習意欲の喚起です。そこに向けた取り組みの一環として、主にN高グループ生を対象として、任意で参加できるオンライン試験を実施しています。この試験は成績評定のためではなく、生徒に自身の学力の伸びを理解してもらうことを目的に設計されています。これにより生徒が持続的学習に取り組んでいけることを目指しています。
特に、学力の伸びを科目ごとに可視化するという部分について,ドワンゴでは項目反応理論 (item response theory, IRT) モデルによる垂直尺度化を活用しています。 IRTとは古くから知られているテスト理論（学力テストや質問紙調査を統計的に分析する学問体系のこと）の統計モデルの一種で、継続的な学力の変化の数値化によく利用されてきました。垂直尺度化とは、対象の学年や難易度が異なる同一教科のテストを比較可能にする分析のことです。 例えば、PISAやTIMSSといった国際学力調査や、全国学力学習状況調査の経年変化分析、埼玉県の学力調査などでもIRTは活用されています。横浜市の学力調査ではIRTによる垂直尺度化が行われています。これらの事例を参考にしながらZEN Studyでも高校生の学力を公平に測るための方法を模索しているのですが、通信制高校とはいえ、単独の高校がオンラインでIRTによって学力を測定することはまだ珍しい試みだと考えています。
ZEN StudyではIRTを用いた垂直尺度化を実装するに当たり、PythonとNumPyroによるIRTモデルのベイズ推定を採用しています。ベイズ推定は、生徒の試験問題への解答データだけでなく、過去の実施結果から分かっている情報を取り込めるため、比較的小規模な試験実施になりがちなケースにおいて力を発揮します。NumPyroはJAXによる自動微分とJITコンパイルを備え、かつ、Pyroのようにわかりやすく確率モデルを記述できるという利点を備えています。IRTモデルのベイズ推定では、やや負荷の高い数値計算が必要になるのですが、このときにNumPyroのハイパフォーマンスな処理が役に立ちます。
今回は人工データによる垂直尺度化の分析例を通して、NumPyroを用いたIRTモデルのベイズ推定の実装方法を解説していきます1。
分析に関するサンプルコードはこちらに公開してあります。
項目反応理論とベイズ推定のおさらい
項目反応理論（IRT）は、受験者の能力と問題（項目）の特性を分離して表現する統計モデルです。もっとも基本的なIRTである1パラメタロジスティックモデル(1PLM)では、受験者の能力（ $\theta$ ）と項目の困難度（ $b$ ）、そしてテスト全体で共通する識別力（ $a$ ）をパラメタとして定義し、受験者が特定の項目に正解する確率をモデル化します。1PLMは、以下のような形で表現されます。
$$ P(X=1|\theta, b) = \frac{1}{1 + e^{-(a(\theta - b))}} $$
ここで、 $X$ は受験者の解答（1: 正解, 0: 不正解）です。IRTモデルでは、受験者の能力と項目の特性を同時に推定できます。1PLMの定義ではテスト全体で共通する識別力 $a$ は省略され、受験者の能力分布の分散（の逆数の平方根）として推定されることが多いですが、このように識別力として推定可能です。
条件付き独立の仮定によりIRTモデルの尤度関数は次のように分解できます。
$$ P(\mathbf{X}|\mathbf{\theta}, \mathbf{b}, a) = \prod_{i=1}^{N} \prod_{j=1}^{J} P(X_{ij}|\theta_i, b_j, a) $$
ベイズ推定では、事前分布と尤度関数を組み合わせて事後分布を求めます。IRTモデルにおいては、受験者の能力パラメタ（ $\theta$ ）や項目パラメタ（ $b, a$ ）に対して適切な事前分布を設定し、観測データ（ $X$ ）から事後分布を推定します。これにより、受験者の能力や項目の特性をより正確に推定できます。各パラメタの事前分布を$P(\theta)$, $P(b)$ , $P(a)$ とすると、ベイズの定理により事後分布は次のように表されます。
$$ P(\mathbf{\theta}, \mathbf{b}, a|\mathbf{X}) = \frac{P(\mathbf{X}|\mathbf{\theta}, \mathbf{b}, a) P(\mathbf{\theta}) P(\mathbf{b}) P(a)}{P(\mathbf{X})} $$
この事後分布の推定方法として計算機的なアプローチ、つまり、コンピューターを用いた疑似乱数生成と、その乱数によって事後分布を近似する方法を用います。このような疑似乱数を使って事後分布を近似する方法のひとつが、MCMC（マルコフ連鎖モンテカルロ法）です。MCMC法の中でも特にNUTS（No-U-Turn Sampler）と呼ばれるアルゴリズムがベイズ推定では広く用いられています。NUTSは、HMC（ハミルトニアンモンテカルロ法）の一種で、複雑なモデルであっても効率的にサンプリングを行うことができます。具体的なアルゴリズムの詳細はNUTSの提案論文2を参照してください。
ここでは実際に手元のPCでIRTモデルの推定プログラムを実装し、NUTSによるMCMC法を用いて受験者の能力や項目の特性を推定する方法を解説します。さらに、得られた事後分布を可視化し、MCMCの収束判定やモデルの適合度を評価する方法についても触れます。NUTSはゼロから自分で実装するには複雑すぎるアルゴリズムですが、NumPyroと呼ばれるPythonのライブラリを使用することで、簡単にNUTSを利用したMCMC法を実行できます。
NumPyroの概要
NumPyro (numpyro) はStanやPyMCのような確率的プログラミング言語（probabilistic programming language, PPL）の一種です。NumPyroはその名の通りPyroと似たような書き方ができるだけでなく、JAXと呼ばれるハイパフォーマンスな数値計算ライブラリと連携することで、JITコンパイルによるHMCの勾配計算の高速化やNumPyroと同じ記法による行列演算といった機能が拡張されています。
PPLの中では、Stan（Rstan, PyStan）やTuring(Julia)も有名だと思います。NumPyroの特長はStanよりも記述が容易であり、Turingと同程度（かそれ以上に）高速で安定したサンプリングが可能なことだと感じています。
例えば、StanはC++のラッパーですから、原則としてユーザーが変数宣言や変数の型指定をStan専用の言語で記述する必要があります。NumPyroではPythonの文法に従って記述すればよく、基本的には型指定せずとも動作する点が魅力的です。
モデル記述が簡単でハイパフォーマンスと謳われるライブラリにJulia言語のTuringがありますが、私の経験上、IRTのような多変数の潜在変数モデルのサンプリングにおいてはNumPyroの方が高速にサンプリングが可能だと思います。
また、ArviZと呼ばれるベイズ推論の評価・可視化ライブラリと連携できるため、目的の分布に到達しているかどうかの確認や事後分布の可視化などを容易に行うことができます。
NumPyroによるIRTモデルの実装
環境構築
まずは分析に必要なライブラリをインストールします。
pip install numpyro arviz matplotlib
次に、インストールした環境でpythonを立ち上げ、以下のコードを実行します。jax は numpyro の依存ライブラリなので自動的にインストールされています。
import numpyro import numpyro.distributions as dist from numpyro.infer import MCMC, NUTS, Predictive import jax.numpy as jnp from jax.random import PRNGKey import matplotlib.pyplot as plt
モデル定義と計算実行
まずはIRTモデルを定義します。以降では、観測データを X 、受験者数を N 、項目数を J とします。観測データは大きさ (N, J) の0と1からなる二値データ配列です。なお、欠測値がある場合は np.nan で表現可能です。
def model(X=None, N=None, J=None): if N is None or J is None: N, J = X.shape with numpyro.plate("students", N): theta = numpyro.sample("theta", dist.Normal(0, 1)) with numpyro.plate("items", J): b = numpyro.sample("b", dist.Normal(0, 1)) a = numpyro.sample("a", dist.LogNormal(0.15, 0.60)) plogits = a*(theta[:, None] - b[None, :]) if X is not None: mask = jnp.isnan(X) else: mask = jnp.ones((N, J), dtype=bool) with numpyro.handlers.mask(mask=~mask): with numpyro.plate_stack("obs", (N, J)): numpyro.sample("obs", dist.BernoulliLogits(plogits), obs=X)
このモデルをプレート表現にすると、下図のようになります。
モデルの内容について上から順に説明していきます。
if N is None or J is None: # Xの形状からNとJを取得 N, J = X.shape
この部分は観測データの次元数、つまり受験者数と項目数を取得しています。モデルの引数で X=None, N=None, J=None としている理由は、観測データを None にし、N, J だけ指定することで、Predictive を用いた人工データの生成（いわゆるprior predictiveによる観測データのサンプリング）を行うためです。
with numpyro.plate("students", N): theta = numpyro.sample("theta", dist.Normal(0, 1)) with numpyro.plate("items", J): b = numpyro.sample("b", dist.Normal(0, 1))
次に、この部分では各パラメタの事前分布を指定しています。numpyro.plate は、同じ分布に従う条件付き独立を表す部分で、コンテキストマネージャーと呼ばれます。この with 句の中ではplateで指定した次元数に沿って自動的に乱数生成がブロードキャストされます。識別力パラメタについては1次元なので、特に numpyro.plate は必要ありません。
numpyro.plate を利用せずとも、多次元配列として numpyro.sample("theta", dist.Normal(0, 1, (N,))) のように指定できます。今回は扱いませんが、確率的変分推論 SVI などの特定の最適化法を用いるためには、このコンテキストマネージャによる明示的な条件付き独立の指定が必要になります。
if X is not None: mask = jnp.isnan(X) else: mask = jnp.ones((N, J), dtype=bool)
続くブロックでは観測データに欠測値があった場合を考慮した処理をしています。mask の変数は次のブロックで使用されます。
with numpyro.handlers.mask(mask=~mask):
変数 mask は欠測値がある場合はその位置を True に、ない場合はすべて False に設定されます。これにより、観測データが欠測している部分を無視してサンプリングを行うことができます。
最後に、反応確率のロジットを計算し、dist.BernoulliLogits を用いて観測データの確率分布を指定しています。通常、識別力は正の値を取ると仮定されるので、ここでは対数正規分布を事前分布に設定しています。これにより非負制約のサンプリングが可能になります。
plogits = a*(theta[:, None] - b[None, :])
with numpyro.plate_stack("obs", (N, J)): numpyro.sample("obs", dist.BernoulliLogits(plogits), obs=X)
やはり、ここでも numpyro.plate_stack という多次元のコンテキストマネージャを用いて、観測変数が条件付き独立であることを明示しています。
人工データの生成
モデルが定義できたので、次に人工データを生成してみましょう。Predictive クラスを用いて、事前分布からのサンプリングを行います。ここでは、受験者数1000人、項目数30問の人工データを生成します。
prior_predictive = Predictive(model, num_samples=1) X = prior_predictive(PRNGKey(0), N=1000, J=30)["obs"][0] print(X)
[[1 0 0 ... 1 0 0] [0 0 0 ... 0 0 0] [1 0 1 ... 1 0 0] ... [1 0 1 ... 0 1 0] [0 0 1 ... 0 0 0] [0 1 1 ... 0 0 1]]
num_samples=1 とすることで、生成するデータを1セットのみに限定しています。包括的なシミュレーションの際には、この数字を任意の整数に変更することで、1回で大量の人工データを生成できます。定義していたモデルには引数が3種類ありましたが、このうち、観測データに対応している X の部分を省略することで、事前分布からの人工データ生成が可能になります。
サンプラーの指定とサンプリングの実行
定義したモデルを用いてMCMCによる事後分布からのサンプリングを実行します。NumPyroではサンプリングまでに3つのステップを記述します。まずサンプリングアルゴリズムを指定します。
kernel = NUTS(model)
デフォルトの設定でNUTSを動かす場合、他に指定する引数はありません。初期値を指定する場合は init_strategy 引数でサンプリングされる変数名（サイト名）をkeyに持つ辞書を渡します。
kernel = NUTS(model, init_strategy=numpyro.infer.init_to_value(values={"a": 0.5, "b": jnp.zeros(30), "theta": jnp.zeros(1000)}))
次に、MCMCのインスタンスを作成します。ここではサンプリング回数やチェイン数などを指定します。
mcmc = MCMC(kernel, num_warmup=500, num_samples=1500, num_chains=3, chain_method="vectorized")
ここではウォームアップ期間を500回、サンプリング回数を1500回、チェイン数を3に設定しています。ウォームアップ期間とはMCMCを実行した初期の反復期間であり、初期値に強く依存したサンプリングになる可能性があるために分析から除外するための期間です。実際には num_samples で指定した数のサンプルが、後続の計算や可視化で用いられます。chain_method="vectorized" は、NumPyroが提供するベクトル化されたチェーン処理を利用するための指定です。chain_method='parallel' の場合、複数のデバイスを起動して並列計算します。'vectorized' は単一デバイス上でベクトル化された計算します。'parallel' は複数デバイスをあらかじめ設定しておく必要があります。numpyro.set_host_device_count(4) を jax を呼び出す前に実行することで4つのデバイスを起動して並列計算を行います。モデルが複雑な場合、コンパイルのオーバーヘッドが大きくなるため parallel だとかえって計算が遅くなることがあります。一方で、単純なモデルでサンプル数が大きい（例えば5000とか）の場合は、vectorized よりも parallel の方が高速かもしれません。
私の環境では、このサンプリング条件下だとvectorized がおよそ30秒に対して、parallel はおよそ20秒でしたので、parallel の方が計算が速かったということになります。一般に、num_samples や num_warmup が大きい場合は、parallel の方が高速になる傾向があると感じています。
最後に、MCMCのサンプリングを実行します。ここで初めてデータをインスタンスに渡します。PRNGKey(0) はNumPyroの乱数生成器のキーを指定しています。NumPyroではJAXの乱数生成器を利用しているため、乱数生成器のキーを明示する必要があります。
mcmc.run(PRNGKey(0), X=X)
サンプリング結果の確認
サンプリングが完了したら、事後分布が目的の分布にきちんと収束しているかどうかを確認します。収束とはMCMCによって得られたサンプルの系列が、モデルで推定したい事後分布からのサンプルとみなせる状態になっていることを指します。一般に、階層モデルやパラメタ間に強い相関があるモデル、裾が厚い分布だったりすると収束が遅くなることが知られています。
収束に関するもっとも簡単な確認方法はNumPyroで提供されている print_summary() メソッドによる要約統計量の出力です。これにより、各パラメータの事後分布の平均、標準偏差、中央値、90%信用（確信）区間、有効サンプルサイズ（n_eff）、Gelman-Rubin統計量（r_hat）などが表示されます。有効サンプルサイズとGelman-Rubin統計量は、MCMCの系列が目的の分布に収束しているかどうかの判断に用いることができます。
有効サンプルサイズとは一個前の値に依存して得られているMCMCによるサンプルの系列について、独立な分布からのサンプルであったら何個分の情報量であるかを表す指標です。今回のサンプリングの設定上では本当は4500個分のサンプルが手元にはあるのですが、連続したサンプル間は似たような数値になっていたら、実質的には4500個よりもずっと小さいサンプルサイズ分の情報しか持っていないかもしれません。有効サンプルサイズが小さい（例えばチェイン数×100を下回る）場合は、うまくサンプルが混合していない可能性があります。
Gelman-Rubin統計量（r_hat, $\hat R$）はチェイン内とチェイン間の分散の比に基づく指標で、1に近いほど系列が収束していることを示します。$\hat R$ がどれくらい小さければよいのかという目安は文献や研究によってまちまちです。経験的に、1.01よりも大きい場合、有効サンプルサイズも小さく、あまり混合していないと考えています。1.1や1.2でも良いとする場合もあるようですが、そのくらいの大きさだとあまり混合していなかったというケースも見かけます。少なくとも今回のような事後分布からのサンプリングが比較的容易なモデルに関して言えば、1.01を超える場合は、サンプリングに関して何らかの問題があると考えた方が良いでしょう。
n_eff がチェイン数×100を下回ったり、$\hat R$ が1.01を超えたりする場合は、パラメタ変換により事後分布からのサンプリングが容易な形式にモデルを変換（reparametrize）したり、事前分布をinformativeなものに変更したり、あるいは間引き(thinning)の数を5や10に増やすなどの工夫をするとよいでしょう。
結果の一部だけを省略して表示しています。このモデルでは（人工データなので当たり前ではありますが）すべてのパラメタで収束していると判断できる結果が得られています。
mcmc.print_summary()
| mean | std | median | 5.0% | 95.0% | n_eff | r_hat | |
|---|---|---|---|---|---|---|---|
| a | 0.54 | 0.02 | 0.54 | 0.51 | 0.58 | 500.54 | 1.00 |
| b[0] | -0.79 | 0.13 | -0.79 | -0.99 | -0.56 | 1183.04 | 1.00 |
| b[1] | -0.26 | 0.13 | -0.27 | -0.47 | -0.04 | 1488.54 | 1.00 |
| b[2] | -0.94 | 0.13 | -0.94 | -1.16 | -0.73 | 1174.92 | 1.00 |
| b[3] | 0.85 | 0.13 | 0.85 | 0.63 | 1.05 | 1218.97 | 1.00 |
| b[4] | -0.83 | 0.14 | -0.83 | -1.06 | -0.61 | 1379.71 | 1.00 |
| b[5] | 0.57 | 0.14 | 0.57 | 0.34 | 0.78 | 1545.37 | 1.00 |
| b[6] | 0.82 | 0.14 | 0.82 | 0.59 | 1.03 | 1578.73 | 1.00 |
| b[7] | -0.20 | 0.13 | -0.21 | -0.45 | -0.01 | 1243.70 | 1.00 |
| b[8] | -1.23 | 0.14 | -1.23 | -1.48 | -1.02 | 1454.54 | 1.00 |
| b[9] | -0.09 | 0.13 | -0.09 | -0.29 | 0.12 | 1201.28 | 1.00 |
| b[10] | -1.16 | 0.14 | -1.16 | -1.38 | -0.93 | 1481.03 | 1.00 |
| b[11] | -0.99 | 0.13 | -0.99 | -1.22 | -0.78 | 1280.52 | 1.00 |
| b[12] | 1.08 | 0.14 | 1.08 | 0.83 | 1.29 | 1772.44 | 1.00 |
| b[13] | 0.19 | 0.13 | 0.19 | -0.03 | 0.40 | 2223.95 | 1.00 |
| b[14] | 0.08 | 0.13 | 0.08 | -0.12 | 0.28 | 1633.67 | 1.00 |
| b[15] | -0.06 | 0.13 | -0.06 | -0.28 | 0.15 | 1444.06 | 1.00 |
| b[16] | -2.25 | 0.17 | -2.24 | -2.54 | -1.98 | 998.65 | 1.00 |
| b[17] | -0.07 | 0.13 | -0.07 | -0.27 | 0.15 | 1430.26 | 1.00 |
| b[18] | 0.60 | 0.13 | 0.60 | 0.39 | 0.80 | 2069.28 | 1.00 |
| b[19] | -0.14 | 0.13 | -0.13 | -0.36 | 0.08 | 2038.62 | 1.00 |
| b[20] | 0.89 | 0.14 | 0.90 | 0.68 | 1.13 | 1542.78 | 1.00 |
| b[21] | -0.08 | 0.13 | -0.08 | -0.29 | 0.12 | 1539.25 | 1.00 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| theta[998] | -0.96 | 0.61 | -0.96 | -1.93 | 0.08 | 1852.54 | 1.00 |
| theta[999] | -0.20 | 0.61 | -0.20 | -1.20 | 0.76 | 2008.64 | 1.00 |
事後サンプリングの可視化
ArviZの利用
ここでは、ArviZを用いて事後分布の可視化と収束判定を行います。az.from_numpyro() 関数を用いて InferenceData オブジェクトを作成します。このオブジェクトは、事後分布の可視化や評価に使用されます。
ArviZに渡す際に、事後分布からのサンプリングも引数に与えることで、後述する事後予測モデルチェックが実行できます。まずはPredictive クラスを用いて事後予測分布からのサンプリングを生成しましょう。ここでは、もとの人工データの行列のサイズと同じ、1000人の受験者に対して30項目のテストを行う条件を指定します（違う数だとエラーになります）。観測データは不要ですので X=None とします。
predictive = Predictive(model, posterior_samples=mcmc.get_samples()) posterior_samples = predictive(PRNGKey(1), X=None, N=1000, J=30)
事後予測分布からのサンプル数はMCMCのサンプルからウォームアップを除いた回数に一致します。ただし、チェイン数の分だけサンプルは独立に行われているため、この場合4500個の予測サンプルが生成できます。
import arviz as az idata = az.from_numpyro(mcmc, posterior_predictive=posterior_samples)
トレースプロット
az.plot_trace 関数はトレースプロットと事後カーネル密度関数を描画します。トレースプロットとはMCMCの事後サンプルの系列を折れ線で結んだグラフで、サンプルが目標の分布に収束しているかどうかを視覚的に確認するために用いられます。事後カーネル密度関数は、サンプルの度数分布をなめらかな線で近似したもので、パラメタの事後分布の形状を視覚的に把握するために用いられます。
多次元のパラメタ（例えば b や theta）の場合、coords 引数を指定することで、項目や受験者ごとに表示します。そうしないと複数のパラメタのトレースが重なってしまい、非常に見にくくなります。'{パラメタ名_dim_0': [次元のインデックス]} と指定します。
以下のコードでは、識別力パラメタ a、項目パラメタ b の最初の項目、受験者能力パラメタ theta の最初の受験者に対するトレースプロットを描画しています。
az.plot_trace(idata, var_names=["a"]) az.plot_trace(idata, var_names=["b"], coords={"b_dim_0": [0]}) az.plot_trace(idata, var_names=["theta"], coords={"theta_dim_0": [0]})
先ほどの有効サンプルサイズの数値上だと a の値がやや小さい結果が確認できていました。a パラメタはロジスティック関数の傾きに比例する値、つまり能力分布の分布の分散に対応するようなパラメタです。正規分布によるギブスサンプリングなどの振る舞いでは有名ですが、分散パラメタのサンプルは自己相関が高く、有効サンプルサイズが小さくなる傾向があります。それと同様に、b や theta に比べて a はこのような結果になったと考えられます。
自己相関
実際に a の自己相関も確認してみましょう。自己相関とは $k$ 個前のサンプルと現在のサンプルの相関のことで、MCMCにおいてはサンプルの順番が $k$ にあたり、これをラグと呼びます。例えば10個前のサンプルとの自己相関はラグ10の自己相関です。一般的に、MCMCではラグが増えるほど事後相関は小さくなるのですが、この減少速度が緩やかだと、サンプルがあまり独立していないことを意味します。
az.plot_autocorr(idata, var_names=["a"], max_lag=20, combined=True)
ラグ3か4あたりまで自己相関がやや高く、10あたりまでうっすら自己相関が残っている様子が確認できます。試しに間引きを10、つまり10回に1回のサンプリング結果だけを保存する設定で自己相関がどのようになるのか確認してみましょう。
kernel = NUTS(model) mcmc_thin10 = MCMC(kernel, num_warmup=500, num_samples=15000, num_chains=3, chain_method="parallel", thinning=10, progress_bar=True) mcmc_thin10.run(PRNGKey(0), X=X) az.plot_autocorr(mcmc_thin10, var_names=["a"], max_lag=20, combined=True, figsize=(10, 5))
先ほどまで残っていた小さいラグでの自己相関が消えています。間引きは自己相関を減らすことに有効な方法と考えられますが、その分サンプリング数を多くとらなければならないという点には注意が必要です。大きなモデル（例えば受験者数が3000人を超えるような場合や、項目数が100を超えるような場合）では計算コストが非常に大きくなる傾向があります。
自己相関は有効サンプルサイズの計算に用いられており、本質的には自己相関が小さければ有効サンプルサイズも大きくなります。自己相関を見ることで有効サンプルサイズが小さいときにどれくらい間引けば良いかの目安を得ることができます。
フォレストプロット
ベイズ推定ではパラメタの事後分布の最高事後密度 (highest posterior density, HPD) 区間を評価できます。HPD区間とは事後分布の密度が高い（パラメタが高い確率で分布すると考えられる）区間のことで、モデルが考えるパラメタの確信度を表現した数値です。そして、これをプロットしたものがフォレストプロットです。ArviZでは az.plot_forest() 関数を用いて、事後分布の平均と94%HPD区間を表示できます。
az.plot_forest(idata, var_names=["a", "b"], combined=True, figsize=(8, 8), hdi_prob=0.94)
事後予測モデルチェック
モデル評価の一環として、事後予測モデルチェック(posterior predictive model checking, PPMC)します。事後予測モデルチェックは、モデルがデータをどの程度うまく説明できているかを評価するための手法で、得られた事後分布から観測データを復元し、その復元データがどの程度、実際の観測データに一致しているかを評価します。IRTのPPMCでは素点分布（周辺分布）などが用いられることがありますが、ArviZには個々の観測変数のデータしか渡されていないので、ここでは項目ごと、あるいは受験者ごとの正答率に関しての事後予測分布を可視化します。
観測データは2次元で、1次元目が受験者、2次元目が項目を表します。例えば最初の項目のPPMC結果は以下のように取得できます。
az.plot_ppc(idata, var_names=["obs"], figsize=(10, 5), coords={"obs_dim_1": [0]})
オレンジの破線が予測平均で、黒の実線が観測データです。見慣れないグラフですが、横軸は観測データに関する軸で、縦軸がその割合です。横軸が非常に見づらいのですが、要するに離散変数が連続変数のように扱われているだけで、ヒストグラムのように解釈できます。中央の1.0の位置での段階関数のように見えますが、1.0よりも左は誤答を、右は正答の割合を表している離散の分布ということです。結局のところ、オレンジと黒の線が重なっていれば事後予測サンプルの平均が観測データの平均と近い、つまりモデル適合は良いことを意味しています。
実際にこのグラフをすべての項目と受験者について1つ1つ確認することは現実的とは言えません。実用的なIRTのPPMCを行うためには、NumPyroのPredictiveメソッドを利用して事後予測サンプルを取得し、そこから周辺分布を計算する必要があります。
ベイズ的 $p$ 値
ベイズ的 $p$ 値(Bayesian p-value)とは、事後予測モデルチェックの一環として、観測データと事後予測分布からのサンプルとの比較するための指標です。ベイズ的 $p$ 値は、観測データが事後予測分布からのサンプルに比べてある統計量$T$がどの程度極端であるかを示します。具体的には、観測データの平均が事後予測分布からのサンプルで求めた平均（の分布）の中でどの程度の位置にあるかを評価します。ベイズ的 $p$ 値が0.5に近い場合、モデルは観測データをうまく説明していると考えられます。一方、0.5から離れるほど、モデルの適合度が低いことを示します。
ここでも最初の項目のベイズ的 $p$ 値を計算してみましょう。az.plot_bpv() 関数を用いて、ベイズ的 $p$ 値を計算し、プロットします。kind="t_stat" は統計量の種類を指定するもので、ここでは平均を指定しています。t_stat="mean" は平均を統計量として使用することを示しています。
az.plot_bpv(idata, var_names=["obs"], coords={"obs_dim_1": [0]}, kind="t_stat", t_stat="mean")
垂直尺度化のためのモデル拡張
ここまではNumPyroでIRTの1PLMのパラメタ推定の方法を説明しました。IRTのパラメタ推定では通常、受験者集団の能力分布が平均0、標準偏差（分散）1になるように尺度が固定されます。その条件は先ほどのモデルでは以下の部分で表現されています。
with numpyro.plate("students", N): theta = numpyro.sample("theta", dist.Normal(0, 1))
さらに実践的なケースを考えましょう。IRTによって2つの異なるテストを分析する場合、典型的には2つのテストの間に同じ試験問題（共通項目）を設けておき、そのパラメタ推定値をアンカーとして2つのテストを対応づけします。この分析を等化（equating）と呼びます。等化の目的は、2つのテストの得点を同じ尺度上で比較できるようにすることです。特に、2つのテストが異なる難易度（例えば高校の履修課程における数学Aと数学Bのような場合）の対応づけの分析を垂直尺度化（vertical scaling）と呼びます。垂直尺度化をするためには2つの試験を受験する集団の学力差に注目することが重要です。一般的に、数学Aは高校1年生が、数学Bは高校2年生が受験することが多いでしょう。したがって、数学Bの受験者の方が数学Aの受験者よりも学力が高いと考えられます。
このまま個別にIRTの分析をしたとしても、数学Aに回答した高校1年生と数学Bに回答した高校2年生の平均と標準偏差がどちらも同じになってしまいます。本当であれば高校1年生よりも高校2年生の方が数学の能力が高いはずです。このように、異なるテストの回答データを独立に推定してしまうと、IRTを用いていても受験者の能力値や項目パラメタを共通尺度上で比較できません。
それでは2つのデータセットを共通項目を鍵に縦に結合して、1つのデータセットにまとめて分析するのはどうでしょうか。それぞれの学年が回答していない部分は欠測値(np.nan)としておけば良いでしょう。
これは半分正解で、半分不正解のやり方です。1回の分析にまとめることで能力尺度を共通化できるため、数学Aと数学Bの項目パラメタを共通尺度化できます。しかし、学年が違えば学力も大きく異なるはずですから、共通の事前分布という仮定は見直した方が良さそうです。
このような考え方に基づく推定方法が、多群（多母集団）IRTモデルの推定です。
多群IRTモデルの実装
多群IRTモデルを実装するためには、欠測値だけでなく集団に関する指示変数をNumPyroで扱えるように拡張する必要があります。欠測値については、先ほどのモデルで mask を用いて処理していましたので、この仕組みをそのまま使うことができます。集団に関する指示変数は、受験者がどの集団に属しているかを示す変数で、例えば数学Aの受験者は0、数学Bの受験者は1といったように指定します。
以下のようにモデルを定義します。
def model_multiplegroup(X=None, N=None, J=None, G=None, design_matrix=None): if N is None or J is None: N, J= X.shape if G is None: Warning("Group indicator vector, G, is not specified.") with numpyro.plate("Hyperparameters", np.unique(G).size-1): _mu_theta = numpyro.sample("_mu_theta", dist.Normal(0.0, 1.0)) _sigma_theta = numpyro.sample("_sigma_theta", dist.HalfNormal(1.0)) mu_theta = jnp.insert(_mu_theta, 0, 0.0) sigma_theta = jnp.insert(_sigma_theta, 0, 1.0) mu_theta = numpyro.deterministic("mu_theta", mu_theta) sigma_theta = numpyro.deterministic("sigma_theta", sigma_theta) with numpyro.plate("students", N): theta = numpyro.sample("theta", dist.Normal(mu_theta[G], sigma_theta[G])) with numpyro.plate("items", J): b = numpyro.sample("b", dist.Normal(0, 1)) a = numpyro.sample("a", dist.LogNormal(0.15, 0.60)) plogits = a*(theta[:, None] - b[None, :]) if X is not None: mask = ~jnp.isnan(X) elif design_matrix is not None: mask = design_matrix else: mask = jnp.ones((N, J), dtype=bool) with numpyro.handlers.mask(mask=mask): with numpyro.plate_stack("obs", (N, J)): numpyro.sample("obs", dist.BernoulliLogits(plogits), obs=X)
だいぶ最初のモデルよりも複雑になってしまいました。主な変更点のあるブロックについて説明していきます。
with numpyro.plate("Hyperparameters", np.unique(G).size-1): _mu_theta = numpyro.sample("_mu_theta", dist.Normal(0.0, 1.0)) _sigma_theta = numpyro.sample("_sigma_theta", dist.HalfNormal(1.0)) mu_theta = jnp.insert(_mu_theta, 0, 0.0) sigma_theta = jnp.insert(_sigma_theta, 0, 1.0) mu_theta = numpyro.deterministic("mu_theta", mu_theta) sigma_theta = numpyro.deterministic("sigma_theta", sigma_theta)
まずこのブロックですが、ここでは基準となる集団以外の集団の受験者能力パラメタ事前分布の平均と標準偏差も推定しています。このように「あるパラメタの事前分布のパラメタに対する事前分布」を超事前分布（hyperprior distribution）と呼びます。基準集団は数学Aの受験者としますので、数学Aの受験者の平均は0.0、標準偏差は1.0に固定されます。np.unique(G).size-1 は、集団数から基準集団を引いた数を指定しています。_mu_theta と _sigma_theta は基準集団以外の集団の受験者能力パラメタの平均と標準偏差を表します。
パラメタを推定するのは2つ目以降の集団だけで良いのですが、このままだと能力値の事前分布を指定する際に場合分けが必要になってしまいます。そこで _mu_theta に固定する部分を挿入した新たな変数 mu_theta を作成し、numpyro.deterministic で宣言しています。
numpyro.deterministic で宣言された変数は、その変数の値がトレースされ、事後分布のサンプリング結果に含まれるようになりますが、特定の一点の値に固定されるというわけではありません。
この例で言えば mu_theta[0] は0.0に固定されますが、mu_theta[1] についてはその直前の _mu_theta = numpyro.sample("_mu_theta", dist.Normal(0.0, 1.0)) 部分で事前分布を指定しているため、事後分布は通常通りサンプルされます。
実際にMCMCを実行し、推定結果を可視化してみましょう。
import numpy as np G = jnp.concatenate([jnp.full(500, i) for i in range(2)]) design_matrix = np.block([ [np.ones((500, 15)), np.ones((500, 15)), np.zeros((500, 15))], [np.zeros((500, 15)), np.ones((500, 15)), np.ones((500, 15))] ]).astype(bool) kernel = NUTS(model_multiplegroup) predict = Predictive(model_multiplegroup, num_samples=1) predictive_samples = predict(PRNGKey(0), X=None, N=1000, J=45, G=G) X = predictive_samples["obs"][0] X = jnp.where(design_matrix, X, jnp.nan) mcmc = MCMC(kernel, num_warmup=500, num_samples=1500, num_chains=3, chain_method="vectorized", progress_bar=True) mcmc.run(PRNGKey(0), X=X, N=1000, J=45, G=G)
フォレストプロットに対して、赤い点で真値を重ね書きしましょう。
idata = az.from_numpyro(mcmc) true_mu_theta = predictive_samples["mu_theta"][0] true_sigma_theta = predictive_samples["sigma_theta"][0] fig, ax = plt.subplots(figsize=(8, 3)) az.plot_forest(idata, var_names=["mu_theta", "sigma_theta"], combined=True, hdi_prob=0.94, ax=ax) y_ticks = ax.get_yticks() for i, true_val in enumerate(reversed(true_sigma_theta)): y_pos = y_ticks[i] ax.plot(true_val, y_pos, 'ro', markersize=8, alpha=0.8) for i, true_val in enumerate(reversed(true_mu_theta)): y_pos = y_ticks[i + len(true_sigma_theta)] ax.plot(true_val, y_pos, 'ro', markersize=8, alpha=0.8) plt.title('Posterior estimates vs True values')
fig, ax = plt.subplots(figsize=(8, 12)) az.plot_forest(idata, var_names=["a", "b"], combined=True, hdi_prob=0.94, ax=ax) true_a = predictive_samples["a"][0] true_b = predictive_samples["b"][0] y_ticks = ax.get_yticks() for i, true_val in enumerate(reversed(true_b)): y_pos = y_ticks[i] ax.plot(true_val, y_pos, 'ro', markersize=6, alpha=0.8) y_pos = y_ticks[len(true_b)] ax.plot(true_a, y_pos, 'ro', markersize=8, alpha=0.8) plt.title('Discrimination and Difficulty: Posterior estimates vs True values')
赤い点が概ねHPD区間に入っているので、良い推定結果が得られていそうです。
固定項目パラメタ法の実装
もうひとつの垂直尺度化の方法として、あらかじめ数学Aのテストの項目パラメタを固定して、数学Bのテストの受験者能力パラメタを推定する方法があります。これを固定項目パラメタ法と呼びます。こちらもNumPyroで実装してみましょう。
def model_fixed(X=None, N=None, J=None, fixed_a=None, fixed_b: None|np.ndarray=None): if N is None or J is None: N, J = X.shape if fixed_b is not None: mu_theta = numpyro.sample("mu_theta", dist.Normal(0.0, 1.0)) sigma_theta = numpyro.sample("sigma_theta", dist.HalfNormal(1.0)) else: mu_theta = numpyro.deterministic("mu_theta", 0.0) sigma_theta = numpyro.deterministic("sigma_theta", 1.0) with numpyro.plate("students", N): theta = numpyro.sample("theta", dist.Normal(mu_theta, sigma_theta)) if fixed_a is not None: # a is fixed to the provided value a = fixed_a numpyro.deterministic("a", a) else: a = numpyro.sample("a", dist.LogNormal(0.15, 0.60)) if fixed_b is not None: nan_mask = np.isnan(fixed_b) n_free = np.sum(nan_mask) if n_free > 0: with numpyro.plate("free", n_free): b_free = numpyro.sample("b_free", dist.Normal(0, 1)) b = jnp.full_like(fixed_b, 0.0) b = b.at[nan_mask].set(b_free) b = b.at[~nan_mask].set(fixed_b[~nan_mask]) else: b = fixed_b numpyro.deterministic("b", b) else: with numpyro.plate("items", J): b = numpyro.sample("b", dist.Normal(0, 1)) if X is not None: mask = ~jnp.isnan(X) else: mask = jnp.ones((N, J), dtype=bool) plogits = a*(theta[:, None] - b[None, :]) with numpyro.handlers.mask(mask=mask): with numpyro.plate_stack("obs", (N, J)): numpyro.sample("obs", dist.BernoulliLogits(plogits), obs=X)
このモデルでは、fixed_a と fixed_b の引数を追加しています。fixed_a は識別力パラメタを、fixed_b は困難度パラメタを固定するための引数です。これらの引数に値を指定してパラメタを固定するのですが、実際にはパラメタ固定しない項目もデータに含まれているはずなので、b については np.nan の部分だけをサンプリングし（つまり、固定せずに）、値が指定されている部分はその値に固定するという風にしています。
余談ですが、以下の部分だけ jax.numpy ではなく、numpy を使っています。
nan_mask = np.isnan(fixed_b) n_free = np.sum(nan_mask)
これは条件分岐に関してのトレースを避けるためです。jax.numpy で指定した変数はすべてトレースされるのですが、トレース中の変数に対して条件分岐を当てはめることができないため TracerBoolConversionError となります。
一個前のモデルでも np.unique(G) の部分で numpy を使っていました。こちらは似たような理由なのですが、jnp.unique としてしまうと違うエラー ConcretizationTypeError が出ます。このエラーはJITコンパイル時に想定したサイズが確定しないことに由来するものです。
JAXを使う上では基本的にはループ構造や、サイズに不確定性が残るような記述は避けるべきなので、ループはベクトル化、サイズは明示的に与えることを意識すると良いです。どうしても必要な場合は、なるべく小さなループにすることと、numpy を使って不要なトレースを避けることがおすすめです。
先ほどの多群IRTモデルでも使ったデータを再利用して固定パラメタ法でMCMCを実行してみましょう。まずは数学Aの受験者であるグループ0単体で項目パラメタを推定し、その点推定値（今回はEAP推定値を用います）を得た上で、数学Bの受験者であるグループ1の受験者能力パラメタと残りの項目パラメタを推定します。
X1 = X[G == 0, :30] kernel = NUTS(model) mcmc_first = MCMC(kernel, num_warmup=500, num_samples=1500, num_chains=3, chain_method="vectorized", progress_bar=True) mcmc_first.run(PRNGKey(0), X=X1, N=X1.shape[0], J=X1.shape[1]) eap_a = mcmc_first.get_samples()["a"].mean(axis=0) eap_b = mcmc_first.get_samples()["b"].mean(axis=0) fixed_b = jnp.concatenate([eap_b, jnp.full(15, jnp.nan)])[15:] X2 = X[G == 1, 15:] kernel_fixed = NUTS(model_fixed) mcmc_fixed = MCMC(kernel_fixed, num_warmup=500, num_samples=1500, num_chains=3, chain_method="vectorized", progress_bar=True) mcmc_fixed.run(PRNGKey(0), X=X2, N=X2.shape[0], J=X2.shape[1], fixed_a=eap_a, fixed_b=fixed_b)
推定結果を、先ほどと同様にフォレストプロットで確認してみましょう。
idata_fixed = az.from_numpyro(mcmc_fixed) fig, ax = plt.subplots(figsize=(8, 12)) az.plot_forest(idata_fixed, var_names=["b", "mu_theta", "sigma_theta"], combined=True, hdi_prob=0.94, ax=ax) true_b_all = predictive_samples["b"][0][15:] true_mu_theta_group1 = predictive_samples["_mu_theta"][0][0] true_sigma_theta_group1 = predictive_samples["_sigma_theta"][0][0] y_ticks = ax.get_yticks() ax.plot(true_sigma_theta_group1, y_ticks[0], 'ro', markersize=8, alpha=0.8) ax.plot(true_mu_theta_group1, y_ticks[1], 'ro', markersize=8, alpha=0.8) for i, true_val in enumerate(reversed(true_b_all)): y_pos = y_ticks[i + 2] ax.plot(true_val, y_pos, 'ro', markersize=6, alpha=0.8) plt.title('Fixed model estimates vs True values (Group 1)\nRed dots show true values')
多群IRTモデルとおおむね同じ結果が得られていそうです。1番目から15番目（インデックスだと0から14まで）の項目パラメタの真値（赤点）はEAP推定値（青丸）とは完全に一致しないため若干のずれがありますが、16番目以降の多くの項目で推定値の94%HPD区間に真値が含まれている様子が確認できているので、推定はうまくいっていると考えられます。
落ち穂拾い【MCMCの実用上の注意点】
ここではMCMC法を実務上で使うに当たって重要になる点をピックアップして説明します。
まず、MCMC法で得られたサンプルがいったい何回目のサンプルから事後分布からのサンプリングと呼んでよいか不明確という点です。マルコフ連鎖は初期値の値に依存していますので、系列の最初の方は事後分布とは関係ない初期分布の影響を強く受けているはずです。そのためサンプリングの最初の方を捨てる必要があるのですが、この最初の方が何回くらいなら良いのかは明確に定まっていません。多くのケースでは全体のサンプルの20~50%程度は捨てるようです。これに関して言えば初期分布（初期値）を適切に設定することで、捨てるサンプルの割合を減らすことができるでしょう。IRTモデルであれば、次のような初期値が考えられます。
$$ \begin{aligned} \theta_{i, init} &= \Phi^{-1}\left(\frac{\sum_{j=1}^{J} X_{ij}}{J}\right) \\ b_{j, init} &= \Phi^{-1}\left(1 - \frac{\sum_{i=1}^{N} X_{ij}}{N}\right) \\ a_{init} &= 1 \end{aligned} $$
ここで、 $\Phi^{-1}$ は標準正規分布の累積分布関数の逆関数（逆正規分布関数）、 $J$ は全項目数を表します。これは、各受験者の正答率と各項目の正答率を基に初期値を設定する方法です。同様に項目パラメタ $b_j$ の初期値は、全受験者の正答率を基に設定しています。識別力 $a$ は1.0を置いています。識別力は尺度の分散に対応するパラメタですので、よほど信頼性が低いテストでない限り1.0で問題ないと考えられます。ただし、固定項目パラメタ法を用いる場合、事前分布の平均や最頻値の値の妥当性は、固定する尺度の原点と単位に依存するため、そのあたりを自動的に調節できるように（能力分布のように）これらの事前分布に超事前分布を仮定することがあります。
次に、マルコフ連鎖に従う系列は自己相関を持っている点です。自己相関とは、ある時点の値が、過去の値にどれだけ依存しているかを表す指標です。自己相関が高いと、サンプル間の情報が重複してしまい、何個サンプリングしても似たような値ばかりが連続して出てしまい、結果として事後分布を上図に近似できるような（独立な）サンプルを得ることが難しくなります。NUTSは自動的にサンプリングに必要なチューニングを行ってくれる賢い方法ですが、事後分布自体が探索しにくい形状だったりすると、NUTSでも自己相関が高くなる傾向があります。特に識別力パラメタについては $b$ パラメタとの関連性が強く出るため自己相関が高くなる傾向があります。関連性が強く出るというのは、例えば $b$ の値の絶対値が極端に大きいケースでは識別力に関する情報をほとんど持っていないため、事後分布も0付近で密度が高くなることが多いです。これは1PLMよりもさらに発展的なモデル（2PLMや3PLM）で顕著になりやすい問題です。私の知る限り、IRTの識別力パラメタの自己相関の高さについては有効な方法はまだ確立されていませんが、対処療法的には間引き(thinning)を5~10程度とることで自己相関を下げることができます。ただしそれにつれて計算時間も長くなりますので、注意が必要です。もしくは、有効な事前情報がある場合にはinformative priorを設定することも有効でしょう。
代表値に何を用いるかも重要な点です。事後分布の代表値としては期待値（Expected A Posteriorの頭文字でEAP）や最頻値（Maximum A Posterioriの頭文字でMAP）などが考えられます。EAPは事後分布の平均値を求める方法で、MAPは事後分布のもっとも高い点を求める方法です。IRTモデルでは、EAPとMAPはほぼ同じ値になることが多いですが、場合によっては大きく異なることもあります。特に、絶対値が大きな値や識別力においては事後分布が非対称になっていることがあるので、どちらを採用するかで値にも大きな違いが出るため、注意が必要です。
最後に、事前分布の設定です。事前分布の主観性はベイズのアキレス腱と呼ばれるほど、理論的な観点からも（実務的にも）重要な要素となっています。項目パラメタ $b_j$ については事前分布に有力な候補がなければ弱情報や無情報事前分布を採用すれば良いのですが、能力パラメタや識別力の場合そうはいきません。能力パラメタの事前分布は尺度の単位と原点を定める意味合いもあるため、基本的には標準正規分布を採用します。どうしても能力パラメタの事前情報を少なくしたいということであれば、先ほど説明した固定項目パラメタ法を利用して、一度、標準正規分布を事前分布としてMCMC法で推定した後に、項目パラメタをEAP推定値で固定した状態で再度、無情報事前分布で $\theta$ の事後分布を推定するという方法が考えられます。識別力には対数正規分布以外を用いることも可能です。特に、対数正規分布は $\log{a}$ が正規分布にしたがうような確率分布であるため、ハイパーパラメタ（事前分布のパラメタ）を直感的に調整しにくいです。半コーシー分布を用いても良いですし、識別力の二乗は尺度の分散パラメタに対応するという点に着目して、 $a^2$ に対して逆ガンマ分布を仮定してもよいでしょう。
おわりに
実際のテストデータの分析場面ではここまで単純なケースはまれですが、比較的実践的なNumPyroの実装例を示すことができたと思います。より複雑なモデルとしては多次元IRTや部分点を扱うようなモデルが存在します。NumPyroを用いればそうした発展的なモデルの実装も比較的容易に行うことができます。 今後は、より発展的なモデルを含むIRTに関する理論的な解説や、より高度なベイズ推定法の実装について解説したいと思います。
私たちドワンゴの教育事業本部では、このような統計モデリングの技術を応用して、一人ひとりの生徒に最適な学びを届けるための研究開発を進めています。
We are hiring!
株式会社ドワンゴの教育事業では、一緒に未来の当たり前の教育をつくるメンバーを募集しています。 カジュアル面談も行っています。 もし、ここで紹介した取り組みに興味を持たれた方がいらっしゃいましたら、お気軽にご連絡ください！
カジュアル面談応募フォームはこちら。
開発チームの取り組み、教育事業の今後については、他の記事や採用資料をご覧ください。
- ベイズ推定に限らず、IRTの理論や実践についてもっと知りたいという方は、テストは何を測るのかという名著がありますので、是非ご一読ください。その他にも、Web上でフリーでアクセスできる2018年の教育心理学会のチュートリアルセミナー資料も、IRTの実践を知る上で参考になります。↩
- NUTSはNo-U-Turn Samplerの略で、HMCの一種です。HMCとはハイブリッド（ハミルトニアン）モンテカルロの略で、高次元の事後分布であっても確率分布の勾配情報を利用して棄却率を抑えつつサンプリングが可能な高性能なMCMCアルゴリズムです。NUTSはHMCの改良版で、HMCで必要とされるハイパーパラメタチューニングを自動的に行ってくれます。↩
