--- title: "段階反応モデル（GRM）をPyMCでベイズ推定する手順｜OrderedLogistic実装と収束診断" url: "https://www.issoh.co.jp/tech/details/12355/" published: 2026-06-03 updated: 2026-10-09 categories: ["Python"] publisher: "株式会社一創" --- # 段階反応モデル（GRM）をPyMCでベイズ推定する手順｜OrderedLogistic実装と収束診断 段階反応モデル（Graded Response Model、GRM）は、5件法のように順序のある選択肢への回答から、回答者の特性θと項目の識別力・境界困難度を推定する項目反応理論（IRT）のモデルです。この記事では、PyMC 6.3.2 でGRMをベイズ推定するコードを、実際に走らせた出力つきで示します。 掲載したコードと数値は、Python 3.12・PyMC 6.3.2・ArviZ 1.3.0 で疑似データ（300人×8項目×5件法）を推定して得たものです。境界の順序制約を外したとき、識別力の符号を縛らなかったときに何が起きるかも、エラー文と診断値で載せています。 ## まとめ：段階反応モデルをPyMCで推定する要点 - GRMは「カテゴリk以上を選ぶ確率」を P(Y≥k|θ)=1/(1+exp(−a(θ−bk))) で表し、隣り合う累積確率の差をカテゴリ確率にします（Samejima, 1969）。 - PyMCでは `pm.OrderedLogistic` に eta=a·θ、cutpoints=a·b を渡せばGRMになります。観測値は0始まり（5件法なら0〜4）に直します。 - 境界bには `transform=pm.distributions.transforms.ordered`、識別力aには LogNormal など正の事前分布を置きます。どちらかを外すと初期化の時点で `SamplingError` になりました。 - 判定は R-hat<1.01、ESS>400、発散0件が目安です（Vehtari et al., 2021）。実測ではchains=4・draws=1000・target\_accept=0.9 で R-hat最大1.004、ESS最小1,100、発散0件でした。 - 大標本で標準的な一次元GRMを当てるだけなら、Rの mirt（`itemtype='graded'`）の周辺最尤推定を先に検討します。PyMCを選ぶ理由は、少標本・疎なカテゴリ・階層化など事前分布が効く場面です。 ## 段階反応モデル（GRM）の定義と累積確率の考え方 GRMは Fumiko Samejima が1969年に Psychometrika のモノグラフ（Psychometric Monograph Supplement No.17、34巻S1号）で示したモデルです。二値の2パラメータロジスティックモデル（2PL）を、順序のある3つ以上のカテゴリへ広げた形になっています。 K件法の項目jでは、カテゴリ境界ごとに「k以上を選ぶ確率」を2PLの曲線で書きます。 ``` P(Y_ij ≥ k | θ_i) = 1 / (1 + exp(−a_j (θ_i − b_jk))) k = 1, …, K−1 P(Y_ij = k | θ_i) = P(Y_ij ≥ k | θ_i) − P(Y_ij ≥ k+1 | θ_i) （P(Y ≥ 0) = 1、P(Y ≥ K) = 0 とおく） ``` パラメータは項目ごとに識別力aが1つ、境界困難度bがK−1個です。5件法なら1項目あたり5個、8項目で40個になり、これに回答者数ぶんのθが加わります。K=2にすると境界が1本になり、2PLと一致します。 カテゴリ確率を差で作るため、bj10 なら、bが昇順である限り a·b も昇順に保たれます。境界をbのまま事前分布に置くので、「境界困難度」をθと同じ尺度で解釈でき、事前分布も Normal(0, 2) のように決めやすくなります。cutpoints を直接パラメータにする傾き切片型も書けますが、その場合は後述する符号の不定性に注意が要ります。 次は疑似データの生成から要約までのコードです。後掲の事後予測チェックと項目情報量のコードは、このファイルのmain()内の末尾に、各行を4スペース字下げして追加してください。両コードはmain()内で定義した変数を使います。 ``` import numpy as np import pymc as pm import arviz as az def main(): # 1. 疑似データ（300人 x 8項目 x 5件法） rng = np.random.default_rng(42) N, J, K = 300, 8, 5 theta_true = rng.normal(0, 1, N) a_true = rng.lognormal(0.2, 0.3, J) b_true = np.sort(rng.normal(0, 1, (J, K - 1)), axis=1) p_ge = 1 / (1 + np.exp(-a_true[None, :, None] * (theta_true[:, None, None] - b_true[None, :, :]))) u = rng.uniform(size=(N, J, 1)) Y = (u < p_ge).sum(axis=2) # 0..K-1 の整数 # 2. ロング形式へ（欠測があれば行ごと落とすだけで済む） person, item = np.indices(Y.shape) person, item, y = person.ravel(), item.ravel(), Y.ravel() coords = {"person": np.arange(N), "item": np.arange(J), "cut": np.arange(K - 1)} with pm.Model(coords=coords) as grm: theta = pm.Normal("theta", 0.0, 1.0, dims="person") a = pm.LogNormal("a", 0.0, 0.5, dims="item") b = pm.Normal( "b", 0.0, 2.0, dims=("item", "cut"), transform=pm.distributions.transforms.ordered, initval=np.tile(np.linspace(-1.5, 1.5, K - 1), (J, 1)), ) eta = a[item] * theta[person] cutpoints = a[item][:, None] * b[item] pm.OrderedLogistic("y", eta=eta, cutpoints=cutpoints, observed=y) idata = pm.sample(1000, tune=1000, chains=4, random_seed=1, nuts={"target_accept": 0.9}) summ = az.summary(idata, var_names=["a", "b"]) print(summ.head(12).to_string()) print("max r_hat", float(summ["r_hat"].max())) print("min ess_bulk", float(summ["ess_bulk"].min()), "min ess_tail", float(summ["ess_tail"].min())) print("divergences", idata["sample_stats"]["diverging"].sum().item()) if __name__ == "__main__": main() ``` θの事前分布を Normal(0, 1) に固定しているのは、尺度の原点と単位を決めるためです。GRMはθとbを同じ量だけ平行移動しても、θとbを同じ正の定数倍にしてaをその定数で割っても尤度が変わりません。尺度を識別する制約がないと、位置と尺度の不定性が推定を難しくします。本記事ではθの母集団分布を標準正規分布に固定しますが、項目パラメータを固定する識別方法もあります。 ordered変換の初期値 `initval` を厳密な昇順で与えるのは、隣接差が0以下になる初期値を避けるためです。等間隔にしているのは、ランダムな初期化で境界が極端な位置から始まるのを避けるためです。 ### NUTSの実行設定とmainガード NUTSの `target_accept` は既定0.8です。上げるほどステップ幅が小さくなり、発散は減る一方で1回の反復は遅くなります。PyMCの公式ドキュメントは、ウォームアップ後に発散が1件でも残る場合は結果が信頼できない可能性があるとして、再パラメータ化するか0.9〜0.95へ上げるよう案内しています。GRMはθの数が多く事後分布の形が場所によって変わるので、最初から0.9で回しました。chains=4 は、Vehtari et al.（2021）が既定として推奨する本数です。 コードを `if __name__ == "__main__":` で包んでいるのには理由があります。chains=4 は別プロセスで並列に走るため、macOSのように fork 以外の方式で子プロセスを起動する環境では、ガードがないと子プロセスがスクリプト全体を再実行します。ガードなしのスクリプトで chains=2・cores=2 のサンプリングを実行すると、`freeze_support()` を案内する `RuntimeError` が出て止まりました。Windowsも fork を使わない起動方式なので、`.py` として実行するなら必須と考えてください。 ## 収束診断と事後予測チェックの判定基準 MCMCの結果は、収束したかを数値で確かめてから読みます。MCMCそのものの仕組みと収束判定の考え方は[MCMC（マルコフ連鎖モンテカルロ法）とは？仕組み・MH法とHMCの違い・収束判定をPythonで解説](/tech/details/10270/)にまとめています。 ### R-hat・ESS・発散の3指標と実測値 Vehtari et al.（2021）は、順位正規化した R-hat が1.01未満のときだけ標本を使うこと、順位正規化した ESS が400を超えることを推奨しています。GelmanとRubin（1992）の1.1より厳しい基準で、1.1未満でも収束していない例が見つかったことが理由として挙げられています。ESSは分布の中心を見る `ess_bulk` と、区間の端を見る `ess_tail` の両方を確認します。 上のコードを実行した結果の先頭12行が次です。 ``` mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat mcse_mean mcse_sd a[0] 1.59 0.243 1.2 2 1462 2073 1.00 0.0064 0.0043 a[1] 0.904 0.154 0.67 1.2 1605 1750 1.00 0.0038 0.0022 a[2] 1.36 0.199 1.1 1.7 1103 1729 1.00 0.006 0.0034 a[3] 1.25 0.202 0.94 1.6 1566 1867 1.00 0.0051 0.0032 a[4] 0.821 0.125 0.63 1 1247 1295 1.00 0.0035 0.0021 a[5] 0.879 0.14 0.67 1.1 1217 1686 1.00 0.004 0.0024 a[6] 0.967 0.169 0.7 1.2 1442 1642 1.00 0.0044 0.0028 a[7] 1.47 0.212 1.1 1.8 1138 1635 1.00 0.0063 0.0039 b[0, 0] 0.577 0.118 0.4 0.78 1577 2173 1.00 0.003 0.0019 b[0, 1] 0.801 0.132 0.61 1 1621 2078 1.00 0.0033 0.0022 b[0, 2] 1.72 0.205 1.4 2.1 1574 1871 1.00 0.0052 0.0037 b[0, 3] 2.2 0.255 1.8 2.6 1638 2226 1.00 0.0064 0.0048 max r_hat 1.0038202526158047 min ess_bulk 1099.8272407270088 min ess_tail 1295.620991042581 divergences 0 ``` a・bの全40パラメータで R-hat は最大1.004、ESSは最小でも1,100、発散は0件で、3指標とも基準を満たしました。θ（300人分）も同じ設定の実行で R-hat最大1.006、ess\_bulk最小3,087でした。個人のθを解釈するときは `az.summary(idata, var_names=["theta"])` で同じ3指標を確かめます。`r_hat` 列は既定の自動表示で小数点以下2桁に丸めて表示されるため、1.01と表示された場合は実際の値が1.01未満かどうかを表示だけでは判断できません。上のように列の最大値を数値で出して確かめます。 同じモデルを draws=500・tune=500・target\_accept 既定（0.8）で回すと、「The rhat statistic is larger than 1.01 for some parameters」という警告が出ました。a・bは R-hat最大1.006で基準内でしたが、θは300人中3人で1.01を超え、最大1.012でした。θは1人あたりの回答が8件しかなく事後分布が広いため、反復を削ると先にθ側が基準を割ります。 ### 疑似データで確かめたパラメータの再現度 真値がわかっている疑似データなので、推定値との一致も確かめました。識別力aの事後平均と真値の相関は0.89、θの事後平均と真値の相関は0.87です。 | 項目 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | | ------- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | | a の真値 | 2.05 | 0.77 | 1.58 | 1.11 | 1.20 | 0.89 | 1.10 | 1.80 | | a の事後平均 | 1.59 | 0.90 | 1.36 | 1.25 | 0.82 | 0.88 | 0.97 | 1.47 | 項目0（真値2.05）と項目4（真値1.20）は、真値が89%区間（1.2〜2.0、0.63〜1.0）の外にありました。300人×8項目の規模では、識別力の点推定が真値から0.4程度ずれることがあります。項目の取捨を識別力の小数第1位で決めるような使い方は、この規模のデータでは避けるべきです。 ### 事後予測チェックでカテゴリ度数を照合する手順 収束しても、モデルがデータに合っているとは限りません。事後予測チェックでは、推定したパラメータから回答を生成し直し、観測した度数がその分布の中に収まるかを見ます。 ``` with grm: pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=5) yrep = idata.posterior_predictive["y"].values # (chain, draw, 回答数) mask = item == 0 obs = np.bincount(y[mask], minlength=K) rep = np.stack([np.bincount(r[mask], minlength=K) for r in yrep.reshape(-1, yrep.shape[-1])]) lo, hi = np.percentile(rep, [3, 97], axis=0) ``` 本番と同じ設定（draws=1000・target\_accept=0.9）で回した項目0の結果は次のとおりで、5カテゴリすべてで観測度数が事後予測の94%区間に収まりました。 ``` obs [200 15 47 14 24] lo [178. 7. 33. 7. 14.] hi [214. 26. 66. 26. 35.] ``` 実データで観測度数が区間から外れる項目があれば、カテゴリの使われ方がGRMの累積構造と合っていない（中間カテゴリだけ選ばれやすい等）か、多次元性が混じっている可能性を疑います。 ## 推定が壊れる失敗パターンと実測した症状 GRMをPyMCで組むとき、制約の入れ忘れはエラーか収束不全のどちらかで表面化します。同じ疑似データで、制約を1つずつ外して実行した結果です。 ### 境界の順序制約を外したときの初期化エラー bから `transform=pm.distributions.transforms.ordered` を外して実行すると、サンプリングが始まる前に止まりました。 ``` pymc.exceptions.SamplingError: Initial evaluation of model at starting point failed! Logp initial evaluation results: {'theta': np.float64(-321.61), 'a': np.float64(-11.74), 'b': np.float64(-53.38), 'y': np.float64(-inf)} ``` 初期値の境界が昇順になっていない項目があると、カテゴリ確率が負になり、観測の対数尤度 `y` が −inf になります。このエラー文で `y` だけが −inf なら、まず境界の順序を疑います。`model.debug()` を呼ぶと、どの値で計算が壊れたかを詳しく出せます。 ### 識別力に負の値を許したときの符号の不定性 aの事前分布を LogNormal から Normal(0, 1) に替えると、本記事のパラメータ化（cutpoints=a·b）では、上と同じ `SamplingError` で止まりました。aが負の項目では a·b が降順になり、OrderedLogistic の前提が崩れるためです。 より厄介なのは、cutpoints を直接パラメータにした傾き切片型（P(Y≥k)=logistic(aθ−ck)）で a〜Normal(0, 1) としたときです。aとθの符号を同時に反転しても尤度が変わらないため、エラーは出ずにサンプリングが完走し、チェーンごとに別の山へ収まりました。 ``` mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat a[0] 0 1.65 -1.9 1.9 6 39 1.73 a[1] 0 1.021 -1.2 1.2 6 40 1.73 ... chain means a[0]: [ 1.64 -1.63 -1.62 1.63] ``` 4本のチェーンのうち2本は a\[0\]≒+1.64、2本は≒−1.63 に張り付き、識別力aのR-hatの表示値は1.73〜1.74、ess\_bulkの表示値は全8項目で6でした。事後平均は0になり、「識別力0の項目」という誤った結論が出ます。R-hatが1.7前後に揃い、チェーン平均が正負で対称になっていたら、符号の不定性です。aを LogNormal や HalfNormal で正に縛れば解消します。 ### 疎なカテゴリと境界の接近 項目0の境界は b\[0,0\]=0.58 と b\[0,1\]=0.80 で、差は0.22しかありません。カテゴリ1の選択確率は2本の累積曲線の差であり、θがこの区間に入る人数とは一致しません。近接した境界によって選択確率が低くなり、この疑似データでの観測度数は300人中15人でした。境界がここまで近い項目は、その選択肢が回答者に区別されていない可能性があります。尺度を改訂するなら、隣のカテゴリとの統合候補になります。 あるカテゴリを誰も選ばなかった項目では、関連する境界の推定が不安定になり、事前分布への依存が強まることがあります。ただし、隣接カテゴリなどの観測も尤度を通じて境界を制約します。推定前に度数表を作り、度数0の原因を確認します。隣接カテゴリとの統合は選択肢の意味と尺度の比較可能性を確認して判断し、項目ごとにKが変わる場合は尤度の実装も変更します。 ## 推定結果の読み方：識別力・境界困難度・項目情報量 ### 識別力aと境界困難度bの解釈 識別力aは、θが変わったときに回答カテゴリがどれだけ鋭く切り替わるかを表します。今回の推定では項目0のaが1.59、項目4が0.82で、項目4はθの違いを回答に反映しにくい項目です。ただし前節のとおり区間が広いので、各項目の89%区間で不確実性を確認し、識別力の大小を判断するときは、同じ事後標本でa同士の差を計算して、その区間や差が正になる確率を確認します。 境界困難度bは、その境界を越える確率が0.5になるθの位置です。項目0は4本の境界が0.58〜2.20とθの高い側に寄っており、事後平均のa・bを代入すると、平均的な回答者（θ=0）が最低カテゴリを選ぶ確率は約71%です。観測度数で300人中200人がカテゴリ0だったのはこのためです。 ### 項目情報量で測定精度の高いθ帯を特定する計算 GRMの項目情報量は I(θ)=Σk Pk′(θ)2/Pk(θ) で、局所独立を仮定すると全項目の和がテスト情報量となり、その逆数の平方根は、項目パラメータを既知とする場合のθの最尤推定の標準誤差の漸近的な目安です。本記事で推定したθの事後標準偏差とは異なります。事後平均のa・bを入れて計算するコードです。 ``` a = idata.posterior["a"].mean(("chain", "draw")).values b = idata.posterior["b"].mean(("chain", "draw")).values grid = np.linspace(-3, 3, 13) def info(aj, bj): ps = 1 / (1 + np.exp(-aj * (grid[:, None] - bj[None]))) # P(Y>=k) Ps = np.hstack([np.ones((len(grid), 1)), ps, np.zeros((len(grid), 1))]) dPs = aj * Ps * (1 - Ps) P = Ps[:, :-1] - Ps[:, 1:] dP = dPs[:, :-1] - dPs[:, 1:] return (dP ** 2 / P).sum(1) test_info = sum(info(a[j], b[j]) for j in range(len(a))) se = 1 / np.sqrt(test_info) ``` | θ | −3 | −2 | −1 | 0 | 0.5 | 1 | 2 | 3 | | ------ | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | | テスト情報量 | 0.73 | 1.53 | 2.51 | 3.18 | 3.36 | 3.33 | 2.62 | 1.39 | | 標準誤差 | 1.17 | 0.81 | 0.63 | 0.56 | 0.55 | 0.55 | 0.62 | 0.85 | この8項目の尺度は θ=0.5〜1 付近で最も精度が高く、θ=−3では標準誤差が1.17と、θの標準偏差（1）を超えます。特性の低い層を見分けたい尺度なら、境界が低い側にある項目を足す必要があります。項目別に見ると、項目0の情報量は θ=1.5 で最大、項目2と項目6は θ=−0.5 で最大でした。短縮版を作るときは、測りたいθ帯に情報量の山を持つ項目を残します。 ## よくある質問 ### GRMとは何の略で、どんなデータに使いますか？ Graded Response Model（段階反応モデル）の略で、Samejima（1969）が提案した項目反応理論のモデルです。「全く当てはまらない」〜「非常に当てはまる」のように選択肢に順序がある質問紙や、論述を段階で採点したデータに使います。選択肢が2つならGRMは2PLと一致します。 ### ベイズ推定は何を推定しているのですか？ パラメータの事前分布と、データが得られる確率（尤度）を掛け合わせた事後分布を推定しています。GRMでは事後分布を式で解けないため、PyMCはNUTS（MCMCの一種）で事後分布から標本を引き、その平均や区間を推定値にします。考え方の土台は[ベイズの定理とは？公式の意味・具体例からナイーブベイズ実装まで](/tech/details/15236/)で解説しています。 ### 弁別力（識別力）とは何ですか？ θの違いを回答の違いとしてどれだけ鋭く反映するかを表すパラメータaです。aが大きい項目は、境界付近でθが少し変わるだけで選ぶカテゴリが切り替わります。GRMでは項目ごとに1つ推定し、本記事の推定では0.82〜1.59でした。 ### PyMCのNUTSサンプラーで発散（divergences）が出たらどうしますか？ まず `target_accept` を既定の0.8から0.9〜0.95へ上げます。それでも残る場合は、事前分布が広すぎないか、aを正に縛っているか、θの尺度を固定しているかを見直し、必要なら非中心化などの再パラメータ化をします。発散が1件でも残る結果は、そのまま解釈しないのが原則です。 ### 境界ごとの曲線とカテゴリごとの曲線は何が違いますか？ 境界ごとの曲線は「k以上を選ぶ確率」P(Y≥k|θ) で、5件法なら4本の右上がりのS字になります。カテゴリごとの曲線は「ちょうどkを選ぶ確率」P(Y=k|θ) で、5本のうち両端は単調、中間カテゴリは山型です。中間カテゴリがどのθでも最頻カテゴリにならない場合は、選択肢の機能を点検する材料になります。ただし、それだけで機能していないとは断定せず、選択肢の意味、観測度数、測定情報も確認します。 ## 関連記事 - [MCMC（マルコフ連鎖モンテカルロ法）とは？仕組み・MH法とHMCの違い・収束判定をPythonで解説](/tech/details/10270/) - [ベイズの定理とは？公式の意味・具体例からナイーブベイズ実装まで](/tech/details/15236/) - [尤度とは？読み方（ゆうど）・確率との違いから尤度関数・最尤推定・尤度比までわかりやすく解説](/column/details/6850/) - [モンテカルロ法とは？乱数で解く仕組みと円周率・リスク評価・強化学習での使い方【2026年版】](/tech/details/15242/) --- 出典: [段階反応モデル（GRM）をPyMCでベイズ推定する手順｜OrderedLogistic実装と収束診断]()（株式会社一創）