このChapterでは項目反応理論 （Item Response Theory [IRT]: Lord & Novick, 1968）のお話をしていきます。
後ほど説明しますが，実はIRTの最もメジャーなモデルの一つは，因子分析と数学的に同じものです。 そのためlavaanでも基本的なIRT分析はできるのですが，IRT専用のパッケージも色々あります1。 この講義ではmirtというパッケージを使用します。
↓本Chapterで使用するファイルのダウンロードはこちらから
8.1 項目反応理論とは
項目反応理論は，古典的テスト理論と同じく，テスト理論と呼ばれるカテゴリに属するものの一つです。 テスト理論の目的をとてもシンプルに言うと，以下の2つです。
- 回答者の潜在的な特性（心理特性の強さ・能力など）を測定する
- 測定に使用した項目（心理尺度の各項目・テストの問題など）の性能を評価する
1点目は直感的にもすぐ分かると思います。 2点目について少し補足しておくと，「項目の性能」とは，単純にはテストの問題の難易度などを指しています。 もしも難しすぎてほとんど誰も正解していない問題があったとしたら，その問題には，いわば床効果が出ているわけなのであまり良い問題とは言えません。 また，因子分析的に見たときに因子負荷が低い項目があったとしたら，その項目への回答は，測定したい「回答者の潜在的な特性」を反映していないので，注意が必要かもしれません。 テスト理論では，そのような形で項目の良し悪しを評価することで，「回答者の潜在的な特性」をより良く測定することを目指します。
そんなテスト理論の中で，項目反応理論が目指しているところを理解するために，仮想的な状況を考えてみましょう。 付録 A では，バーンアウトを例に心理尺度のお話をしました。 その時には，現在ではMBI(-GS)とUWESという2種類の尺度が用いられることが多いということを紹介しました。 そもそもの構成概念の定義の違いや項目数・回答方向の違いはありますが，いずれの尺度も0点から6点の7件法で構成されています。 ここでは，一旦構成概念が完全に同一とみなして，尺度得点として「平均点」を使う状況を考えてみます。 どちらの尺度についても0点が「最もやる気に満ち溢れた状態」，6点が「完全に燃え尽きた状態」だとみなす，ということです。 このような仮定のもとで，以下のAさんとBさんはどちらのほうがよりバーンアウト状態になっていると言えるでしょうか。
- Aさんは，MBI-GSによって測定されたバーンアウト得点が3.4点でした。
- Bさんは，UWESによって測定されたバーンアウト得点が4.2点でした。
ぱっと見では，Bさんのほうが得点が高いのでより燃え尽きかけているように見えますが，実はこの情報だけではそうとは言い切れません。 ポイントは2つの尺度（テスト）の平均値が同じとは限らないという点です。 もしもMBI-GSとUWESの回答者平均点がそれぞれ3.0点と4.5点だとしたら，Aさんは平均点以上，Bさんは平均点以下なので，Aさんのほうがよりバーンアウト状態に近いかも，という気がします。 このように，心理尺度によって得られる得点は尺度および項目に依存するもの（項目依存性）です。
項目依存性を避けて2つの尺度得点を比較可能にするためのナイーブなアイデアは，標準化得点を使うことです。 平均値を引いてから標準偏差で割った標準化得点ならば，平均値の違いなどを気にせず比較できそうです。 ということで先程の2人の標準化得点を計算してみたところ，それぞれAさんが0.2点，Bさんが-0.25点だったとしましょう。 2人はそれぞれ平均点以上・以下だったので，もちろん標準化得点でもAさんのほうが高い得点になっています。 ただ，実はまだこの情報だけではAさんのほうがよりバーンアウトに近いとは言い切れません。 それは，2人の所属するグループの平均値が同じとは限らないためです。
標準化得点は，当然ながらそのグループの得点分布に依存します。 もしかしたら，UWESの回答者平均点のほうがMBI-GSよりも高い理由は，UWESのほうが高得点になりやすい項目なのではなく，単にAさんの所属するグループよりもBさんの所属するグループの方が平均的にバーンアウト傾向が強いためかもしれません。 このように，心理尺度得点に基づく評価は項目だけでなく集団にも依存します（集団依存性）。
IRTでは，項目依存性と集団依存性をクリアした得点を算出するために，まず「集団によらない項目パラメータ」と「項目によらない個人特性値」の存在を考えています。 もちろんこれらのパラメータは，単純な和得点・平均値からはわかりません。 が，仮にそのようなパラメータがあるとして，項目パラメータが\(\beta\)，個人特性値が\(\theta\)で表せるとしたら，項目\(i\)に対する回答者\(p\)の反応（回答）は \[ y_{pi} = f(\theta_p,\beta_i), \tag{8.1}\] つまり\(\beta_i\)と\(\theta_p\)による何らかの関数（項目反応関数: item response function [IRF]）になるはずです。 項目反応関数の形はいろいろと考えられますが，例えば多くの心理尺度に対しては，\(\theta_p\)の値が大きいほど「当てはまる」寄りの回答をする（＝\(y_{pi}\)の値が大きくなる）と考えられます。 そのようなデータをIRTで分析する場合には，項目反応関数として「\(\theta_p\)に対して単調増加の関数」を持ってきたら良いわけです。 …という感じでIRTでは，観測されたデータにおいて想定される\(y_{pi},\theta_p,\beta_i\)の関係性をうまく表現できるような項目反応関数が色々と提案されてきました。 後ほどIRTの代表的なモデルを紹介しますが，IRTの本質はこのように，回答を(集団によらない)項目パラメータと(項目によらない)個人特性値に分離して分析することにあります2。
ということで，一般的なIRTの利用場面では，
- 特定の関数形（IRF）を仮定する
- そのモデルのもとで項目パラメータを推定する
- 推定した項目パラメータをもとに，項目の性能を評価する
- 合わせて個人特性値を推定する
- 推定した個人特性値を使ってあとはご自由に
といったことが行われます。 それでは，まずは「とりあえずIRTで分析」ができるようになるために必要な最低限の知識を学んでいきましょう。
8.2 基本的なモデル
最もシンプルなIRTモデルを考えるため，一旦ここからは項目への反応を二値（当てはまらない・あてはまる）に限定して考えてみます。 というのも，もともとIRTは「正解・不正解」の二値で考えるような学力テストなどを想定して作られているもののためです。
8.2.1 正規累積モデル
以前カテゴリカルデータの相関（ポリコリック・ポリシリアル相関）の計算時には，背後にある連続量を考えるという説明をしました。 ここでも，二値反応\(y_{pi}\)の背後に連続量\(\theta_p\)が存在し，\(\theta_p\)が閾値を超えていたら1，超えていなければ0と回答すると考えます。 そして項目の閾値を\(\beta_i\)と表すと， \[ y_{pi} = \left\{ \begin{aligned} 0 & \quad \theta_p < \beta_i\\ 1 & \quad \theta_p \geq \beta_i \end{aligned} \right. \tag{8.2}\] となります（ 図 8.1 ）。 \(\beta_i\)の値が大きい項目ほど，\(y_{pi}=1\)になる\(\theta_p\)の範囲が狭くなるため，「当てはまらない」と回答する人の割合が高くなるだろう，ということです。
ただし，このままでは能力\(\theta_p\)を持つ\(p\)さんは難易度が\(\beta_i\)以下の問題には，どんなジャンル・内容の問題でも100%正解することになってしまいます（\(\theta_p \geq \beta_i\)なので）。 しかし，もちろん客観的に見た難易度が全く同じであっても，内容によって\(p\)さんが答えられる・答えられない問題があるはずです。 さらに言えば，\(p\)さんの体調や時の運（たまたま昨日一夜漬けした内容に入っていた，など）によっても答えられるかどうかは変わるはずです。 というように，心理尺度やテストへの回答では，常に「その他の要因」による変動を考えてきました。 言い換えると，古典的テスト理論でも回帰分析でも因子分析でも，観測値は必ず何らかの誤差を伴って発生するものだと考えているわけです。 そこで先程の(8.2)式に誤差項\(\varepsilon_{pi}\)を追加して \[ y_{pi} = \left\{ \begin{aligned} 0 & \quad (\theta_p - \varepsilon_{pi}) < \beta_i\\ 1 & \quad (\theta_p - \varepsilon_{pi}) \geq \beta_i \end{aligned} \right. \quad \varepsilon_{pi} \sim N(0,\sigma_{\varepsilon}^2) \tag{8.3}\] とします。なお誤差項は，回帰分析などと同じく，最も一般的な「平均0の正規分布に従う」と仮定してあります。
誤差項を追加することによって，「特性値が\(\theta_p\)の人が項目\(j\)に『当てはまる』と回答する確率3」を考えることが出来ます。 図 8.1 では「全ての人の\(\theta_p\)」の分布を示していた一方で， 図 8.2 では「\(\theta_p\)が特定の値の人」における分布になっている点に注意してください。
さらに\(\beta_i\)と\(\theta_p\)を一つの関数\(f(\theta_p,\beta_i)\)として見るために\(\beta_i\)を移項して， \[ y_{pi} = \left\{ \begin{aligned} 0 & \quad (\theta_p- \beta_i - \varepsilon_{pi}) < 0\\ 1 & \quad (\theta_p-\beta_i - \varepsilon_{pi}) \geq 0 \end{aligned} \right. \quad \varepsilon_{pi} \sim N(0,\sigma_{\varepsilon}^2) \tag{8.4}\] とします（ 図 8.3 ）。 当然ながら， 図 8.2 と 図 8.3 は，閾値と分布が平行移動しただけなので\(P(y_{pi}=1)\)は変わりません。
図 8.3 の色つきの確率分布は，正規分布\(N(\theta_p - \beta_i,\sigma_{\varepsilon}^2)\)と表すことができます。 そして式から，\(y_{pi}=1\)となるのは誤差\(\varepsilon_{pi}\)が\((\theta_p - \beta_i)\)以下のときだとわかります。 したがって，仮に\(\varepsilon_{pi}\)が標準正規分布に従う（\(\sigma_{\varepsilon}^2=1\)）としてみると，「特性値が\(\theta_p\)の人が項目\(j\)に当てはまると回答する確率」（ 図 8.4 の青い部分）は \[ \begin{aligned} P(y_{pi}=1)&= P(\varepsilon_{pi} \leq \theta_p - \beta_i)\\ &= \int_{-\infty}^{(\theta_p - \beta_i)}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{2}\right)dz \end{aligned} \tag{8.5}\] と表すことができます。教科書などで見たことがある形かもしれませんが，この積分の中身は標準正規分布の確率密度関数です。 このように，\(P(y_{pi}=1)\)を標準正規分布の累積確率で表すモデルのことを正規累積モデル(normal ogive model)と呼びます。 この後項目パラメータを1つ増やすので，それに対してこの正規累積モデルは，特に1パラメータ正規累積モデルと呼ばれます。
ここまでは，途中で話を簡単にするために勝手に\(\sigma_{\varepsilon}=1\)として進めていましたが，ここからはもう少し一般化して\(\sigma_{\varepsilon}^2\)が項目ごとに異なる(\(\sigma_{\varepsilon}^2=1/\alpha_i\))とします。 このとき，誤差の確率分布は\(\varepsilon_{pi} \sim N\left(0,\frac{1}{\alpha_i}\right)\)となるため，両辺を\(\alpha_i\)倍して\(\alpha_i\varepsilon_{pi} \sim N\left(0,1\right)\)と表すことができます。 \(\alpha_i\varepsilon_{pi}\)が標準正規分布に従うならば，(8.5)式の正規累積モデルは \[ \begin{aligned} P(y_{pi}=1)&= P\left[\varepsilon_{pi} \leq (\theta_p - \beta_i)\right]\\ &= P\left[\alpha_i\varepsilon_{pi} \leq \alpha_{i}(\theta_p - \beta_i)\right]\\ &= \int_{-\infty}^{\alpha_{i}(\theta_p - \beta_i)}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{2}\right)dz \end{aligned} \tag{8.6}\] という形で書き換えることができます。また正規累積モデルが出てきました。(8.5)式との違いは積分の上限が\(\alpha_{i}(\theta_p - \beta_i)\)に変わった点だけです。 いま，突拍子もなく\(\alpha_{i}\)というパラメータを導入しましたが，とりあえずこれによって項目\(i\)が持つパラメータは\((\beta_i, \alpha_i)\)の2つになりました。 ということで，(8.6)式の正規累積モデルは2パラメータ正規累積モデルと呼ばれます。
8.2.2 項目特性曲線
ここで，2つの項目パラメータ\((\beta_i, \alpha_i)\)が何を表しているのかを明らかにしていくため， 図 8.3 を別の形で表してみます。
図 8.5 の下半分は，横軸に\(\theta_p - \beta_i\)を，縦軸に\(P(y_{pi}=1)\)をとったグラフです。 図の上半分にあるように，\(\theta_p - \beta_i\)の値が大きくなるほど分布の青い部分が増える，つまり\(P(y_{pi}=1)\)の値が大きくなることを表しています。
これを踏まえて，パラメータ\(\beta_i\)の挙動を見るために，横軸に\(\theta_p\)を，縦軸に\(P(y_{pi}=1)\)をとったグラフを，\(\beta_i\)の値を変えながらいくつか重ねてみます（ 図 8.6 ）。 誤差は平均0の正規分布に従うため，\(P(y_{pi}=1)=0.5\)になるのは，\(\theta_p = \beta_i\)のときです。 したがって\(\beta_i\)の値が大きい項目ほど，\(\theta_p\)の値が大きくないと「当てはまる」と回答する確率が高くなりません。 そのような意味で，\(\beta_i\)は項目困難度(item difficulty)と呼ばれます4。
次は\(\alpha_i\)の役割を見るために， 図 8.6 よりも誤差分散(\(\alpha_i\))が小さい場合や大きい場合を見てみましょう。 図 8.7 の左側に示されているように，誤差分散が小さいほど分布の広がりがなくなるため，\(\theta_p - \beta_i\)の値の変化に対する\(P(y_{pi}=1)\)の変化が急激になっています。
これを踏まえて，\(\beta_i\)の時と同じように横軸に\(\theta_p\)をとり，縦軸に\(P(y_{pi}=1)\)をとったグラフを，今度は\(\alpha_i\)の値を変えながらいくつかプロットしたものが 図 8.8 です（\(\beta_i=0\)）。 誤差分散の逆数である\(\alpha_i\)の値が大きい項目ほど，\(\theta_p\)が低い人と高い人の間で\(P(y_{pi}=1)\)が大きく変化しています。 言い換えると\(\alpha_i\)の値が大きい項目ほど，その項目への回答によって\(\theta_p\)の高低を識別する能力が高いといえます。 ということで\(\alpha_i\)のことを項目識別力(item discrimination)と呼びます。
図 8.6 や 図 8.8 のように，横軸に\(\theta_p\)をとり，縦軸に\(P(y_{pi}=1)\)をとったグラフは，いわば回答者の特性値に対する項目の特性を表しています。 ということで，これらのグラフのことを項目特性曲線(item characteristic curve [ICC])と呼びます。
8.2.3 ロジスティックモデル
2パラメータ正規累積モデルの式を再掲しておきましょう。 \[ \begin{aligned} P(y_{pi}=1)&= \int_{-\infty}^{\alpha_{i}(\theta_p - \beta_i)}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{2}\right)dz \end{aligned} \] 正規累積モデルには，積分計算が含まれています。 ある程度想像がつくかもしれませんが，コンピュータは積分計算が比較的得意ではありません。 ということで，積分計算がいらない別のモデルを考えてみましょう。
正規累積モデルがやっていることは，プロビット回帰と同じです。 そんな正規累積モデルの代わりになるものということは，正規分布のように
- 累積分布関数がS字になっていて
- 確率密度関数は左右対称の山型
な分布があれば良いわけですが…そのような確率分布として，ロジスティック分布というものがありました。 ロジスティック分布の確率密度および累積分布関数は \[ f(x) = \frac{\exp(-x)}{\left[1+\exp(-x)\right]^2}, \ \ \ \ \ F(x) = \frac{1}{1+\exp(-x)} \tag{8.7}\] と表され，確率密度関数を正規分布と重ねて見ると 図 8.9 のようになります。 こうしてみると，裾の重さは違いますが確かにロジスティック分布は正規分布のような左右対称形になっているようです。
また，ロジスティック分布の\(x\)をおよそ1.7倍すると，累積分布関数が正規分布とかなり近くなることが知られています( 図 8.10 )5。
ということで，正規分布の代わりにロジスティック分布を使った(2パラメータ)ロジスティックモデルは \[ \begin{aligned} P(y_{pi}=1)&= \frac{1}{1+\exp\left[-1.7\alpha_{i}(\theta_p - \beta_i)\right]} \end{aligned} \tag{8.8}\] と表されます。 いま説明したように，この1.7はあくまでも「正規累積モデルに近づけるための定数」なので，特に正規累積モデルとの比較をする予定が無いならば1.7は無くても問題ありません。 ただ，結果を報告する際には「正規累積モデルとロジスティックモデルのどちらを使用したのか」「ロジスティックモデルならば1.7倍した値なのか」がハッキリわかるようにしておきましょう。 論文などでは，モデル式を書いておけば一発です。
（以後の説明では，特に正規累積モデルと比較することもないので1.7は省略していきます）
また，1パラメータロジスティックモデルはラッシュモデル(Rasch model)と呼ばれることもあります。 これは，IRTとは独立にRaschという人が理論的に「1パラメータロジスティックモデルと同じ関数」を導出した （Rasch, 1960/1980） ことに由来します。 ということで「ラッシュモデル」にはやや高尚な理念があったりしますが，私もよく理解できていないのでそこの説明はしません。 とりあえず「ラッシュモデル」という名前が出てきたら「困難度だけのロジスティックモデルなんだな」と理解できていれば良いでしょう。 厳密な使い分けとしては，ラッシュモデルは「正規累積モデルの近似」を目的としていないので，1.7倍しない1パラメータモデルを指す，という場合もあるようです。
8.2.4 モデルの使い分け
正規累積モデルとロジスティックモデルのどちらを使用するかは，どちらでも問題ありません（ほぼ同じ関数なので）6。 コンピュータが今ほど発達していなかった昔では，計算にかかる時間がクリティカルな問題であったなどの理由でロジスティックモデルのほうが良かったかもしれませんが，現代のコンピュータの性能であればその差は無視してよいレベルでしょう。 とはいえ，今でも計算の容易さからロジスティックモデルを使用するケースのほうが多いような気がします。
パラメータ数に関しては，考え方の違いによるところが大きい気がします。 ざっくりと1パラメータモデルと2パラメータモデルを比べると，
- 1パラメータモデルのほうがサンプルサイズが少なくても推定が安定する
- 2パラメータモデルのほうが推定の正確さ（モデル適合度）が高くなりやすい
といえます。 また，項目特性曲線的に考えてみると，識別力の異なる2つの項目間に関して，「特性値が低い人では項目Aのほうが当てはまりやすい一方で，特性値が高い人では項目Bのほうが当てはまりやすい」といった逆転現象のようなものが起こります（ 図 8.8）。 困難度パラメータは「その項目の難しさ」を表すと言ったものの，識別力が異なる2パラメータモデルではシンプルな「難しさ」として考えるのは案外難しいのかもしれません。
それに対して1パラメータモデルでは，項目特性曲線どうしが交差することは決してない（ 図 8.6 ）ので，どの項目間でも当てはまりやすさの大小関係が変わることはありません。 したがって，困難度パラメータを文字通りの「難しさ」として捉えることができます。
…という感じでいろいろ御託を並べてみましたが，基本的には迷ったら2パラメータ（ロジスティック）モデルを選んでおけばとりあえずは安泰です。 後述しますが，2パラメータモデルは因子分析モデルと数学的に等価である一方で，1パラメータモデルは「因子負荷をすべて1に固定した因子分析モデル」に相当します。 そう考えると「2パラメータモデルでいっか」という気持ちになりますね。
実は世の中にはもっとパラメータ数の多いモデルが存在します。
- 【3パラメータ】
- 3パラメータモデルでは，2パラメータモデルに加えて「当て推量」のパラメータ\(c_i\)が追加されます。 \[ P(y_{pi}=1)= c_i + \frac{1-c_i}{1-\exp\left[\alpha_i(\theta_p-\beta_i)\right]} \tag{8.9}\] \(c_i\)は項目特性曲線の下限を操作します。\(\theta_p\)がどんなに低い人でも，必ず\(P(y_{pi}=1)\)は\(c_i\)以上の値になります。 例えば4択問題であれば，どんなに能力が低くても，勘で答えれば正答率は1/4になります。これが\(c_i\)です。
- 【4パラメータ】
- 4パラメータモデルでは，3パラメータモデルに加えて「上限」のパラメータ\(d_i\)が追加されます。 \[ P(y_{pi}=1)= c_i + \frac{d_i-c_i}{1-\exp\left[\alpha_i(\theta_p-\beta_i)\right]} \tag{8.10}\] 4パラメータモデルでは，\(\theta_p\)がどんなに高い人でも，必ず\(P(y_{pi}=1)\)は\(d_i\)以下の値になります。 どんなに能力がある人でも100%にはならないもの，例えばバスケのフリースローなどをイメージすると良いかもしれません。
- 【5パラメータ】
- 5パラメータモデルでは，4パラメータモデルに加えて「非対称性」のパラメータ\(e_i\)が追加されます。 \[ P(y_{pi}=1)= c_i + \frac{d_i-c_i}{\left(1-\exp\left[\alpha_i(\theta_p-\beta_i)\right]\right)^{e_i}} \tag{8.11}\] 通常のIRTモデルでは，\(P(y_{pi}=1)\)の動き方は0.5を軸に対称になっています。 これに対して，例えば「最初は急激に成功率が上がるが，上に行くほどきつくなる」的なことを表現したい場合には非対称性を考慮する必要があるでしょう。
ただ，これらのモデルは安定した推定に必要なサンプルサイズが莫大であったり，そもそも推定量に一致性がないなどの問題を抱えていることから，実務場面で使われることはほとんどありません（強いて言えば，当て推量パラメータを「選択肢数分の1」に固定した3パラメータモデルなどはありえる）。 ので，「そういうのもあるんだなぁ」くらいに思っておけば良いと思います。
8.3 因子分析との関係
IRTの基本的なモデルが分かったところで，IRTモデルと因子分析の特定のモデルが数学的には等価であること （Takane & de Leeuw, 1987） を紹介しておきたいと思います。 このことは，心理尺度の分析において因子分析だけでなくときにはIRTも使用できることを意味しています。 それぞれの分析手法は出自が異なる以上，分析の焦点や技法にも違いがあるのですが，適切な設定のもとでは，因子分析とIRTのいいとこ取り（あるいは併用）ができるということは理解しておくと良いかもしれません。
多母集団同時分析のところ （ セクション 7.9 ） で説明したように，標準化していない（＝切片が0でない）因子分析モデル（1因子モデル）は \[ y_{pi} = \tau_i + f_pb_i - \varepsilon_{pi} \qquad \varepsilon_{pi}\sim N(0,\sigma_{\varepsilon}) \tag{8.12}\] という形になっています7。 そしてカテゴリカル因子分析では，観測変数\(y_{pi}\)はその背後にある連続量によって決まる，という考え方をしていました。 心理尺度の項目が2件法だとしたら，観測変数\(y_{pi}\)とその背後の連続量の関係は 図 8.3 と同じようにして 図 8.11 のように表すことができます。
したがって，標準正規分布に従う誤差\(\varepsilon_{pi}\)が\(\tau_i + f_pb_i\)より小さいときに\(y_{pi} = 1\)となるので，その確率は \[ \begin{aligned} P(y_{pi}=1)&= P(\varepsilon_{pi} \leq \tau_i + f_pb_i)\\ &= \int_{-\infty}^{(\tau_i + f_pb_i)}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{2}\right)dz \end{aligned} \tag{8.13}\] となります。これを2パラメータ正規累積モデル（ 8.6 式）に対応させて見ると， \[ \begin{aligned} \tau_i + f_pb_i &= b_i\left(f_p + \frac{\tau_i}{b_i}\right) \\ &= \alpha_{i}(\theta_p - \beta_i) \end{aligned} \tag{8.14}\] ということで，\((f_p, b_i, \tau_i) = (\theta_p, \alpha_{i}, \alpha_{i}\beta_i)\)としてみればこれら2つのモデルは完全に同一だということがわかります。
IRTでは，(8.14)式に基づく（因子分析的な）パラメータ化を行うことがあります。 つまり2パラメータロジスティックモデルで言えば \[ \begin{aligned} P(y_{pi}=1)&= \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_i)\right]} \\ &= \frac{1}{1+\exp\left[-(\alpha_{i}\theta_p - \tau_i)\right]} \end{aligned} \tag{8.15}\] という形で，困難度パラメータの代わりに切片パラメータを求めるということです。 これもどちらの設定を使用しても良いのですが，ソフトウェアによってはこの切片パラメータを使用する方法がデフォルトになっていることもあるので，出力を見る際は気をつけましょう。 なおIRTの枠組みでは，どちらかというと困難度パラメータ\(\beta_i\)を用いる解釈のほうがよく使われます。 推定結果で切片パラメータが出てきた場合には，自分で\(\beta_i = \frac{\tau_i}{\alpha_i}\)と変換する必要があるかもしれません。
8.4 RでIRT
IRTの基本を学んだところで，いよいよ実際にRでIRTを実行してみます。 これまでの説明に合わせて，まずは二値で試してみるため，それ用のデータを作成します。 元々のdatは6件法への回答データでしたが，各選択肢は1から3が「当てはまらない」寄り，4から6が「当てはまる」寄りでした。 そこで，「どちらかと言えば当てはまると回答したかどうか」の二値データに変換したものを使用していきます。
そしてIRTでの分析には，mirtパッケージの中にあるmirt()関数を使用します。 基本的な引数は以下の3つです。他にも色々ありますが，後ほどもう少し紹介します。
data-
推定に用いるデータです。
lavaanの関数と異なり，推定に使う変数のみが入ったデータを用意する必要があります。 model-
lavaanのように各変数と因子の関係を記述したモデル式です。先程説明したように，IRTは因子分析と同じものなので，モデル式を記述することで多次元のIRTモデルも実行可能になる，というわけです（詳細は後ほど）。変数名が連番になっている場合，後述のようにハイフンでまとめて指定できるのでおすすめです。 itemtype- モデルの指定です。たぶん正規累積モデルは選べず，すべてロジスティック関数ベースのモデルです。これまでに紹介したモデルでは
'Rasch'：ラッシュモデル，つまり1パラメータモデルです。後述のものに合わせて'1PL'と指定することはできないので気をつけてください。'2PL'：2パラメータモデルです。同じように'3PL', '4PL'も指定可能です。
得られた推定値を見る場合には，coef()関数を使います。
$Q1_1
a1 d g u
par 0.416 1.262 0 1
$Q1_2
a1 d g u
par 1.023 2.359 0 1
$Q1_3
a1 d g u
par 1.027 1.911 0 1
$Q1_4
a1 d g u
par 1.033 1.737 0 1
$Q1_5
a1 d g u
par 1.01 1.791 0 1
$Q1_6
a1 d g u
par 1.203 1.911 0 1
$Q1_7
a1 d g u
par 1.367 1.667 0 1
$Q1_8
a1 d g u
par 1.225 1.582 0 1
$Q1_9
a1 d g u
par 1.339 1.365 0 1
$Q1_10
a1 d g u
par 1.27 0.031 0 1
$GroupPars
MEAN_1 COV_11
par 0 1デフォルトでは，1項目ごとのlist型で返すため少し見にくく，また今後の取り扱いがやや面倒になりがちです。 引数simplify = TRUEとすることで，データフレームの形で返してくれるようになります。
最初の$item部分に表示されているものが項目パラメータです。左から
a1- 識別力パラメータの推定値
d- 切片パラメータの推定値
g- 当て推量パラメータ (guessing) の推定値（2パラメータモデルなら全て0）
u- 上限パラメータ (upper bound) の推定値（2パラメータモデルなら全て1）
となっています。このように，mirtはデフォルトでは切片パラメータを出力します。困難度パラメータを見たい場合にはd列をa1列で割るか，引数IRTpars = TRUEをつけて出力しましょう。
これでd列の代わりにb（困難度パラメータ）列が表示されるようになりました。 出力を見ると，例えばQ1_1では，\(1.262 / 0.416 \simeq 3.031\)ということで，確かにd列をa1列で割った値（の符号を入れ替えたもの）が表示されています。
$itemの下に表示されている$meansと$covはそれぞれ特性値（因子得点）\(\theta\)の平均値と分散です。 この時点では一般的な因子分析の制約（平均0，分散1）に沿った値のみが表示されているだけなので特に注目する必要はありません。 後ほど紹介する多次元モデルや多母集団モデルになると，ここに表示される値も重要な意味を持つことになります。
plot()関数を使うと，引数typeに応じて，推定結果に基づく様々なプロットを全項目まとめて出すことができます。 図 8.6 のようなICCを出したい場合にはtype = "trace"としてあげてください。
デフォルトではすべての項目のICCを縦横に並べてくれるのですが，項目数が多くなると見るのが大変そうです。 この場合，引数facet_items = FALSEとすると，すべての項目のプロットを重ねて表示してくれます。 また引数which.itemsを指定すると，特定の項目だけのICCを出すことも出来ます。
ということで 図 8.12 や 図 8.13 を見ると，識別力がダントツで低い項目1では，ICCの傾きが緩やかになっていることがわかります。
IRTによる項目分析では，まずこのように推定されたパラメータをもとにICCを描くのが重要です。これは 図 4.16 と同じように解釈することができるので，例えば
- \(\theta\)の広い領域でICCがほぼ上か下に貼り付いている（\(\beta_i\)がとても小さいか大きい）
- ICCの傾きが緩やかである（\(\alpha_i\)がとても小さい）
ような項目は，\(\theta\)の測定という面からはあまり良くない項目であろう，ということがわかります（除外するかは要検討）。 この他にも，複数の項目を組み合わせたときにどうなるかという分析などもあるわけですが，それはおいおい紹介していきます。
話は変わって，項目パラメータを固定した上での特性値\(\theta_p\)の推定値を出す場合は，fscores()関数を使います。 このあたりは因子分析やCFAでの「項目パラメータを推定」→「項目パラメータを固定して個人特性値を推定」と全く同じです。
8.5 IRTの前提条件
IRT，なかでも2パラメータモデルは基本的には因子分析と同じようなものなのですが，理論的背景の違いなどから，因子分析を行う際とは留意すべき点も少し異なってきます。 ここでは，IRTを適用するにあたって重要となる2つの前提条件について見ていきます。
8.5.1 局所独立性
局所独立性(local independence)の仮定は，因子分析のときと概ね同じです。 探索的因子分析では，独自因子の間が完全に無相関であるという仮定をおき，その中でデータの相関行列を最もよく復元できる因子負荷行列と独自因子の分散(\(\symbf{B}^\top\symbf{B} + \symbf{\Psi}\))を求めていました。 独自因子の間に相関があるということは，その当該項目の間に何か別の共通因子があるということです。
局所独立性をもう少し固い言い方で説明すると，「\(\theta_p\)で条件づけたとき，項目\(i\)と項目\(j~(i\neq j)\)への回答は完全に独立」となります。 項目への回答は二値なので，\(2\times2\)クロス表をイメージしてみましょう。 \(\theta_p\)が特定の値の人を1000人ほど集めて（つまり\(\theta_p\)で条件づけて） 表 8.1 のように2つの項目\(i,j\)の回答のクロス表を作ると，もし局所独立性が満たされているならば連関はゼロになります。 項目\(i\)に「当てはまる」と回答した人を見ても，「当てはまらない」と回答した人を見ても，項目\(j\)に「当てはまる」と回答した人の割合は全体の割合と同じ6:4になっています。
| 項目\(i\) | 当てはまらない | 当てはまる | 計 |
|---|---|---|---|
| 当てはまらない | 180 | 420 | 600 |
| 当てはまる | 120 | 280 | 400 |
| 計 | 300 | 700 | 1000 |
一方で，\(\theta_p\)以外の別の共通因子がある場合，クロス表に連関が現れます。 例えば「他者への助けを求める傾向」の尺度において，宗教観の影響を受けると思われる項目（e.g., 困ったときは神に祈る）どうしのクロス表を見てみると， 表 8.2 のように
- 項目\(i\)に「当てはまる」と答えた人では\(340/400=85\%\)の人が項目\(j\)でも「当てはまる」と回答
- 項目\(i\)に「当てはまらない」と答えた人では\(360/600=60\%\)の人が項目\(j\)でも「当てはまる」と回答
しているため，項目\(i\)の回答と項目\(j\)の回答が関連を持っていることになります。 項目\(i\)に「当てはまる」と答えた人のほうが項目\(j\)でも「当てはまる」と回答する確率が高い，ということです。
| 項目\(i\) | 当てはまらない | 当てはまる | 計 |
|---|---|---|---|
| 当てはまらない | 240 | 60 | 600 |
| 当てはまる | 360 | 340 | 400 |
| 計 | 300 | 700 | 1000 |
局所独立の仮定とは，このようにクロス表を作成したときに，どの項目ペア間でも，どの特性値の値でも連関がない，ということを指しています8。
局所独立性がなぜ重要なのでしょうか。その理由の一つは計算上の問題です。 IRTでは基本的に最尤法によってパラメータを推定します。 もしも局所独立性が満たされていれば，\(\theta_p\)で条件づけたときの各項目に対する反応確率は独立なので，回答者\(p\)に対する尤度関数は単純にすべての項目をかけ合わせた \[ P(\symbf{y}_p|\theta_p) = \prod_{i=1}^{I}P(y_{pi}=1|\theta_p)^{y_{pi}}P(y_{pi}=0|\theta_p)^{1-y_{pi}} \tag{8.16}\] と表せます9。 これに対してもしも局所独立が満たされていない場合（ 表 8.2 ），\(P(y_{pj}=1|\theta_p)\)の値が\(y_{pi}\)によって変わることになります。 そのため単純に掛け算ができなくなってしまい，もっと複雑な計算が必要になってしまうのです。
ということで，IRTモデルは「局所独立性が満たされている」という仮定にもとづいて定式化されています。 したがって，データを収集する時点で局所独立性が満たされていることをしっかりと確認しておく必要があります。 なお後半では，局所独立性が満たされているかを事後的に評価する指標を紹介します。
8.5.2 一次元性
これまで紹介してきたIRTモデルはいわば一因子の因子分析モデルです。 因子分析では「項目をグルーピングする」ことに重きを置きがちだったので，因子数は「いい感じにグループができる数」にすることが一般的でした。 これに対してIRTはもともと単一のテストに対する単一の能力を測定することが主たる目的と言えるため，基本的には因子数は1となります10。 そういう意味では，心理尺度とは相性が良い側面もあると言えそうです。 そんなわけで，IRTの重要な前提条件としてすべての項目が同じ能力・特性のみを反映していること（一次元性）が必要となります。
一次元性を検証するためのフレームワークは，多くの場合（探索的・確認的）因子分析と共通です。 つまり， セクション 6.11 で紹介したスクリープロットや平行分析，MAPなどの基準によって「因子数1」が提案されたら良い，ということです。 ですが，因子分析的な次元性の確認方法では，項目数が増えるとなかなか「因子数1が最適」とはならないケースが出てきます。 したがってIRTでは，「一次元とみなしても問題ないか」といった別の視点から一次元性の評価を行うことがあります。 これは，例えば一因子のCFAを適用した際の適合度指標（RMSEAやCFA, AGFIなど）が許容範囲にあるかによって確認することが出来ます。
またIRTでは，局所独立性の観点から一次元性を確認する指標もいくつか提案されています。 一次元性が満たされているということは，回答行動に影響する共通因子が一つ（\(\theta_p\)）だけということです。 したがって一次元性が満たされていない場合には， セクション 6.4.1 （ 図 6.9 ） で少し説明したように，共通因子で説明出来ない部分（独自因子）の間に共分散が発生することになります。 そのように考えると，一次元性は局所独立と同じようなものだとみなすことができるわけです。 厳密には局所独立性が満たされている場合，一次元性も満たされていると言って良いだろう，ということになります11。 ということで，より洗練された一次元性確認の方法として，例えばDIMTEST（Stout et al., 1996）やDETECT（Zhang & Stout, 1999）という名前の方法が提案されていたりします。
一次元性が満たされているだけでは局所独立性が満たされるとは限りません。 例えば
- 変数\(x\)の平均値を求めよ
- 変数\(x\)の分散を求めよ
という2つの問題があった場合，この2つはいずれも「統計の知識」を問うているという意味では同じ潜在特性に関する項目であり，一次元性を満たしているといえます。 ですが，分散を計算する過程では必ず平均値を計算する必要があるため，項目1に間違えた人はほぼ確実に項目2に間違えます。 したがって，この2問は一次元性は満たされているが局所独立性は満たされていない状態だといえます。 このように，ある項目への回答が他の項目への回答に直接的に影響を与える場合を実験的独立性の欠如と呼びます。 心理尺度でいえば，
といった指示を出した場合には，実験的独立性が欠如していると言えそうです。 他には，回答時間が足りない場合の後ろの方の項目も，「前の項目に回答出来ていない人は次の項目にも回答できていない」ため，実験的独立性が欠如しているとみなすことが出来ます。
ということで局所独立性を正確に表すと，一次元性に加えてこの実験的独立性が必要といえます （南風原，2000; Lord & Novick, 1968） 。
8.6 多値型モデル
ここまで，二値の項目反応データに対する基本的なモデルを紹介し，Rでの実行方法を示しました。 しかし実際の場面では，リッカート式尺度ならばたいてい5件法から7件法の多値項目です。 そしてテストの場面でも，部分点や組問12を考えると多値型のモデルが必要となるでしょう。 ということで，代表的な2つの多値型モデルを紹介したいと思います。 どちらを使っても本質的には違いは無いので，説明を読んでしっくり来たほうを使えば良いと思います。
8.6.1 段階反応モデル
段階反応モデル （Graded response model [GRM]: Samejima, 1969）では，複数の二値IRTモデルを組み合わせて多値反応を表現します。 項目\(i\)に対する回答者\(p\)の回答\(y_{pi}=k~(k=1,2,\cdots,K)\)について「\(k\)以上のカテゴリを選ぶ確率」を考えると，これはまだ二値（\(k\)以上 or \(k\)より小さい）なので，2パラメータロジスティックモデルならば \[ P(y_{pi}\geq k) = \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_{ik})\right]} \tag{8.17}\] と表せます。なお，困難度パラメータは「項目\(i\)のカテゴリ\(k\)」ごとに用意されるため\(\beta_{ik}\)としています。
異なるカテゴリを選ぶ確率は当然排反な事象なので，二値モデルを組み合わせると「ちょうど\(k\)番目のカテゴリを選ぶ確率」は \[ P(y_{pi} = k) = P(y_{pi}\geq k) - P(y_{pi}\geq k+1) \tag{8.18}\] と表すことができます。 ただし，端のカテゴリに関する確率はそれぞれ \[ \begin{aligned} P(y_{pi}\geq 1) &= 1 \\ P(y_{pi}\geq K+1) &= 0 \end{aligned} \tag{8.19}\] とします。 「1以上のカテゴリを選ぶ確率」は当然1になるため，これに対する困難度パラメータ\(\beta_{i1}=-\infty\)となります。 同様に，「\(K+1\)以上のカテゴリを選ぶ確率」は0になります。\(P(y_{pi} = K|\theta_p)\)を考えるときだけ仮想的な「さらに上のカテゴリ」を想定する，ということです。 このようにGRMは二値型モデルの素直な拡張のため，数学的には正規累積型のGRMはカテゴリカル因子分析と同じものだとみなすことができます （Takane & de Leeuw, 1987） 。
図 8.14 に，4件法の項目に対する段階反応モデル（\((\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(1,[-1.5, 0, 2])\)）のイメージを示しました。 左右で色の塗り方や分布の位置は変えていません。点線の位置だけです。
左は，\(\theta_p\)の位置に応じて各カテゴリの反応確率が変化する状況を表しています。 反応確率は全カテゴリ合計で1となっているため，点線の位置での各カテゴリの長さの比率がそのままカテゴリ選択確率となっています。
右は，\(\beta_{ik}\)の位置に点線を引いています。 カテゴリの困難度\(\beta_{ik}\)は二値モデルの時と同じように\(P(y_{pi}\geq k)=0.5\)になる\(\theta_p\)の値を表しています。 そのため，上の確率密度分布を見ると，色が変わるポイントはすべての分布で同じになっていることがわかります。 したがって，\(\theta_p\)の値が大きい（正規分布が右に動く）ほど，また\(\beta_{ik}\)の値が小さい（閾値の点線が左に動く）ほど，上位カテゴリを選択する確率が高くなっていくわけです。
続いて， 図 8.15 には，いくつかの項目パラメータのセットにおける項目特性曲線(ICC)を示しました。 図 8.14 の下半分では\(P(y_{pi}\geq k)\)のプロットを示しましたが，通常「項目特性曲線」というと\(P(y_{pi} = k)\)のプロットです。 二値項目の場合はどちらでも同じですが，多値項目の場合は異なるため改めてICCを示しています。 なお，多値型モデルの場合は 図 8.15 の一本一本の線のことを項目反応カテゴリ特性曲線(item response category characteristic curve [IRCCC])と呼ぶこともあります。 が，そこまで細かい名前を覚える必要はないでしょう。
左上のプロットは， 図 8.15 と同じ項目パラメータ（\((\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(1,[-1.5, 0, 2])\)）でのIRCCCです。 基本的には，このように\(\theta_p\)の値に応じて左から順に各カテゴリの選択確率のピークが訪れます。 点線は中間カテゴリの選択確率が最も高くなる箇所ですが，これは具体的にカテゴリ困難度\(\beta_{ik}\)と\(\beta_{i,k+1}\)のちょうど中間になります。
右上のプロットは，左上からすべてのカテゴリ困難度の値を1ずつ大きくしました（\((\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(1,[-0.5,1,3])\)）。 二値モデルの時と同じように，すべてのカテゴリ困難度が同じ値だけ変化するときにはIRCCCは平行移動します。
左下のプロットは，左上から識別力を半分にしました（\((\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(0.5,[-1.5, 0, 2])\)）。 また右下のプロットは，左上から識別力を2倍にしています（\((\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4]})=(2,[-1.5, 0, 2])\)）。 識別力は各IRCCCの山のきつさを調整しています。 識別力が低い場合，\(\theta_p\)の大小に対して各カテゴリの選択確率があまり変化しなくなるため，反対に「あるカテゴリを選択した」という情報が\(\theta_p\)の推定に及ぼす影響も小さくなってしまいます。
また，左下のプロットのように，どの\(\theta_p\)の値においてもあるカテゴリの選択確率が最大とならないこともありえます。 識別力が低い場合や，カテゴリ困難度が近すぎる場合，あるいはそもそもカテゴリ数が多い場合にはそのような現象が起こりえます。 ただ何も問題ではないので，そのようなカテゴリが出現した場合には「選ばれにくいんだなぁ」くらいに思っておいてください。
8.6.2 一般化部分採点モデル
有名な多値型モデルのもう一つが，一般化部分採点モデル（generalized partial credit model [GPCM]: Muraki, 1992）です。 GPCMでは，隣のカテゴリとの関係を別の視点から考えます。 ということで，まずは「隣のカテゴリとの関係」を二値（2パラメータロジスティック）モデルで考えてみましょう。 二値モデルでは，「カテゴリ0を選ぶ確率」と「カテゴリ1を選ぶ確率」の和は当然1になります。 そして，「カテゴリ1を選ぶ（e.g., 『当てはまる』を選ぶ，正解する）」確率は \[ \begin{aligned} P(y_{pi}=1) &= \frac{P(y_{pi}=1)}{P(y_{pi}=1)+P(y_{pi}=0)} \\ &= \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_{i})\right]} \\ &= \frac{\exp\left[\alpha_{i}(\theta_p - \beta_{i})\right]}{1+\exp\left[\alpha_{i}(\theta_p - \beta_{i})\right]} \end{aligned} \tag{8.20}\] となります。 つまりこの式は，隣のカテゴリとの二択において，当該カテゴリを選ぶ確率という見方ができるわけです。 これを多値\((k=1,2,\cdots,K)\)に置き換えると，「カテゴリ\(k-1\)と\(k\)の二択だったら\(k\)の方を選ぶ確率」に拡張することができそうです。 ということで，GPCMではこの確率\(P(y^*_{pi}=k)_{k-1,k}\)を \[ \begin{aligned} P(y^*_{pi}=k)_{k-1,k} &= \frac{P(y_{pi}=k)}{P(y_{pi}=k)+P(y_{pi}=k-1)} \\ &= \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_{ik})\right]} \\ &= \frac{\exp\left[\alpha_{i}(\theta_p - \beta_{ik})\right]}{1+\exp\left[\alpha_{i}(\theta_p - \beta_{ik})\right]} \\ &= \frac{\exp\left[\pi_{pik}\right]}{1+\exp\left[\pi_{pik}\right]} \end{aligned} \tag{8.21}\] と表します。なお，この後の式展開を考えて\(\pi_{pik}=\alpha_{i}(\theta_p - \beta_{ik})\)と置き換えています。
ここから，全カテゴリの中でカテゴリ\(k\)を選ぶ確率を導出しましょう。 (8.21)式を\(P(y_{pi}=k)=\)の形に変形すると， \[ P(y_{pi}=k) = P(y_{pi}=k-1)\frac{P(y^*_{pi}=k)_{k-1,k}}{1-P(y^*_{pi}=k)_{k-1,k}} \tag{8.22}\] となります。最後の右辺の\(\frac{P}{1-P}\)の部分は2つのカテゴリ間のオッズを表しているため， \[ \frac{P(y^*_{pi}=k)_{k-1,k}}{1-P(y^*_{pi}=k)_{k-1,k}} = \exp(\pi_{pik}) \tag{8.23}\] を用いて \[ P(y_{pi}=k) = P(y_{pi}=k-1)\exp(\pi_{pik}) \tag{8.24}\] という漸化式が得られ，特定のカテゴリ\(k\)を選ぶ確率が表現できるようになります。 例えば4カテゴリの場合には， \[ \begin{aligned} P(y_{pi}=1) &= P(y_{pi}=1) \\ P(y_{pi}=2) &= P(y_{pi}=1)\exp(\pi_{pi2}) \\ P(y_{pi}=3) &= P(y_{pi}=1)\exp(\pi_{pi2})\exp(\pi_{pi3}) \\ P(y_{pi}=4) &= P(y_{pi}=1)\exp(\pi_{pi2})\exp(\pi_{pi3})\exp(\pi_{pi4}) \end{aligned} \tag{8.25}\] という形で\(P(y_{pi}=1)\)を起点に「一個上のカテゴリに上がるときのオッズ」の積によってすべてのカテゴリの反応確率が表されるわけです。 ということで，この式を一般化して \[ \begin{aligned} P(y_{pi}=k) &= P(y_{pi}=k-1)\exp(\pi_{pik}) \\ &= P(y_{pi}=k-2)\exp(\pi_{pi,k-1})\exp(\pi_{pik}) \\ &\qquad\vdots \\ &= P(y_{pi}=1)\prod_{c=2}^{k} \exp(\pi_{pic}) \\ &= P(y_{pi}=1)\exp\left(\sum_{c=2}^{k}\pi_{pic}\right) \end{aligned} \tag{8.26}\] と書き下します。 ここではまだ\(P(y_{pi}=1)\)が残っていますが，GPCMの考え方ではこれ以上計算できません13。 改めて(8.26)式を見ると，すべてのカテゴリの反応確率は「\(P(y_{pi}=1)\)の何倍」という形になっていることがわかります。 そこで一旦\(\pi_{pi1}=0\)として得られる\(P(y_{pi}=1)=\exp(\pi_{pi1})=1\)を(8.26)式に代入して， \[ P(y_{pi}=k) = \exp\left(\sum_{c=1}^{k}\pi_{pic}\right) \tag{8.27}\] とまとめてしまいます。 これは，各カテゴリの反応確率を「一番下のカテゴリが選択される確率の何倍か」という形で表すことを意味しています。 これによって，ひとまず\(P(y_{pi}=k)\)が定式化されましたが，\(P(y_{pi}=1)=1\)と設定している以上\(P(y_{pi}=k)\)の全カテゴリの総和が1になるかわかりません。
ということで最後に，\(P(y_{pi}=1)\)の値が何であっても全カテゴリでの総和が1になるようにするため，総和で割ることを考えます。 実際に \[ \begin{aligned} \frac{P(y_{pi}=k)}{\sum_{l=1}^{K}P(y_{pi}=l)} &=\frac{P(y_{pi}=k)}{P(y_{pi}=1)+P(y_{pi}=2)+\cdots +P(y_{pi}=K)} \\ &= \frac{\exp\left(\sum_{c=1}^{k}\pi_{pic}\right)}{\exp\left(\sum_{l=1}^{1}\pi_{pil}\right) + \exp\left(\sum_{l=1}^{2}\pi_{pil}\right) + \cdots + \exp\left(\sum_{l=1}^{K}\pi_{pil}\right)} \\ &= \frac{\exp\left(\sum_{c=1}^{k}\pi_{pic}\right)}{\sum_{l=1}^{K}\exp\left(\sum_{c=1}^{l}\pi_{pil}\right)}\\ \end{aligned} \tag{8.28}\] とすると，分子も分母もすべての項に共通の\(\exp(\pi_{pi1})\)がかかっているため，この項の値が何であっても全体の割合に占める各カテゴリの反応確率の割合は変わらずに済みます。 ということで，ようやくGPCMの式が導出できました。 図 8.16 に，ここまでの考え方をまとめてみました。 特定のカテゴリの反応確率\(P(y_{pi}=k)\)が直接計算できないため，いったん一番下のカテゴリの反応確率\(P(y_{pi}=1)=1\)とおいて，他のカテゴリは「前のカテゴリの何倍」という形で表します。 最後に全体の総和が1になるように調整をしたものがGPCMにおける反応確率というわけです。
図 8.17 は，GPCMでの項目特性曲線をプロットしたものです。 GPCMでは，カテゴリ困難度\(\beta_{ik}\)は「カテゴリ\(k-1\)と\(k\)の二択」における選択確率がちょうど50%になるポイントを表しています。 したがって，隣接するカテゴリの反応確率はちょうど\(\beta_{ik}\)のところで大小関係が入れ替わります。 \(\beta_{ik}\)はそのように解釈できるわけです。
GRMでは，カテゴリ困難度は必ず単調増加になっている必要がありますが，GPCMでは単調増加で無くても問題ありません。 図 8.18 にカテゴリ困難度が単調増加にならないケースをプロットしてみました。 この場合，「カテゴリ2がカテゴリ1を上回る」よりも先に「カテゴリ3がカテゴリ2を上回る」ため，結果的にカテゴリ2が最大になることがありません。
8.6.3 Rで多値型モデル
それでは，mirtで2つの多値型モデルを実行してみましょう。 …といっても引数itemtypeにgraded(GRM)かgpcmを指定したら良いだけなので難しいことはありません。
$items
a b1 b2 b3 b4 b5
Q1_1 0.572 -6.382 -3.936 -2.323 -1.040 1.269
Q1_2 1.181 -4.051 -2.765 -2.087 -0.840 0.821
Q1_3 1.238 -3.245 -2.243 -1.643 -0.564 1.019
Q1_4 1.074 -3.293 -2.174 -1.648 -0.706 0.415
Q1_5 1.116 -3.951 -2.472 -1.664 -0.460 1.216
Q1_6 1.062 -3.988 -2.698 -1.726 -0.440 1.429
Q1_7 1.177 -3.454 -2.084 -1.314 -0.178 1.485
Q1_8 1.097 -3.671 -2.186 -1.386 -0.023 1.744
Q1_9 1.365 -3.306 -2.016 -1.039 -0.312 0.921
Q1_10 1.200 -2.226 -1.048 -0.056 0.461 1.586
$means
f_1
0
$cov
f_1
f_1 1datの各項目は6件法なので，カテゴリ困難度bが5つ出力されました。
二値型モデルの時と同じように，項目特性曲線を出力することももちろん可能です。
また，itemplot(., type = "threshold")関数を使うと， 図 8.14 の下半分のような「カテゴリ\(k\)以上の反応をする確率」のプロットを出すこともできます。
同じようにGPCMでもやってみましょう。 やり方も結果の見方も同じなので詳しい説明は省略します。
$items
a b1 b2 b3 b4 b5
Q1_1 0.222 -5.012 -2.300 -1.027 -3.288 -0.347
Q1_2 0.681 -2.878 -1.254 -2.419 -1.103 0.536
Q1_3 0.672 -2.103 -1.046 -2.050 -0.935 0.791
Q1_4 0.484 -2.141 -0.254 -2.357 -0.891 -0.851
Q1_5 0.582 -3.195 -1.219 -1.986 -0.831 0.978
Q1_6 0.421 -3.016 -1.949 -2.352 -1.175 1.566
Q1_7 0.491 -3.108 -1.031 -1.960 -0.787 1.528
Q1_8 0.467 -3.297 -0.963 -2.341 -0.371 1.841
Q1_9 0.553 -3.276 -1.887 -0.518 -0.992 0.491
Q1_10 0.438 -1.844 -0.932 1.298 -1.047 0.940
$means
f_1
0
$cov
f_1
f_1 1
8.7 多次元モデル
これまでのIRTモデルでは一次元性の仮定に基づく一次元モデルのみを扱ってきました。 ですが多くの心理尺度やテストでは，複数の次元の特性を同時に測定したいというニーズがあります。 先程紹介したようにIRTはカテゴリカル因子分析と数学的には同値なので，IRTでも多次元モデルを考えていきましょう。
多次元因子分析モデルでは，項目反応を \[ y_{pi} = \tau_i + \sum_{t=1}^{T}{f_{pt} b_{ti}} - \varepsilon_{pi} \tag{8.29}\] と表しました。これと同じように考えるならば，多次元IRTモデル（2パラメータロジスティックモデル）での項目反応関数は \[ P(y_{pi}=1) = \frac{1}{1+\exp\left[-\sum_{t=1}^{T}\alpha_{ti}\theta_{pt} - \tau_i\right]} \tag{8.30}\] と表せそうです。
(8.30)式のモデルにおける項目パラメータは，識別力が次元の数（\(T\)個）と，切片が一つです。 また，項目の全体的な識別力を評価する際には，次元ごとの識別力を一つにまとめた多次元識別力（multidimensional discrimination: 豊田，2013） \[ \mathrm{MDISC}_i=\sqrt{\sum_{t=1}^{T}\alpha_{ti}^2} \tag{8.31}\] を用いることもあります。これを用いると，切片パラメータ\(\tau_i\)を困難度に変換することもできるようになります： \[ \beta_i=\frac{-\tau_i}{\mathrm{MDISC}_i}. \tag{8.32}\]
二次元までならば，ICCと同じように\(\symbf{\theta}_p\)と\(P(y_{pi}=1)\)の関係をプロットすることも可能です。 図 8.22 の左は，とある二次元項目における項目反応曲面 (item response surface [IRS])です。 X軸，Y軸がそれぞれの\(\theta_t\)の値を表しており，Z軸の値が反応確率\(P(y_{pi}=1)\)を表しています。 この項目では，\(([\alpha_{i1},\alpha_{i2}],\tau_i)=([1.5,0.5],0.7)\)となっていることから，\(\theta_{p2}\)の値はさほど反応確率に影響を及ぼしていないことが図からもわかります。 なお 図 8.22 の右は，2つの\(\theta_{pt}\)の値に対して\(P(y_{pi}=1)\)が等しくなるラインを表しています。
8.7.1 補償モデルと非補償モデル
ここまで紹介した多次元モデルでは 図 8.22 からも分かるように，反応確率に対して\(\theta_{pt}\)が及ぼす影響は加法的といえます。 例えば\(([\alpha_{i1},\alpha_{i2}],\tau_i)=([1.5,1],0)\)の項目の場合，\(\symbf{\theta}_p=\begin{bmatrix}-1 & 1\end{bmatrix}\)の人と\(\symbf{\theta}_p=\begin{bmatrix}0 & -0.5\end{bmatrix}\)の人では反応確率は同じになります。 また，ある次元の特性値がとても低い場合でも，別の次元の特性値がとても高い場合には反応確率が高くなる可能性があります。 このような次元間の関係を反映した先程の多次元IRTモデルは補償型 (compensatory) モデルと呼ばれます。
一方で，複数の次元の要素がすべて揃って初めて反応確率が高くなるようなケースも考えられます。 例えば数学のテストが日本語ではなく英語で実施される場合，これは「英語」と「数学」の2つの能力が試される二次元構造ではありますが，補償型モデルのように「英語が苦手でも数学力が高ければ解答できる」とは考えにくいです。 英語で書かれた問題文を理解するだけの英語力と，計算を実行できるだけの数学力の両方が揃って初めて正答できるでしょう。 このような状況では，一方の特性値の低さをもう一方の特性値の高さでカバーすることはできません。 こうした次元間の関係を反映したモデルを非補償 (noncompensatory) モデルあるいは部分補償 (partially compensatory) モデルと呼びます14。
最もシンプルな非補償モデルの項目反応関数は \[ P(y_{pi}=1) = \prod_{t=1}^{T}\frac{1}{1+\exp\left[-\alpha_{it}(\theta_{pt} - \beta_{it})\right]} \tag{8.33}\] と表されます。 各次元に関して独立に算出された「その次元の閾値を超える確率」を計算し，最後にすべての積をとっています。 そのため，一つでも特性値の低い次元があれば反応確率は一気に小さくなってしまう，というわけです。 非補償モデルの項目反応曲面は 図 8.23 の左のようになります。 図 8.22 と比べると，両方の次元の\(\theta_{p}\)が高くないと\(P(y_{pi}=1)\)が高くならないことが見て取れます。
補償モデルと非補償モデルを比べると，補償モデルのほうがよく用いられているような印象です。 その理由の一つとして，非補償モデルのほうがパラメータの推定が難しいという点が挙げられると思います。 上の式を見るとわかるように，非補償モデルでは困難度パラメータ\(\beta_{it}\)が次元の数だけあるためそもそものパラメータ数が補償型よりも多くなっています。 また，非補償モデルは項目反応関数の全体の形が積になっています。 コンピュータは足し算よりも掛け算のほうが苦手なので，そういった意味でも計算の難易度が高いのです。 また，心理尺度的には補償型のほうが従来の因子分析モデルと同じ形になっているため，扱いやすいというのもあると思います。 とはいえ，多次元の考え方として非補償型もあるんだということを頭の片隅に置いておくと，いつかどこかで何かの参考になるかもしれません15。
8.7.2 Rでやってみる
ここからは補償型の多次元IRTモデルをmirtで実行する方法を紹介します16。 補償型モデルはカテゴリカル因子分析と同値なので，因子分析と同じように「探索的モデル」と「検証的モデル」の2種類があります。 ただ探索的モデルでは因子の回転に応じて\(\theta_{pt}\)の意味が変化してしまうため実用上はなかなか使いにくいような気もします。
探索的補償型モデル
探索的モデルを実行する場合には，引数modelに「次元数」を与えるだけです。 itemtypeはこれまでと同じように自由に選んでください。二値だったら2PLとか，多値だったらgradedとかgpcmとか。
$items
a1 a2 d1 d2 d3 d4 d5
Q1_1 -0.251 -0.877 3.879 2.433 1.457 0.665 -0.785
Q1_2 -0.817 -1.700 5.685 3.977 3.034 1.234 -1.208
Q1_3 -0.990 -2.403 5.789 4.099 3.032 1.048 -1.891
Q1_4 -0.775 -0.903 3.646 2.414 1.831 0.781 -0.465
Q1_5 -0.757 -1.523 5.105 3.271 2.221 0.612 -1.628
Q1_6 -1.449 0.151 4.667 3.195 2.067 0.537 -1.707
Q1_7 -1.645 0.183 4.602 2.822 1.798 0.252 -2.011
Q1_8 -1.329 0.047 4.263 2.559 1.625 0.028 -2.047
Q1_9 -1.918 0.131 5.227 3.255 1.695 0.508 -1.501
Q1_10 -1.436 0.000 2.864 1.361 0.075 -0.595 -2.040
$means
F1 F2
0 0
$cov
F1 F2
F1 1 0
F2 0 1結果を見ると，識別力パラメータが次元の数(a1, a2)と，切片パラメータが（カテゴリ数-1）個(d1-d5)推定されています。 その下の$meanと$covには，それぞれ\(\theta_{pt}\)の平均と共分散行列が出ています。 これを見る限り，まだ因子間相関が0なので回転前の因子負荷（識別力）が出ているようです。 さらに詳しく見ると，最後の項目のa2が0.000となっています。 実はmirtでの探索的モデルの推定の際には，このように一部の因子負荷を0に固定した状態で推定を行っています17。 「探索的」と紹介していましたが，実際には「モデルが識別できる最小限の制約を置いて推定する」ということを行っているのです。 IRTに限らずCFAもそうなのですが，複数の因子がある場合に識別性を持たせるには，最低でも\(\frac{T(T-1)}{2}\)個の因子負荷を固定する必要があることが知られています。 そのため，2因子の場合は第2因子の因子負荷を一つ0にしているわけです。 同様に，3因子の場合にはさらに第3因子の因子負荷を2つ0にする必要があり，4因子の場合にはさらにさらに第4因子の因子負荷を3つ0にする必要があり…という要領です。
このままでは，いくら探索的とは言え各次元の解釈が出来ません。 mirt()で得られるパラメータの推定値はいわば「初期解」なので，ここからさらに回転させてあげれば良さそうです。 回転させる場合は，psych::fa()の時と同じように引数rotateを指定してあげます。
Rotation: oblimin
$items
a1 a2 d1 d2 d3 d4 d5
Q1_1 -0.146 0.953 3.879 2.433 1.457 0.665 -0.785
Q1_2 0.037 1.872 5.685 3.977 3.034 1.234 -1.208
Q1_3 -0.107 2.634 5.789 4.099 3.032 1.048 -1.891
Q1_4 0.351 1.021 3.646 2.414 1.831 0.781 -0.465
Q1_5 0.057 1.680 5.105 3.271 2.221 0.612 -1.628
Q1_6 1.474 -0.050 4.667 3.195 2.067 0.537 -1.707
Q1_7 1.678 -0.070 4.602 2.822 1.798 0.252 -2.011
Q1_8 1.311 0.051 4.263 2.559 1.625 0.028 -2.047
Q1_9 1.920 0.007 5.227 3.255 1.695 0.508 -1.501
Q1_10 1.394 0.110 2.864 1.361 0.075 -0.595 -2.040
$means
F1 F2
0 0
$cov
F1 F2
F1 1.00 0.35
F2 0.35 1.00回転後の識別力は，いい感じに最初の5項目と最後の5項目に分かれてくれました。
検証的補償型モデル
検証的モデルを行う際には，lavaanに似た要領で，引数modelに各次元と項目の関係を記述してあげる必要があります。
$items
a1 a2 d1 d2 d3 d4 d5
Q1_1 0.895 0.000 3.865 2.425 1.453 0.665 -0.779
Q1_2 1.882 0.000 5.693 3.979 3.032 1.228 -1.212
Q1_3 2.590 0.000 5.789 4.095 3.026 1.042 -1.888
Q1_4 1.109 0.000 3.575 2.353 1.780 0.750 -0.465
Q1_5 1.701 0.000 5.118 3.275 2.219 0.605 -1.632
Q1_6 0.000 1.463 4.670 3.197 2.070 0.537 -1.712
Q1_7 0.000 1.606 4.538 2.779 1.770 0.249 -1.980
Q1_8 0.000 1.330 4.263 2.557 1.621 0.025 -2.047
Q1_9 0.000 1.979 5.311 3.310 1.719 0.509 -1.532
Q1_10 0.000 1.421 2.853 1.352 0.071 -0.596 -2.032
$means
f_1 f_2
0 0
$cov
f_1 f_2
f_1 1 0
f_2 0 1きちんと項目1-5の次元2の識別力（a2）と項目6-10の次元1の識別力（a1）が0に固定されています。 先程の探索的モデルと見比べても，因子の順序こそ違えど概ね同じような推定値が得られていますね。
ちなみに，多次元モデルにおけるICCも，一次元モデルの時と同じようにplot()やitemplot()によって出力可能です。
また，引数rotをうまく指定すると，プロットを回すことができます。 うまく回せば，多少は見やすくなるかもしれません。 例えば確認的モデルでは，各項目の識別力を一つだけ非ゼロにしていました。 そのため，識別力がゼロである次元の方向は潰すように回転させると良いかもしれません（見やすくなるとは言っていない）。
8.8 次元性の検証
セクション 8.5.2 で少し触れたように，IRTでは基本的にそのテストが測定しようとする次元数は理論的に決まっていることが多く，そのため「何因子が最適か」を探索的に決める方法よりは「\(n\)因子で問題ないか」を検証する方法が採用されることが多いと思います。 ここでは，そのようなIRTの文脈で用いられることの多い（と思われる）いくつかの次元性検証法を紹介していきます。
8.8.1 パラメトリックな方法
パラメトリックな方法では，（多次元・一次元の）IRTや因子分析モデルに基づいてパラメータの推定を行い，その結果を用いて次元性を評価していきます。 例えば確認的因子分析によって推定を行い，SEM的なモデル適合度をチェックする方法はパラメトリックな検証法の一つと言えるでしょう。 その他にも，IRTの文脈で用いられている検証法がいくつか存在します。 代表的と思われる方法の一つがNOHARM（Gessaroli & De Champlain, 1996）です18。 この方法では，項目の各ペアでの残差共分散を，推定されたパラメータおよび観測されたデータ（各項目の正答率）から計算していきます。 そして，全ペアに対して計算された残差共分散（を\(z\)変換したもの）の二乗和（を何倍かしたもの）が\(\chi^2\)分布に従うことを利用して，「残差共分散が全てゼロである」という帰無仮説に対する検定を行います。 そして帰無仮説が棄却されなかった場合，「その次元数で全ての項目ペア間の相関関係を表現できていないとは言えない」ということで，その次元数でOK，という判断がくだされます。 考え方としてはBartlett検定に近いかもしれません。
8.8.2 ノンパラメトリックな方法
IRTの文脈で提案されたノンパラメトリックな次元性評価法の中でも代表的なものが，ここで紹介するDIMTEST（Stout et al., 1996）やDETECT（Zhang & Stout, 1999）です。 これらの方法では，「その次元（共通因子）数があれば局所独立性が満たされるか」という視点から評価を行います。 セクション 8.5.1 では，局所独立性を「どの項目ペア間でも，どの特性値の値でも連関がない」ことと表現していましたが，これを数式で表すと \[ \mathrm{Cov}(\symbf{y}_i, \symbf{y}_j\vert\theta) = 0 ~ \mathrm{for~all}~i,j,\theta \tag{8.34}\] となります19。したがって，\(\theta\)で条件づけた各項目ペアの（条件付き）共分散を計算していけば，局所独立性を満たすために必要な次元数を考えることができるわけです。
具体的な方法の説明に移る前に，次元性と条件付き共分散のイメージを確認しておきます。 図 8.26 には，2因子モデルにおいて，なるべく単純構造になるように回転した後の各項目の因子負荷をベクトルとして表したものです。また\(\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}\)と表記されているベクトルは，各クラスタ内の項目の因子負荷の平均を取ったようなものを表しています。同様に\(\symbf{\Theta}_{TT}\)と表記されているベクトルは，クラスタをまたいだ全項目のベクトルの平均を取ったもので，端的に言えばこれが「テスト全体で測定しようとしている特性のベクトル」というわけです。 よくある例としては，\(\symbf{\Theta}_{TT}\)が「数学の学力」，\(\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}\)がそれぞれ「代数」「幾何」といった感じです。 図 8.26 では，\(\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}\)は完全に無相関というわけではないですが，さほど強くもないので，一次元性が満たされているとは言いづらい感じがします。
このテストが一次元性を満たしているかを確認するためには，「\(\symbf{\Theta}_{TT}\)で条件付けられた」共分散を考えることになります。 そして条件付き共分散は， 図 8.26 においてはベクトル\(\symbf{\Theta}_{TT}\)に向かって伸ばした足の長さの積のような形で表されます。 これを理解するために， 図 8.26 に少し描き足してみました。 図 8.27 に追加されたオレンジの線は，\(\symbf{\Theta}_{TT}=\frac{\symbf{\Theta}_{C1}+\symbf{\Theta}_{C2}}{2}\)となる（次元間で重みが同じ）場合に\(\symbf{\Theta}_{TT}\)の値が同じになる\(\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}\)の組みあわせを表しています。
この図をもとに，項目ペアに対して計算できる条件付き共分散を考えてみましょう。 まず，クラスター\(C1\)の中の2つの項目の条件付き共分散は正の値になることが期待されます。 \(\symbf{\Theta}_{TT}\)の値が同じ2人（AさんとBさん）について考えてみると，Aさんはクラスター\(C1\)の項目に正解するのに必要な能力（\(\symbf{\Theta}_2\)）の値が高いためクラスター\(C1\)の項目には正解する確率が高いと予想される一方で，Bさんは\(\symbf{\Theta}_2\)の値が低いので，正解する確率は低そうです。 これをまとめると，「\(\symbf{\Theta}_{TT}\)で条件づけたとき，クラスター\(C1\)内のある項目に正解する人ほど，同じクラスター内の別の項目にも正解する確率が高い」わけなので，条件付き共分散は正の値になります。
同様に，クラスター\(C1\)の中のある項目と，別のクラスター\(C2\)の中のある項目の条件付き共分散を考えてみると，Aさんはクラスター\(C1\)の項目には正解するがクラスター\(C2\)の項目には正解できないと予想されます。 そしてBさんは反対に，クラスター\(C2\)の項目にのみ正解できる可能性が高いでしょう。 したがって，「\(\symbf{\Theta}_{TT}\)で条件づけたとき，クラスター\(C1\)内のある項目に正解する人ほど，別のクラスター内のある項目に正解する確率は低い」ことになり，条件付き共分散は負の値になります。
以上をまとめると， 図 8.27 における条件付き共分散の符号は，2つの項目の因子負荷ベクトルが「ベクトル\(\symbf{\Theta}_{TT}\)から見て」同じ方向にあるかどうかに対応するわけです。 そして，仮に全ての項目が近い方向のベクトルで表される場合には，当然\(\symbf{\Theta}_{TT}\)ベクトルも同じような方向になるため，条件付き共分散も0に近くなっているはずです。 という感じで，(8.34)式が局所独立性，そして一次元性を表しています。
DIMTEST
DIMTEST（Stout et al., 1996）では，まずテストを「（次元性の）評価用 (assessment subtest: AT)」と「データの分割用 (partitioning subtest: PT)」の2つに分割します。 この分割では，理論的に想定される分かれ方（ 図 8.26 の\(\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}\)に相当）でわけたり，探索的因子分析をもとに分けたりします。 つまり，一旦「このテストが分かれる（2次元だ）としたらこうなるだろう」という形で分けてみるのです。 もしこのテストが実際には一次元だとしたら，そのように分けたとしても結局全ての項目のベクトルの向きは近いので，片方のテストの得点によって条件付けたとしても，もう一方のテストの中の共分散は0に近いはずです。 一方で，テストが（\(\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}\)のように）分かれる場合には，例えば\(\symbf{\Theta}_{C1}\)で条件づけた際の\(\symbf{\Theta}_{C2}\)内の条件付き共分散は全て正になってしまいます。
以上の考え方に基づいて，DIMTESTでは以下の統計量を計算します。 \[ T = \sum_k \sum_{(i,j \in AT)} \widehat{\mathrm{Cov}}(\symbf{y}_{i},\symbf{y}_{j}|S_{PT}=k). \tag{8.35}\] ここで\(S_{PT}\)はPTにおける正答数を表しています。 したがって，(8.35)式は「PTにおける正答数で（＝\(\symbf{\Theta}_{PT}\)で）条件づけたときのAT内の共分散」を求めているわけです。 統計量\(T\)が0に近いほど局所独立が満たされている（そのテストは一次元である）と言えるので，あとはこれを標準化した値を用いて標準正規分布に基づく\(z\)検定 （Stout, 1987） を行います。
DIMTESTの成功のカギは，ATとPTの分割にあるわけですが，どのように分けたら良いかについては Stout et al. （1996） などを参照してください。
DETECT
DETECTも条件付き共分散をベースにした方法です。 以下では，二値データに対して考案されたもの（Zhang & Stout, 1999）を紹介します。 多値データの場合にはこれを発展させた方法（Zhang, 2007）が提案されています。
DETECTでも，まずはテストをいくつかのクラスター（\(C1, C2, \cdots, Ck\)）に分割します。 そして以下のようにテスト全体で推定された特性値\(\symbf{\Theta}_{TT}\)によって条件づけた共分散の統計量を計算します。 \[ D = \frac{2}{I(I-1)}\sum_{1\leq i\leq j\leq I}\delta_{i,j}E[\mathrm{Cov}(\symbf{y}_i,\symbf{y}_j|\symbf{\Theta}_{TT})]. \tag{8.36}\] ここで重要なのは指示変数\(\delta_{i,j}\)です。 これは，
- 2つの項目\(i,j\)が同じクラスターにある場合には1
- 2つの項目\(i,j\)が異なるクラスターにある場合には-1
になります。
DETECTが計算しようとしているものを 図 8.26 をもとに考えてみます。 図 8.26 は，あるテストの真のクラスターを表しているとします。 つまりあるテストはきれいに2つのクラスターに分かれている，ということです。 このときに，この真のクラスターと一致するようにテストを分けた上でDETECTの統計量\(D\)を計算してみましょう。
- 2つの項目\(i,j\)が同じクラスターにある場合には，条件付き共分散（の期待値）\(E[\mathrm{Cov}(\symbf{y}_i,\symbf{y}_j|\symbf{\Theta}_{TT})]\)も\(\delta_{i,j}\)も正の値
- 2つの項目\(i,j\)が異なるクラスターにある場合には，条件付き共分散（の期待値）\(E[\mathrm{Cov}(\symbf{y}_i,\symbf{y}_j|\symbf{\Theta}_{TT})]\)も\(\delta_{i,j}\)も負の値なので，積は正の値になる
ということで，全ての項目ペア間で計算される値（\(\sum\)の中身）が全て正の値になります。 このように，統計量\(D\)は（クラスター数も含めて）最も適切なクラスター基づいて計算したときに最大値を取る値です。 そして，そもそも全ての項目のベクトルの方向が近ければ（次元性が満たされていれば），条件付き共分散（の期待値）は全て0に近いので，その平均である統計量\(D\)も0に近い値になっているはずです。
ということで，この統計量を用いて一次元性を評価する場合には，「統計量\(D\)が最大でいくつになるか」を見たら良いと言えます。 例えば「理論的に分かれるとしたらこう」で分けてみたり，探索的因子分析の結果をもとに分けてみたり，あるいはとにかく手当たり次第にランダムに様々な分け方をひたすら試してみたりして，統計量\(D\)が高々どれくらいの値に収まるかを確認しましょう。 一次元性の検証という意味では，もちろん\(D\)の値は小さいほどよいのですが，二値データの場合例えば Roussos & Ozbek （2006） では最大値が0.2以下なら一次元とみなしてよいのではないか，と言われていたりします。
そんなDETECTについては，Rで計算するための関数としてsirtパッケージにconf.detect()およびexpl.detect()というものが用意されています20。 これらは二値・多値のどちらにも対応していますが，ここでは二値データへの適用例を見てみます。
まずはconf.detect()です。confは”confirmatory”のことで，その名の通り分析者が指定したクラスター構造のもとでDETECTを計算します。 そのため，事前に想定されるクラスターがある（数学のテストの中に複数の単元がある，心理尺度の中に複数の下位概念がある）ような場合に使える方法です。
-----------------------------------------------------------
Confirmatory DETECT Analysis
Conditioning on 1 Scores
Bandwidth Scale: 1.1
DETECT Calculation Score 1
-----------------------------------------------------------
NScores Mean SD Min Max
DETECT Unweighted 1 1.064 NA 1.064 1.064
DETECT Weighted 1 1.064 NA 1.064 1.064
ASSI Unweighted 1 0.511 NA 0.511 0.511
ASSI Weighted 1 0.511 NA 0.511 0.511
RATIO Unweighted 1 0.664 NA 0.664 0.664
RATIO Weighted 1 0.664 NA 0.664 0.664結果の中のDETECT Unweightedにかかれている値が(8.36)式で計算される\(D\)の値です。この値が0.2より小さければ，そのテストは一次元であるとみなせそうなのですが…今回使用した10項目は元々2次元だったため，「一次元ではない」と正しく判定されています（一般に1を超えると「かなり強く多次元である」と言えるようです）。
続いて，本来一次元であることが想定されるケースとして，項目間の相関が高い（＝内的一貫性の高い）Q1_6からQ1_10の5項目のみを使ってconf.detect()してみましょう。
-----------------------------------------------------------
Confirmatory DETECT Analysis
Conditioning on 1 Score
Bandwidth Scale: 1.1
-----------------------------------------------------------
unweighted weighted
DETECT 0.425 0.425
ASSI 0.200 0.200
RATIO 0.140 0.140
MADCOV100 3.041 3.041
MCOV100 -3.041 -3.041DETECT Unweightedの値は先程よりは小さくなっていますが，それでも0.2よりは大きくなっています…理由はよく分からないのですが，もしかしたら項目数が少なすぎる場合には小さくなりにくいのかもしれません（要出典）。
ここまでは，引数itemclusterによって項目のクラスタリングを与えていました。 一次元性を検証するためには，本当はこれをすべてのパターン（\(I\)項目を\(K\)個のクラスタに分けるとしたらおよそ\(K^I\)通り）で計算して，その最大値を出す必要があります。 ですがこれはどう考えても非効率な（というかすぐ無理になる）話です。 そこで，各クラスタ数における「\(D\)が最大になる分割」を探索的に求めることを考えましょう。
expl.detect()では，階層的クラスター分析を用いて，指定したクラスタ数の分割の中で\(D\)が最大になるものを探索的に決定してくれます。 考え方はシンプルで，条件付き共分散（の期待値）\(E[\mathrm{Cov}(\symbf{y}_i,\symbf{y}_j\vert\symbf{\Theta}_{TT})]\)が大きいペアは近くにある＝同じクラスターに属すると考えて，距離の近い項目から順に，最終的に\(K\)個のクラスターにまとまるまでくっつけていき，そのクラスター構造の元で計算した\(D\)を出力してくれます。 ということで早速試してみましょう。
Pairwise Estimation of Conditional Covariances
...........................................................
Nonparametric ICC estimation
...........................................................
Nonparametric Estimation of conditional covariances
5% 10% 15% 20% 25% 30% 35% 40% 45% 50%
55% 60% 65% 70% 75% 80% 85% 90% 95%
Pairwise Estimation of Conditional Covariances
...........................................................
Nonparametric ICC estimation
...........................................................
Nonparametric Estimation of conditional covariances
5% 10% 15% 20% 25% 30% 35% 40% 45% 50%
55% 60% 65% 70% 75% 80% 85% 90% 95%
DETECT (unweighted)
Optimal Cluster Size is 4 (Maximum of DETECT Index)
N.Cluster N.items N.est N.val size.cluster DETECT.est ASSI.est RATIO.est
1 2 10 1216 1216 5-5 1.133 0.511 0.674
2 3 10 1216 1216 5-3-2 1.590 0.778 0.947
3 4 10 1216 1216 5-2-1-2 1.659 0.867 0.988
4 5 10 1216 1216 4-1-2-1-2 1.644 0.867 0.979
MADCOV100.est MCOV100.est DETECT.val ASSI.val RATIO.val MADCOV100.val
1 1.68 -1.284 1.136 0.556 0.680 1.669
2 1.68 -1.284 1.543 0.822 0.924 1.669
3 1.68 -1.284 1.612 0.822 0.966 1.669
4 1.68 -1.284 1.583 0.733 0.949 1.669
MCOV100.val
1 -1.283
2 -1.283
3 -1.283
4 -1.283
結果を見ると，このデータ（Q1_1からQ1_10）に対してはクラスター数4のときが最大値になるようです。 図 8.28 にその時の分かれ方が示されています。 （やっぱり項目数が少ないとクラスター数が多めになるのか？）
一点注意として，expl.detect()の内部では，デフォルトではデータを半分に分割して，
- 前半のデータで条件付き共分散を計算→クラスター構造を決定
- 後半のデータで\(D\)を計算
という流れで計算を行っています。これは多分オーバーフィッティングを避けるためだと思います。 そのため，N.Clusterが2のところのDETECT.est（手順1で使用したデータで計算された\(D\)）もDETECTR.val（手順2で使用したデータで計算された\(D\)）も，conf.detect()での結果1.064とは異なる値になっています（N.estおよびN.valがいずれも元のdatの行数の半分になっている）。
…ということで，今回のデータではきれいな一次元性を確認することはできなかったのですが，このような方法を用いると，「一次元とみなして分析しても問題ないか」を評価することができるのです。
8.9 適合度
IRTでもSEMと同じように適合度を考えることができます。 特に「個人」と「項目」の両方に関心があることの多いIRTでは，それぞれに対しての適合度を考える必要があります。 もちろんこれから紹介する適合度の考え方は，そのまま因子分析やSEMにも応用しようと思えばできるので，そういった意味でも参考にしてみてください。
8.9.1 局所独立性の確認
まずはじめに，データとモデルの適合度の一つとして，データが局所独立の仮定に適合しているかをチェックします。 セクション 8.5.1 で説明したように，局所独立性が満たされている場合には，\(\theta_p\)で条件づけたときの2項目の回答の間に相関が無いはずです。
\(X^2\)統計量
\(X^2\)統計量 （Chen & Thissen, 1997） では，カテゴリ変数のクロス表に対して用いられる\(\chi^2\)検定の枠組みを利用して局所独立性を評価します。 ある二値項目\(i\)の項目パラメータ\((\alpha_i,\beta_i)\)が分かっている場合，その項目の母集団全体での期待正答率は \[ \begin{aligned} E[P(y_{pi}=1)] = \int_{-\infty}^{\infty}P(y_{pi}=1|\theta)f(\theta)d\theta \end{aligned} \tag{8.37}\] という形で求めることができます。 ここで\(f(\theta)\)は特性値の確率分布（ふつうは標準正規分布）における\(x=\theta\)での確率密度を表しています。 同様に，その項目の誤答率の期待値を求めるならば\(P(y_{pi}=1|\theta)\)の部分を\(1-P(y_{pi}=1|\theta)\)にしたら良いだけです。
この考え方を2項目に広げると，母集団全体において項目\((i,j)\)に両方とも正解する確率はもし2項目が独立ならば \[ \begin{aligned} E[P(y_{pi}=1 \land y_{pj}=1)] = \int_{-\infty}^{\infty}P(y_{pi}=1|\theta)P(y_{pj}=1|\theta)f(\theta)d\theta \end{aligned} \tag{8.38}\] と表せます。 ということは，サンプルサイズが\(N\)の場合，項目\((i,j)\)に両方とも正解する人数の期待値\(E_{11}\)は\(N\times P(y_{pi}=1 \land y_{pj}=1)\)です。 同様にして，一方の項目にだけ正解する人数および両方の項目に誤答する人数の期待値も求めると， 表 8.3 のようになります。
| 項目\(i\) | 当てはまらない | 当てはまる | 計 |
|---|---|---|---|
| 当てはまらない | \(\begin{aligned}E_{00}= \\ N\times P(y_{pi}=0 \land y_{pj}=0)\end{aligned}\) | \(\begin{aligned}E_{01}= \\ N\times P(y_{pi}=0 \land y_{pj}=1)\end{aligned}\) | – |
| 当てはまる | \(\begin{aligned}E_{10}= \\ N\times P(y_{pi}=1 \land y_{pj}=0)\end{aligned}\) | \(\begin{aligned}E_{11}= \\ N\times P(y_{pi}=1 \land y_{pj}=1)\end{aligned}\) | – |
| 計 | – | – | \(N\) |
この期待度数分布はもし2項目が局所独立ならばこうなるだろうという状態を表しており， 表 8.1 に相当するものを\(\theta\)について周辺化して作成したものといえます。 つまり局所独立ではなくなるほど，実際のデータでのクロス表は（ 表 8.2 のように）この期待度数から離れるだろう，ということです。
クロス表での連関の検定は（学部の統計でやっているかもしれない）\(\chi^2\)検定です。 その検定統計量は \[ \chi^2=\sum_{s=0}^{1}\sum_{t=0}^{1}\frac{(O_{st}-E_{st})^2}{E_{st}} \tag{8.39}\] で計算されます。ここで\(O_{st}\)は，実際のデータで\(y_{pi}=s \land y_{pj}=t\)だった人数を表しています。
こうして計算された\(\chi^2\)統計量は，自由が(行数-1) \(\times\) (列数-1)の\(\chi^2\)分布に従うことが知られています。 二値データの場合は\(1\times1=1\)ということです。 そしてここまでのプロセスはそのまま多値項目に対しても拡張可能です。 例えばdatの項目のように6件法どうしの場合は，自由度は\(5\times5=25\)となるわけです。 ということで，項目のペア毎に\(\chi^2\)統計量を計算して統計的仮説検定を行い，有意だったペアの間は局所独立ではない可能性が疑われる，ということになります。
\(Q_3\)統計量
\(Q_3\)統計量 （Yen, 1984） では，回帰分析でいうところの偏相関係数にあたるものを利用します。 偏相関係数は「変数\(z\)の影響を取り除いた時の\(x\)と\(y\)の相関係数」で，疑似相関の影響を検討する際などに利用されるものです。 これをIRTの文脈に当てはめると「変数\(\theta_p\)の影響を取り除いた時の\(y_{pi}\)と\(y_{pj}\)の相関係数」ということで，これが0であれば局所独立だろうと言えそうです。
回帰分析では，\(x\)と\(y\)をそれぞれ\(z\)で回帰した時の予測値\(\hat{x},\hat{y}\)に対して，残差の相関\(r_{(x-\hat{x},y-\hat{y})}\)のことを偏相関係数と呼んでいました。 IRT（ロジスティックモデル）もその中身はロジスティック回帰なので，同じようにして\(y_{pi}\)と\(y_{pj}\)をそれぞれ\(\theta_p\)で予測した時の予測値（期待正答率）を \[ \hat{P}(y_{pi}=1)=\frac{1}{1+\exp\left[-\hat{\alpha}_{i}(\hat{\theta}_p - \hat{\beta}_i)\right]} \tag{8.40}\] というように，\(\alpha_i,\beta_i,\theta_p\)にそれぞれ推定値\(\hat{\alpha}_i,\hat{\beta}_i,\hat{\theta}_p\)を入れることで計算が出来ます。 すると，残差を\(d_{pi}=y_{pi}-\hat{P}(y_{pi}=1)\)と求めることができます。 図 8.29 に\(d_{pi}\)を表してみました。
この項目にはAさんもBさんも正解しているため，ふたりとも\(y\)軸の値は1です。 ですが他の全項目から推定された2人の\(\hat{\theta}_p\)はそれぞれ\(-2,2\)だったとします。 この場合，Aさんは他の項目の出来から考えるとこの問題にはほぼ間違えるはずなのに正解しています。 ということで\(d_{pi}\)の値が大きくなっています。 このようになる理由としては，やはり\(\theta_p\)以外の別の要因によって正解できた，と考えるのが妥当でしょう。 これを全体に広げてみた時に，\(d_{pi}\)が高い人ほど別の項目でも同様に\(d_{pj}\)の値が大きくなっているとしたら，この2つの項目には何か正解に寄与する\(\theta_p\)以外の別の要因がある，と考えることができるわけです。
ということで\(Q_3\)統計量は \[ Q_3=r_{\symbf{d}_i,\symbf{d}_j} \tag{8.41}\] という形になります。 \(Q_3\)統計量は相関係数なので，サンプルサイズに依存しない効果量の指標として見ることができ，絶対値が0.2を超えてくると怪しいと判断できるようです。
Rでやってみる
mirtではresidualsという関数でいま紹介した指標を出してくれます。 引数typeで，出してほしい統計量を指定できます。 ちょっと厄介なのですが，\(\chi^2\)統計量を出してほしい場合にはtype = "LD"と指定してください。 LDはlocal dependenceの頭文字ですが，そういったら\(Q_3\)だってLDだろ，と思われるかもしれません。 これは Chen & Thissen （1997） の書き方に則っているのだと思います。我慢してください。
Degrees of freedom (lower triangle) and p-values:
Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10
Q1_1 NA 0 0 0.019 0 0.000 0.000 0.093 0.824 0.004
Q1_2 1 NA 0 0.000 0 0.172 0.000 0.208 0.000 0.000
Q1_3 1 1 NA 0.000 0 0.006 0.003 0.000 0.000 0.000
Q1_4 1 1 1 NA 0 0.001 0.193 0.000 0.001 0.485
Q1_5 1 1 1 1.000 NA 0.003 0.000 0.004 0.000 0.019
Q1_6 1 1 1 1.000 1 NA 0.000 0.002 0.099 0.477
Q1_7 1 1 1 1.000 1 1.000 NA 0.032 0.019 0.550
Q1_8 1 1 1 1.000 1 1.000 1.000 NA 0.045 0.121
Q1_9 1 1 1 1.000 1 1.000 1.000 1.000 NA 0.000
Q1_10 1 1 1 1.000 1 1.000 1.000 1.000 1.000 NA
LD matrix (lower triangle) and standardized values.
Upper triangle summary:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.109 -0.067 -0.026 0.002 0.063 0.219
Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10
Q1_1 NA 0.144 0.111 0.048 0.082 -0.098 -0.090 -0.034 0.005 -0.059
Q1_2 50.489 NA 0.171 0.077 0.130 -0.028 -0.071 -0.026 -0.100 -0.103
Q1_3 29.999 71.332 NA 0.102 0.219 -0.055 -0.061 -0.078 -0.109 -0.083
Q1_4 5.504 14.503 25.105 NA 0.094 -0.067 -0.026 -0.079 -0.065 0.014
Q1_5 16.380 41.085 116.151 21.292 NA -0.059 -0.080 -0.058 -0.104 -0.048
Q1_6 23.577 1.863 7.479 10.768 8.539 NA 0.101 0.063 0.033 -0.014
Q1_7 19.698 12.128 8.934 1.695 15.588 24.903 NA 0.044 0.048 0.012
Q1_8 2.820 1.583 14.878 15.199 8.112 9.698 4.621 NA 0.041 0.031
Q1_9 0.049 24.351 29.035 10.195 26.099 2.723 5.519 4.012 NA 0.123
Q1_10 8.358 25.794 16.877 0.487 5.500 0.507 0.357 2.411 37.068 NA上の行列の上三角が\(p\)値です。したがってこれが0.05を下回っている場合は「局所独立ではない」という判断になるのですが，見た感じかなり多くのペアで有意になっています。 これは，統計的仮説検定なのでサンプルサイズの影響を受けているためです。 また上の行列の下三角は自由度です。 今回は二値に変換したデータに対して行っているので自由度は1になっています。 もしも6件法のまま(GRMやGPCMで)推定した結果に対して出した場合には，自由度は先程説明した通り\(5\times5=25\)となります。
そして下の行列は実際に計算された\(\chi^2\)統計量（下三角）および標準化した値（上三角；たぶんクラメールの連関係数）です。 ということで，この行列の上三角について大きいペアが，局所独立からより強く離れているといえます。
続いて\(Q_3\)統計量です。
Q3 summary statistics:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.235 -0.166 -0.101 -0.067 0.033 0.244
Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10
Q1_1 1.000 0.150 0.109 0.033 0.074 -0.149 -0.148 -0.072 -0.027 -0.115
Q1_2 0.150 1.000 0.188 0.056 0.131 -0.101 -0.173 -0.097 -0.209 -0.193
Q1_3 0.109 0.188 1.000 0.082 0.244 -0.149 -0.170 -0.182 -0.235 -0.187
Q1_4 0.033 0.056 0.082 1.000 0.071 -0.166 -0.122 -0.186 -0.175 -0.068
Q1_5 0.074 0.131 0.244 0.071 1.000 -0.152 -0.197 -0.152 -0.226 -0.145
Q1_6 -0.149 -0.101 -0.149 -0.166 -0.152 1.000 0.054 0.005 -0.048 -0.114
Q1_7 -0.148 -0.173 -0.170 -0.122 -0.197 0.054 1.000 -0.041 -0.047 -0.105
Q1_8 -0.072 -0.097 -0.182 -0.186 -0.152 0.005 -0.041 1.000 -0.045 -0.069
Q1_9 -0.027 -0.209 -0.235 -0.175 -0.226 -0.048 -0.047 -0.045 1.000 0.040
Q1_10 -0.115 -0.193 -0.187 -0.068 -0.145 -0.114 -0.105 -0.069 0.040 1.000ところどころ高めの相関が出ており，完全な局所独立ではないといえそうです。 というのも今回のデータは，本来2因子に分かれるはずの10項目を無理やり1因子とみなして分析しています。 そもそもデータが一次元性を持っていないはずなので，局所独立性も保たれてはいない，ということになります。 （本来こういう時には一旦立ち止まり，項目の内容を吟味したりモデルを変更したりといった手段を考えるべきですが，ここでは気にせず進めます。）
ちなみに，引数suppressを与えると，絶対値がsuppress以下の値がNAに置き換わるため，項目がたくさんある場合でも「閾値を超えているペア」を探しやすくなるかもしれません。
8.9.2 個人適合度
続いて個人の適合度です。 IRTにおける個人適合度の基本的な考え方としては「\(\theta_p\)が高いのに困難度が低い項目に間違えるのはおかしい」「\(\theta_p\)が低いのに困難度が高い項目に正解するのはおかしい」というものです。 するとこれは， 図 8.29 的な考え方になります。 Aさんのほうが予測値と実測値の間の差\(d_{pi}\)が大きいことは，前述の通り「\(\theta_p\)以外に回答に影響を与える要因がある（のでモデルが間違っている）」という可能性の他に，「（もしモデルは正しいと仮定したら）回答者の方に，何か項目\(i\)の正解に導いた要因があるのでは」と考える，というわけです。
回答者の方の原因としては，心理尺度でいえば以下のようなものが考えられます（Meijer et al., 2016）。
- 回答者の言語能力が不十分なために項目の内容を十分に理解できていない可能性。低学齢への実施や第一言語以外での実施で起こる可能性がある。
- 心理尺度が測定しようとしている心理特性に対する認識がない状態（Traitedness）。例えば「そんなこと考えたこともない」といった場合，項目の意図を他の人とは異なるニュアンスで捉えてしまい，結果的に「大多数が当てはまると回答する項目に当てはまらないと回答する」といった反応につながる。
- 回答のモチベーションが低い可能性。やる気が無いので適当に回答している？
- 特定の心理特性を持っている人は平均的な人より一貫性の低い回答をする可能性がある。例えば精神的な問題を抱えていたり，希死念慮を持ち合わせている人はそうでない人と比べて適合度が低いという先行研究もあるらしい。
ということで，個人適合度が低いイコール適当な回答者（だから除外してもよい）とも言い切れないのですが，特徴的な回答者をあぶり出すなどの目的でも個人適合度は使い道がありそうです。
実際に個人適合度の指標として提案されているものは相当数ある （Karabatsos, 2003; Meijer & Sijtsma, 2001） のですが，わかりやすい＆統計的仮説検定のフレームワークが整備されていると言った理由から，\(z_h\)（またの名を\(l_z\)）という統計量（Drasgow et al., 1984）が用いられることが多いようです （Meijer et al., 2016）。
\(z_h(l_z)\)統計量
\(z_h(l_z)\)統計量では，ある回答者の実際の反応パタン\(\symbf{y}_p=\begin{bmatrix}y_{p1}&y_{p2} &\cdots&y_{pI}\end{bmatrix}\)と「考えられる全反応パタン」での尤度を比較することにより，その回答者の反応パタンの「普通さ」を評価します。
まず，回答者\(p\)の特性値の最尤推定値\(\hat{\theta}_p\)に関して，その時の尤度を計算します。 局所独立性が満たされているならば，尤度は単純に全ての項目をかけ合わせたら良いので， \[ l_0=\prod_{i=1}^{I}P(y_{pi}=1)^{y_{pi}}P(y_{pi}=0)^{1-y_{pi}} \tag{8.42}\] となります。 ただ掛け算は（コンピュータ的に）大変なので，対数尤度 \[ ll_0 = \sum_{i=1}^{I}y_{pi} \ln P(y_{pi}=1) + (1-y_{pi})\ln P(y_{pi}=0) \tag{8.43}\] を使います。 この尤度\(ll_0\)は，「特性値が\(\hat{\theta}_p\)の人が\(\symbf{y}_p\)という回答をする確率」のようなもの21です。 ここで\(\symbf{y}_p\)以外の回答パタンに目を向けてみると，もしかしたら「\(\symbf{y}_p\)よりも\(\hat{\theta}_p\)と整合的な回答パタン」があるかもしれません。 もちろん反対に「\(\hat{\theta}_p\)と全く合わない回答パタン」というものもあるかもしれません。
そこで，考えられる全回答パタン（2件法10項目ならば\(2^{10}\)通り）の中で\(\symbf{y}_p\)という回答パタンが\(\hat{\theta}_p\)とどの程度整合的かを相対的にチェックすることを考えます。 全回答パタン内で\(\symbf{y}_p\)という回答パタンと\(\hat{\theta}_p\)の整合性（尤度）が平均くらいであれば，その回答パタンはまあ起こりうるものだろうとみなせる一方で，もし整合性が低い場合には，\(\symbf{y}_p\)という回答パタンは起こりにくい（ズレが大きい）だろうと判断できる，ということです。
実際に「全回答パタンでの対数尤度の期待値」は（各項目への回答は離散変数なので） \[ E_h = \sum_{i=1}^{I}P(y_{pi}=0)\ln P(y_{pi}=0) + P(y_{pi}=1)\ln P(y_{pi}=1) \tag{8.44}\] として求めることができます。 つまり\(ll_0-E_h\)が「平均的な回答パタンでの尤度と実際の回答パタン\(\symbf{y}_p\)での尤度の差」であり，これがマイナス方向に大きいほどモデルとデータの適合度が悪いということを意味します。
ここで全回答パタンにおける対数尤度の分散を \[ \sigma_h^2 = \sum_{i=1}^{I}\left\{P(y_{pi}=1)P(y_{pi}=0)\left[\ln\frac{P(y_{pi}=1)}{P(y_{pi}=0)}\right]^2\right\} \tag{8.45}\] として求めることができるのですが，これを用いると，\(ll_0\)を全回答パタン内で標準化した値\(z_h=\frac{ll_0-E_h}{\sigma_h}\)が経験的に標準正規分布に従うということが知られています。 ということで，\(|z_h| > 1.96\)となった回答者はどうやら「特性値が\(\hat{\theta}_p\)のもとでは平均的ではない回答パタンである」と判断することができるわけです。
図 8.30 に，\(z_h\)統計量の考え方を表しました。 ヒストグラムは母集団からサンプリングされた様々な\(\hat{\theta}_p\)の人と，その人の（ランダムサンプリングされた）回答ベクトル\(\symbf{y}_p\)のもとで計算した\(z_h\)を集めたものというイメージです。 これが近似的に正規分布に従っているということで，この外側5%を棄却域とみなして検定を行うことができます。 Bさんの\(z_h\)は，この分布の中で割と平均的なところにあるため，Bさんの回答は「（Bさんの推定値\(\hat{\theta}_{\mathrm{B}}\)における）ふつう」に近いとみなせます。 一方Aさんの\(z_h\)は，この分布の中でもかなり小さい方にあります。 ということは，Aさんの回答パタンはAさんの推定値\(\hat{\theta}_{\mathrm{A}}\)からすると当てはまりが悪く（つまり困難度の低い項目に「当てはまらない」と回答し，困難度の高い項目で「当てはまる」と回答している可能性が高く），なんなら母集団全体の中でも相当悪いほうだと判断でき，何らかの問題があるのではないかと疑いをかけられるわけです。
ちなみに，正規分布的には上側の外れ値も「ふつうではない」ということになりますが，こちらはどう考えたら良いかよくわかりません。 理論的には「特性値\(\hat{\theta}_p\)に対してあまりにも適合しすぎている」という状態で，これは例えば困難度が\(\hat{\theta}\)以下の項目では全て「当てはまる」と回答し，\(\hat{\theta}\)以上の項目ではすべて「当てはまらない」と回答しているような感じなのですが，だからといって即座に問題とは言えなさそうです。 ということで，基本的には\(z_h < -1.96\)の回答者に着目しておけば良いと思います。
Rでやってみる
※個人適合度と回答パタンの関係がわかりやすいので，ここはGRMで推定した結果を使用します。
mirtでは，personfit()という関数によって\(z_h\)統計量を出すことが出来ます。
outfit z.outfit infit z.infit Zh
1 0.3043410 -2.51418387 0.2901436 -2.630507534 1.11739286
2 0.4319025 -1.53744653 0.4991522 -1.313386687 1.54458648
3 0.5417319 -1.17103185 0.6082584 -0.964496321 1.30478354
4 0.6961604 -0.71567303 0.7021298 -0.707746942 0.38159507
5 0.5979698 -0.95977207 0.5796844 -1.045595986 0.06132109
6 1.6611227 0.99612493 2.2816175 1.585552919 -0.57702558
7 0.2629479 -2.06106332 0.2970130 -1.974928887 1.15924498
8 1.2820317 0.82619048 1.2682515 0.803602305 -1.12899488
9 0.4296002 -0.79327083 0.4483505 -0.806286109 0.88382209
10 0.6018210 -0.83208882 0.6646485 -0.682369844 0.51748259
11 0.7508775 -0.39482320 0.8851335 -0.095554151 0.66704225
12 0.9285353 0.04354098 1.1524578 0.453682964 0.30674141
13 1.9537886 1.83531480 1.9490471 1.865766458 -1.36629763
14 1.0385754 0.23709522 1.0081517 0.168761287 -0.60996496
15 0.8936027 -0.09915012 0.9407077 0.005928737 0.05273411
16 0.6633546 -0.62999115 0.7489752 -0.432539405 0.73371137
17 0.5941524 -1.07615247 0.6027536 -1.058618799 0.98045292
18 1.7516559 1.78554390 1.7275073 1.788745281 -1.75844255
19 1.4892405 1.27332749 1.5781245 1.468740378 -1.18728452
20 1.3794907 0.81105690 1.5238636 1.050442554 -0.07950879
[ reached 'max' / getOption("max.print") -- omitted 2412 rows ]outfitやinfitというものは，\(z_h\)よりもシンプルに予測値と実測値の間の差\(d_{pi}\)の二乗和をもとにした適合度指標（Smith et al., 1995）です22。 ということでこれらについても大きいほど適合度が悪いと判断できます。 実際にこれらの指標の相関を見てみると，いずれもかなり高い相関になっていることがわかります。
outfit z.outfit infit z.infit Zh
outfit 1.0000000 0.9029153 0.9767846 0.8895739 -0.8593484
z.outfit 0.9029153 1.0000000 0.8962045 0.9898224 -0.8803762
infit 0.9767846 0.8962045 1.0000000 0.9083434 -0.8420095
z.infit 0.8895739 0.9898224 0.9083434 1.0000000 -0.8617816
Zh -0.8593484 -0.8803762 -0.8420095 -0.8617816 1.0000000では実際に，当てはまりが悪い人の回答パタンを見てみましょう。
いま扱っているモデルは10項目1因子のモデルです。 そして10項目の識別力は，いずれも正の値になっています。 つまり，基本的にはすべての項目間には正の相関があるため，この回答者のように「ある項目では1，別の項目では6」という回答の付け方には違和感があります。
Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 outfit z.outfit
45 6 5 3 2 3 3 6 3 6 4 1.362101 0.8901956
75 6 6 3 3 1 6 6 5 6 1 2.777980 2.4882677
96 6 6 4 1 1 6 4 5 5 6 2.782816 2.5103934
98 1 1 1 1 2 3 6 3 4 3 2.006426 2.2242279
122 6 6 6 6 6 6 6 3 1 2 3.089210 2.3332034
123 6 6 6 5 6 6 1 3 3 2 2.065644 1.7906200
126 4 6 6 4 6 2 3 2 2 2 1.459869 1.1442225
172 5 2 1 6 3 6 6 6 6 6 4.695997 3.3697866
186 6 5 2 2 2 2 2 6 5 4 1.817030 1.8103247
226 1 1 1 4 6 5 6 5 6 6 4.784881 3.9067830
infit z.infit Zh
45 1.401457 0.9775542 -2.193202
75 2.751967 2.5542871 -3.279082
96 2.855421 2.6749604 -2.107685
98 2.099394 2.4577322 -3.172528
122 3.034305 2.3832158 -2.030530
123 1.955780 1.7068100 -1.971380
126 1.314690 0.8576219 -1.980124
172 3.651095 2.8513846 -3.188176
186 1.791534 1.7810377 -2.201609
226 4.276645 3.7176248 -4.057956
[ reached 'max' / getOption("max.print") -- omitted 121 rows ]同様にして，\(z_h < -1.96\)の人たちを見てみると，項目ごとにかなりバラバラな回答をしている人たちがあぶり出されています。 ということで，\(z_h\)統計量によって怪しい回答者をあぶり出すことができそうだということがわかりました。
ちなみに反対に\(z_h > 1.96\)の人たちを見たところですべての項目で同じような値をつけており，やはり「適合度が悪い」という判断はなかなか難しそうです（ストレートライン気味な感じはありますが）。
Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 outfit z.outfit
119 5 4 4 4 4 4 4 4 3 2 0.1059937 -3.933211
136 6 5 5 6 5 5 5 5 5 5 0.1749942 -2.032175
221 2 2 1 1 2 2 1 2 2 1 0.3892825 -1.326634
674 5 5 5 4 5 5 4 4 5 2 0.2410515 -2.270639
684 4 4 5 5 4 4 4 4 3 2 0.1893306 -3.039619
694 5 5 5 4 5 5 4 4 3 2 0.2993084 -2.123527
739 5 5 5 6 4 4 4 4 3 3 0.2791157 -2.191891
804 5 5 5 5 4 5 5 5 5 3 0.1527831 -2.557928
826 5 5 4 4 4 4 4 3 3 1 0.2959759 -2.472784
834 5 5 4 5 5 4 4 4 3 3 0.1765065 -2.919233
infit z.infit Zh
119 0.1122620 -3.899088 2.556666
136 0.2246921 -1.890473 2.065611
221 0.4313157 -1.407809 2.117805
674 0.2776608 -2.150842 1.970796
684 0.2021914 -2.982029 2.001422
694 0.3096588 -2.123376 2.017472
739 0.3059683 -2.111215 1.971686
804 0.1559106 -2.653563 1.963533
826 0.3014073 -2.463963 2.015447
834 0.1741386 -3.000314 2.105461
[ reached 'max' / getOption("max.print") -- omitted 9 rows ]ここで紹介した\(z_h\)統計量のような個人適合度は，モデルベースで適当回答者を除外するのに使うこともできます。 ただし，個人適合度はあくまでも（ロジスティックなどのIRT）モデルに基づいて推定されたパラメータをもとに算出される点には少し注意が必要です。 そもそもモデルの式や因子構造の設定が間違っている場合には，正しくない検出をしてしまうかもしれません。 そんなわけで，もちろん適合度の値だけで機械的に除外するのは危ないですが，例えば項目パラメータの推定値がなんだかおかしい場合には，適合度指標を参考にしたデータクリーニングによって，推定をより安定させられるかもしれません。
8.9.3 項目適合度
IRTでは項目一つ一つの性能（識別力や困難度）に関心があることが多いため，項目レベルでも適合度を考えます23。 基本的な考え方としては，推定されたパラメータに基づくICCが，データとどの程度一致するかを確認していきます。 適合度が悪い項目が見つかった場合には，項目の内容を吟味して削除するかを検討したり，場合によってはモデル自体を修正（因子数や因子構造の修正など）する必要があるかもしれません。 項目適合度は，具体的には以下のような手順によって算出していきます （Orlando & Thissen, 2000） 。
- 項目パラメータと\(\theta_p\)を推定する
- 回答者を\(\theta_p\)の値で並び替える
- 回答者を\(\theta_p\)の値によっていくつかのグループに分ける
- 各グループでの反応確率を計算する
- 4.で求めた反応確率と「モデル上の期待反応確率」を\(\chi^2\)的な統計量や視覚的方法を用いて比較する
視覚的なチェック
ということで，まずは視覚的に比較してみましょう。 mirtでは，項目適合度を見るための関数としてitemfit()が用意されています。 これに引数empirical.plotを与えることでプロットを見ることができます。
図 8.31 を見ると，各カテゴリごとにICCが表示されています。これに重なっている小さいマルのY座標は，手順3で分割した各グループ（ここでは10グループ）における\(P(y_{pi}=1)\)を表しています（多値型の場合はカテゴリ選択確率が表示される）。 また，X座標はそのグループにおける\(\theta_p\)の平均値を表しています。 ということで，各カテゴリにおいて，小さいマルと実線が近いほどデータとモデルの当てはまりが良い，と判断することができるわけです。
ただカテゴリ数や項目数が多くなってくるとすべてを目視で確認するのはかなり大変です。 また視覚的なチェックだけではどの項目が最も乖離しているのかを判断するのも大変でしょう。 ということで，怪しい項目にアタリをつけるための統計量を使っていきましょう。
カイ二乗的な統計量
\(\chi^2\)統計量は独立性の検証のところ（ セクション 8.9.1 ）で登場したものです。 IRTモデルに基づいて項目および回答者のパラメータが推定されると，各項目について，指定されたサンプルサイズのうち何人がどの選択肢を選ぶと期待されるか，を計算することができます。 これを期待度数\(E\)，また実際にその項目に対する反応の分布を観測度数\(O\)として \[ \chi^2=\sum\frac{(O-E)^2}{E} \tag{8.46}\] と計算すると，「帰無仮説：モデルがデータに適合している」が正しいならば小さい値におさまるはずだと考えられます。
Yen （1981） は，回答者を\(\theta_p\)の値の高低によってグループ分けし，グループごとに\(\chi^2\)の計算をすることにしました。 特にグループ数を10に固定して算出した指標は\(Q_1\)と呼ばれています。 グループ\(g~(g=1,2,\cdots,10)\)における特性値\(\hat{\theta}_p\)の平均値を\(\bar{\theta}_g\)とおくと，期待正答率は（2パラメータロジスティックモデルならば） \[ E_{ig} = \frac{1}{1+\exp\left[-\hat{\alpha}_{i}\left(\bar{\theta}_g - \hat{\beta}_i\right)\right]} \tag{8.47}\] と求めることができます。 すると，二値項目における\(Q_1\)統計量は \[ \begin{aligned} Q_1 &= \sum_{g=1}^{10}N_g\left(\frac{(O_{ig}-E_{ig})^2}{E_{ig}} + \frac{\left[(1-O_{ig})-(1-E_{ig})\right]^2}{1-E_{ig}}\right) \\ &= \sum_{g=1}^{10}N_g\left(\frac{(1-E_{ig})(O_{ig}-E_{ig})^2 + E_{ig}(O_{ig}-E_{ig})^2}{E_{ig}(1-E_{ig})}\right) \\ &= \sum_{g=1}^{10}N_g\left(\frac{(O_{ig}-E_{ig})^2}{E_{ig}(1-E_{ig})}\right) \end{aligned} \tag{8.48}\] となります。 なお1行目の右辺は，第1項が正答セルにおける乖離度を表しています。 同様に，第2項は\((1-O_{ig}),(1-E_{ig})\)に置き換わっていることから，誤答セルにおける乖離度を表していることがわかります。 したがって，\(Q_1\)統計量は「グループ」\(\times\)「項目\(i\)の正誤」という\(10\times2\)クロス表における\(\chi^2\)統計量を表しているのです。 そしてこの\(Q_1\)統計量は自由度が\(10-\)パラメータ数の\(\chi^2\)分布に従うため，これを用いて検定ができるわけです。
同様に，差ではなく比を用いた統計量として\(G^2\)統計量 （McKinley & Mills, 1985） というものもあります。 \[ G^2 = \sum_{g=1}^{10}N_g\left[O_{ig}\ln\frac{O_{ig}}{E_{ig}} + (1-O_{ig})\ln\frac{1-O_{ig}}{1-E_{ig}}\right] \tag{8.49}\]
ほかにもいくつか適合度指標は提案されているのですが，基本的には前述の通り「期待反応確率と実際の反応確率の差」を何らかの形で統計量に落とし込んでいます。 細かな設定としては
- グループの数
- \(\theta_p\)の推定方法（基本的には最尤法，ただベイズ推定などでも可）
- 各グループの特性値の代表値の決め方（平均値 or 中央値）
- 期待反応確率の計算方法（個人の\(\hat{\theta}_p\)における期待反応確率\(E_{ip}\)のグループ平均という手もあり）
といったところで様々なオプションが考えられるようです（加藤 他，2014）。
Rでやってみる
項目適合度をを出す場合はitemfit()関数を使います。 引数fit_statsに，どの統計量を出してほしいかを指定します。なお，X2を指定すると内部では\(Q_1\)統計量（i.e., グループ数10）を出してきます。
item X2 df.X2 RMSEA.X2 p.X2 G2 df.G2 RMSEA.G2 p.G2
1 Q1_1 375.251 8 0.137 0 458.736 6 0.176 0
2 Q1_2 73.800 8 0.058 0 99.817 5 0.088 0
3 Q1_3 92.664 8 0.066 0 135.273 5 0.104 0
4 Q1_4 102.569 8 0.070 0 153.430 5 0.111 0
5 Q1_5 126.263 8 0.078 0 169.111 5 0.116 0
6 Q1_6 157.177 8 0.088 0 191.087 5 0.124 0
7 Q1_7 151.358 8 0.086 0 219.612 4 0.149 0
8 Q1_8 125.809 8 0.078 0 154.268 6 0.101 0
9 Q1_9 151.495 8 0.086 0 207.812 5 0.129 0
10 Q1_10 572.082 8 0.170 0 754.816 5 0.248 0G2の自由度(df.G2)がいろいろな値になっていますがあまり気にしなくてOKです。 今回も例によって検定の結果(p.*)は全て有意となりました。 このように，サンプルサイズが大きい場合これらの統計量はあまり参考にならない可能性が高いです。
…じゃあどうするんだ，という話になるわけですが，一つの方法としては効果量的な考え方をとるために「実際の差\((O_{ig}-E_{ig})\)や比\(\frac{O_{ig}}{E_{ig}}\)を計算する」という方法が考えられます。 例えば\(\frac{1}{10}\sum_{g=1}^{10}|O_{ig}-E_{ig}|\)が大きい項目は，基本的に視覚的にプロットを見ても差があるでしょう。 ただ\((O_{ig}-E_{ig})\)などを直接計算してくれる関数は用意されていないので，多少頑張る必要があります。
8.10 項目情報量・テスト情報量
IRTにかぎらず，統計的推測では推測の精度が問われます。 例えば正規分布の母平均の推測の場合，標準誤差は母分散\(\sigma^2\)およびサンプルサイズ\(n\)を用いて\(\sqrt{\frac{\sigma^2}{n}}\)と表されました。 この式は，データの数が多いほど推測の精度が高くなるという自然な考えを表しています。
テストや心理尺度では，標準誤差を真値と誤差の関係から考えます。 セクション 4.6 で紹介しましたが，潜在特性の測定（古典的テスト理論）では，観測得点\(x\)と真値\(t\)と誤差\(e\)の関係を \[ x = t + e \tag{8.50}\] と表し，信頼性係数を \[ \rho = 1-\frac{\sigma_e^2}{\sigma_x^2} \tag{8.51}\] と定義していました。 標準誤差は「真値が\(t\)の人が測定を繰り返した場合に，（平均的に）どの程度観測得点\(x\)がばらつくか」を表しています。 したがって，\(x\sim N(t,\sigma_e^2)\)とすると，心理測定における標準誤差は\(\sigma_e\)にあたります。 信頼性係数\(\rho\)に関しては，上の式を変形させると \[ \sigma_e = \sigma_x^2\sqrt{1-\rho^2} \tag{8.52}\] となることから，これが間接的に測定の精度を表す指標だということが言えます。 なお実際に標準誤差を計算する場合には，\(\rho\)の推定値として\(\alpha\)係数などを代入します。
さて，信頼性係数に基づく上記の標準誤差の式には真値\(t\)は含まれていません。 つまり真値の値が何であれ標準誤差は同じということを意味しています。 ですがこの考え方は割と不自然なものです。 特性値の真値が平均くらいの人では，尺度の項目はほとんどすべてがその人の特性値を測定するのに意味のある（正解であるほど真値\(t\)は高いだろうという判断につながる）項目です。 一方で，特性値の真値が非常に高い人では，実はたいていの項目はその人の特性値を測定するのにほぼ無意味です。 というのも，ある程度難易度の高い項目に正解する人は，難易度の低い項目にはほぼ確実に正解するだろうと言えるためです。 つまり，「ある程度難易度の高い項目」の情報さえあれば「難易度の低い項目」の情報は無くても特に問題ない（真値の推定には影響しない）と言えます。 先ほど，標準誤差はサンプルサイズが多くなるほど小さくなる，ということを説明しました。 本来，特性値の真値によって「その人の特性値を測定するのに意味のある」項目の数は変わるはずだと考えると，標準誤差も個人（の特性値の値）ごとに変わると考えるのが自然な気がしてきます。
IRTでは，この問題を克服するために項目情報量 (item information)およびテスト情報量 (test information)というものを用います。 この項目情報量・テスト情報量は，後述するテストの設計などにおいても大きな力を発揮してくれます。
例えば大学受験生の数学力を測定するテスト（大学入学共通テストなど）を考えてみると，「東大を狙うレベルの受験生」と「大学のレベルは問わない受験生」では，どう考えても「その人の能力を測定する上で意味のある問題」は異なるはずです。 そのように考えると，特性値の真値（受験者の能力のレベル）によって標準誤差が異なるはずだ，という点は結構重要なことだと思われるかもしれません。
実際にこのようなテストで標準誤差の不均一性をあまり気にしないのは，テストの設計時に細心の注意を払っているためと言えます。 例えば大学入学共通テストやTOEICなどでは，簡単な問題から難しい問題までをバランスよく組み合わせることで，少なくともそのテストの目的である範囲の能力の測定においては十分な精度で測定ができるようになっています（たぶん）。 一方で，もちろんそのテストが意図した範囲外の回答者に対して同様に使えるテストであることは保証されません。 例えば大学入学共通テストは，中学生の数学力を測るのにも，また超天才の数学力を測るのにも不向きでしょう。
まとめると，古典的テスト理論の範囲では，ある特定の範囲においては十分な精度で特性値を測定可能であることを求めています。 そして，実際には特性値のレベルごとに異なる精度（標準誤差）を，平均化（周辺化）することで一つの値として表しているのです（これが，(-式 8.51)式中の観測値の分散\(\sigma_x^2\)に表れています）。
反対に言えば，このあと紹介する「情報量」の考え方を利用すると，上記の課題をクリアすることできます。 つまり，もともとある特定の範囲を対象に設計したテストを，その範囲外に適用せざるを得なくなったときにどのような精度になるかが分かったり，あるいは同じ特性値の異なる範囲を測定する異なるテスト間で測定精度を統一的に評価したり，といったことが可能になります。
8.10.1 情報量とは
項目情報量の話に入る前に，少しだけ「情報量」とは何かについてごくごく簡単に説明しておきましょう。 我々は，一般名詞として「情報」という言葉を使うことが多々あります。 この時の「情報」は，何かしらの事柄について「知らなかった」状態から「知っている」状態に遷移させるものという見方ができます。
この「情報」を確率の世界に持ち込むため，以下の例え話における「情報」について考えてみましょう。
この例え話からは，情報量に関する以下の性質がわかります。
- 情報は予測の精度を高める： 情報量が多くなるほど，正確に予測ができるようになります。
- より不確実性を下げる情報のほうが情報量が多い： マークは4通りである一方で，数字は13通りあるため，数字の情報のほうが不確実性をより大きく下げている＝情報量が多いと解釈できます。
- 独立な情報の情報量は加法的である： マークと数字は互いに独立なので，先に「エースを持っている」という情報を得たとしてもカードを当てる確率はやはり13倍（\(1/52\)から\(1/4\)）になっていました。
ここまでの話をIRTに置き換えてみます。 \(\theta_p\)の予測に際して，項目がもつ情報量を考えてみましょう。 図 8.32 には困難度がそれぞれ\((0,3)\)の2つの項目のICCが描かれています。
ある回答者\(p\)は，それまでの項目への解答状況から\(\hat{\theta}_p=0\)だと予測されているとします。 すると，この回答者\(p\)にとって，緑の項目（\(\beta_i=0\)）は正解するかどうかがかなり不確実（50%）です。 そしてこの項目への反応によって\(\hat{\theta}_p=0\)の予測は変化し，多少なりともより正確なものになるでしょう。 そう考えると，緑の項目は\(\hat{\theta}_p=0\)の人にとってそれなりの情報量を持っていると考えられます。
一方で赤い項目（\(\beta_i=3\)）を見ると，\(\hat{\theta}_p=0\)の人はほぼ確実に正解しないだろうという予想がつきます。 つまり，出題時点で赤い項目に対する不確実性はほぼゼロ，ということです。 実際にこの回答者\(p\)が赤い項目に不正解だったとしても，「そうでしょうね」という話なので\(\theta_p\)の予測に対して特に新規の情報をもたらすものとは言えません24。 同様に，緑の項目でも\(\hat{\theta}_p=3\)の人にとってはほぼ情報量のない項目，ということになります。 このように，ある項目は，その項目の困難度が回答者の\(\theta\)に対してちょうどよいほど，\(\theta\)の推定に対して多くの情報量を持っているのです。 ということである項目が持つ項目情報量は，実際には\(\theta\)の値によって変動する関数になっています。 この関数のことを項目情報関数 (item information function [IIF])と呼びます25。 この節のはじめにお話したように，IRTでは真値によって標準誤差が変わるようにすることを考えます。 この後具体的な関数の形を示しますが，項目情報関数を\(\theta\)の関数として表すことによってその目的を達成しようというわけです。
8.10.2 項目情報関数
IRTにおける情報量の考え方のイメージを確認したので，具体的な項目情報関数の式を見てみましょう。 二値型モデルの場合，項目情報関数\(I_i(\theta)\)は \[ I_i(\theta)=\frac{P'(y_{i}=1)^2}{P(y_{i}=1)(1-P(y_{i}=1))} \tag{8.53}\] という式になります。ここで\(P'(y_{i}=1)\)は，\(P(y_{i}=1)\)を\(\theta\)で微分した導関数です。 …といってもよくわからないと思うので，実際に見てみましょう。 図 8.33 は， 図 8.32 に示したICCを持つ2項目のロジスティックモデルにおける項目情報関数です。 どちらの線も先程説明したように，正解するかが最も不確実な\(\theta_p=\beta_i\)のところで情報量が最大になっています。
IIFの式 \[ I_i(\theta)=\frac{P'(y_{i}=1)^2}{P(y_{i}=1)(1-P(y_{i}=1))} \tag{8.54}\] を求めるためには，\(P(y_{i}=1)\)の導関数を求める必要があります。 正規累積モデルは項目反応関数の中に積分が含まれているためこれを求めるのはかなり大変なのですが，それと比べるとロジスティックモデルでは解析的にIIFの式を求めることができます。 例えば2パラメータロジスティックモデルの場合，IIFは \[ I_i(\theta)=\alpha_i^2P(y_{i}=1)(1-P(y_{i}=1)) \tag{8.55}\] となります。 この式からも，IIFは\(P(y_{i}=1)=0.5\)となる\(\theta=\beta_i\)の点において最大値をとることがわかります。 同時に，識別力\(\alpha_i\)が大きいほど項目情報量も多くなることがわかります。
図 8.34 に，識別力が異なる2つの項目(\(\beta_i=0\); ICCは 図 8.8 ) のロジスティックモデルでのIIFを示しました。 識別力が高い項目ほど，項目特性関数の傾きが大きくなっていました。 そのため\(\theta_p\)が高い人は正解し\(\theta_p\)が低い人は不正解する，というメリハリがついていました。 \(\theta_p\)の推定という視点から言えば，\(\alpha_i\)が大きい項目ほどその項目への正誤が特性値の推定により確かな情報を与えるという意味で，\(\alpha\)と項目情報量には密接な関係があるわけです。
しかし改めて 図 8.34 をよく見ると，\(\theta_p=0\)付近では識別力の高い緑の項目のほうが高い情報量を持っている一方で，\(\theta_p\)が0から離れたところではむしろ識別力の低い青い項目のほうが高い情報量を持っています。 そして対応する 図 8.8 を見ると，識別力\(\alpha_i\)が高くなるほど，\(\theta_p\)が\(\beta_i\)から離れると急激に\(P(y_{i}=1)\)が0.5から離れていることが分かります。 つまり「正解するかが不確実」な（≒情報量を持ちうる）\(\theta_p\)の範囲は識別力の高さと反比例しており，その結果\(\theta_p\)が\(\beta_i\)から離れたところでの項目情報量も少なくなってしまうのです。
この事実は「識別力は高ければ高いほど良いわけではない」ことを教えてくれています。 極端な例として「\(\theta_p<\beta_i\)の人は確実に間違えるが，\(\theta_p\geq\beta_i\)の人は確実に正解する」項目を考えてみましょう。 その項目は識別力\(\alpha_i\rightarrow\infty\)であり， 図 8.35 のようなICCをもつ項目です。
この項目は「\(\theta_p\)が\(\beta_i\)より上か下か」に関しては完全な情報を教えてくれます(\(I_i(\theta_p=\beta_i)=\infty\))。 その一方で\(\theta_p>\beta_i\)の範囲では，\(\theta_p\)が\(\beta_i\)よりどの程度高かろうと\(P(y_{i}=1)\)は100%（逆もまた然り）なので，「\(\theta_p\)が\(\beta_i\)よりどの程度上か下か」に関しては全くわかりません(\(I_i(\theta_p\neq\beta_i)\simeq0\))。 このように，識別力が高すぎる項目は有効な\(\theta_p\)の範囲が狭すぎるため，多くの回答者にとって意味のない項目になってしまう可能性があるわけです。
（いつか書きたいという気持ちだけは常に持っている）
8.10.3 テスト情報関数
ここまで，項目情報関数の定義およびその性質を見てきました。 この「情報」は，\(\theta_p\)を推定するための「情報」を表しているわけですが，\(\theta_p\)の推定は尺度・テスト内のすべての項目への回答を総合して行われるものです。 そのため，次は項目情報関数を拡張して尺度全体での情報関数を考えてみます。
情報量には「独立な情報の情報量は加法的である」という性質がありました。 IRTでもこの性質は健在です。 もしテスト全体が局所独立の仮定を満たしているならば，テスト情報関数(test information function [TIF]) は単純に \[ I(\theta) = \sum_{i=1}^{I}I_{i}(\theta) \tag{8.56}\] と表せます。 そしてテスト情報量の重要な性質として，真の特性値が\(\theta_p\)の人における最尤推定量\(\hat{\theta}_p\)の誤差分散は漸近的に \[ \sigma^2_{(\hat{\theta}_p|\theta_p)}=\frac{1}{I(\theta_p)} \tag{8.57}\] となることが知られています。 ということで，標準誤差は \[ \sigma_{(\hat{\theta}_p|\theta_p)}=\frac{1}{\sqrt{I(\theta_p)}} \tag{8.58}\] となります。 これは古典的テスト理論の表記に合わせると，最尤推定値\(\hat{\theta}_p=\theta_p+\varepsilon_p\)とした場合の\(\sigma_{\varepsilon}\)に相当するものです。 このように，IRTでは測定の精度を表す標準誤差を，\(\theta_p\)の関数として求めることができるわけです。
8.10.4 Rで情報量を確認する
項目・テスト情報関数の考え方がわかったところで，実際の項目での関数を確認してみましょう。 これまでICCを出すために使っていたplot()関数によって，IIFも出すことができます。 項目情報関数を出す場合にはtype = "infotrace"とします。
先ほど説明した通り，低識別力の項目1は他の項目と比べて圧倒的に情報量が少ないことがわかります。 また、\(\theta_p>2\)のエリアでは、Q1_10以外のほぼどの項目も情報量が無いため、この10項目では特性値が高い人についてはまともな推定はできなさそうだということがわかります。 この項目情報関数を見れば，例えばどうしても1項目削除しないといけないとしたら，基本的にはQ1_1を削除したら良いが，もしも\(\theta_p=-6\)付近がメインターゲットだとしたらむしろQ1_1は残したほうが良いとか、\(\theta_p>2\)がメインターゲットならば現状の項目セットではどうしようもないので新規項目を追加する必要がある、というように状況に応じた決定をサポートしてくれるでしょう。
具体的に特定の\(\theta_p\)の値における項目情報量を出したい場合にはiteminfo()関数を使います。
続いてテスト情報関数です。 こちらは引数をtype = "info"とすることで表示できます。
尺度（10項目）全体を見ても，項目困難度がまんべんなく存在している-3から0付近の情報量が高い一方で，すべての項目の困難度が0以上なので\(\theta_p\)が正のエリアでは情報量がさほど高くないことがわかります。
プロットする際に引数typeを"infoSE"にすると，標準誤差を合わせて出してくれるようになります。
もしかしたら，以下のようなエラーが出てプロットが表示されないかもしれません。
Error in eval(expr, envir, enclos): latticeExtra package is not available. Please install.その場合は，指示に従ってinstall.packages("latticeExtra")をしてから再度実行してください。
最も情報量が高い\(\theta_p=-1.4\)付近では，標準誤差はだいたい\(\frac{1}{\sqrt{2.82}}\simeq0.60\)くらいになっています。 仮に\(\theta_p\)の真値が\(-1.4\)の人が大量にいた場合，その人達の推定値は平均で0.6くらいは上下するだろう，ということです。 裏を返せばその程度の精度でしか推定できていないならば，例えば\(\hat{\theta}_p=-1.3\)の人と\(\hat{\theta}_p=-1.5\)の人がいたとしても，前者のほうが真の特性値が高いとはなかなか言い切れないわけですね。
具体的に特定の\(\theta_p\)の値におけるテスト情報量を出したい場合にはtestinfo()関数を使います。
\(\theta_p=4\)の人では，テスト情報量が0.5程度しかなく，標準誤差が5.28と非常に大きくなっています。 つまりこの10項目で推定を行っても，\(\theta\)の真値が4の人に対する推定値は平均で\(\pm 5\)くらいには変動するわけです。 そう考えると，この10項目ではどうあがいても\(\theta_p\)が高い層においてまともな推定値は得られないということが分かります。
8.11 （おまけ）多母集団モデルと特異項目機能
多母集団同時分析のところ （ セクション 7.9 ） では，多母集団同時分析という枠組みを紹介しました。 当然ながらIRTにおいても，多母集団同時分析を実行することができます。 そもそもSEMやIRTにおいて多母集団モデルを行う理由が何だったのかを思い出してみると，その一つには「集団ごとに異なる回答傾向の原因を探る」という点がありました。 SEMにおける多母集団モデルでは，例えばある（心理）尺度への回答において，特定のグループ（e.g., 男性・公立高校生・日本人）と別のグループ（e.g., 女性・私立高校生・米国人）で平均点に差があったとしたら，
- まず測定不変性（因子負荷および切片が同じ）を確認した上で
- 因子得点の平均値を集団ごとに自由推定することで差を見る
という手続きを取るのが一般的です。 ここで重要になってくるのが，そもそも回答傾向に違いが生じる原因は2種類あるという点です。 先程の手続きの2つのステップに対応する形で
- そもそも項目・尺度が測っている構成概念の意味がグループ間で異なる
- 構成概念の意味は同じだが，その特性の強さがグループ間で異なる
という2つの可能性が考えられます。 SEMの多母集団モデルでは，多くの場合2番目に強い関心があることに加えて，ひとつひとつの項目よりは尺度全体としての測定不変性に関心があるためか，1番目についてはあくまでも前提条件扱いであまり重要視されていない気がします。
これに対してIRTの場合，観測されるデータは回答者と項目の交互作用として考えられる側面が強いため，「具体的にどの項目がグループ間で異なる挙動をしているか」というような項目側の要因にも関心があることが多いです。 この「異なるグループ間で項目が異なる挙動をする」ことを，一般的には特異項目機能 (Differential Item Functioning [DIF])と呼びます。 そんなわけで，IRTにおける多母集団モデルは主に以下の2種類の用途で用いられることがあります。
- 特性値\(\theta\)の分布はグループ間で共通として，グループごとに推定される項目パラメータ\(\alpha, \beta\)の差異（=DIF）を見る
- 項目パラメータ\(\alpha, \beta\)はグループ間で共通だとして，特性値\(\theta\)の分布のグループごとの差異を見る（SEMで言うところの強測定不変モデルと同じ）
mirtパッケージで多母集団IRTモデルを実行するためには，mirt()関数の代わりにmultipleGroup()関数を使用します。 multipleGroup()関数では，mirt()関数に引数に加えてgroupという引数が登場します。 lavaanパッケージのときにも多母集団モデルでは引数groupを指定しましたが，その時とは指定の仕方が異なるのでご注意ください。
結果を見るための関数は概ねmirt()のときと同じように使えます。 とりあえず推定された項目パラメータを見てみましょう。
$`1`
$items
a b g u
Q1_1 0.447 -2.061 0 1
Q1_2 1.057 -1.751 0 1
Q1_3 0.971 -1.545 0 1
Q1_4 1.169 -1.271 0 1
Q1_5 1.073 -1.385 0 1
Q1_6 1.249 -1.455 0 1
Q1_7 1.281 -1.086 0 1
Q1_8 1.273 -1.111 0 1
Q1_9 1.045 -0.930 0 1
Q1_10 1.156 0.203 0 1
$means
F1
0
$cov
F1
F1 1
$`2`
$items
a b g u
Q1_1 0.294 -4.870 0 1
Q1_2 0.853 -3.029 0 1
Q1_3 0.962 -2.178 0 1
Q1_4 0.858 -2.118 0 1
Q1_5 0.871 -2.184 0 1
Q1_6 1.234 -1.613 0 1
Q1_7 1.449 -1.268 0 1
Q1_8 1.228 -1.369 0 1
Q1_9 1.574 -1.041 0 1
Q1_10 1.322 -0.127 0 1
$means
F1
0
$cov
F1
F1 1結果を見ると，$`1`と$`2`という2つに分かれて，それぞれ項目パラメータの推定値が出ています。 これは，引数groupで指定したグループの値ごとに項目パラメータをそれぞれ推定した結果です。 （2つのグループで\(\theta\)の分布が同じという前提のもとで）もしも推定された項目パラメータが大きく異なる場合には，DIFが発生していると考えられるわけです。 また，各グループにおいて$meansおよび$covはそれぞれ0，1となっています。 multipleGroup()関数は，デフォルトではこのように「各グループの\(\theta\)の母集団分布がそれぞれ標準正規分布である」という設定のもとで項目パラメータを推定します。
特性値\(\theta\)の分布のグループごとの差異を見たい場合には，SEMの多母集団モデルと同じ要領で等値制約を表す引数を与える必要があります。 lavaanでの多母集団モデルを実行したときには，因子負荷の等値制約などを引数group.equalで指定しました。 multipleGroup()では，invarianceという引数が用意されています。 この引数には， 表 8.4 に示す値を入れることが出来ます。 free_meansおよびfree_varsという値があるということは，これらを設定しない限り，全てのグループの\(\theta\)は標準正規分布に固定されてしまうということなのでご注意ください。
invarianceに指定できるもの
| 指定 | 制約されるもの |
|---|---|
free_means |
グループ1以外の\(\theta\)の平均値を自由推定する |
free_vars |
グループ1以外の\(\theta\)の分散を自由推定する |
slopes |
項目の識別力パラメータをグループ間で同じとする |
intercepts |
項目の切片パラメータをグループ間で同じとする |
改めて推定されたパラメータを見てみましょう。
$`1`
$items
a b g u
Q1_1 0.445 -2.584 0 1
Q1_2 1.072 -1.978 0 1
Q1_3 1.067 -1.555 0 1
Q1_4 1.041 -1.414 0 1
Q1_5 1.038 -1.482 0 1
Q1_6 1.133 -1.396 0 1
Q1_7 1.294 -1.003 0 1
Q1_8 1.162 -1.078 0 1
Q1_9 1.272 -0.792 0 1
Q1_10 1.208 0.242 0 1
$means
F1
0
$cov
F1
F1 1
$`2`
$items
a b g u
Q1_1 0.445 -2.584 0 1
Q1_2 1.072 -1.978 0 1
Q1_3 1.067 -1.555 0 1
Q1_4 1.041 -1.414 0 1
Q1_5 1.038 -1.482 0 1
Q1_6 1.133 -1.396 0 1
Q1_7 1.294 -1.003 0 1
Q1_8 1.162 -1.078 0 1
Q1_9 1.272 -0.792 0 1
Q1_10 1.208 0.242 0 1
$means
F1
0.4
$cov
F1
F1 1.013確かに項目パラメータの推定値がグループ間で同じになり，$`2`の$meansおよび$covの値が変わりました。 この結果は（2つのグループでDIFが発生している項目が無いという前提のもとで）グループ2のほうが\(\theta\)の平均値が高いことを示しているわけです。
8.11.1 DIFの評価
ここからは，DIFが発生しているかを評価するいくつかの方法を紹介します。 とりあえず二値IRTモデルについての方法を紹介していきますが，多くの方法はすでに多値モデルに対する拡張も提案されているはずなので，興味があれば探してみてください。
ノンパラメトリックな方法
DIFの検出法は，大きく分けるとIRTモデルに基づく，つまり項目パラメータの推定値によって行われるパラメトリックな方法と，IRTを使用しないノンパラメトリックな方法に分けられます。 ノンパラメトリックな方法として最も有名なのが，\(2\times2\)クロス表に対する\(\chi^2\)検定を拡張したMantel-Haenszel検定 （マンテル・ヘンツェル検定: Mantel & Haenszel, 1959）に基づく方法です。 手法自体は非常に古いものですが，その使いやすさから近年でもDIF検出の方法として非常に多く用いられているようです （Berrío et al., 2020）。
MH検定では，まず各集団を何らかの変数の値によって\(K\)層に層化します。 DIF検出の場合には，IRTで推定された特性値\(\theta\)や合計点などによって層化するのが一般的です。 表 8.5 は，そうして層化されたうちの第\(k\)層について，ある二値項目に対する回答をグループごとに集計したものです。 MH検定では， 表 8.5 のような\(2\times2\)クロス表が\(K\)個あるときに，それらを全部ひっくるめて独立性の検定を行う手法と言えます。
| 項目\(i\) | 当てはまらない | 当てはまる | 計 |
|---|---|---|---|
| グループA | \(x_{A0}^{(k)}\) | \(x_{A1}^{(k)}\) | \(n_A^{(k)}\) |
| グループB | \(x_{B0}^{(k)}\) | \(x_{B1}^{(k)}\) | \(n_B^{(k)}\) |
| 計 | \(n_0^{(k)}\) | \(n_1^{(k)}\) | \(N\) |
まずは一つの層から考えることで，具体的な方法を紹介していきましょう。 例えば\(\theta\)が最も低い人達が集められた第1層において，グループAの人たちのほうが「当てはまる」と回答した割合が高かったとします。 このように，特定の層に絞ったときに回答傾向に違いが見られるというのも立派なDIFです。 そしてDIFが発生している場合には，2つのグループにおける回答の割合つまりオッズ比 \[ \mathrm{OR}^{(k)} = \frac{x_{A1}^{(k)}/n_A^{(k)}}{x_{B1}^{(k)}/n_B^{(k)}} \tag{8.59}\] が1ではない，という見方ができます。 MH検定では，このオッズ比を全ての層について統合したもの \[ \mathrm{OR}^{(\mathrm{MH})} = \frac{\sum_{k=1}^{K}x_{A1}^{(k)}/n_A^{(k)}}{\sum_{k=1}^{K}x_{B1}^{(k)}/n_B^{(k)}} \tag{8.60}\] が1であるかどうかを検定します。
MH検定の結果が有意であれば，データ全体として一方のグループのほうが「当てはまる」と回答した割合が高かったと言えます。 ICC的に言えば 図 8.39 のような状況であると言えるわけです。
一方で，MH検定では検出できないDIFも存在します。 図 8.40 に示されているDIFでは，\(\theta_p\)の低い層ではグループAのほうが「当てはまる」割合が高い一方で，\(\theta_p\)の高い層では逆転が起こっています。 \(\mathrm{OR}^{(k)}\)を一つ一つ見れば1ではないわけですが，これを統合した\(\mathrm{OR}^{(\mathrm{MH})}\)は1になってしまうわけです。 MH検定による方法は， 図 8.39 に示した均一DIFは検出可能な一方で， 図 8.40 のような不均一DIFは検出できないかもしれない方法であることに気をつけましょう。
項目パラメータの推定値を比べる方法
IRTに基づくパラメトリックな方法として，まずは項目パラメータの推定値の差を直接検定する方法を紹介します。 Wald検定は，パラメータの推定値に対して用いられる検定の一種です。 一般的に用いられる一変量Wald検定では，パラメータ\(\xi\)の推定値\(\hat{\xi}\)に対して \[ W = \left(\frac{\hat{\xi}-\xi_0}{\sigma_{\hat{\xi}}}\right)^2 \tag{8.61}\] という検定統計量を考えます。ただし\(\xi_0\)は帰無仮説における値，\(\sigma_{\hat{\xi}}\)は\(\hat{\xi}\)の標準誤差です。 検定統計量\(W\)は自由度1の\(\chi^2\)分布に従うことが知られているため，これを用いて「\(\hat{\xi}\)が\(\xi_0\)である」という帰無仮説を検定します。
Lord （1980） はこのWald検定によってDIFを検出する方法を提案しました。 DIFが発生している場合には，グループごとに推定された項目パラメータに差があると考えられるため，(8.61)式の中身を「2つのグループでの推定値」に置き換えることにより，例えば項目\(i\)の困難度\(\beta_i\)であれば \[ \begin{aligned} \begin{split} W &= \left(\frac{\beta_i^{(A)}-\beta_i^{(B)}}{\sigma_{\beta_i^{(A)}-\beta_i^{(B)}}}\right)^2 \\ &= \left(\frac{\beta_i^{(A)}-\beta_i^{(B)}}{\sigma_{\beta_i^{(A)}} + \sigma_{\beta_i^{(B)}}}\right)^2 \end{split} \end{aligned} \tag{8.62}\] という形で，同じようにWald検定を実行することができます。
Wald検定で用いるWald統計量\(W\)は「推定値の差」を「標準誤差」で割った形をしているため，パラメータの数が複数であっても同じように定義することができます。 例えば項目\(i\)の識別力と困難度を並べたベクトル\(\symbf{\xi}_i = \begin{bmatrix}\alpha_i & \beta_i\end{bmatrix}\)に対して \[ W = \symbf{\xi}_i\symbf{\Sigma}^{-1}\symbf{\xi}_i^\top \tag{8.63}\] という統計量（\(\symbf{\Sigma}\)は標準誤差行列）を考えてやれば，これが自由度2（＝パラメータの数）の\(\chi^2\)分布に従うため，同じようにして検定ができるのです。
項目特性曲線を比べる方法
IRTでは，推定された項目パラメータをもとに項目特性曲線（ICC）を描くことが出来ました。 ICCは，いわば「得点の期待値」を表したグラフです。 したがって「グループ間でICCのずれが大きい」場合にはDIFが発生していると考えることが出来そうです （Raju, 1988）。 ICCのずれとは， 図 8.41 でオレンジ色に塗りつぶされた部分の面積のことです。 これが大きい場合には，確かにDIFが発生していると言えそうです。
グループAにおけるICC（正答確率）は，2パラメータロジスティックモデルならば \[ P_i^{(A)}(\theta) = P(y_{pi}=1|\theta) = \frac{1}{1+\exp[\alpha_i^{(A)}(\theta-\beta_i^{(A)})]} \tag{8.64}\] と表せるため，ICCのズレの面積に基づくDIFの指標として，例えば \[ \mathrm{DIF}_i = \int_{-\infty}^{\infty}\left\vert P_i^{(A)}(\theta) - P_i^{(B)}(\theta)\right\vert f(\theta)d\theta \tag{8.65}\] というものを考えることができるわけです（e.g., Chalmers, 2018; Raju et al., 1995）。ただし\(f(\theta)\)は\(\theta\)の母集団分布（普通は標準正規分布）を指しています。
Rでやってみる
DIFの検出法は相当たくさんあります （レビューとして Berrío et al., 2020; Kleinman & Teresi, 2016 など）。 ということで検出法の紹介はこのあたりでおしまいにして，mirtパッケージに実装されているDIF検出法を試してみましょう26。 これまでに紹介した方法と異なるのですが，DIF()関数を使うと，モデル比較および尤度比検定の観点からDIFの有無を評価してくれます。 DIF()関数のデフォルトの挙動では，第一引数で与えたモデルをベースラインに特定の1項目のみに等値制約をかけた（＝その項目にはDIFが無いと仮定した）ときの変化を計算します。 ここで等値制約がかかるのは引数which.parで指定したパラメータです。 したがって，2パラメータモデルの場合はwhich.par = c("a1", "d")と指定することで「その項目の全ての項目パラメータのDIFを丸ごと評価する」ことを意味します。 （あまり無いケースですが）例えば識別力には興味がなく，困難度（切片）のDIFだけに関心がある場合にはwhich.par = "d"としたら良いわけです。
groups converged AIC SABIC HQ BIC X2 df p
Q1_1 1,2 TRUE -24.578 -19.339 -20.363 -12.985 28.578 2 0
Q1_2 1,2 TRUE -32.846 -27.608 -28.632 -21.253 36.846 2 0
Q1_3 1,2 TRUE -18.162 -12.924 -13.948 -6.569 22.162 2 0
Q1_4 1,2 TRUE -13.108 -7.870 -8.894 -1.515 17.108 2 0
Q1_5 1,2 TRUE -13.712 -8.473 -9.497 -2.119 17.712 2 0
Q1_6 1,2 TRUE 2.217 7.456 6.432 13.810 1.783 2 0.41
Q1_7 1,2 TRUE -4.187 1.052 0.028 7.406 8.187 2 0.017
Q1_8 1,2 TRUE -1.149 4.090 3.066 10.444 5.149 2 0.076
Q1_9 1,2 TRUE -15.366 -10.128 -11.151 -3.773 19.366 2 0
Q1_10 1,2 TRUE -9.558 -4.320 -5.344 2.035 13.558 2 0.001結果を見てみると，SEMの多母集団モデルのときにも登場したAICやBICなどの情報量規準およびX2（\(\chi^2\)）統計量とこれに関連した自由度・\(p\)値が表示されています。 例えばQ1_1行のAICの値-24.578というのは，「Q1_1のみグループごとに項目パラメータが同じ＝DIFが無いという制約をかけることで，AICが24.578悪化した」ということを意味します。 したがって，AICのみで判断するならばQ1_1にはDIFが生じている，と言えるわけです。 同様に，\(\chi^2\)検定の結果も「その項目にはDIFが発生していないとする（グループ間で等値制約を置く）ことでモデル適合度がどの程度悪化するか」という観点に基づく尤度比検定が行われています27。
実際にDIFがありそうと判断された項目についてICCを見るには，引数plotdif = TRUEを与えてください。
groups converged AIC SABIC HQ BIC X2 df p
Q1_1 1,2 TRUE -24.578 -19.339 -20.363 -12.985 28.578 2 0
Q1_2 1,2 TRUE -32.846 -27.608 -28.632 -21.253 36.846 2 0
Q1_3 1,2 TRUE -18.162 -12.924 -13.948 -6.569 22.162 2 0
Q1_4 1,2 TRUE -13.108 -7.870 -8.894 -1.515 17.108 2 0
Q1_5 1,2 TRUE -13.712 -8.473 -9.497 -2.119 17.712 2 0
Q1_6 1,2 TRUE 2.217 7.456 6.432 13.810 1.783 2 0.41
Q1_7 1,2 TRUE -4.187 1.052 0.028 7.406 8.187 2 0.017
Q1_8 1,2 TRUE -1.149 4.090 3.066 10.444 5.149 2 0.076
Q1_9 1,2 TRUE -15.366 -10.128 -11.151 -3.773 19.366 2 0
Q1_10 1,2 TRUE -9.558 -4.320 -5.344 2.035 13.558 2 0.001
確かにICCにはそこそこのずれがありそうです。
「項目1つ1つについてDIFの有無をチェックする」という意味では，全く逆方向からのアプローチを取ることも可能です。 パラメータ推定時にinvariance = c("free_means","free_vars","slopes","intercepts")を与えると，全ての項目パラメータが同じ，すなわち全ての項目にDIFが無いことを仮定した上での推定を行うことが出来ました。 このときの結果（result_group2）をベースラインとして考えると，DIF()関数には特定の1項目のみから等値制約を外した（＝その項目だけDIFがあると仮定した）ときの変化を計算してもらいたいところです。 これを実現するためには，DIF()関数の引数schemeを指定する必要があります。
groups converged AIC SABIC HQ BIC X2 df p
Q1_1 1,2 TRUE -10.691 -5.452 -6.476 0.902 14.691 2 0.001
Q1_2 1,2 TRUE -8.581 -3.343 -4.367 3.012 12.581 2 0.002
Q1_3 1,2 TRUE -0.172 5.066 4.042 11.420 4.172 2 0.124
Q1_4 1,2 TRUE 1.756 6.995 5.971 13.349 2.244 2 0.326
Q1_5 1,2 TRUE 2.772 8.010 6.986 14.365 1.228 2 0.541
Q1_6 1,2 TRUE -2.267 2.972 1.948 9.326 6.267 2 0.044
Q1_7 1,2 TRUE 1.468 6.706 5.683 13.061 2.532 2 0.282
Q1_8 1,2 TRUE 0.590 5.829 4.805 12.183 3.41 2 0.182
Q1_9 1,2 TRUE -1.618 3.620 2.596 9.975 5.618 2 0.06
Q1_10 1,2 TRUE 3.055 8.293 7.270 14.648 0.945 2 0.623結果の見方は先程と同じです。情報量規準の値がマイナスになっている・検定が有意になっている（\(p < .05\)）項目は等値制約を外したほうが適合度が良くなっている＝DIFの疑いあり，と判断することができます。
続いてはLordのWald検定の結果を出してみます。 ただ，Wald検定を行うためには「推定値の標準誤差行列」が必要になるため，まずは標準誤差を推定するように引数SE = TRUEを付け加えてやり直す必要があります。
推定された標準誤差行列は一応result_group@vcovで見ることができますが，見たところでよくわからないと思うのでそっとしておいてあげてください。
無事標準誤差が推定できたので，これを用いてWald検定に望みます。 Wald検定は，DIF()関数で引数Wald = TRUEを与えることで実行可能です。
結果の見方はこれまでと基本的に同じです。 つまり検定結果が有意になった（\(p < .05\)）項目はグループ間でパラメータの推定値が有意に異なる＝DIFの疑いあり，ということになります。
最後はICCのズレに基づく評価方法です。 こちらは Chalmers （2018） による方法がDRF()という関数に実装されています。 DRF()関数は，デフォルトでは項目特性曲線（ICC）の代わりにテスト特性曲線（TCC）の差について評価を行います28。 項目ごとのDIFを見たい場合には，引数DIF = TRUEを与える必要があるのでご注意ください。 また引数plot = TRUEを与えると，「ICCの差のプロット」を出してくれるようになります。
結果を見ると，3種類のDIFが表示されています。
sDIF- 符合つき（signed）のズレの面積（単純に差を積分したもの）を表しています。 図 8.41 のようにズレが上側と下側にある場合，これらが打ち消し合ってsDIFの値は小さくなります。つまりsDIFは「総合的に見てどちらのグループのほうが平均点が高いか」を表していると言えるでしょう。
uDIF- 符号なし（unsigned）のズレの面積（差の絶対値を積分したもの: 8.65 式）を表しています。したがって，sDIFよりも大きな値が表示されているはずです。
dDIF- (8.65)式とは別のやり方で符号なしのズレを求めるために，絶対値の代わりに差を二乗することを考えたDIFです。二乗した分を戻すため，最後に全体にルートをかけます。「二乗平均のルートをとる」という考え方は，いわゆる逸脱度（deviance）的な考え方なのでdDIFと名付けられているのだと思います。多分。uDIFと比べると，ズレの面積が同じだとしても，大きくズレている箇所があるほど，二乗されている分だけdDIFのほうが大きな値になると思われます。
いずれにしても，この(s|u|d)DIFの値が大きい項目はDIFの疑いあり，と言えるわけです。 ICCのズレに基づく検定を行うためには「ICCのズレの標準誤差」のようなものが必要になりますが，これはパラメータの標準誤差だけでは対応できないようです。 そのため，DRF()関数では代わりにブートストラップ標準誤差を用いて検定を行います。 これを行うためには，引数drawsに適当な数を指定してください。
$sDIF
groups item sDIF CI_2.5 CI_97.5 X2 df p
1 1,2 Q1_1 0.096 0.060 0.132 27.476 1 0
2 1,2 Q1_2 0.086 0.054 0.118 29.952 1 0
3 1,2 Q1_3 0.077 0.044 0.117 20.188 1 0
4 1,2 Q1_4 0.066 0.031 0.100 14.482 1 0
5 1,2 Q1_5 0.069 0.037 0.102 17.731 1 0
6 1,2 Q1_6 0.022 -0.009 0.056 1.789 1 0.181
7 1,2 Q1_7 0.050 0.015 0.086 7.794 1 0.005
8 1,2 Q1_8 0.041 0.006 0.080 4.808 1 0.028
9 1,2 Q1_9 0.070 0.034 0.108 14.702 1 0
10 1,2 Q1_10 0.078 0.032 0.123 12.502 1 0
$uDIF
groups item uDIF CI_2.5 CI_97.5 X2 df p
1 1,2 Q1_1 0.096 0.061 0.132 28.194 2 0
2 1,2 Q1_2 0.086 0.055 0.118 27.086 2 0
3 1,2 Q1_3 0.077 0.045 0.117 16.579 1 0
4 1,2 Q1_4 0.067 0.038 0.101 14.895 2 0.001
5 1,2 Q1_5 0.069 0.038 0.103 15.664 2 0
6 1,2 Q1_6 0.022 0.005 0.058 2.18 1 0.14
7 1,2 Q1_7 0.050 0.019 0.087 7.796 2 0.02
8 1,2 Q1_8 0.041 0.015 0.081 5.256 2 0.072
9 1,2 Q1_9 0.079 0.049 0.125 18.671 2 0
10 1,2 Q1_10 0.078 0.040 0.124 13.037 2 0.001
$dDIF
groups item dDIF CI_2.5 CI_97.5
1 1,2 Q1_1 0.105 0.069 0.148
2 1,2 Q1_2 0.112 0.073 0.156
3 1,2 Q1_3 0.087 0.053 0.137
4 1,2 Q1_4 0.094 0.043 0.142
5 1,2 Q1_5 0.089 0.049 0.137
6 1,2 Q1_6 0.027 0.006 0.076
7 1,2 Q1_7 0.054 0.022 0.097
8 1,2 Q1_8 0.048 0.018 0.101
9 1,2 Q1_9 0.084 0.053 0.135
10 1,2 Q1_10 0.084 0.045 0.137
結果の見方はやはり同じです。 つまり検定結果が有意になった（\(p < .05\)）項目はIRFがグループ間でズレている＝DIFの疑いあり，ということになります。
DIF検出の実際
ここまでDIFを検出するいくつかの方法を紹介しました。 色々な観点からの評価がある，という意味では，やはり複数の指標を組み合わせて総合的に判断するのが良いのだろうと思います。
それはそれとして，DIFの有無を評価しようとする場合，そもそも関心のあるグループの間で特性値\(\theta\)の分布に差が無いということはあまり考えられません（無かったら非常にラクなのですが…）。 したがって，そもそも「\(\theta\)の分布が同じだとしてDIFを検出する」というスキームは現実的に結構厳しいのです。 この問題を克服するためには，「グループ間で挙動が変わらない＝DIFが無いことが事前にわかっている項目」であるアンカー項目 (anchor items) を加えておくことが考えられます。 アンカー項目の項目パラメータを基準とすることで，\(\theta\)の分布がグループで異なる場合であっても特定の項目にDIFが生じているかを評価することが出来るわけです。 具体的にアンカー項目をどのように決定したら良いかはまた難しい話（レビューとして Kopf et al., 2015）なので，ここではとりあえずアンカー項目がある場合のやり方だけ紹介しておきます。
multipleGroup()の引数invarianceには，実は項目の名前（データの列名）を与えることも出来ます。
あとの流れはほぼ同じなので，紹介はここまでとします。
$`1`
$items
a b g u
Q1_1 0.444 -2.360 0 1
Q1_2 1.013 -1.831 0 1
Q1_3 0.962 -1.456 0 1
Q1_4 1.165 -1.274 0 1
Q1_5 1.061 -1.396 0 1
Q1_6 1.255 -1.450 0 1
Q1_7 1.291 -1.082 0 1
Q1_8 1.279 -1.108 0 1
Q1_9 1.056 -0.924 0 1
Q1_10 1.167 0.202 0 1
$means
F1
0
$cov
F1
F1 1
$`2`
$items
a b g u
Q1_1 0.444 -2.360 0 1
Q1_2 1.013 -1.831 0 1
Q1_3 0.962 -1.456 0 1
Q1_4 0.949 -1.177 0 1
Q1_5 0.957 -1.247 0 1
Q1_6 1.344 -0.734 0 1
Q1_7 1.578 -0.418 0 1
Q1_8 1.353 -0.502 0 1
Q1_9 1.735 -0.204 0 1
Q1_10 1.441 0.625 0 1
$means
F1
0.741
$cov
F1
F1 0.828 groups converged AIC SABIC HQ BIC X2 df p
Q1_1 1,2 TRUE 0.000 0.000 0.000 0.000 0 0 NaN
Q1_2 1,2 TRUE 0.000 0.000 0.000 0.000 0 0 NaN
Q1_3 1,2 TRUE 0.000 0.000 0.000 0.000 0 0 NaN
Q1_4 1,2 TRUE -1.404 3.835 2.811 10.189 5.404 2 0.067
Q1_5 1,2 TRUE 0.570 5.809 4.785 12.163 3.43 2 0.18
Q1_6 1,2 TRUE -18.972 -13.733 -14.757 -7.379 22.972 2 0
Q1_7 1,2 TRUE -15.059 -9.821 -10.845 -3.467 19.059 2 0
Q1_8 1,2 TRUE -14.729 -9.491 -10.515 -3.136 18.729 2 0
Q1_9 1,2 TRUE -14.361 -9.122 -10.146 -2.768 18.361 2 0
Q1_10 1,2 TRUE -10.458 -5.220 -6.243 1.135 14.458 2 0.001結果を見て分かるように，アンカー項目を適切に設定しないと，その他の項目のDIFの評価も適切に行うことができません。 そういった意味では，DIFの検出というのは結構職人技だったりするのです。
8.12 いろいろなIRTモデル
IRTの本質は，(8.1)式に集約されていると個人的には思っています。 それは項目への反応（データ）は回答者と項目の交互作用によって確率的に決定するという点です。 そのような挙動をとるモデルの形の一つとして，因子分析と同じようなモデル（正規累積・ロジスティック）や多値反応に対するモデル（GRM・GPCMなど）を紹介してきました。 しかし世の中にはもっと様々なIRTモデルが存在しています。 全てを紹介するのは不可能ですが，ここではいくつかのモデルのさわりだけを（個人的な関心に基づいて）紹介していきたいと思います29。
8.12.1 多値反応に対するその他のIRTモデル
多値反応に対するIRTモデルとして， セクション 8.6 ではGRMやGPCMといったモデルを紹介してきました。 多値反応に対するIRTモデルには他にも様々なものが提案されており，近年注目を集めている（気がする）のがIRTree （Böckenholt, 2017） と呼ばれるモデル群です。 これらのモデルでは，回答決定のプロセスとして名前通り「ツリー構造」を仮定します。 例えば「とても嫌い(1)」「嫌い(2)」「どちらでもない(3)」「好き(4)」「とても好き(5)」という5件法の項目があったとしましょう。 この項目への回答決定プロセスの考え方の一例としては， 図 8.45 のように
- まず「どちらでもない(3)」を選択（態度保留）するかを二択で考える
- 態度を表明する場合には「好き(4-5)」か「嫌い(1-2)」かを二択で考える
- 「とても(1,5)」の選択肢を選ぶほど好き or 嫌いかを二択で考える
という感じで３つの二択を順番に行う，というものがあるでしょう。
するとこの多値項目の背後には，実は３つの二値項目（query）が存在していると考えることができ，それぞれが
- 「どちらでもない」を選ぶ傾向（中心傾向）：\(\theta_p^{(M)}\)が強い人かどうか
- 実際にその対象物が好きかどうか：\(\theta_p^{(A)}\)
- 極端な選択肢を選ぶ傾向（極端傾向）：\(\theta_p^{(E)}\)が強い人かどうか
という異なる個人特性によって動いていると仮定することができそうです（もちろんプロセス単位で完全に独立しているとは思いませんが，ある程度は）30。 各queryに対する反応確率は，普通の正規累積・ロジスティックモデルなどで表現できます。 例えば中心傾向についてのqueryだけで表現できる「どちらでもない」を選ぶ確率は \[ \begin{aligned} \begin{split} P(y_{pi}=3|\theta_p^{(M)}) &= \int_{-\infty}^{\alpha_{i}^{(M)}\left(\theta_p^{(M)} - \beta_i^{(M)}\right)}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{2}\right)dz \\ &= \symbf{\Phi}\left[\alpha_{i}^{(M)}\left(\theta_p^{(M)} - \beta_i^{(M)}\right)\right] \end{split} \end{aligned} \tag{8.66}\] などとなるわけです。IRTreeモデルでは，項目パラメータにも上付き添字がついています。 したがって，項目ごとに「どちらでもない」および極端な選択肢の選ばれやすさが異なる，という設定がされています。
同様に，「好き(4)」を選ぶ確率は，「『どちらでもない』を選ばない」かつ「好き」かつ「『とても』ではない」の積なので \[ \begin{aligned} \begin{split} &P(y_{pi}=4|\theta_p^{(M)},\theta_p^{(A)}\theta_p^{(E)}) = \\ &\left(1-\symbf{\Phi}\left[\alpha_{i}^{(M)}\left(\theta_p^{(M)} - \beta_i^{(M)}\right)\right]\right)\symbf{\Phi}\left[\alpha_{i}^{(A)}\left(\theta_p^{(A)} - \beta_i^{(A)}\right)\right]\left(1-\symbf{\Phi}\left[\alpha_{i}^{(E)}\left(\theta_p^{(E)} - \beta_i^{(E)}\right)\right]\right) \end{split} \end{aligned} \tag{8.67}\] という形で表現できるわけです。
IRTreeの目的あるいは利点は，「測定したい回答者の特性値」と「回答傾向」を分離して考えられる点です。 一般的に，リッカート尺度への回答には個人のクセのようなものが混ざっていることが知られています。 そのせいで，「どちらでもない(3)」を選んだ人が「好きと嫌いのちょうど中間」なのか「本当は『どちらかといえば好き』くらいなんだけど，その程度では『好き』は選びづらい」のかは区別がつきにくいものです。 IRTreeモデルならば，上の例で言えば中心傾向と極端傾向の強さを，個人特性\(\theta_p\)の一種として考えることができます。 そのため，残った\(\theta_p^{(A)}\)はより純粋な「測定したい回答者の特性値」として解釈できるかもしれないというわけです。 当然ながら，IRTreeモデルでは 図 8.45 に示したような回答決定プロセスの仮定がかなり重要になってきます。 逆に言えば，複数のツリー構造間でモデル比較をすることで，どのプロセスが最も当てはまりそうか，ということを評価したりも出来るかもしれません…。
8.12.2 強制選択形式に対するIRTモデル
IRTreeモデルはいわば，得られた回答データから回答傾向を統計的に分離しようと頑張っているモデルの一種です。 そもそもリッカート尺度では回答傾向の影響を逃れることができないわけですが，その対抗策としては他にも「そもそも中間選択肢をなくそう」という方向を考えることができます。 また，入社試験などの重要な局面では「ウソをついて本当の自分より良い人に見せよう」という思惑が働いて回答が歪められる可能性があったりします （フェイキング: e.g., Jackson et al., 2000）。 これを防ぐためには，例えば「友だちが多い」「真面目である」といった複数の短文を同時に提示し，その中から一つだけを選択させたり，ランキングを付けさせるという回答方法が考えられます。 これならば，前述の「中間選択肢を選びがち」といった問題も防げますし，フェイキングの問題もかなり抑制できるでしょう。
複数の選択肢から一つの選択肢を選ばせる回答方法は，強制選択 (Forced Choice) 型項目などと呼ばれます。 また，特に選択肢数が２個の場合は一対比較 (Paired Comparison) と呼ばれることもあります。 そしてこれらの出題形式に対しては，やはり普通のIRTとは少し異なるモデルが考えられています。 最も代表的なモデルが，Thurstonian IRT (TIRT) モデル （Brown & Maydeu-Olivares, 2011）でしょう。 TIRTモデルでは，各選択肢がそれぞれ独立に「効用」を持っていると考えます。 この「効用」とは，経済学などで出てくるものと同じようなものだと考えてOKです。 あるいは単純に「各選択肢が自分に当てはまっていると感じる程度」というイメージでも構いません。 とにかく強制選択型項目への回答では，回答者は効用を最大化する選択肢を選ぶのだ，と考えているわけです。 ここでは，回答者\(p\)が項目\(i\)の２つの選択肢\(A, B\)に対してそれぞれ効用\(u_{pi}^{(A)},u_{pi}^{(B)}\)を持っているとしましょう。 すると，回答者\(p\)は効用を最大化する選択肢を選ぶため，観測されるデータとしては \[\begin{aligned} y_{pi} = \begin{cases} A & \mathrm{if}~u_{pi}^{(A)} - u_{pi}^{(B)} \geq 0 \\ B & \mathrm{if}~u_{pi}^{(A)} - u_{pi}^{(B)} < 0 \end{cases} \end{aligned} \tag{8.68}\] と表わせます。
TIRTモデルでは，選択肢の効用\(u_{pi}^{(A)}\)を，よくある因子分析モデルの形 \[ u_{pi}^{(A)} = \tau_{i}^{(A)} + f_{p}^{(A)}b_i^{(A)} + \varepsilon_{pi}^{(A)}, ~~ \varepsilon_{pi}^{(A)} \sim N(0, \sigma_{\varepsilon_{pi}^{(A)}}) \tag{8.69}\] と定義します。つまり，平均的な効用\(\tau_{i}^{(A)}\)と，その選択肢が反映している因子の得点\(f_{i}^{(A)}\)と因子負荷\(b_p^{(A)}\)の和によって表される効用の期待値に，そこに正規分布に従う誤差（独自因子得点）\(\varepsilon_{pi}^{(A)}\)が加わって効用の実現値が確率的に決定すると考えているわけです。 すると，２つの選択肢\(A, B\)の中から\(A\)を選ぶ確率は最終的に \[ \begin{aligned} \begin{split} P(y_{pi}=A) &= P\left(u_{pi}^{(A)} - u_{pi}^{(B)} \geq 0\right) \\ &= \symbf{\Phi}\left[\alpha_{i}^{(AB)} + \left(f_{i}^{(A)}\beta_p^{(A)} - f_{i}^{(B)}\beta_p^{(B)}\right)\right] \end{split} \end{aligned} \tag{8.70}\] と表されます31。ただし表記を簡略化するため \[ \begin{aligned} \alpha_{i}^{(AB)} &= \frac{\tau_{i}^{(A)} - \tau_{i}^{(A)}}{\sqrt{\sigma_{\varepsilon_{pi}^{(A)}} + \sigma_{\varepsilon_{pi}^{(B)}}}}, \\ \beta_i^{(A)} &= \frac{b_i^{(A)}}{\sqrt{\sigma_{\varepsilon_{pi}^{(A)}} + \sigma_{\varepsilon_{pi}^{(B)}}}}, \\ \beta_i^{(B)} &= \frac{b_i^{(B)}}{\sqrt{\sigma_{\varepsilon_{pi}^{(B)}} + \sigma_{\varepsilon_{pi}^{(B)}}}} \end{aligned} \tag{8.71}\] という形で別の記号を用いて表現しています。 つまりTIRTモデルは，一対比較の場合には，各項目が２つの因子のみに負荷量を持っていると設定したCFAと実質的には同じことをしています。 そして一つの項目に３つ以上の選択肢がある場合は，対ごとに異なる項目を構成しているとみなして分析を行います。 例えば項目\(i\)が\(A, B, C, D\)の４択だとすると，この項目は「AとB」「AとC」「AとD」「BとC」「BとD」「CとD」という６つの二値項目があると考えます。 そして回答者\(p\)が選択肢\(A\)を選んだとすると，「AとB」「AとC」「AとD」については（いずれもAを選んだという）結果が観測され，残りの「BとC」「BとD」「CとD」については欠測扱いにします。 そして観測された３項目についてはそれぞれ \[ \begin{aligned} \begin{split} P(y_{pi}=A|A,B) &= \symbf{\Phi}\left[\alpha_{i}^{(AB)} + \left(f_{i}^{(A)}\beta_p^{(A)} - f_{i}^{(B)}\beta_p^{(B)}\right)\right]\\ P(y_{pi}=A|A,C) &= \symbf{\Phi}\left[\alpha_{i}^{(AC)} + \left(f_{i}^{(A)}\beta_p^{(A)} - f_{i}^{(C)}\beta_p^{(C)}\right)\right]\\ P(y_{pi}=A|A,D) &= \symbf{\Phi}\left[\alpha_{i}^{(AD)} + \left(f_{i}^{(A)}\beta_p^{(A)} - f_{i}^{(D)}\beta_p^{(D)}\right)\right] \end{split} \end{aligned} \tag{8.72}\] という形でそれぞれ選択肢Aに関する部分が共通しているという設定のもとで定式化されるわけです。
8.12.3 回答時間を利用するIRTモデル
コンピュータによる測定 (Computer Based Testing [CBT]) では，項目反応以外の情報を得ることが容易になりました。 もっとも代表的な情報として，回答時間 (Response Time [RT]) が挙げられるでしょう。 セクション 3.2.2 では，回答時間を「適当な回答をしている人を検出する」といった目的で使用しましたが，その一方で回答者の特性値推定に利用するという方向性も相当に研究されています。 そういうわけで，回答時間を利用するIRTモデルを挙げていくとキリがないのですが，ここではシンプルなモデルとして，「当てはまる」「当てはまらない」の二択の項目に対するモデル（Ferrando & Lorenzo-Seva, 2007）を紹介します。
回答時間を利用するIRTモデルでは，そもそも回答者の特性値と回答時間にはどのような関係があるかを考えるところから始まります。 性格特性を尋ねる項目の場合には，「特性の強さが中途半端な人ほど回答時間が長い」という逆U字の関係になることが経験的に知られています。 例えば「友だちが多い」という項目に対して，友だちが多い／少ない自信がある人，つまり「外向性」といった特性が極端な人は，きっとすぐに選択肢を選ぶことでしょう。 これに対して，『友だちはそこそこ居るけど』という人は，きっと「当てはまる」「当てはまらない」のどちらを選ぶかに迷いが生じるはずです。 したがって，同じ「当てはまる」を選んだ人同士であっても，素早く「当てはまる」を選んだ人のほうがより強い特性を持っていると考えるのが自然です。
ということで，この関係を尤度関数に組み込んでIRTモデルを構築していきます。 と言っても考え方は簡単で，回答時間\(t_{pi}\)についても項目反応\(y_{pi}\)の(8.1)式と同じように \[ t_{pi} = g(\omega_p, \gamma_i) \tag{8.73}\] という形で，回答者の要因\(\omega_p\)と項目の要因\(\gamma_i\)の交互作用として考えます。 Ferrando & Lorenzo-Seva （2007） では，これを具体的に \[ \begin{aligned} \begin{split} \ln(t_{pi}) &= \mu + \omega_p + \gamma_i + b\delta_{pi} + \varepsilon_{pi}, ~~ \varepsilon_{pi} \sim N(0, \sigma_\varepsilon^2) \\ \delta_{pi} &= \sqrt{\alpha_i^2(\theta_p - \beta_i)^2} \end{split} \end{aligned} \tag{8.74}\] と置きました。 ここで\(\mu\)は全体的な切片項，\(\omega_p\)は回答者\(p\)の回答スピード，\(\gamma_i\)は項目\(i\)の所要時間を表現したパラメータです。 また，\(\delta_{pi}\)の中身には\((\theta_p - \beta_i)^2\)というものが入っています。この項は，\(\theta_p = \beta_i\)のとき，つまり反応確率が五分五分のときに最小値を取ります。 係数\(b\)が負の値であれば，回答者の特性値\(\theta_p\)が項目困難度\(\beta_i\)に近いほど回答時間の期待値が長くなることで，逆U字の関係を表現しているわけです。
回答時間などの追加情報を用いて特性値を推定するIRTモデルでは，うまく行けばより精度の高い，あるいは妥当性の高い推定を行うことが出来ると期待されています。 ちなみに，IRTに限らず社会調査などの分野では，このように回答内容以外の情報（位置情報や回答デバイスの情報などを含む）をパラデータ （2013） と呼んでいます。 パラデータは，やはり適当な回答の検出に使われるだけでなく，よりノイズの少ないデータを収集するための調査改善など様々な用途で利用されていたりします。
8.12.4 展開型IRTモデル
これまでに紹介したIRTモデルはすべて，\(\theta\)の値が高ければ高いほど高い点数をとる（i.e., 「当てはまる」寄りの回答をする，正解する）という仮定がされていました。 項目特性曲線にしても，（\(\alpha_i>0\)ならば）右上がり（単調増加）でした。 ですが質問のタイプによっては「過ぎたるはなお及ばざるが如し」的なものが存在しています。 よくあるのは「態度」に関する項目です。 例えば「消費税を5%にするべきだ」という項目は，「消費税は完全に撤廃してほしい」人も「消費税はこのままで良い」人も「そう思わない」を選択するでしょう。 同様に，心理特性に関しても「カフェで友人と静かに会話をするのは楽しい（因子：外向性）」（Stark et al., 2006）という項目に対しては「そもそも他人と話したくない人」も「もっとワイワイ騒ぎたい人」も「当てはまらない」を選ぶと思われます。 つまりこうした項目における項目特性曲線は， 図 8.46 のように，あるところを境に単調減少に転じると考えられるわけです。
ロジスティック・正規累積モデルでは，「項目(\(\beta_i\))」と「回答者(\(\theta_p\))」は共通の軸の上でマッピングされていると考えることが出来ます。 すごく簡単にいえば「困難度が\(\beta_i\)の項目は\(\theta_p=\beta_i\)の回答者にとって（易しすぎず難しすぎず）ちょうどよい難易度（50%）の項目である」ということです。 これを先程の項目に置き換えて考えると，「\(\beta_i\)と\(\theta_p\)の距離が近いほど，その人によってちょうどよい程度の（理想の）意見・態度を表している」ということができそうです。 このような「理想点(ideal point)」の考え方に基づくIRTモデルを展開型 (unfolding) IRTモデルと呼びます32。
展開型IRTモデルでは，反応確率\(P(y_{pi}=1)\)が「\(\beta_i\)と\(\theta_p\)の距離」によって規定されます。 したがって，最も基本的な一次元・二値モデルの場合は \[ P(y_{pi}=1) = \exp\left[-\frac{1}{2}(\beta_i-\alpha_i\theta_p)^2\right] \tag{8.75}\] と表されます（Maydeu-olivares et al., 2006）。 ここでも2つの項目パラメータの意味は概ね通常の場合と同じです。 \(\beta_i\)はその項目の「位置」を表しており，\(\theta_p\)がこの値に近い人ほど「当てはまる」と回答する確率が高くなることがわかります。 \(\alpha_i\)は項目識別力を表していますが，展開型モデルの場合は，正規分布の分散パラメータのように山の狭さを規定します（ 図 8.47 ）。 ですが結局「\(\theta_p\)が理想点\(\beta_i\)付近か否かを区別する強さ」という意味では，これまでに出てきた項目識別力と同じような役割を果たしているわけです。
mirtでは，itemtype="ideal"と指定することで展開型モデルを実行可能です。
通常のIRTモデルと同じようにcoef()でパラメータの推定値を見ることもできるのですが，どうやらIRTpars = TRUEには対応していないようなのでご注意ください。
また，ICCのプロットも同じように実行可能です。
例えば1番目の（一番山が緩い）項目は，\((\alpha_i,\beta_i)=(0.118,-0.716)\)ということで，\(\theta=0.716/0.118\simeq 6.068\)のあたりで最大値\(P(y_{pi}=1)=1\)となります。そのため，\(\theta\)が6までしか表示されていない 図 8.48 では右上がりのように見えています。
多値の展開型IRTモデルの代表的なものは，一般化段階展開型モデル （Generalized graded unfolding model［GGUM］: Roberts, 2019）です。 このモデルは，GPCMをベースに展開型IRTモデルを構成しています。
例えば4カテゴリの項目がある場合には，「全く当てはまらない（下すぎて）」「あまり当てはまらない（下すぎて）」～「あまり当てはまらない（上すぎて）」「全く当てはまらない（上すぎて）」という要領で8個のカテゴリがあると考えます。 こうすることで，隣同士のカテゴリの比較は今まで通りできるのでGPCM的な方法が適用できる，というわけです。 そして「全く当てはまらない（下すぎて）」と「全く当てはまらない（上すぎて）」の確率を合わせたものを「全く当てはまらない」の選択確率として表します。 このままだとパラメータ数が多すぎる（4カテゴリに対してカテゴリ困難度が7個必要になってしまう）ので，理想点を軸にして各カテゴリ閾値間の距離は左右で等しいという仮定をおいた上で定式化しています。
なおmirtでは，itemtype = "ggum"と指定することで実行可能です。
8.13 IRTの使い道
最後に，IRTの考え方によって出来るようになることをいくつか紹介します33。
8.13.1 等化
セクション 8.1 では，IRTを項目依存性と集団依存性をクリアするための理論だと紹介しました。 ですが実際には，IRTモデル（2パラメータのGRM）はカテゴリカル因子分析と同じなので\(\theta_p\)はデータをとった集団ごとに標準化されています。 つまり最初に例示したような，異なる集団に対して実施した異なる尺度の比較はこのままではできません。 もちろん反対に，異なる尺度を用いて実施した異なる集団同士の比較もこのままではできません。 こうした比較を可能にするために，IRTでは等化 (equating) またはリンキング (linking) という手続き （Kolen & Brennan, 2014）を行います34。 等化の技術は，TOEICのような大規模テストでも実際に利用されています。 そのおかげで，毎回異なるはずの難易度によらず「TOEIC ***点」というように点数を履歴書に書けるわけです。
（正直言うと，等化＝異なる集団や尺度の比較を考えない限りは，特に心理尺度においてわざわざ因子分析やSEMではなくIRTを選ぶ理由はほぼ無い気すらします）
等化の基本的な考え方を説明しましょう。 ここからはIRTで推定されるパラメータが集団によらない項目パラメータと項目によらない個人特性だという仮定のもとで話を進めます。
共通回答者による等化
まず，個人特性値が「項目によらない」のであれば，どんな尺度で推定した場合でも\(\theta_p\)の値は同じになるはずです。 いま，\(p\)さんが同じ構成概念を測定する尺度\(1\)と尺度\(2\)による2回の調査の両方に回答したとします。 この2回のデータをもとにパラメータを推定する場合には，いままでの手順にのっとってそれぞれ\(\theta_p\sim N(0,1)\)と仮定して推定を行います。 その結果
- 尺度1での推定値\(\hat{\theta}_p^{(1)}=0.5\)
- 尺度2での推定値\(\hat{\theta}_p^{(2)}=-0.5\)
だったとしましょう。 仮に尺度1による推定が「項目によらない個人特性」の値とすると，尺度2での推定値は本来の値より1小さい値になっているようです。 \(p\)さんの特性値が尺度1回答時と尺度2回答時で変わっていないとしたら，この1のズレは「尺度2は尺度1よりも特性値を小さめに推定する尺度である」ことが原因と考えられます。 そこで，\(\hat{\theta}_p^{(2)}=-0.5\)がどうにかして\(\hat{\theta}_p^{(1)}\)と同じ\(0.5\)になるように調整する（ゲタを履かせる）必要がありそうです。
もう少し規模を広げて，同じようにこの尺度1と尺度2による2回の調査の両方に回答した人が100人いたとします。 その結果，2回の推定値の散布図が 図 8.49 のようになったとします。 そしてこの2回の推定値の回帰直線を求めたら \[ \hat{\theta}_p^{(2)} = A\hat{\theta}_p^{(1)} + B \tag{8.76}\] であったとします。 この場合，全員の\(\hat{\theta}_p^{(2)}\)をなるべく\(\hat{\theta}_p^{(1)}\)に近づくように調整した値\(\hat{\theta}_p^{(2-1)}\)は， \[ \hat{\theta}_p^{(2-1)} = \frac{\hat{\theta}_p^{(2)}-B}{A} \tag{8.77}\] という形で求められそうだとわかります。
2パラメータロジスティックモデルでは，項目反応関数は \[ \begin{aligned} P(y_{pi}=1)&= \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_i)\right]} \end{aligned} \tag{8.78}\] でした。ここで， \[ \begin{aligned} \theta_p^{*} &= \frac{\theta_p - B}{A} \\ \beta_i^{*} &= \frac{\beta_i - B}{A} \\ \alpha_i^{*} &= A\alpha_i \end{aligned} \tag{8.79}\] とおくと， \[ \begin{aligned} P(y_{pi}=1)&= \frac{1}{1+\exp\left[-\alpha_{i}^{*}(\theta_p^{*} - \beta_i^{*})\right]} \\ &= \frac{1}{1+\exp\left[-A\alpha_{i}\left(\left[\frac{\theta_p - B}{A}\right] - \left[\frac{\beta_i - B}{A}\right]\right)\right]} \\ &= \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_i)\right]} \end{aligned} \tag{8.80}\] となり，もとの項目反応関数と同じになります。 探索的因子分析では，潜在変数のスケールの不定性の問題から，全ての観測変数を標準化していました。 等化では，このスケールの不定性をある種逆手に取ることで\(\alpha, \beta\)や\(\theta\)を適当なスケールに変換しているわけです。 具体的には，2つの尺度に両方とも回答した全員の\(\hat{\theta}_p^{(2)}\)がなるべく\(\hat{\theta}_p^{(1)}\)に近づくようにスケールを変換し，これに応じて項目パラメータにも変換をかけます。 こうして得られた変換後の項目パラメータならば，尺度Aと尺度Bでも問題なく比較ができるようになります。
共通項目による等化
先程説明した等化の流れは，共通の項目がある場合でも同じように考えられます。 つまり，2つの集団に対して実施する際に，それぞれの尺度に共通の項目をいくつか入れておけば， \[ \begin{aligned} \beta_i^{(2-1)} &= \frac{\beta_i^{(2)} - B}{A} \\ \alpha_i^{(2-1)} &= A\alpha_i^{(2)} \end{aligned} \tag{8.81}\] となるような等化係数\(A,B\)を求めることができ，それに応じた変換が可能となります。 等化係数\(A,B\)を求める方法としては，推定値の要約統計量を用いる方法(mean-mean法, mean-sigma法) や，「変換後のICCがなるべく重なるようにする」方法 （Stocking-Lord法: Stocking & Lord, 1983; Haebara法: Haebara, 1980）が一般的です。
等化はIRTの可能性を大きく広げてくれる重要な技術ですが，その適用には様々なハードルがあります。 何より2つの尺度が測定している\(\theta\)が等質であるという大前提が重要です。 もしも尺度1と尺度2がそもそも異なる構成概念を測定しているならば，パラメータを調整したところで同じモノサシの上で比較することはナンセンスです。 本来は全部まとめてデータを取った上で一次元性の検証ができれば良いのですが，等化をしたい場合にはそんなことができないので，理論的側面などから測定内容の等質性を担保してあげる必要があります。 また，等化係数はパラメータの推定値をもとに計算されます。 つまり，等化の係数には少なからず誤差がある，ということです。 厄介なのは，尺度2での項目パラメータの推定値それ自体にも誤差があるという点です。 つまり等化後のパラメータには「等化前の誤差＋等化係数の誤差」という2つの誤差が重なっているので，等化をしない場合と比べるとどうしても推定値は不安定になりやすいです。 理論的には，等化は繰り返し行うことも可能です（尺度4→尺度3→尺度2→尺度1）が，実際に等化を繰り返すとなかなか厳しいことになるので注意が必要です。
8.13.2 テストの設計
テストの一番の目的は，回答者の特性値\(\theta\)を精度良く推定することです。 ですがテストの目的によって「どの範囲の\(\theta\)をどの程度の精度で推定できたら良いか」は変わってきます。 例えばふつうの心理尺度では，比較的広い範囲の特性値をそれなりの精度で推定したいと思うでしょう。 この場合には，困難度が高い項目から低い項目までをまんべんなく用意するのが良さそうです。 一方で，心療内科が扱うような「診断」を目的とした尺度や合格／不合格を決めるテストなどの場合，ある特性の値が具体的にどの程度かよりも「基準値以上か以下か」に強い関心があったりします。 このような場合には，困難度が基準値周辺の，識別力ができるだけ高い項目ばかりを用意するほうが良さそうです。 このように最適な項目の特性はテストの目的に応じて異なります。
IRTでは項目・テスト情報関数が\(\theta\)の関数になっているおかげで，様々な\(\theta\)の値に合わせてオーダーメイドで尺度を構成することができます。 幸いなことにテスト情報関数は局所独立の仮定が満たされていれば項目情報関数の和なので，他の項目との関係を気にすることなく，一つ一つの項目に対して「その項目を追加したらテスト情報量がどの程度増加するか」のみを考えたらよいわけです。 特に自動テスト構成 (automated test assembly) では，コンピュータの力を使って，例えば「\(\theta=(-4,-2,0,2,4)\)の各点における標準誤差が0.2以下＝テスト情報量が25以上になる最小項目数のセット」などの条件にあった最適な項目セットを探し出してくれます35。
ただこれだけのテストの設計を行うためにはかなりの事前準備が必要となります。 というのも項目情報関数は項目パラメータが既知でないとわかりません。 また，目的に応じて最適な項目セットを作成するために，なるべく候補となる項目の数は多いほうが良いです。 もちろんすべての項目のパラメータを一気に推定するわけにはいかないので，多くの場合は複数回に分けて集めたデータをもとにパラメータを等化したりして項目プールを作成する必要があるでしょう。 ということで心理尺度の場合にはなかなかそこまで大量の項目を用意すること自体が難しいので，ここまで理論的な尺度構成を行っている先行事例はたぶんほとんど存在しません。 一応用途としては，100-200項目くらい用意した候補の中から「最高の20項目」を選んで心理尺度を作成する，といった方法はありそうな気がします36。
8.13.3 適応型テスト
テストや尺度の中には通常様々な困難度の項目が混ざっているため，ほぼ確実に「その人にはほぼ意味のない項目」が含まれることになります。 この問題は，先ほど紹介したテストの設計の考え方を持ってしても避けられません。 そこで自動テスト構成をもっと局所的に行うと「個人に応じて最適な項目セットを選び出す」という考えになります。 つまり視力検査のように\(\theta_p\)が高そうな人にはより困難度の高い項目を，\(\theta_p\)が低そうな人にはより困難度の低い項目を出していけば，効率的にテスト情報量を稼ぐことができそうです。 コンピュータ上で行われるコンピュータ適応型テスト (computerized adaptive testing [CAT]) 37では，この考えに基づいて個人ごとに次の項目を出し分けていくことができます。
ここではIRTに基づく最もシンプルな項目選択法を紹介しましょう。 推定の標準誤差は \[ \sigma_{(\hat{\theta}_p|\theta_p)}=\frac{1}{\sqrt{I(\theta_p)}} \] で表されていました。 したがって推定精度をなるべく上げるためには，\(I(\theta_p)\)を大きくする＝項目情報量\(I_i(\theta_p)\)が大きい項目を選んでいけばよいわけです。 ただし\(I_i(\theta_p)\)は「真値\(\theta_p\)における項目情報量」のため，真値\(\theta_p\)がわからない状況では\(I_i(\theta_p)\)の値もわかりません。 そこで，1問解答するたびに推定を行い，その時点での暫定的な推定値\(\hat{\theta}_p^*\)における項目情報量\(I_i(\hat{\theta}_p^*)\)が最大となる項目をとりあえず選び続けたらなんだかうまく行きそうです38。
CATはうまく行けば従来式と同程度の推定精度を半分以下の項目数で達成できたりもする，かなり夢のある技術です。 ただすべての回答者にとって最適な項目を用意できないと効果が薄れてしまうため，単純に自動テスト構成するときよりももっと大量の項目を用意しておく必要があります。 また，テストを実施しリアルタイムで次の項目を選択するためのシステムの構築など，実装に向けたハードルはかなり高いものです。 とりあえずは「そういう技術もIRTによって実現できるんだなぁ」くらいに思っておいていただければ良いと思います39。
Choi & Asilkalkan （2019） によれば，IRTができるRパッケージは45個もあるらしいです（今はもっと多い？）。細かく目的や用途が異なると思うので，特定の場面ではもっと使い勝手の良いパッケージがあったりするかもしれません。あとは計算の精度や速度の違いなど？↩︎
「項目反応理論」というと，ほとんどの場合この後紹介する代表的なモデルを指すのですが，モデル式（項目反応関数）がそれだけしか無いのではなく，実際には様々なモデル式による形が考えられるというのは重要なポイントです。「項目反応理論」という言葉はあくまでも項目依存性や集団依存性を解消したパラメータを得るための理論的な方法を考える土台としての理論である，という認識です。↩︎
頻度論的に言えば，ここでの「確率」は「ある\(\theta_p\)の一個人が何回も回答した際の当てはまるの割合」や「母集団内での\(\theta_p\)の人たちにおける当てはまるの割合」と解釈できます。これはどちらでもOKです。↩︎
もとが学力テストなので「困難度」という呼び方をしています。正解・不正解がない心理尺度でも特に気にせず「困難度」と呼びます。↩︎
もっと近づけるためには1.704倍のほうがよかったり，ズレを極限まで減らすためには\(0.07056x^3+1.5976x\)にすると良い（Bowling et al., 2009）ということは分かっていますが，そこまでの精度は必要ない，ということで，わかりやすい1.7が用いられています。↩︎
正規累積モデルはプロビットモデル，ロジスティックモデルはロジットモデルに対応します。…と考えると「本当は理論的にはどちらを使ったほうが良いか」が明確なケースもあるのかもしれません。実用上はどちらでも良いとされていますが，本当にそのまま1.7倍の変換をするだけで交換可能なのか，というところについてはやや議論があるようです （Cho, 2023）。↩︎
IRTの説明に合わせて誤差の\(\varepsilon_{pi}\)の符号を変えています。↩︎
厳密には，「独立」と「無相関」は別の概念です。理想は独立なのですが，実用上は無相関であれば良いだろう，という感じになっています。厳密な独立とは異なるという意味で，弱局所独立性と呼ぶこともあります。↩︎
「\(\theta_p\)で条件づけた」ということが明確になるようにここだけ\(P(.|\theta_p)\)としていますが，これまでの\(P(.)\)と意味は同じです。↩︎
後半では少しだけ多次元モデルの紹介をしますが，これは一次元モデルの拡張として展開されるものであり，IRTの基本は一次元だといえます。↩︎
もちろん，多因子モデルの場合には，一次元とは限りません。ただ同様に，局所独立性が満たされているならば\(n\)次元性が満たされている，ということはできると思います。↩︎
問1の答えを使わないと問2が解けないような場合，局所独立性が満たされなくなってしまいます。そこで「問1が解けたら1点，問2も解けたら2点」というように得点化してこれを一つの大きな問題として考えることで局所独立性を確保しようという考えに至るわけです。↩︎
\(P(y_{pi}=1)=\exp(\pi_{pi1})P(y_{pi}=0)\)と表せますが，カテゴリ\(k=0\)は存在しない，つまり\(P(y_{pi}=0)=0\)です。ということは\(P(y^*_{pi}=1)_{0,1}=\infty\)となるため，\(\exp(\pi_{pi1})=\infty\)となってしまい計算ができないのです。↩︎
完全には非補償ではない（ある程度なら別の次元の特性値でリカバーできる）という意味で「部分補償モデル」のほうがより正しい名前らしいですが，そこまでこだわるほどの違いはありません。↩︎
あるかはわかりませんが，例えば「2つの特性を両方とも持っている人だけ当てはまると回答する尺度・項目」のようなものを考えた場合には，普通の因子分析よりも非補償型のモデルを適用したほうが良い（回答メカニズムに則している）可能性があります。↩︎
非補償型モデルは，二値ならば
itemtype='PC2PL'またはitemtype='PC3PL'で実行可能です（PCはpartially compensatoryの頭文字）。ただパラメータの推定は結構安定しないので利用する際は気をつけてください。↩︎psych::fa()のように固有値分解など使えばいいじゃないと思われるかもしれませんが，mirt()関数として他の様々なIRTモデルと同じ文法で実装するためにはこうするのが最も自然なのです。↩︎NOHARMという名のソフトウェアで出力される推定結果をもとに計算する方法なので，「NOHARMに基づく方法」などと言ったほうが正確かもしれません。↩︎
これは厳密には弱局所独立性 (weak liocl independence)と呼ばれるものです。一方で，ある反応ベクトル\(\symbf{y}_p\)が得られる確率が各項目の反応確率\(y_{pi}\)の積に等しい，という形でも局所独立性は表されるのですが，こちらは強局所独立性 (strong liocl independence)と呼ばれています。本来はこの強局所独立性が必要なのですが，現実世界の多くの場面では，弱局所独立性が満たされている場合には強局所独立性も満たされていることが多いということで，ここでは弱局所独立性の検証を行っていきます。↩︎
正確に言えば，「\(\symbf{y}_p\)という回答をした人の特性値として\(\hat{\theta}_p\)という値がどの程度もっともらしいか」という感じですが，あえて曲解させています。↩︎
outfitは重み付けない二乗和，infitは情報量（後述）で重み付けた二乗和のようです。↩︎もちろんモデル全体での適合度を考えることも可能です。例えば1パラメータと2パラメータの両モデルを
lavaanで実行して，compareFit()することでどちらのモデルが良さそうかを判断したりというプロセスが考えられます。↩︎もちろんこの回答者\(p\)が赤い項目に正解した場合には，\(\theta_p\)の予測を大きく動かすため，情報量が多くなります。ですがそもそも「回答者\(p\)が赤い項目に正解する」という事象の発生確率がほぼゼロなので，「この項目への回答が\(\theta_p\)の予測にもたらす情報量の期待値」という観点ではやはりほぼゼロなのです。↩︎
項目情報「量」という場合には特定の\(\theta\)における値を指す一方で，項目情報「関数」という場合には\(\theta\)全域を対象とすることが多い気がするので，使い分けに気をつけてください。↩︎
MH検定は
mirtパッケージには用意されていないようです。RではdifR::difMH()などの関数によってMH検定を行うことができます。↩︎項目の数だけ検定を繰り返しているので，実はここで表示されている結果には多重検定の問題が生じています。これを解消するためには，引数
p.adjustを指定すると良いです。例えばp.adjust = "holm"とするとHolmの方法によって調整された\(p\)値を表示してくれるようになります。↩︎この「テスト全体でのグループ間の挙動の違い」はDIFと同じように特異テスト機能 (Differentian Test Functioning [DTF])と呼ばれたりします。↩︎
ここで紹介するモデルを始め，多くの（特殊な）IRTモデルは
mirtパッケージに限らずRのパッケージとしては実装されていなかったりします。そういったモデルについては，自分で尤度関数をプログラムしてoptim()関数による数値最適化を利用する，あるいはstanなどを利用してマルコフ連鎖モンテカルロ法（MCMC）によって推定する，といった方法を取る必要があります。stanを使うためにはまずベイズ統計学を理解する必要があるのがちょっと大変ですが，慣れるとほぼなんでもできるようになります。↩︎\(\theta_p\)の上付き添字のM, A, EはそれぞれMiddle/Mid-point, Attitude/Agree, Extremeなどの頭文字です。↩︎
２つの選択肢の比較のときに，その背後にある効用を比較していて，しかも効用が確率的に変動する，という考え方をサーストンの比較判断の原則 （law of comparative judgment: Thurstone, 1927）と呼ぶため，このモデルはThurstonian IRTと名付けられています。↩︎
これに対して通常のIRTモデルはdominance modelと呼びます（日本語訳がわからない…）。↩︎
Rでの実践例までやると膨大になってしまうので，概念の紹介だけにとどめます。気になる人は自分で調べたり直接質問してください。↩︎
Rでは
equateIRTやplinkといったパッケージによって等化のメソッドが提供されています。↩︎RではlpSolveやglpkAPI，Rsymphonyなどのパッケージを使うと自動テスト構成ができるようになります。↩︎100-200項目すべてに真面目に回答してくれる人たちを数百人用意する必要があるのですが，それさえクリアできれば…↩︎
RではcatRやmirtCATといったパッケージが役に立つかもしれません。↩︎もちろん他にも項目の選び方は色々あります。ざっくりいうと「推定精度をなるべく高める」方向性と「なるべく人によって異なる項目が出るようにする（項目情報が外部に漏れるのを防ぐため）」という方向性があります。↩︎
IRTによらないCATもあることにはあります。ただIRTに基づくCATのほうが圧倒的マジョリティです。↩︎