8  項目反応理論

著者
所属

神戸大学大学院経営学研究科

初版公開日

July 4, 2022

最終更新日

January 22, 2025

abstract
データから「回答者の特性値」と「項目の特性値」を同時に推定する分析法である項目反応理論について,前半では理論的な考え方とmirtパッケージでのやり方を説明し,後半では実践上のいくつかのトピック(多値,多次元モデル・テスト情報量・様々な拡張など)を解説しています。
Keywords

R, 多変量解析, 項目反応理論, mirt

本資料は,神戸大学経営学研究科で2022年度より担当している「統計的方法論特殊研究(多変量解析)」の講義資料です。CC BY-NC 4.0ライセンスの下に提供されています。


このChapterでは項目反応理論 (Item Response Theory [IRT]: Lord & Novick, 1968のお話をしていきます。

後ほど説明しますが,実はIRTの最もメジャーなモデルの一つは,因子分析と数学的に同じものです。 そのためlavaanでも基本的なIRT分析はできるのですが,IRT専用のパッケージも色々あります1。 この講義ではmirtというパッケージを使用します。

コード 8.1: 事前準備
# install.packages("mirt")
library(mirt)
dat <- readRDS("chapter04.rds")
# 毎回列名を指定するのが面倒なので先に作っておく
cols <- paste0("Q1_", 1:10)

8.1 項目反応理論とは

項目反応理論は,古典的テスト理論と同じく,テスト理論と呼ばれるカテゴリに属するものの一つです。 テスト理論の目的をとてもシンプルに言うと,以下の2つです。

  • 回答者の潜在的な特性(心理特性の強さ・能力など)を測定する
  • 測定に使用した項目(心理尺度の各項目・テストの問題など)の性能を評価する

1点目は直感的にもすぐ分かると思います。 2点目について少し補足しておくと,「項目の性能」とは,単純にはテストの問題の難易度などを指しています。 もしも難しすぎてほとんど誰も正解していない問題があったとしたら,その問題には,いわば床効果が出ているわけなのであまり良い問題とは言えません。 また,因子分析的に見たときに因子負荷が低い項目があったとしたら,その項目への回答は,測定したい「回答者の潜在的な特性」を反映していないので,注意が必要かもしれません。 テスト理論では,そのような形で項目の良し悪しを評価することで,「回答者の潜在的な特性」をより良く測定することを目指します

そんなテスト理論の中で,項目反応理論が目指しているところを理解するために,仮想的な状況を考えてみましょう。 Appendix 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で表せるとしたら,項目iiに対する回答者ppの反応(回答)は ypi=f(θp,βi),(8.1) y_{pi} = f(\theta_p,\beta_i), \qquad(8.1) つまりβi\beta_iθp\theta_pによる何らかの関数(項目反応関数: item response function [IRF])になるはずです。 項目反応関数の形はいろいろと考えられますが,例えば多くの心理尺度に対しては,θp\theta_pの値が大きいほど「当てはまる」寄りの回答をする(=ypiy_{pi}の値が大きくなる)と考えられます。 そのようなデータをIRTで分析する場合には,項目反応関数として「θp\theta_pに対して単調増加の関数」を持ってきたら良いわけです。 …という感じでIRTでは,観測されたデータにおいて想定されるypi,θp,βiy_{pi},\theta_p,\beta_iの関係性をうまく表現できるような項目反応関数が色々と提案されてきました。 後ほどIRTの代表的なモデルを紹介しますが,IRTの本質はこのように,回答を(集団によらない)項目パラメータと(項目によらない)個人特性値に分離して分析することにあります2

ということで,一般的なIRTの利用場面では,

  1. 特定の関数形(IRF)を仮定する
  2. そのモデルのもとで項目パラメータを推定する
  3. 推定した項目パラメータをもとに,項目の性能を評価する
  4. 合わせて個人特性値を推定する
  5. 推定した個人特性値を使ってあとはご自由に

といったことが行われます。 それでは,まずは「とりあえずIRTで分析」ができるようになるために必要な最低限の知識を学んでいきましょう。

8.2 基本的なモデル

最もシンプルなIRTモデルを考えるため,一旦ここからは項目への反応を二値(当てはまらない・あてはまる)に限定して考えてみます。 というのも,もともとIRTは「正解・不正解」の二値で考えるような学力テストなどを想定して作られているもののためです。

8.2.1 正規累積モデル

以前カテゴリカルデータの相関(ポリコリック・ポリシリアル相関)の計算時には,背後にある連続量を考えるという説明をしました。 ここでも,二値反応ypiy_{pi}の背後に連続量θp\theta_pが存在し,θp\theta_pが閾値を超えていたら1,超えていなければ0と回答すると考えます。 そして項目の閾値をβi\beta_iと表すと, ypi={0θp<βi1θpβi(8.2) y_{pi} = \left\{ \begin{aligned} 0 & \quad \theta_p < \beta_i\\ 1 & \quad \theta_p \geq \beta_i \end{aligned} \right. \qquad(8.2) となります( 図 8.1 )。 βi\beta_iの値が大きい項目ほど,ypi=1y_{pi}=1になるθp\theta_pの範囲が狭くなるため,「当てはまらない」と回答する人の割合が高くなるだろう,ということです。

図 8.1: 二値反応と潜在的な連続量の関係

ただし,このままでは能力θp\theta_pを持つppさんは難易度がβi\beta_i以下の問題には,どんなジャンル・内容の問題でも100%正解することになってしまいます(θpβi\theta_p \geq \beta_iなので)。 しかし,もちろん客観的に見た難易度が全く同じであっても,内容によってppさんが答えられる・答えられない問題があるはずです。 さらに言えば,ppさんの体調や時の運(たまたま昨日一夜漬けした内容に入っていた,など)によっても答えられるかどうかは変わるはずです。 というように,心理尺度やテストへの回答では,常に「その他の要因」による変動を考えてきました。 言い換えると,古典的テスト理論でも回帰分析でも因子分析でも,観測値は必ず何らかの誤差を伴って発生するものだと考えているわけです。 そこで先程の(8.2)式に誤差項εpi\varepsilon_{pi}を追加して ypi={0(θpεpi)<βi1(θpεpi)βiεpiN(0,σε2)(8.3) 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) \qquad(8.3) とします。なお誤差項は,回帰分析などと同じく,最も一般的な「平均0の正規分布に従う」と仮定してあります。

誤差項を追加することによって,「特性値がθp\theta_pの人が項目jjに『当てはまる』と回答する確率3」を考えることが出来ます。 図 8.1 では「全ての人のθp\theta_p」の分布を示していた一方で, 図 8.2 では「θp\theta_pが特定の値の人」における分布になっている点に注意してください。

図 8.2: 誤差を追加すると

さらにβi\beta_iθp\theta_pを一つの関数f(θp,βi)f(\theta_p,\beta_i)として見るためにβi\beta_iを移項して, ypi={0(θpβiεpi)<01(θpβiεpi)0εpiN(0,σε2)(8.4) 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) \qquad(8.4) とします( 図 8.3 )。 当然ながら, 図 8.2図 8.3 は,閾値と分布が平行移動しただけなのでP(ypi=1)P(y_{pi}=1)は変わりません。

図 8.3: βi\beta_iを移項すると

図 8.3 の色つきの確率分布は,正規分布N(θpβi,σε2)N(\theta_p - \beta_i,\sigma_{\varepsilon}^2)と表すことができます。 そして式から,ypi=1y_{pi}=1となるのは誤差εpi\varepsilon_{pi}(θpβi)(\theta_p - \beta_i)以下のときだとわかります。 したがって,仮にεpi\varepsilon_{pi}が標準正規分布に従う(σε2=1\sigma_{\varepsilon}^2=1)としてみると,「特性値がθp\theta_pの人が項目jjに当てはまると回答する確率」( 図 8.4 の青い部分)は P(ypi=1)=P(εpiθpβi)=(θpβi)12πexp(z22)dz(8.5) \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} \qquad(8.5) と表すことができます。教科書などで見たことがある形かもしれませんが,この積分の中身は標準正規分布の確率密度関数です。 このように,P(ypi=1)P(y_{pi}=1)を標準正規分布の累積確率で表すモデルのことを正規累積モデル(normal ogive model)と呼びます。 この後項目パラメータを1つ増やすので,それに対してこの正規累積モデルは,特に1パラメータ正規累積モデルと呼ばれます。

図 8.4: 最終的な反応確率の分布

ここまでは,途中で話を簡単にするために勝手にσε=1\sigma_{\varepsilon}=1として進めていましたが,ここからはもう少し一般化してσε2\sigma_{\varepsilon}^2が項目ごとに異なる(σε2=1/αi\sigma_{\varepsilon}^2=1/\alpha_i)とします。 このとき,誤差の確率分布はεpiN(0,1αi)\varepsilon_{pi} \sim N\left(0,\frac{1}{\alpha_i}\right)となるため,両辺をαi\alpha_i倍してαiεpiN(0,1)\alpha_i\varepsilon_{pi} \sim N\left(0,1\right)と表すことができます。 αiεpi\alpha_i\varepsilon_{pi}が標準正規分布に従うならば,(8.5)式の正規累積モデルは P(ypi=1)=P[εpi(θpβi)]=P[αiεpiαi(θpβi)]=αi(θpβi)12πexp(z22)dz(8.6) \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} \qquad(8.6) という形で書き換えることができます。また正規累積モデルが出てきました。(8.5)式との違いは積分の上限がαi(θpβi)\alpha_{i}(\theta_p - \beta_i)に変わった点だけです。 いま,突拍子もなくαi\alpha_{i}というパラメータを導入しましたが,とりあえずこれによって項目iiが持つパラメータは(βi,αi)(\beta_i, \alpha_i)の2つになりました。 ということで,(8.6)式の正規累積モデルは2パラメータ正規累積モデルと呼ばれます。

8.2.2 項目特性曲線

ここで,2つの項目パラメータ(βi,αi)(\beta_i, \alpha_i)が何を表しているのかを明らかにしていくため, 図 8.3 を別の形で表してみます。

図 8.5: パラメータと正答率の関係

図 8.5 の下半分は,横軸にθpβi\theta_p - \beta_iを,縦軸にP(ypi=1)P(y_{pi}=1)をとったグラフです。 図の上半分にあるように,θpβi\theta_p - \beta_iの値が大きくなるほど分布の青い部分が増える,つまりP(ypi=1)P(y_{pi}=1)の値が大きくなることを表しています。

これを踏まえて,パラメータβi\beta_iの挙動を見るために,横軸にθp\theta_pを,縦軸にP(ypi=1)P(y_{pi}=1)をとったグラフを,βi\beta_iの値を変えながらいくつか重ねてみます( 図 8.6 )。 誤差は平均0の正規分布に従うため,P(ypi=1)=0.5P(y_{pi}=1)=0.5になるのは,θp=βi\theta_p = \beta_iのときです。 したがってβi\beta_iの値が大きい項目ほど,θp\theta_pの値が大きくないと「当てはまる」と回答する確率が高くなりません。 そのような意味で,βi\beta_i項目困難度(item difficulty)と呼ばれます4

図 8.6: βi\beta_iの役割

次はαi\alpha_iの役割を見るために, 図 8.6 よりも誤差分散(αi\alpha_i)が小さい場合や大きい場合を見てみましょう。 図 8.7 の左側に示されているように,誤差分散が小さいほど分布の広がりがなくなるため,θpβi\theta_p - \beta_iの値の変化に対するP(ypi=1)P(y_{pi}=1)の変化が急激になっています。

図 8.7: 左:誤差分散が小さい場合,右:誤差分散が大きい場合

これを踏まえて,βi\beta_iの時と同じように横軸にθp\theta_pをとり,縦軸にP(ypi=1)P(y_{pi}=1)をとったグラフを,今度はαi\alpha_iの値を変えながらいくつかプロットしたものが 図 8.8 です(βi=0\beta_i=0)。 誤差分散の逆数であるαi\alpha_iの値が大きい項目ほど,θp\theta_pが低い人と高い人の間でP(ypi=1)P(y_{pi}=1)が大きく変化しています。 言い換えるとαi\alpha_iの値が大きい項目ほど,その項目への回答によってθp\theta_pの高低を識別する能力が高いといえます。 ということでαi\alpha_iのことを項目識別力(item discrimination)と呼びます。

図 8.8: αi\alpha_iの役割

図 8.6図 8.8 のように,横軸にθp\theta_pをとり,縦軸にP(ypi=1)P(y_{pi}=1)をとったグラフは,いわば回答者の特性値に対する項目の特性を表しています。 ということで,これらのグラフのことを項目特性曲線(item characteristic curve [ICC])と呼びます。

8.2.3 ロジスティックモデル

2パラメータ正規累積モデルの式を再掲しておきましょう。 P(ypi=1)=αi(θpβi)12πexp(z22)dz \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)=exp(x)[1+exp(x)]2,F(x)=11+exp(x)(8.7) f(x) = \frac{\exp(-x)}{\left[1+\exp(-x)\right]^2}, \ \ \ \ \ F(x) = \frac{1}{1+\exp(-x)} \qquad(8.7) と表され,確率密度関数を正規分布と重ねて見ると 図 8.9 のようになります。 こうしてみると,裾の重さは違いますが確かにロジスティック分布は正規分布のような左右対称形になっているようです。

図 8.9: 正規分布とロジスティック分布

また,ロジスティック分布のxxをおよそ1.7倍すると,累積分布関数が正規分布とかなり近くなることが知られています( 図 8.10 )5

図 8.10: 1.7倍ロジスティック分布

ということで,正規分布の代わりにロジスティック分布を使った(2パラメータ)ロジスティックモデルP(ypi=1)=11+exp[1.7αi(θpβi)](8.8) \begin{aligned} P(y_{pi}=1)&= \frac{1}{1+\exp\left[-1.7\alpha_{i}(\theta_p - \beta_i)\right]} \end{aligned} \qquad(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パラメータモデルに加えて「当て推量」のパラメータcic_iが追加されます。 P(ypi=1)=ci+1ci1exp[αi(θpβi)](8.9) P(y_{pi}=1)= c_i + \frac{1-c_i}{1-\exp\left[\alpha_i(\theta_p-\beta_i)\right]} \qquad(8.9) cic_iは項目特性曲線の下限を操作します。θp\theta_pがどんなに低い人でも,必ずP(ypi=1)P(y_{pi}=1)cic_i以上の値になります。 例えば4択問題であれば,どんなに能力が低くても,勘で答えれば正答率は1/4になります。これがcic_iです。
【4パラメータ】
4パラメータモデルでは,3パラメータモデルに加えて「上限」のパラメータdid_iが追加されます。 P(ypi=1)=ci+dici1exp[αi(θpβi)](8.10) P(y_{pi}=1)= c_i + \frac{d_i-c_i}{1-\exp\left[\alpha_i(\theta_p-\beta_i)\right]} \qquad(8.10) 4パラメータモデルでは,θp\theta_pがどんなに高い人でも,必ずP(ypi=1)P(y_{pi}=1)did_i以下の値になります。 どんなに能力がある人でも100%にはならないもの,例えばバスケのフリースローなどをイメージすると良いかもしれません。
【5パラメータ】
5パラメータモデルでは,4パラメータモデルに加えて「非対称性」のパラメータeie_iが追加されます。 P(ypi=1)=ci+dici(1exp[αi(θpβi)])ei(8.11) 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}} \qquad(8.11) 通常のIRTモデルでは,P(ypi=1)P(y_{pi}=1)の動き方は0.5を軸に対称になっています。 これに対して,例えば「最初は急激に成功率が上がるが,上に行くほどきつくなる」的なことを表現したい場合には非対称性を考慮する必要があるでしょう。

ただ,これらのモデルは安定した推定に必要なサンプルサイズが莫大であったり,そもそも推定量に一致性がないなどの問題を抱えていることから,実務場面で使われることはほとんどありません(強いて言えば,当て推量パラメータを「選択肢数分の1」に固定した3パラメータモデルなどはありえる)。 ので,「そういうのもあるんだなぁ」くらいに思っておけば良いと思います。

8.3 因子分析との関係

IRTの基本的なモデルが分かったところで,IRTモデルと因子分析の特定のモデルが数学的には等価であること Takane & de Leeuw, 1987 を紹介しておきたいと思います。 このことは,心理尺度の分析において因子分析だけでなくときにはIRTも使用できることを意味しています。 それぞれの分析手法は出自が異なる以上,分析の焦点や技法にも違いがあるのですが,適切な設定のもとでは,因子分析とIRTのいいとこ取り(あるいは併用)ができるということは理解しておくと良いかもしれません。

多母集団同時分析のところ Section 7.9 で説明したように,標準化していない(=切片が0でない)因子分析モデル(1因子モデル)は ypi=τi+fpbiεpiεpiN(0,σε)(8.12) y_{pi} = \tau_i + f_pb_i - \varepsilon_{pi} \qquad \varepsilon_{pi}\sim N(0,\sigma_{\varepsilon}) \qquad(8.12) という形になっています7。 そしてカテゴリカル因子分析では,観測変数ypiy_{pi}はその背後にある連続量によって決まる,という考え方をしていました。 心理尺度の項目が2件法だとしたら,観測変数ypiy_{pi}とその背後の連続量の関係は 図 8.3 と同じようにして 図 8.11 のように表すことができます。

図 8.11: 因子分析モデルの背後に

したがって,標準正規分布に従う誤差εpi\varepsilon_{pi}τi+fpbi\tau_i + f_pb_iより小さいときにypi=1y_{pi} = 1となるので,その確率は P(ypi=1)=P(εpiτi+fpbi)=(τi+fpbi)12πexp(z22)dz(8.13) \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} \qquad(8.13) となります。これを2パラメータ正規累積モデル( 8.6 式)に対応させて見ると, τi+fpbi=bi(fp+τibi)=αi(θpβi)(8.14) \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} \qquad(8.14) ということで,(fp,bi,τi)=(θp,αi,αiβi)(f_p, b_i, \tau_i) = (\theta_p, \alpha_{i}, \alpha_{i}\beta_i)としてみればこれら2つのモデルは完全に同一だということがわかります。

IRTでは,(8.14)式に基づく(因子分析的な)パラメータ化を行うことがあります。 つまり2パラメータロジスティックモデルで言えば P(ypi=1)=11+exp[αi(θpβi)]=11+exp[(αiθpτi)](8.15) \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} \qquad(8.15) という形で,困難度パラメータの代わりに切片パラメータを求めるということです。 これもどちらの設定を使用しても良いのですが,ソフトウェアによってはこの切片パラメータを使用する方法がデフォルトになっていることもあるので,出力を見る際は気をつけましょう。 なおIRTの枠組みでは,どちらかというと困難度パラメータβi\beta_iを用いる解釈のほうがよく使われます。 推定結果で切片パラメータが出てきた場合には,自分でβi=τiαi\beta_i = \frac{\tau_i}{\alpha_i}と変換する必要があるかもしれません。

8.4 RでIRT

IRTの基本を学んだところで,いよいよ実際にRでIRTを実行してみます。 これまでの説明に合わせて,まずは二値で試してみるため,それ用のデータを作成します。 元々のdatは6件法への回答データでしたが,各選択肢は1から3が「当てはまらない」寄り,4から6が「当てはまる」寄りでした。 そこで,「どちらかと言えば当てはまると回答したかどうか」の二値データに変換したものを使用していきます。

コード 8.2: 二値データの作成
# 使用する部分だけ取り出す
# 細かいことは考えずに10項目をそのまま使います
dat_bin <- dat[, cols]
# 1から3は0,4から6は1に変換して二値化
dat_bin <- (dat_bin >= 4)
# as.numericで数字に変換
dat_bin <- apply(dat_bin, 2, as.numeric)

そしてIRTでの分析には,mirtパッケージの中にあるmirt()関数を使用します。 基本的な引数は以下の3つです。他にも色々ありますが,後ほどもう少し紹介します。

data
推定に用いるデータです。lavaanの関数と異なり,推定に使う変数のみが入ったデータを用意する必要があります。
model
lavaanのように各変数と因子の関係を記述したモデル式です。先程説明したように,IRTは因子分析と同じものなので,モデル式を記述することで多次元のIRTモデルも実行可能になる,というわけです(詳細は後ほど)。変数名が連番になっている場合,後述のようにハイフンでまとめて指定できるのでおすすめです。
itemtype
モデルの指定です。たぶん正規累積モデルは選べず,すべてロジスティック関数ベースのモデルです。これまでに紹介したモデルでは
  • 'Rasch':ラッシュモデル,つまり1パラメータモデルです。後述のものに合わせて'1PL'と指定することはできないので気をつけてください。
  • '2PL':2パラメータモデルです。同じように'3PL', '4PL'も指定可能です。
コード 8.3: Rではじめての項目反応理論
# モデル式の指定
# lavaanとは違い,ふつうにイコールで結びます
model <- "f_1 = Q1_1-Q1_10"
# ハイフンの指定は以下の書き方の簡略化
# model <- 'f_1 = Q1_1 + Q1_2 + Q1_3 + (中略) + Q1_10'

result_2plm <- mirt(dat_bin, model = model, itemtype = "2PL")

得られた推定値を見る場合には,coef()関数を使います。

コード 8.4: 推定値を見る
coef(result_2plm)
$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とすることで,データフレームの形で返してくれるようになります。

コード 8.5: 推定値をシンプルに見る
coef(result_2plm, simplify = TRUE)
$items
         a1     d g u
Q1_1  0.416 1.262 0 1
Q1_2  1.023 2.359 0 1
Q1_3  1.027 1.911 0 1
Q1_4  1.033 1.737 0 1
Q1_5  1.010 1.791 0 1
Q1_6  1.203 1.911 0 1
Q1_7  1.367 1.667 0 1
Q1_8  1.225 1.582 0 1
Q1_9  1.339 1.365 0 1
Q1_10 1.270 0.031 0 1

$means
f_1 
  0 

$cov
    f_1
f_1   1

最初の$item部分に表示されているものが項目パラメータです。左から

a1
識別力パラメータの推定値
d
切片パラメータの推定値
g
当て推量パラメータ (guessing) の推定値(2パラメータモデルなら全て0)
u
上限パラメータ (upper bound) の推定値(2パラメータモデルなら全て1)

となっています。このように,mirtはデフォルトでは切片パラメータを出力します。困難度パラメータを見たい場合にはd列をa1列で割るか,引数IRTpars = TRUEをつけて出力しましょう。

コード 8.6: 困難度パラメータを見る
coef(result_2plm, simplify = TRUE, IRTpars = TRUE)
$items
          a      b g u
Q1_1  0.416 -3.031 0 1
Q1_2  1.023 -2.307 0 1
Q1_3  1.027 -1.861 0 1
Q1_4  1.033 -1.681 0 1
Q1_5  1.010 -1.773 0 1
Q1_6  1.203 -1.588 0 1
Q1_7  1.367 -1.220 0 1
Q1_8  1.225 -1.292 0 1
Q1_9  1.339 -1.019 0 1
Q1_10 1.270 -0.025 0 1

$means
f_1 
  0 

$cov
    f_1
f_1   1

これでd列の代わりにb(困難度パラメータ)列が表示されるようになりました。 出力を見ると,例えばQ1_1では,1.262/0.4163.0311.262 / 0.416 \simeq 3.031ということで,確かにd列をa1列で割った値(の符号を入れ替えたもの)が表示されています。

$itemの下に表示されている$means$covはそれぞれ特性値(因子得点)θ\thetaの平均値と分散です。 この時点では一般的な因子分析の制約(平均0,分散1)に沿った値のみが表示されているだけなので特に注目する必要はありません。 後ほど紹介する多次元モデルや多母集団モデルになると,ここに表示される値も重要な意味を持つことになります。

plot()関数を使うと,引数typeに応じて,推定結果に基づく様々なプロットを全項目まとめて出すことができます。 図 8.6 のようなICCを出したい場合にはtype = "trace"としてあげてください。

コード 8.7: 項目特性曲線
plot(result_2plm, type = "trace")
図 8.12: 項目特性曲線

デフォルトではすべての項目のICCを縦横に並べてくれるのですが,項目数が多くなると見るのが大変そうです。 この場合,引数facet_items = FALSEとすると,すべての項目のプロットを重ねて表示してくれます。 また引数which.itemsを指定すると,特定の項目だけのICCを出すことも出来ます。

コード 8.8: 特定の項目のICCだけを重ねて表示
plot(result_2plm, type = "trace", which.items = 1:5, facet_items = FALSE)
図 8.13: 最初の5項目の項目特性曲線

ということで 図 8.12図 8.13 を見ると,識別力がダントツで低い項目1では,ICCの傾きが緩やかになっていることがわかります。

IRTによる項目分析では,まずこのように推定されたパラメータをもとにICCを描くのが重要です。これは 図 4.16 と同じように解釈することができるので,例えば

  • θ\thetaの広い領域でICCがほぼ上か下に貼り付いている(βi\beta_iがとても小さいか大きい)
  • ICCの傾きが緩やかである(αi\alpha_iがとても小さい)

ような項目は,θ\thetaの測定という面からはあまり良くない項目であろう,ということがわかります(除外するかは要検討)。 この他にも,複数の項目を組み合わせたときにどうなるかという分析などもあるわけですが,それはおいおい紹介していきます。

話は変わって,項目パラメータを固定した上での特性値θp\theta_pの推定値を出す場合は,fscores()関数を使います。 このあたりは因子分析やCFAでの「項目パラメータを推定」→「項目パラメータを固定して個人特性値を推定」と全く同じです。

コード 8.9: $ heta_p$の推定値を出す
head(fscores(result_2plm))
              f_1
[1,] -1.537468218
[2,] -0.007789283
[3,]  0.226340769
[4,] -0.643698451
[5,]  0.085994618
[6,]  0.788036256

8.5 IRTの前提条件

IRT,なかでも2パラメータモデルは基本的には因子分析と同じようなものなのですが,理論的背景の違いなどから,因子分析を行う際とは留意すべき点も少し異なってきます。 ここでは,IRTを適用するにあたって重要となる2つの前提条件について見ていきます。

8.5.1 局所独立性

局所独立性(local independence)の仮定は,因子分析のときと概ね同じです。 探索的因子分析では,独自因子の間が完全に無相関であるという仮定をおき,その中でデータの相関行列を最もよく復元できる因子負荷行列と独自因子の分散(𝐁𝐁+𝚿\symbf{B}^\top\symbf{B} + \symbf{\Psi})を求めていました。 独自因子の間に相関があるということは,その当該項目の間に何か別の共通因子があるということです。

局所独立性をもう少し固い言い方で説明すると,「θp\theta_pで条件づけたとき,項目iiと項目j(ij)j~(i\neq j)への回答は完全に独立」となります。 項目への回答は二値なので,2×22\times2クロス表をイメージしてみましょう。 θp\theta_pが特定の値の人を1000人ほど集めて(つまりθp\theta_pで条件づけて) 表 8.1 のように2つの項目i,ji,jの回答のクロス表を作ると,もし局所独立性が満たされているならば連関はゼロになります。 項目iiに「当てはまる」と回答した人を見ても,「当てはまらない」と回答した人を見ても,項目jjに「当てはまる」と回答した人の割合は全体の割合と同じ6:4になっています。

表 8.1: 局所独立な2項目
項目jj
項目ii 当てはまらない 当てはまる
当てはまらない 180 420 600
当てはまる 120 280 400
300 700 1000

一方で,θp\theta_p以外の別の共通因子がある場合,クロス表に連関が現れます。 例えば「他者への助けを求める傾向」の尺度において,宗教観の影響を受けると思われる項目(e.g., 困ったときは神に祈る)どうしのクロス表を見てみると, 表 8.2 のように

  • 項目iiに「当てはまる」と答えた人では340/400=85%340/400=85\%の人が項目jjでも「当てはまる」と回答
  • 項目iiに「当てはまらない」と答えた人では360/600=60%360/600=60\%の人が項目jjでも「当てはまる」と回答

しているため,項目iiの回答と項目jjの回答が関連を持っていることになります。 項目iiに「当てはまる」と答えた人のほうが項目jjでも「当てはまる」と回答する確率が高い,ということです。

表 8.2: 局所独立ではない2項目
項目jj
項目ii 当てはまらない 当てはまる
当てはまらない 240 60 600
当てはまる 360 340 400
300 700 1000

局所独立の仮定とは,このようにクロス表を作成したときに,どの項目ペア間でも,どの特性値の値でも連関がない,ということを指しています8

局所独立性がなぜ重要なのでしょうか。その理由の一つは計算上の問題です。 IRTでは基本的に最尤法によってパラメータを推定します。 もしも局所独立性が満たされていれば,θp\theta_pで条件づけたときの各項目に対する反応確率は独立なので,回答者ppに対する尤度関数は単純にすべての項目をかけ合わせた P(𝐲p|θp)=i=1IP(ypi=1|θp)ypiP(ypi=0|θp)1ypi(8.16) 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}} \qquad(8.16) と表せます9。 これに対してもしも局所独立が満たされていない場合( 表 8.2 ),P(ypj=1|θp)P(y_{pj}=1|\theta_p)の値がypiy_{pi}によって変わることになります。 そのため単純に掛け算ができなくなってしまい,もっと複雑な計算が必要になってしまうのです。

ということで,IRTモデルは「局所独立性が満たされている」という仮定にもとづいて定式化されています。 したがって,データを収集する時点で局所独立性が満たされていることをしっかりと確認しておく必要があります。 なお後半では,局所独立性が満たされているかを事後的に評価する指標を紹介します。

8.5.2 一次元性

これまで紹介してきたIRTモデルはいわば一因子の因子分析モデルです。 因子分析では「項目をグルーピングする」ことに重きを置きがちだったので,因子数は「いい感じにグループができる数」にすることが一般的でした。 これに対してIRTはもともと単一のテストに対する単一の能力を測定することが主たる目的と言えるため,基本的には因子数は1となります10。 そういう意味では,心理尺度とは相性が良い側面もあると言えそうです。 そんなわけで,IRTの重要な前提条件としてすべての項目が同じ能力・特性のみを反映していること(一次元性)が必要となります。

一次元性を検証するためのフレームワークは,多くの場合(探索的・確認的)因子分析と共通です。 つまり, Section 6.11 で紹介したスクリープロットや平行分析,MAPなどの基準によって「因子数1」が提案されたら良い,ということです。 ですが,因子分析的な次元性の確認方法では,項目数が増えるとなかなか「因子数1が最適」とはならないケースが出てきます。 したがってIRTでは,「一次元とみなしても問題ないか」といった別の視点から一次元性の評価を行うことがあります。 これは,例えば一因子のCFAを適用した際の適合度指標(RMSEAやCFA, AGFIなど)が許容範囲にあるかによって確認することが出来ます。

またIRTでは,局所独立性の観点から一次元性を確認する指標もいくつか提案されています。 一次元性が満たされているということは,回答行動に影響する共通因子が一つ(θp\theta_p)だけということです。 したがって一次元性が満たされていない場合には, Section 6.4.1図 6.9 で少し説明したように,共通因子で説明出来ない部分(独自因子)の間に共分散が発生することになります。 そのように考えると,一次元性は局所独立と同じようなものだとみなすことができるわけです。 厳密には局所独立性が満たされている場合,一次元性も満たされていると言って良いだろう,ということになります11。 ということで,より洗練された一次元性確認の方法として,例えばDIMTESTStout et al., 1996やDETECTZhang & Stout, 1999という名前の方法が提案されていたりします。

局所独立と一次元性

一次元性が満たされているだけでは局所独立性が満たされるとは限りません。 例えば

  • 変数xxの平均値を求めよ
  • 変数xxの分散を求めよ

という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モデルを組み合わせて多値反応を表現します。 項目iiに対する回答者ppの回答ypi=k(k=1,2,,K)y_{pi}=k~(k=1,2,\cdots,K)について「kk以上のカテゴリを選ぶ確率」を考えると,これはまだ二値(kk以上 or kkより小さい)なので,2パラメータロジスティックモデルならば P(ypik)=11+exp[αi(θpβik)](8.17) P(y_{pi}\geq k) = \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_{ik})\right]} \qquad(8.17) と表せます。なお,困難度パラメータは「項目iiのカテゴリkk」ごとに用意されるためβik\beta_{ik}としています。

異なるカテゴリを選ぶ確率は当然排反な事象なので,二値モデルを組み合わせると「ちょうどkk番目のカテゴリを選ぶ確率」は P(ypi=k)=P(ypik)P(ypik+1)(8.18) P(y_{pi} = k) = P(y_{pi}\geq k) - P(y_{pi}\geq k+1) \qquad(8.18) と表すことができます。 ただし,端のカテゴリに関する確率はそれぞれ P(ypi1)=1P(ypiK+1)=0(8.19) \begin{aligned} P(y_{pi}\geq 1) &= 1 \\ P(y_{pi}\geq K+1) &= 0 \end{aligned} \qquad(8.19) とします。 「1以上のカテゴリを選ぶ確率」は当然1になるため,これに対する困難度パラメータβi1=\beta_{i1}=-\inftyとなります。 同様に,「K+1K+1以上のカテゴリを選ぶ確率」は0になります。P(ypi=K|θp)P(y_{pi} = K|\theta_p)を考えるときだけ仮想的な「さらに上のカテゴリ」を想定する,ということです。 このようにGRMは二値型モデルの素直な拡張のため,数学的には正規累積型のGRMはカテゴリカル因子分析と同じものだとみなすことができます Takane & de Leeuw, 1987

図 8.14 に,4件法の項目に対する段階反応モデル((αi,[βi2,βi3,βi4])=(1,[1.5,0,2])(\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(1,[-1.5, 0, 2]))のイメージを示しました。 左右で色の塗り方や分布の位置は変えていません。点線の位置だけです。

図 8.14: 段階反応モデル

左は,θp\theta_pの位置に応じて各カテゴリの反応確率が変化する状況を表しています。 反応確率は全カテゴリ合計で1となっているため,点線の位置での各カテゴリの長さの比率がそのままカテゴリ選択確率となっています。

右は,βik\beta_{ik}の位置に点線を引いています。 カテゴリの困難度βik\beta_{ik}は二値モデルの時と同じようにP(ypik)=0.5P(y_{pi}\geq k)=0.5になるθp\theta_pの値を表しています。 そのため,上の確率密度分布を見ると,色が変わるポイントはすべての分布で同じになっていることがわかります。 したがって,θp\theta_pの値が大きい(正規分布が右に動く)ほど,またβik\beta_{ik}の値が小さい(閾値の点線が左に動く)ほど,上位カテゴリを選択する確率が高くなっていくわけです。

続いて, 図 8.15 には,いくつかの項目パラメータのセットにおける項目特性曲線(ICC)を示しました。 図 8.14 の下半分ではP(ypik)P(y_{pi}\geq k)のプロットを示しましたが,通常「項目特性曲線」というとP(ypi=k)P(y_{pi} = k)のプロットです。 二値項目の場合はどちらでも同じですが,多値項目の場合は異なるため改めてICCを示しています。 なお,多値型モデルの場合は 図 8.15 の一本一本の線のことを項目反応カテゴリ特性曲線(item response category characteristic curve [IRCCC])と呼ぶこともあります。 が,そこまで細かい名前を覚える必要はないでしょう。

図 8.15: GRMの項目特性曲線

左上のプロットは, 図 8.15 と同じ項目パラメータ((αi,[βi2,βi3,βi4])=(1,[1.5,0,2])(\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(1,[-1.5, 0, 2]))でのIRCCCです。 基本的には,このようにθp\theta_pの値に応じて左から順に各カテゴリの選択確率のピークが訪れます。 点線は中間カテゴリの選択確率が最も高くなる箇所ですが,これは具体的にカテゴリ困難度βik\beta_{ik}βi,k+1\beta_{i,k+1}のちょうど中間になります。

右上のプロットは,左上からすべてのカテゴリ困難度の値を1ずつ大きくしました((αi,[βi2,βi3,βi4])=(1,[0.5,1,3])(\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(1,[-0.5,1,3]))。 二値モデルの時と同じように,すべてのカテゴリ困難度が同じ値だけ変化するときにはIRCCCは平行移動します。

左下のプロットは,左上から識別力を半分にしました((αi,[βi2,βi3,βi4])=(0.5,[1.5,0,2])(\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4}])=(0.5,[-1.5, 0, 2]))。 また右下のプロットは,左上から識別力を2倍にしています((αi,[βi2,βi3,βi4])=(2,[1.5,0,2])(\alpha_i,[\beta_{i2},\beta_{i3},\beta_{i4]})=(2,[-1.5, 0, 2]))。 識別力は各IRCCCの山のきつさを調整しています。 識別力が低い場合,θp\theta_pの大小に対して各カテゴリの選択確率があまり変化しなくなるため,反対に「あるカテゴリを選択した」という情報がθp\theta_pの推定に及ぼす影響も小さくなってしまいます。

また,左下のプロットのように,どのθp\theta_pの値においてもあるカテゴリの選択確率が最大とならないこともありえます。 識別力が低い場合や,カテゴリ困難度が近すぎる場合,あるいはそもそもカテゴリ数が多い場合にはそのような現象が起こりえます。 ただ何も問題ではないので,そのようなカテゴリが出現した場合には「選ばれにくいんだなぁ」くらいに思っておいてください。

8.6.2 一般化部分採点モデル

有名な多値型モデルのもう一つが,一般化部分採点モデル(generalized partial credit model [GPCM]: Muraki, 1992です。 GPCMでは,隣のカテゴリとの関係を別の視点から考えます。 ということで,まずは「隣のカテゴリとの関係」を二値(2パラメータロジスティック)モデルで考えてみましょう。 二値モデルでは,「カテゴリ0を選ぶ確率」と「カテゴリ1を選ぶ確率」の和は当然1になります。 そして,「カテゴリ1を選ぶ(e.g., 『当てはまる』を選ぶ,正解する)」確率は P(ypi=1)=P(ypi=1)P(ypi=1)+P(ypi=0)=11+exp[αi(θpβi)]=exp[αi(θpβi)]1+exp[αi(θpβi)](8.20) \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} \qquad(8.20) となります。 つまりこの式は,隣のカテゴリとの二択において,当該カテゴリを選ぶ確率という見方ができるわけです。 これを多値(k=1,2,,K)(k=1,2,\cdots,K)に置き換えると,「カテゴリk1k-1kkの二択だったらkkの方を選ぶ確率」に拡張することができそうです。 ということで,GPCMではこの確率P(ypi*=k)k1,kP(y^*_{pi}=k)_{k-1,k}P(ypi*=k)k1,k=P(ypi=k)P(ypi=k)+P(ypi=k1)=11+exp[αi(θpβik)]=exp[αi(θpβik)]1+exp[αi(θpβik)]=exp[πpik]1+exp[πpik](8.21) \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} \qquad(8.21) と表します。なお,この後の式展開を考えてπpik=αi(θpβik)\pi_{pik}=\alpha_{i}(\theta_p - \beta_{ik})と置き換えています。

ここから,全カテゴリの中でカテゴリkkを選ぶ確率を導出しましょう。 (8.21)式をP(ypi=k)=P(y_{pi}=k)=の形に変形すると, P(ypi=k)=P(ypi=k1)P(ypi*=k)k1,k1P(ypi*=k)k1,k(8.22) 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}} \qquad(8.22) となります。最後の右辺のP1P\frac{P}{1-P}の部分は2つのカテゴリ間のオッズを表しているため, P(ypi*=k)k1,k1P(ypi*=k)k1,k=exp(πpik)(8.23) \frac{P(y^*_{pi}=k)_{k-1,k}}{1-P(y^*_{pi}=k)_{k-1,k}} = \exp(\pi_{pik}) \qquad(8.23) を用いて P(ypi=k)=P(ypi=k1)exp(πpik)(8.24) P(y_{pi}=k) = P(y_{pi}=k-1)\exp(\pi_{pik}) \qquad(8.24) という漸化式が得られ,特定のカテゴリkkを選ぶ確率が表現できるようになります。 例えば4カテゴリの場合には, P(ypi=1)=P(ypi=1)P(ypi=2)=P(ypi=1)exp(πpi2)P(ypi=3)=P(ypi=1)exp(πpi2)exp(πpi3)P(ypi=4)=P(ypi=1)exp(πpi2)exp(πpi3)exp(πpi4)(8.25) \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} \qquad(8.25) という形でP(ypi=1)P(y_{pi}=1)を起点に「一個上のカテゴリに上がるときのオッズ」の積によってすべてのカテゴリの反応確率が表されるわけです。 ということで,この式を一般化して P(ypi=k)=P(ypi=k1)exp(πpik)=P(ypi=k2)exp(πpi,k1)exp(πpik)=P(ypi=1)c=2kexp(πpic)=P(ypi=1)exp(c=2kπpic)(8.26) \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} \qquad(8.26) と書き下します。 ここではまだP(ypi=1)P(y_{pi}=1)が残っていますが,GPCMの考え方ではこれ以上計算できません13。 改めて(8.26)式を見ると,すべてのカテゴリの反応確率は「P(ypi=1)P(y_{pi}=1)の何倍」という形になっていることがわかります。 そこで一旦πpi1=0\pi_{pi1}=0として得られるP(ypi=1)=exp(πpi1)=1P(y_{pi}=1)=\exp(\pi_{pi1})=1を(8.26)式に代入して, P(ypi=k)=exp(c=1kπpic)(8.27) P(y_{pi}=k) = \exp\left(\sum_{c=1}^{k}\pi_{pic}\right) \qquad(8.27) とまとめてしまいます。 これは,各カテゴリの反応確率を「一番下のカテゴリが選択される確率の何倍か」という形で表すことを意味しています。 これによって,ひとまずP(ypi=k)P(y_{pi}=k)が定式化されましたが,P(ypi=1)=1P(y_{pi}=1)=1と設定している以上P(ypi=k)P(y_{pi}=k)の全カテゴリの総和が1になるかわかりません。

ということで最後に,P(ypi=1)P(y_{pi}=1)の値が何であっても全カテゴリでの総和が1になるようにするため,総和で割ることを考えます。 実際に P(ypi=k)l=1KP(ypi=l)=P(ypi=k)P(ypi=1)+P(ypi=2)++P(ypi=K)=exp(c=1kπpic)exp(l=11πpil)+exp(l=12πpil)++exp(l=1Kπpil)=exp(c=1kπpic)l=1Kexp(c=1lπpil)(8.28) \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} \qquad(8.28) とすると,分子も分母もすべての項に共通のexp(πpi1)\exp(\pi_{pi1})がかかっているため,この項の値が何であっても全体の割合に占める各カテゴリの反応確率の割合は変わらずに済みます。 ということで,ようやくGPCMの式が導出できました。 図 8.16 に,ここまでの考え方をまとめてみました。 特定のカテゴリの反応確率P(ypi=k)P(y_{pi}=k)が直接計算できないため,いったん一番下のカテゴリの反応確率P(ypi=1)=1P(y_{pi}=1)=1とおいて,他のカテゴリは「前のカテゴリの何倍」という形で表します。 最後に全体の総和が1になるように調整をしたものがGPCMにおける反応確率というわけです。

図 8.16: GPCMのイメージ

図 8.17 は,GPCMでの項目特性曲線をプロットしたものです。 GPCMでは,カテゴリ困難度βik\beta_{ik}は「カテゴリk1k-1kkの二択」における選択確率がちょうど50%になるポイントを表しています。 したがって,隣接するカテゴリの反応確率はちょうどβik\beta_{ik}のところで大小関係が入れ替わります。 βik\beta_{ik}はそのように解釈できるわけです。

図 8.17: GPCMの項目特性曲線

GRMでは,カテゴリ困難度は必ず単調増加になっている必要がありますが,GPCMでは単調増加で無くても問題ありません。 図 8.18 にカテゴリ困難度が単調増加にならないケースをプロットしてみました。 この場合,「カテゴリ2がカテゴリ1を上回る」よりも先に「カテゴリ3がカテゴリ2を上回る」ため,結果的にカテゴリ2が最大になることがありません。

図 8.18: カテゴリ困難度が単調増加にならないケース

8.6.3 Rで多値型モデル

それでは,mirtで2つの多値型モデルを実行してみましょう。 …といっても引数itemtypegraded(GRM)かgpcmを指定したら良いだけなので難しいことはありません。

コード 8.10: 多値型モデル (GRM)
# 多値型なのでdatをそのまま突っ込む
result_grm <- mirt(dat[, cols], model = "f_1 = Q1_1-Q1_10", itemtype = "graded")
コード 8.11: 係数をチェック
coef(result_grm, IRTpars = TRUE, simplify = TRUE)
$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   1

datの各項目は6件法なので,カテゴリ困難度bが5つ出力されました。

二値型モデルの時と同じように,項目特性曲線を出力することももちろん可能です。

コード 8.12: GRMで項目特性曲線
# カテゴリ数のぶんだけ線を引いてくれる
plot(result_grm, type = "trace")
図 8.19: GRMで項目特性曲線

また,itemplot(., type = "threshold")関数を使うと, 図 8.14 の下半分のような「カテゴリkk以上の反応をする確率」のプロットを出すこともできます。

コード 8.13: GRMで項目特性曲線2
# itemplot()は1項目しかプロットできない
# 何番目の項目をプロットするかを引数itemで指定
itemplot(result_grm, item = 6, type = "threshold")
図 8.20: GRMで項目特性曲線2

同じようにGPCMでもやってみましょう。 やり方も結果の見方も同じなので詳しい説明は省略します。

コード 8.14: 多値型モデル (GPCM)
# 多値型なのでdatをそのまま突っ込む
result_gpcm <- mirt(dat[, cols], model = "f_1 = Q1_1-Q1_10", itemtype = "gpcm")
コード 8.15: 係数をチェック
coef(result_gpcm, IRTpars = TRUE, simplify = TRUE)
$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.16: GPCMで項目特性曲線
# カテゴリ数のぶんだけ線を引いてくれる
plot(result_gpcm, type = "trace")
図 8.21: GPCMで項目特性曲線

8.7 多次元モデル

これまでのIRTモデルでは一次元性の仮定に基づく一次元モデルのみを扱ってきました。 ですが多くの心理尺度やテストでは,複数の次元の特性を同時に測定したいというニーズがあります。 先程紹介したようにIRTはカテゴリカル因子分析と数学的には同値なので,IRTでも多次元モデルを考えていきましょう。

多次元因子分析モデルでは,項目反応を ypi=τi+t=1Tfptbtiεpi(8.29) y_{pi} = \tau_i + \sum_{t=1}^{T}{f_{pt} b_{ti}} - \varepsilon_{pi} \qquad(8.29) と表しました。これと同じように考えるならば,多次元IRTモデル(2パラメータロジスティックモデル)での項目反応関数は P(ypi=1)=11+exp[t=1Tαtiθptτi](8.30) P(y_{pi}=1) = \frac{1}{1+\exp\left[-\sum_{t=1}^{T}\alpha_{ti}\theta_{pt} - \tau_i\right]} \qquad(8.30) と表せそうです。

(8.30)式のモデルにおける項目パラメータは,識別力が次元の数(TT個)と,切片が一つです。 また,項目の全体的な識別力を評価する際には,次元ごとの識別力を一つにまとめた多次元識別力(multidimensional discrimination: 豊田,2013 MDISCi=t=1Tαti2(8.31) \mathrm{MDISC}_i=\sqrt{\sum_{t=1}^{T}\alpha_{ti}^2} \qquad(8.31) を用いることもあります。これを用いると,切片パラメータτi\tau_iを困難度に変換することもできるようになります: βi=τiMDISCi.(8.32) \beta_i=\frac{-\tau_i}{\mathrm{MDISC}_i}. \qquad(8.32)

二次元までならば,ICCと同じように𝛉p\symbf{\theta}_pP(ypi=1)P(y_{pi}=1)の関係をプロットすることも可能です。 図 8.22 の左は,とある二次元項目における項目反応曲面 (item response surface [IRS])です。 X軸,Y軸がそれぞれのθt\theta_tの値を表しており,Z軸の値が反応確率P(ypi=1)P(y_{pi}=1)を表しています。 この項目では,([αi1,αi2],τi)=([1.5,0.5],0.7)([\alpha_{i1},\alpha_{i2}],\tau_i)=([1.5,0.5],0.7)となっていることから,θp2\theta_{p2}の値はさほど反応確率に影響を及ぼしていないことが図からもわかります。 なお 図 8.22 の右は,2つのθpt\theta_{pt}の値に対してP(ypi=1)P(y_{pi}=1)が等しくなるラインを表しています。

図 8.22: 項目反応曲面Reckase, 2019, fig. 12.1

8.7.1 補償モデルと非補償モデル

ここまで紹介した多次元モデルでは 図 8.22 からも分かるように,反応確率に対してθpt\theta_{pt}が及ぼす影響は加法的といえます。 例えば([αi1,αi2],τi)=([1.5,1],0)([\alpha_{i1},\alpha_{i2}],\tau_i)=([1.5,1],0)の項目の場合,𝛉p=[11]\symbf{\theta}_p=\begin{bmatrix}-1 & 1\end{bmatrix}の人と𝛉p=[00.5]\symbf{\theta}_p=\begin{bmatrix}0 & -0.5\end{bmatrix}の人では反応確率は同じになります。 また,ある次元の特性値がとても低い場合でも,別の次元の特性値がとても高い場合には反応確率が高くなる可能性があります。 このような次元間の関係を反映した先程の多次元IRTモデルは補償型 (compensatory) モデルと呼ばれます。

一方で,複数の次元の要素がすべて揃って初めて反応確率が高くなるようなケースも考えられます。 例えば数学のテストが日本語ではなく英語で実施される場合,これは「英語」と「数学」の2つの能力が試される二次元構造ではありますが,補償型モデルのように「英語が苦手でも数学力が高ければ解答できる」とは考えにくいです。 英語で書かれた問題文を理解するだけの英語力と,計算を実行できるだけの数学力の両方が揃って初めて正答できるでしょう。 このような状況では,一方の特性値の低さをもう一方の特性値の高さでカバーすることはできません。 こうした次元間の関係を反映したモデルを非補償 (noncompensatory) モデルあるいは部分補償 (partially compensatory) モデルと呼びます14

最もシンプルな非補償モデルの項目反応関数は P(ypi=1)=t=1T11+exp[αit(θptβit)](8.33) P(y_{pi}=1) = \prod_{t=1}^{T}\frac{1}{1+\exp\left[-\alpha_{it}(\theta_{pt} - \beta_{it})\right]} \qquad(8.33) と表されます。 各次元に関して独立に算出された「その次元の閾値を超える確率」を計算し,最後にすべての積をとっています。 そのため,一つでも特性値の低い次元があれば反応確率は一気に小さくなってしまう,というわけです。 非補償モデルの項目反応曲面は 図 8.23 の左のようになります。 図 8.22 と比べると,両方の次元のθp\theta_{p}が高くないとP(ypi=1)P(y_{pi}=1)が高くならないことが見て取れます。

図 8.23: 非補償モデルの項目反応曲面Reckase, 2019, fig. 12.2

補償モデルと非補償モデルを比べると,補償モデルのほうがよく用いられているような印象です。 その理由の一つとして,非補償モデルのほうがパラメータの推定が難しいという点が挙げられると思います。 上の式を見るとわかるように,非補償モデルでは困難度パラメータβit\beta_{it}が次元の数だけあるためそもそものパラメータ数が補償型よりも多くなっています。 また,非補償モデルは項目反応関数の全体の形が積になっています。 コンピュータは足し算よりも掛け算のほうが苦手なので,そういった意味でも計算の難易度が高いのです。 また,心理尺度的には補償型のほうが従来の因子分析モデルと同じ形になっているため,扱いやすいというのもあると思います。 とはいえ,多次元の考え方として非補償型もあるんだということを頭の片隅に置いておくと,いつかどこかで何かの参考になるかもしれません15

8.7.2 Rでやってみる

ここからは補償型の多次元IRTモデルをmirtで実行する方法を紹介します16。 補償型モデルはカテゴリカル因子分析と同値なので,因子分析と同じように「探索的モデル」と「検証的モデル」の2種類があります。 ただ探索的モデルでは因子の回転に応じてθpt\theta_{pt}の意味が変化してしまうため実用上はなかなか使いにくいような気もします。

探索的補償型モデル

探索的モデルを実行する場合には,引数modelに「次元数」を与えるだけです。 itemtypeはこれまでと同じように自由に選んでください。二値だったら2PLとか,多値だったらgradedとかgpcmとか。

コード 8.17: 探索的補償型モデル
# GRMでやってみましょう
# このあたりまで来ると推定にもそこそこ時間がかかるかもしれません。
result_expl_comp <- mirt(dat[, cols], model = 2, itemtype = "graded")
コード 8.18: 推定値のチェック
# 多次元の場合引数IRTparsは使えないので注意
coef(result_expl_comp, simplify = T)
$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には,それぞれθpt\theta_{pt}の平均と共分散行列が出ています。 これを見る限り,まだ因子間相関が0なので回転前の因子負荷(識別力)が出ているようです。 さらに詳しく見ると,最後の項目のa20.000となっています。 実はmirtでの探索的モデルの推定の際には,このように一部の因子負荷を0に固定した状態で推定を行っています17。 「探索的」と紹介していましたが,実際には「モデルが識別できる最小限の制約を置いて推定する」ということを行っているのです。 IRTに限らずCFAもそうなのですが,複数の因子がある場合に識別性を持たせるには,最低でもT(T1)2\frac{T(T-1)}{2}個の因子負荷を固定する必要があることが知られています。 そのため,2因子の場合は第2因子の因子負荷を一つ0にしているわけです。 同様に,3因子の場合にはさらに第3因子の因子負荷を2つ0にする必要があり,4因子の場合にはさらにさらに第4因子の因子負荷を3つ0にする必要があり…という要領です。

このままでは,いくら探索的とは言え各次元の解釈が出来ません。 mirt()で得られるパラメータの推定値はいわば「初期解」なので,ここからさらに回転させてあげれば良さそうです。 回転させる場合は,psych::fa()の時と同じように引数rotateを指定してあげます。

コード 8.19: 識別力の回転
# psych::fa()のデフォルトであるoblimin回転をしてみる
coef(result_expl_comp, rotate = "oblimin", simplify = TRUE)

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に各次元と項目の関係を記述してあげる必要があります。

コード 8.20: 検証的補償型モデル
# model式の右辺は列番号だけでもOKです
model <- "
f_1 = 1-5
f_2 = 6-10
"

result_conf_comp <- mirt(dat[, cols], model = model, itemtype = "graded")
coef(result_conf_comp, simplify = TRUE)
$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()によって出力可能です。

コード 8.21: 多次元でも項目反応カテゴリ特性曲線
# もちろん非補償型でも探索的でも可能です
# 多値型項目だとすべてのカテゴリのIRCCCが重なるのでめっちゃ見にくい
# 全項目やると重いので,which.itemで少しずつやるのがおすすめ
plot(result_conf_comp, type = "trace", which.item = 4)
図 8.24: 項目反応カテゴリ特性曲線

また,引数rotをうまく指定すると,プロットを回すことができます。 うまく回せば,多少は見やすくなるかもしれません。 例えば確認的モデルでは,各項目の識別力を一つだけ非ゼロにしていました。 そのため,識別力がゼロである次元の方向は潰すように回転させると良いかもしれません(見やすくなるとは言っていない)。

コード 8.22: 多次元でもitemplot()
# 項目4はa2=0だったので,a1方向に見るとIRCCCが見やすいかも?
plot(result_conf_comp, type = "trace", which.item = 4, rot = list(xaxis = -90, yaxis = 0, zaxis = 0))
図 8.25: グラフを回してみる

8.8 次元性の検証

Section 8.5.2 で少し触れたように,IRTでは基本的にそのテストが測定しようとする次元数は理論的に決まっていることが多く,そのため「何因子が最適化」を探索的に決める方法よりは「nn因子で問題ないか」を検証する方法が採用されることが多いと思います。 ここでは,そのようなIRTの文脈で用いられることの多い(と思われる)いくつかの次元性検証法を紹介していきます。

8.8.1 パラメトリックな方法

パラメトリックな方法では,(多次元・一次元の)IRTや因子分析モデルに基づいてパラメータの推定を行い,その結果を用いて次元性を評価していきます。 例えば確認的因子分析によって推定を行い,SEM的なモデル適合度をチェックする方法はパラメトリックな検証法の一つと言えるでしょう。 その他にも,IRTの文脈で用いられている検証法がいくつか存在します。 代表的と思われる方法の一つがNOHARMGessaroli & De Champlain, 1996です18。 この方法では,項目の各ペアでの残差共分散を,推定されたパラメータおよび観測されたデータ(各項目の正答率)から計算していきます。 そして,全ペアに対して計算された残差共分散(をzz変換したもの)の二乗和(を何倍かしたもの)がχ2\chi^2分布に従うことを利用して,「残差共分散が全てゼロである」という帰無仮説に対する検定を行います。 そして帰無仮説が棄却されなかった場合,「その次元数で全ての項目ペア間の相関関係を表現できていないとは言えない」ということで,その次元数でOK,という判断がくだされます。 考え方としてはBartlett検定に近いかもしれません。

8.8.2 ノンパラメトリックな方法

IRTの文脈で提案されたノンパラメトリックな次元性評価法の中でも代表的なものが,ここで紹介するDIMTESTStout et al., 1996やDETECTZhang & Stout, 1999です。 これらの方法では,「その次元(共通因子)数があれば局所独立性が満たされるか」という視点から評価を行います。 Section 8.5.1 では,局所独立性を「どの項目ペア間でも,どの特性値の値でも連関がない」ことと表現していましたが,これを数式で表すと Cov(𝐲i,𝐲j|θ)=0foralli,j,θ(8.34) \mathrm{Cov}(\symbf{y}_i, \symbf{y}_j\vert\theta) = 0 ~ \mathrm{for~all}~i,j,\theta \qquad(8.34) となります19。したがって,θ\thetaで条件づけた各項目ペアの(条件付き)共分散を計算していけば,局所独立性を満たすために必要な次元数を考えることができるわけです。

具体的な方法の説明に移る前に,次元性と条件付き共分散のイメージを確認しておきます。 図 8.26 には,2因子モデルにおいて,なるべく単純構造になるように回転した後の各項目の因子負荷をベクトルとして表したものです。また𝚯C1,𝚯C2\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}と表記されているベクトルは,各クラスタ内の項目の因子負荷の平均を取ったようなものを表しています。同様に𝚯TT\symbf{\Theta}_{TT}と表記されているベクトルは,クラスタをまたいだ全項目のベクトルの平均を取ったもので,端的に言えばこれが「テスト全体で測定しようとしている特性のベクトル」というわけです。 よくある例としては,𝚯TT\symbf{\Theta}_{TT}が「数学の学力」,𝚯C1,𝚯C2\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}がそれぞれ「代数」「幾何」といった感じです。 図 8.26 では,𝚯C1,𝚯C2\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}は完全に無相関というわけではないですが,さほど強くもないので,一次元性が満たされているとは言いづらい感じがします。

図 8.26: 2次元のテストの因子負荷プロットStout et al., 1996, fig. 2

このテストが一次元性を満たしているかを確認するためには,「𝚯TT\symbf{\Theta}_{TT}で条件付けられた」共分散を考えることになります。 そして条件付き共分散は, 図 8.26 においてはベクトル𝚯TT\symbf{\Theta}_{TT}に向かって伸ばした足の長さの積のような形で表されます。 これを理解するために, 図 8.26 に少し描き足してみました。 図 8.27 に追加されたオレンジの線は,𝚯TT=𝚯C1+𝚯C22\symbf{\Theta}_{TT}=\frac{\symbf{\Theta}_{C1}+\symbf{\Theta}_{C2}}{2}となる(次元間で重みが同じ)場合に𝚯TT\symbf{\Theta}_{TT}の値が同じになる𝚯C1,𝚯C2\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}の組みあわせを表しています。

図 8.27: 条件付き共分散Stout et al., 1996, fig. 2 に少し加筆)

この図をもとに,項目ペアに対して計算できる条件付き共分散を考えてみましょう。 まず,クラスターC1C1の中の2つの項目の条件付き共分散は正の値になることが期待されます。 𝚯TT\symbf{\Theta}_{TT}の値が同じ2人(AさんとBさん)について考えてみると,AさんはクラスターC1C1の項目に正解するのに必要な能力(𝚯2\symbf{\Theta}_2)の値が高いためクラスターC1C1の項目には正解する確率が高いと予想される一方で,Bさんは𝚯2\symbf{\Theta}_2の値が低いので,正解する確率は低そうです。 これをまとめると,「𝚯TT\symbf{\Theta}_{TT}で条件づけたとき,クラスターC1C1内のある項目に正解する人ほど,同じクラスター内の別の項目にも正解する確率が高い」わけなので,条件付き共分散は正の値になります。

同様に,クラスターC1C1の中のある項目と,別のクラスターC2C2の中のある項目の条件付き共分散を考えてみると,AさんはクラスターC1C1の項目には正解するがクラスターC2C2の項目には正解できないと予想されます。 そしてBさんは反対に,クラスターC2C2の項目にのみ正解できる可能性が高いでしょう。 したがって,「𝚯TT\symbf{\Theta}_{TT}で条件づけたとき,クラスターC1C1内のある項目に正解する人ほど,別のクラスター内のある項目に正解する確率は低い」ことになり,条件付き共分散は負の値になります。

以上をまとめると, 図 8.27 における条件付き共分散の符号は,2つの項目の因子負荷ベクトルが「ベクトル𝚯TT\symbf{\Theta}_{TT}から見て」同じ方向にあるかどうかに対応するわけです。 そして,仮に全ての項目が近い方向のベクトルで表される場合には,当然𝚯TT\symbf{\Theta}_{TT}ベクトルも同じような方向になるため,条件付き共分散も0に近くなっているはずです。 という感じで,(8.34)式が局所独立性,そして一次元性を表しています。

DIMTEST

DIMTESTStout et al., 1996では,まずテストを「(次元性の)評価用 (assessment subtest: AT)」と「データの分割用 (partitioning subtest: PT)」の2つに分割します。 この分割では,理論的に想定される分かれ方( 図 8.26𝚯C1,𝚯C2\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}に相当)でわけたり,探索的因子分析をもとに分けたりします。 つまり,一旦「このテストが分かれる(2次元だ)としたらこうなるだろう」という形で分けてみるのです。 もしこのテストが実際には一次元だとしたら,そのように分けたとしても結局全ての項目のベクトルの向きは近いので,片方のテストの得点によって条件付けたとしても,もう一方のテストの中の共分散は0に近いはずです。 一方で,テストが(𝚯C1,𝚯C2\symbf{\Theta}_{C1},\symbf{\Theta}_{C2}のように)分かれる場合には,例えば𝚯C1\symbf{\Theta}_{C1}で条件づけた際の𝚯C2\symbf{\Theta}_{C2}内の条件付き共分散は全て正になってしまいます。

以上の考え方に基づいて,DIMTESTでは以下の統計量を計算します。 T=k(i,jAT)Cov̂(𝐲i,𝐲j|SPT=k).(8.35) T = \sum_k \sum_{(i,j \in AT)} \widehat{\mathrm{Cov}}(\symbf{y}_{i},\symbf{y}_{j}|S_{PT}=k). \qquad(8.35) ここでSPTS_{PT}はPTにおける正答数を表しています。 したがって,(8.35)式は「PTにおける正答数で(=𝚯PT\symbf{\Theta}_{PT}で)条件づけたときのAT内の共分散」を求めているわけです。 統計量TTが0に近いほど局所独立が満たされている(そのテストは一次元である)と言えるので,あとはこれを標準化した値を用いて標準正規分布に基づくzz検定 Stout, 1987 を行います。

DIMTESTの成功のカギは,ATとPTの分割にあるわけですが,どのように分けたら良いかについては Stout et al. (1996 などを参照してください。

DETECT

DETECTも条件付き共分散をベースにした方法です。 以下では,二値データに対して考案されたものZhang & Stout, 1999を紹介します。 多値データの場合にはこれを発展させた方法Zhang, 2007が提案されています。

DETECTでも,まずはテストをいくつかのクラスター(C1,C2,,CkC1, C2, \cdots, Ck)に分割します。 そして以下のようにテスト全体で推定された特性値𝚯TT\symbf{\Theta}_{TT}によって条件づけた共分散の統計量を計算します。 D=2I(I1)1ijIδi,jE[Cov(𝐲i,𝐲j|𝚯TT)].(8.36) 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})]. \qquad(8.36) ここで重要なのは指示変数δi,j\delta_{i,j}です。 これは,

  • 2つの項目i,ji,jが同じクラスターにある場合には1
  • 2つの項目i,ji,jが異なるクラスターにある場合には-1

になります。

DETECTが計算しようとしているものを 図 8.26 をもとに考えてみます。 図 8.26 は,あるテストの真のクラスターを表しているとします。 つまりあるテストはきれいに2つのクラスターに分かれている,ということです。 このときに,この真のクラスターと一致するようにテストを分けた上でDETECTの統計量DDを計算してみましょう。

  • 2つの項目i,ji,jが同じクラスターにある場合には,条件付き共分散(の期待値)E[Cov(𝐲i,𝐲j|𝚯TT)]E[\mathrm{Cov}(\symbf{y}_i,\symbf{y}_j|\symbf{\Theta}_{TT})]δi,j\delta_{i,j}も正の値
  • 2つの項目i,ji,jが異なるクラスターにある場合には,条件付き共分散(の期待値)E[Cov(𝐲i,𝐲j|𝚯TT)]E[\mathrm{Cov}(\symbf{y}_i,\symbf{y}_j|\symbf{\Theta}_{TT})]δi,j\delta_{i,j}も負の値なので,積は正の値になる

ということで,全ての項目ペア間で計算される値(\sumの中身)が全て正の値になります。 このように,統計量DD(クラスター数も含めて)最も適切なクラスター基づいて計算したときに最大値を取る値です。 そして,そもそも全ての項目のベクトルの方向が近ければ(次元性が満たされていれば),条件付き共分散(の期待値)は全て0に近いので,その平均である統計量DDも0に近い値になっているはずです。

ということで,この統計量を用いて一次元性を評価する場合には,「統計量DDが最大でいくつになるか」を見たら良いと言えます。 例えば「理論的に分かれるとしたらこう」で分けてみたり,探索的因子分析の結果をもとに分けてみたり,あるいはとにかく手当たり次第にランダムに様々な分け方をひたすら試してみたりして,統計量DDが高々どれくらいの値に収まるかを確認しましょう。 一次元性の検証という意味では,もちろんDDの値は小さいほどよいのですが,二値データの場合例えば Roussos & Ozbek (2006 では最大値が0.2以下なら一次元とみなしてよいのではないか,と言われていたりします。

そんなDETECTについては,Rで計算するための関数としてsirtパッケージにconf.detect()およびexpl.detect()というものが用意されています20。 これらは二値・多値のどちらにも対応していますが,ここでは二値データへの適用例を見てみます。

まずはconf.detect()です。confは”confirmatory”のことで,その名の通り分析者が指定したクラスター構造のもとでDETECTを計算します。 そのため,事前に想定されるクラスターがある(数学のテストの中に複数の単元がある,心理尺度の中に複数の下位概念がある)ような場合に使える方法です。

コード 8.23: DETECTで次元性の検証
# install.packages("sirt")
library(sirt)

# テスト全体で測定している特性値
# シンプルに和得点などを使用しても良い
fs <- fscores(result_2plm)
# 各項目がどのクラスタに属するか
# 今回は元のデータの想定に合わせて設定しておきます
cluster <- rep(1:2, each = 5)
res_detect <- conf.detect(dat_bin, score = fs, itemcluster = cluster)
# conf.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)式で計算されるDDの値です。この値が0.2より小さければ,そのテストは一次元であるとみなせそうなのですが…今回使用した10項目は元々2次元だったため,「一次元ではない」と正しく判定されています(一般に1を超えると「かなり強く多次元である」と言えるようです)。

続いて,本来一次元であることが想定されるケースとして,項目間の相関が高い(=内的一貫性の高い)Q1_6からQ1_10の5項目のみを使ってconf.detect()してみましょう。

コード 8.24: 一次元性が示されそうな例
# 今度は和得点でやってみます
sum_s <- apply(dat_bin[, 6:10], 1, sum)
# 各項目がどのクラスタに属するか
# 今回は適当に奇数・偶数で分けてみます
cluster <- c(1, 2, 1, 2, 1)
res_detect <- conf.detect(dat_bin[, 6:10], score = sum_s, itemcluster = cluster)
-----------------------------------------------------------
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.041

DETECT Unweightedの値は先程よりは小さくなっていますが,それでも0.2よりは大きくなっています…理由はよく分からないのですが,もしかしたら項目数が少なすぎる場合には小さくなりにくいのかもしれません(要出典)。

ここまでは,引数itemclusterによって項目のクラスタリングを与えていました。 一次元性を検証するためには,本当はこれをすべてのパターン(II項目をKK個のクラスタに分けるとしたらおよそKIK^I通り)で計算して,その最大値を出す必要があります。 ですがこれはどう考えても非効率な(というかすぐ無理になる)話です。 そこで,各クラスタ数における「DDが最大になる分割」を探索的に求めることを考えましょう。

expl.detect()では,階層的クラスター分析を用いて,指定したクラスタ数の分割の中でDDが最大になるものを探索的に決定してくれます。 考え方はシンプルで,条件付き共分散(の期待値)E[Cov(𝐲i,𝐲j|𝚯TT)]E[\mathrm{Cov}(\symbf{y}_i,\symbf{y}_j\vert\symbf{\Theta}_{TT})]が大きいペアは近くにある=同じクラスターに属すると考えて,距離の近い項目から順に,最終的にKK個のクラスターにまとまるまでくっつけていき,そのクラスター構造の元で計算したDDを出力してくれます。 ということで早速試してみましょう。

コード 8.25: 探索的DETECT
fs <- fscores(result_2plm)
# nclustersは想定される最大クラスター数
# 想定よりやや大きめの値にしておくのが無難
res_detect_expl <- expl.detect(dat_bin, score = fs, nclusters = 5)
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
図 8.28: DDが最大になるクラスターのデンドログラム

結果を見ると,このデータ(Q1_1からQ1_10)に対してはクラスター数4のときが最大値になるようです。 図 8.28 にその時の分かれ方が示されています。 (やっぱり項目数が少ないとクラスター数が多めになるのか?)

一点注意として,expl.detect()の内部では,デフォルトではデータを半分に分割して,

  1. 前半のデータで条件付き共分散を計算→クラスター構造を決定
  2. 後半のデータでDDを計算

という流れで計算を行っています。これは多分オーバーフィッティングを避けるためだと思います。 そのため,N.Clusterが2のところのDETECT.est(手順1で使用したデータで計算されたDD)もDETECTR.val(手順2で使用したデータで計算されたDD)も,conf.detect()での結果1.064とは異なる値になっています(N.estおよびN.valがいずれも元のdatの行数の半分になっている)。

…ということで,今回のデータではきれいな一次元性を確認することはできなかったのですが,このような方法を用いると,「一次元とみなして分析しても問題ないか」を評価することができるのです。

8.9 適合度

IRTでもSEMと同じように適合度を考えることができます。 特に「個人」と「項目」の両方に関心があることの多いIRTでは,それぞれに対しての適合度を考える必要があります。 もちろんこれから紹介する適合度の考え方は,そのまま因子分析やSEMにも応用しようと思えばできるので,そういった意味でも参考にしてみてください。

8.9.1 局所独立性の確認

まずはじめに,データとモデルの適合度の一つとして,データが局所独立の仮定に適合しているかをチェックします。 Section 8.5.1 で説明したように,局所独立性が満たされている場合には,θp\theta_pで条件づけたときの2項目の回答の間に相関が無いはずです。

X2X^2統計量

X2X^2統計量 Chen & Thissen, 1997 では,カテゴリ変数のクロス表に対して用いられるχ2\chi^2検定の枠組みを利用して局所独立性を評価します。 ある二値項目iiの項目パラメータ(αi,βi)(\alpha_i,\beta_i)が分かっている場合,その項目の母集団全体での期待正答率は E[P(ypi=1)]=P(ypi=1|θ)f(θ)dθ(8.37) \begin{aligned} E[P(y_{pi}=1)] = \int_{-\infty}^{\infty}P(y_{pi}=1|\theta)f(\theta)d\theta \end{aligned} \qquad(8.37) という形で求めることができます。 ここでf(θ)f(\theta)は特性値の確率分布(ふつうは標準正規分布)におけるx=θx=\thetaでの確率密度を表しています。 同様に,その項目の誤答率の期待値を求めるならばP(ypi=1|θ)P(y_{pi}=1|\theta)の部分を1P(ypi=1|θ)1-P(y_{pi}=1|\theta)にしたら良いだけです。

この考え方を2項目に広げると,母集団全体において項目(i,j)(i,j)に両方とも正解する確率はもし2項目が独立ならば E[P(ypi=1ypj=1)]=P(ypi=1|θ)P(ypj=1|θ)f(θ)dθ(8.38) \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} \qquad(8.38) と表せます。 ということは,サンプルサイズがNNの場合,項目(i,j)(i,j)に両方とも正解する人数の期待値E11E_{11}N×P(ypi=1ypj=1)N\times P(y_{pi}=1 \land y_{pj}=1)です。 同様にして,一方の項目にだけ正解する人数および両方の項目に誤答する人数の期待値も求めると, 表 8.3 のようになります。

表 8.3: 期待度数(独立な2項目の分布)
項目jj
項目ii 当てはまらない 当てはまる
当てはまらない E00=N×P(ypi=0ypj=0)\begin{aligned}E_{00}= \\ N\times P(y_{pi}=0 \land y_{pj}=0)\end{aligned} E01=N×P(ypi=0ypj=1)\begin{aligned}E_{01}= \\ N\times P(y_{pi}=0 \land y_{pj}=1)\end{aligned}
当てはまる E10=N×P(ypi=1ypj=0)\begin{aligned}E_{10}= \\ N\times P(y_{pi}=1 \land y_{pj}=0)\end{aligned} E11=N×P(ypi=1ypj=1)\begin{aligned}E_{11}= \\ N\times P(y_{pi}=1 \land y_{pj}=1)\end{aligned}
NN

この期待度数分布はもし2項目が局所独立ならばこうなるだろうという状態を表しており, 表 8.1 に相当するものをθ\thetaについて周辺化して作成したものといえます。 つまり局所独立ではなくなるほど,実際のデータでのクロス表は( 表 8.2 のように)この期待度数から離れるだろう,ということです。

クロス表での連関の検定は(学部の統計でやっているかもしれない)χ2\chi^2検定です。 その検定統計量は χ2=s=01t=01(OstEst)2Est(8.39) \chi^2=\sum_{s=0}^{1}\sum_{t=0}^{1}\frac{(O_{st}-E_{st})^2}{E_{st}} \qquad(8.39) で計算されます。ここでOstO_{st}は,実際のデータでypi=sypj=ty_{pi}=s \land y_{pj}=tだった人数を表しています。

こうして計算されたχ2\chi^2統計量は,自由が(行数-1) ×\times (列数-1)のχ2\chi^2分布に従うことが知られています。 二値データの場合は1×1=11\times1=1ということです。 そしてここまでのプロセスはそのまま多値項目に対しても拡張可能です。 例えばdatの項目のように6件法どうしの場合は,自由度は5×5=255\times5=25となるわけです。 ということで,項目のペア毎にχ2\chi^2統計量を計算して統計的仮説検定を行い,有意だったペアの間は局所独立ではない可能性が疑われる,ということになります。

Q3Q_3統計量

Q3Q_3統計量 Yen, 1984 では,回帰分析でいうところの偏相関係数にあたるものを利用します。 偏相関係数は「変数zzの影響を取り除いた時のxxyyの相関係数」で,疑似相関の影響を検討する際などに利用されるものです。 これをIRTの文脈に当てはめると「変数θp\theta_pの影響を取り除いた時のypiy_{pi}ypjy_{pj}の相関係数」ということで,これが0であれば局所独立だろうと言えそうです。

回帰分析では,xxyyをそれぞれzzで回帰した時の予測値x̂,ŷ\hat{x},\hat{y}に対して,残差の相関r(xx̂,yŷ)r_{(x-\hat{x},y-\hat{y})}のことを偏相関係数と呼んでいました。 IRT(ロジスティックモデル)もその中身はロジスティック回帰なので,同じようにしてypiy_{pi}ypjy_{pj}をそれぞれθp\theta_pで予測した時の予測値(期待正答率)を P̂(ypi=1)=11+exp[α̂i(θ̂pβ̂i)](8.40) \hat{P}(y_{pi}=1)=\frac{1}{1+\exp\left[-\hat{\alpha}_{i}(\hat{\theta}_p - \hat{\beta}_i)\right]} \qquad(8.40) というように,αi,βi,θp\alpha_i,\beta_i,\theta_pにそれぞれ推定値α̂i,β̂i,θ̂p\hat{\alpha}_i,\hat{\beta}_i,\hat{\theta}_pを入れることで計算が出来ます。 すると,残差をdpi=ypiP̂(ypi=1)d_{pi}=y_{pi}-\hat{P}(y_{pi}=1)と求めることができます。 図 8.29dpid_{pi}を表してみました。

図 8.29: 残差

この項目にはAさんもBさんも正解しているため,ふたりともyy軸の値は1です。 ですが他の全項目から推定された2人のθ̂p\hat{\theta}_pはそれぞれ2,2-2,2だったとします。 この場合,Aさんは他の項目の出来から考えるとこの問題にはほぼ間違えるはずなのに正解しています。 ということでdpid_{pi}の値が大きくなっています。 このようになる理由としては,やはりθp\theta_p以外の別の要因によって正解できた,と考えるのが妥当でしょう。 これを全体に広げてみた時に,dpid_{pi}が高い人ほど別の項目でも同様にdpjd_{pj}の値が大きくなっているとしたら,この2つの項目には何か正解に寄与するθp\theta_p以外の別の要因がある,と考えることができるわけです。

ということでQ3Q_3統計量は Q3=r𝐝i,𝐝j(8.41) Q_3=r_{\symbf{d}_i,\symbf{d}_j} \qquad(8.41) という形になります。 Q3Q_3統計量は相関係数なので,サンプルサイズに依存しない効果量の指標として見ることができ,絶対値が0.2を超えてくると怪しいと判断できるようです。

Rでやってみる

mirtではresidualsという関数でいま紹介した指標を出してくれます。 引数typeで,出してほしい統計量を指定できます。 ちょっと厄介なのですが,χ2\chi^2統計量を出してほしい場合にはtype = "LD"と指定してください。 LDはlocal dependenceの頭文字ですが,そういったらQ3Q_3だってLDだろ,と思われるかもしれません。 これは Chen & Thissen (1997 の書き方に則っているのだと思います。我慢してください。

コード 8.26: χ2\chi^2統計量を出す
# 引数df.p = TRUEとすると,自由度とp値を出してくれる
residuals(result_2plm, type = "LD", df.p = TRUE)
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

上の行列の上三角がpp値です。したがってこれが0.05を下回っている場合は「局所独立ではない」という判断になるのですが,見た感じかなり多くのペアで有意になっています。 これは,統計的仮説検定なのでサンプルサイズの影響を受けているためです。 また上の行列の下三角は自由度です。 今回は二値に変換したデータに対して行っているので自由度は1になっています。 もしも6件法のまま(GRMやGPCMで)推定した結果に対して出した場合には,自由度は先程説明した通り5×5=255\times5=25となります。

そして下の行列は実際に計算されたχ2\chi^2統計量(下三角)および標準化した値(上三角;たぶんクラメールの連関係数)です。 ということで,この行列の上三角について大きいペアが,局所独立からより強く離れているといえます。

続いてQ3Q_3統計量です。

コード 8.27: Q3Q_3統計量を出す
# 仮説検定は無いので引数df.pは書いても無視される
residuals(result_2plm, type = "Q3")
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.28: 表示する閾値を設定
residuals(result_2plm, type = "Q3", suppress = 0.2)
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                                                      
Q1_2        1.000                                   -0.209      
Q1_3               1.000       0.244                -0.235      
Q1_4                        1                                   
Q1_5               0.244       1.000                -0.226      
Q1_6                                    1                       
Q1_7                                         1                  
Q1_8                                              1             
Q1_9       -0.209 -0.235      -0.226                 1.000      
Q1_10                                                          1

8.9.2 個人適合度

続いて個人の適合度です。 IRTにおける個人適合度の基本的な考え方としては「θp\theta_pが高いのに困難度が低い項目に間違えるのはおかしい」「θp\theta_pが低いのに困難度が高い項目に正解するのはおかしい」というものです。 するとこれは, 図 8.29 的な考え方になります。 Aさんのほうが予測値と実測値の間の差dpid_{pi}が大きいことは,前述の通り「θp\theta_p以外に回答に影響を与える要因がある(のでモデルが間違っている)」という可能性の他に,「(もしモデルは正しいと仮定したら)回答者の方に,何か項目iiの正解に導いた要因があるのでは」と考える,というわけです。

回答者の方の原因としては,心理尺度でいえば以下のようなものが考えられますMeijer et al., 2016

  • 回答者の言語能力が不十分なために項目の内容を十分に理解できていない可能性。低学齢への実施や第一言語以外での実施で起こる可能性がある。
  • 心理尺度が測定しようとしている心理特性に対する認識がない状態(Traitedness)。例えば「そんなこと考えたこともない」といった場合,項目の意図を他の人とは異なるニュアンスで捉えてしまい,結果的に「大多数が当てはまると回答する項目に当てはまらないと回答する」といった反応につながる。
  • 回答のモチベーションが低い可能性。やる気が無いので適当に回答している?
  • 特定の心理特性を持っている人は平均的な人より一貫性の低い回答をする可能性がある。例えば精神的な問題を抱えていたり,希死念慮を持ち合わせている人はそうでない人と比べて適合度が低いという先行研究もあるらしい。

ということで,個人適合度が低いイコール適当な回答者(だから除外してもよい)とも言い切れないのですが,特徴的な回答者をあぶり出すなどの目的でも個人適合度は使い道がありそうです。

実際に個人適合度の指標として提案されているものは相当数ある Karabatsos, 2003; Meijer & Sijtsma, 2001 のですが,わかりやすい&統計的仮説検定のフレームワークが整備されていると言った理由から,zhz_h(またの名をlzl_z)という統計量Drasgow et al., 1984が用いられることが多いようです Meijer et al., 2016

zh(lz)z_h(l_z)統計量

zh(lz)z_h(l_z)統計量では,ある回答者の実際の反応パタン𝐲p=[yp1yp2ypI]\symbf{y}_p=\begin{bmatrix}y_{p1}&y_{p2} &\cdots&y_{pI}\end{bmatrix}と「考えられる全反応パタン」での尤度を比較することにより,その回答者の反応パタンの「普通さ」を評価します。

まず,回答者ppの特性値の最尤推定値θ̂p\hat{\theta}_pに関して,その時の尤度を計算します。 局所独立性が満たされているならば,尤度は単純に全ての項目をかけ合わせたら良いので, l0=i=1IP(ypi=1)ypiP(ypi=0)1ypi(8.42) l_0=\prod_{i=1}^{I}P(y_{pi}=1)^{y_{pi}}P(y_{pi}=0)^{1-y_{pi}} \qquad(8.42) となります。 ただ掛け算は(コンピュータ的に)大変なので,対数尤度 ll0=i=1IypilnP(ypi=1)+(1ypi)lnP(ypi=0)(8.43) ll_0 = \sum_{i=1}^{I}y_{pi} \ln P(y_{pi}=1) + (1-y_{pi})\ln P(y_{pi}=0) \qquad(8.43) を使います。 この尤度ll0ll_0は,「特性値がθ̂p\hat{\theta}_pの人が𝐲p\symbf{y}_pという回答をする確率」のようなもの21です。 ここで𝐲p\symbf{y}_p以外の回答パタンに目を向けてみると,もしかしたら「𝐲p\symbf{y}_pよりもθ̂p\hat{\theta}_pと整合的な回答パタン」があるかもしれません。 もちろん反対に「θ̂p\hat{\theta}_pと全く合わない回答パタン」というものもあるかもしれません。

そこで,考えられる全回答パタン(2件法10項目ならば2102^{10}通り)の中で𝐲p\symbf{y}_pという回答パタンがθ̂p\hat{\theta}_pとどの程度整合的かを相対的にチェックすることを考えます。 全回答パタン内で𝐲p\symbf{y}_pという回答パタンとθ̂p\hat{\theta}_pの整合性(尤度)が平均くらいであれば,その回答パタンはまあ起こりうるものだろうとみなせる一方で,もし整合性が低い場合には,𝐲p\symbf{y}_pという回答パタンは起こりにくい(ズレが大きい)だろうと判断できる,ということです。

実際に「全回答パタンでの対数尤度の期待値」は(各項目への回答は離散変数なので) Eh=i=1IP(ypi=0)lnP(ypi=0)+P(ypi=1)lnP(ypi=1)(8.44) E_h = \sum_{i=1}^{I}P(y_{pi}=0)\ln P(y_{pi}=0) + P(y_{pi}=1)\ln P(y_{pi}=1) \qquad(8.44) として求めることができます。 つまりll0Ehll_0-E_hが「平均的な回答パタンでの尤度と実際の回答パタン𝐲p\symbf{y}_pでの尤度の差」であり,これがマイナス方向に大きいほどモデルとデータの適合度が悪いということを意味します。

ここで全回答パタンにおける対数尤度の分散を σh2=i=1I{P(ypi=1)P(ypi=0)[lnP(ypi=1)P(ypi=0)]2}(8.45) \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\} \qquad(8.45) として求めることができるのですが,これを用いると,ll0ll_0を全回答パタン内で標準化した値zh=ll0Ehσhz_h=\frac{ll_0-E_h}{\sigma_h}が経験的に標準正規分布に従うということが知られています。 ということで,|zh|>1.96|z_h| > 1.96となった回答者はどうやら「特性値がθ̂p\hat{\theta}_pのもとでは平均的ではない回答パタンである」と判断することができるわけです。

図 8.30 に,zhz_h統計量の考え方を表しました。 ヒストグラムは母集団からサンプリングされた様々なθ̂p\hat{\theta}_pの人と,その人の(ランダムサンプリングされた)回答ベクトル𝐲p\symbf{y}_pのもとで計算したzhz_hを集めたものというイメージです。 これが近似的に正規分布に従っているということで,この外側5%を棄却域とみなして検定を行うことができます。 Bさんのzhz_hは,この分布の中で割と平均的なところにあるため,Bさんの回答は「(Bさんの推定値θ̂B\hat{\theta}_{\mathrm{B}}における)ふつう」に近いとみなせます。 一方Aさんのzhz_hは,この分布の中でもかなり小さい方にあります。 ということは,Aさんの回答パタンはAさんの推定値θ̂A\hat{\theta}_{\mathrm{A}}からすると当てはまりが悪く(つまり困難度の低い項目に「当てはまらない」と回答し,困難度の高い項目で「当てはまる」と回答している可能性が高く),なんなら母集団全体の中でも相当悪いほうだと判断でき,何らかの問題があるのではないかと疑いをかけられるわけです。

図 8.30: 個人適合度の考え方

ちなみに,正規分布的には上側の外れ値も「ふつうではない」ということになりますが,こちらはどう考えたら良いかよくわかりません。 理論的には「特性値θ̂p\hat{\theta}_pに対してあまりにも適合しすぎている」という状態で,これは例えば困難度がθ̂\hat{\theta}以下の項目では全て「当てはまる」と回答し,θ̂\hat{\theta}以上の項目ではすべて「当てはまらない」と回答しているような感じなのですが,だからといって即座に問題とは言えなさそうです。 ということで,基本的にはzh<1.96z_h < -1.96の回答者に着目しておけば良いと思います。

Rでやってみる

※個人適合度と回答パタンの関係がわかりやすいので,ここはGRMで推定した結果を使用します。

mirtでは,personfit()という関数によってzhz_h統計量を出すことが出来ます。

コード 8.29: 個人適合度を出す
personfit(result_grm)
      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 ]

outfitinfitというものは,zhz_hよりもシンプルに予測値と実測値の間の差dpid_{pi}の二乗和をもとにした適合度指標Smith et al., 1995です22。 ということでこれらについても大きいほど適合度が悪いと判断できます。 実際にこれらの指標の相関を見てみると,いずれもかなり高い相関になっていることがわかります。

コード 8.30: 指標間の相関
# z_hだけは方向が違う(小さいほど悪い)ので他の指標と負の相関
cor(personfit(result_grm))
             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

では実際に,当てはまりが悪い人の回答パタンを見てみましょう。

コード 8.31: 個人適合度が最も悪い人
# 引数stats.only = FALSEにすると,統計量に加えて回答パタンも残してくれる
pf <- personfit(result_grm, stats.only = FALSE)
# Zhの値が最小の人を見てみる
pf[which.min(pf[, "Zh"]), ]
    Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10   outfit z.outfit
581    1    1    1    1    1    6    6    6    6     6 5.804625 5.127515
       infit  z.infit        Zh
581 5.552154 5.122718 -7.275079

いま扱っているモデルは10項目1因子のモデルです。 そして10項目の識別力は,いずれも正の値になっています。 つまり,基本的にはすべての項目間には正の相関があるため,この回答者のように「ある項目では1,別の項目では6」という回答の付け方には違和感があります。

コード 8.32: 個人適合度が悪い人
subset(pf, Zh < -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
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 ]

同様にして,zh<1.96z_h < -1.96の人たちを見てみると,項目ごとにかなりバラバラな回答をしている人たちがあぶり出されています。 ということで,zhz_h統計量によって怪しい回答者をあぶり出すことができそうだということがわかりました。

ちなみに反対にzh>1.96z_h > 1.96の人たちを見たところですべての項目で同じような値をつけており,やはり「適合度が悪い」という判断はなかなか難しそうです(ストレートライン気味な感じはありますが)。

コード 8.33: 個人適合度が良すぎる?人
subset(pf, Zh > 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 ]

ここで紹介したzhz_h統計量のような個人適合度は,モデルベースで適当回答者を除外するのに使うこともできます。 ただし,個人適合度はあくまでも(ロジスティックなどのIRT)モデルに基づいて推定されたパラメータをもとに算出される点には少し注意が必要です。 そもそもモデルの式や因子構造の設定が間違っている場合には,正しくない検出をしてしまうかもしれません。 そんなわけで,もちろん適合度の値だけで機械的に除外するのは危ないですが,例えば項目パラメータの推定値がなんだかおかしい場合には,適合度指標を参考にしたデータクリーニングによって,推定をより安定させられるかもしれません。

8.9.3 項目適合度

IRTでは項目一つ一つの性能(識別力や困難度)に関心があることが多いため,項目レベルでも適合度を考えます23。 基本的な考え方としては,推定されたパラメータに基づくICCが,データとどの程度一致するかを確認していきます。 適合度が悪い項目が見つかった場合には,項目の内容を吟味して削除するかを検討したり,場合によってはモデル自体を修正(因子数や因子構造の修正など)する必要があるかもしれません。 項目適合度は,具体的には以下のような手順によって算出していきます Orlando & Thissen, 2000

  1. 項目パラメータとθp\theta_pを推定する
  2. 回答者をθp\theta_pの値で並び替える
  3. 回答者をθp\theta_pの値によっていくつかのグループに分ける
  4. 各グループでの反応確率を計算する
  5. 4.で求めた反応確率と「モデル上の期待反応確率」をχ2\chi^2的な統計量や視覚的方法を用いて比較する

視覚的なチェック

ということで,まずは視覚的に比較してみましょう。 mirtでは,項目適合度を見るための関数としてitemfit()が用意されています。 これに引数empirical.plotを与えることでプロットを見ることができます。

コード 8.34: 項目適合度を視覚的にチェック
# 引数empirical.plotに項目番号か列名を指定します
itemfit(result_2plm, empirical.plot = 8)
図 8.31: 反応確率の比較プロット

図 8.31 を見ると,各カテゴリごとにICCが表示されています。これに重なっている小さいマルのY座標は,手順3で分割した各グループ(ここでは10グループ)におけるP(ypi=1)P(y_{pi}=1)を表しています(多値型の場合はカテゴリ選択確率が表示される)。 また,X座標はそのグループにおけるθp\theta_pの平均値を表しています。 ということで,各カテゴリにおいて,小さいマルと実線が近いほどデータとモデルの当てはまりが良い,と判断することができるわけです。

ただカテゴリ数や項目数が多くなってくるとすべてを目視で確認するのはかなり大変です。 また視覚的なチェックだけではどの項目が最も乖離しているのかを判断するのも大変でしょう。 ということで,怪しい項目にアタリをつけるための統計量を使っていきましょう。

カイ二乗的な統計量

χ2\chi^2統計量は独立性の検証のところ( Section 8.9.1 )で登場したものです。 IRTモデルに基づいて項目および回答者のパラメータが推定されると,各項目について,指定されたサンプルサイズのうち何人がどの選択肢を選ぶと期待されるか,を計算することができます。 これを期待度数EE,また実際にその項目に対する反応の分布を観測度数OOとして χ2=(OE)2E(8.46) \chi^2=\sum\frac{(O-E)^2}{E} \qquad(8.46) と計算すると,「帰無仮説:モデルがデータに適合している」が正しいならば小さい値におさまるはずだと考えられます。

Yen (1981 は,回答者をθp\theta_pの値の高低によってグループ分けし,グループごとにχ2\chi^2の計算をすることにしました。 特にグループ数を10に固定して算出した指標はQ1Q_1と呼ばれています。 グループg(g=1,2,,10)g~(g=1,2,\cdots,10)における特性値θ̂p\hat{\theta}_pの平均値をθg\bar{\theta}_gとおくと,期待正答率は(2パラメータロジスティックモデルならば) Eig=11+exp[α̂i(θgβ̂i)](8.47) E_{ig} = \frac{1}{1+\exp\left[-\hat{\alpha}_{i}\left(\bar{\theta}_g - \hat{\beta}_i\right)\right]} \qquad(8.47) と求めることができます。 すると,二値項目におけるQ1Q_1統計量は Q1=g=110Ng((OigEig)2Eig+[(1Oig)(1Eig)]21Eig)=g=110Ng((1Eig)(OigEig)2+Eig(OigEig)2Eig(1Eig))=g=110Ng((OigEig)2Eig(1Eig))(8.48) \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} \qquad(8.48) となります。 なお1行目の右辺は,第1項が正答セルにおける乖離度を表しています。 同様に,第2項は(1Oig),(1Eig)(1-O_{ig}),(1-E_{ig})に置き換わっていることから,誤答セルにおける乖離度を表していることがわかります。 したがって,Q1Q_1統計量は「グループ」×\times「項目iiの正誤」という10×210\times2クロス表におけるχ2\chi^2統計量を表しているのです。 そしてこのQ1Q_1統計量は自由度が1010-パラメータ数のχ2\chi^2分布に従うため,これを用いて検定ができるわけです。

同様に,差ではなく比を用いた統計量としてG2G^2統計量 McKinley & Mills, 1985 というものもあります。 G2=g=110Ng[OiglnOigEig+(1Oig)ln1Oig1Eig](8.49) 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] \qquad(8.49)

ほかにもいくつか適合度指標は提案されているのですが,基本的には前述の通り「期待反応確率と実際の反応確率の差」を何らかの形で統計量に落とし込んでいます。 細かな設定としては

  1. グループの数
  2. θp\theta_pの推定方法(基本的には最尤法,ただベイズ推定などでも可)
  3. 各グループの特性値の代表値の決め方(平均値 or 中央値)
  4. 期待反応確率の計算方法(個人のθ̂p\hat{\theta}_pにおける期待反応確率EipE_{ip}のグループ平均という手もあり)

といったところで様々なオプションが考えられるようです加藤 他,2014

Rでやってみる

項目適合度をを出す場合はitemfit()関数を使います。 引数fit_statsに,どの統計量を出してほしいかを指定します。なお,X2を指定すると内部ではQ1Q_1統計量(i.e., グループ数10)を出してきます。

コード 8.35: 項目適合度をチェック
# まとめて複数の統計量を出すこともできます
itemfit(result_2plm, fit_stats = c("X2", "G2"))
    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    0

G2の自由度(df.G2)がいろいろな値になっていますがあまり気にしなくてOKです。 今回も例によって検定の結果(p.*)は全て有意となりました。 このように,サンプルサイズが大きい場合これらの統計量はあまり参考にならない可能性が高いです。

…じゃあどうするんだ,という話になるわけですが,一つの方法としては効果量的な考え方をとるために「実際の差(OigEig)(O_{ig}-E_{ig})や比OigEig\frac{O_{ig}}{E_{ig}}を計算する」という方法が考えられます。 例えば110g=110|OigEig|\frac{1}{10}\sum_{g=1}^{10}|O_{ig}-E_{ig}|が大きい項目は,基本的に視覚的にプロットを見ても差があるでしょう。 ただ(OigEig)(O_{ig}-E_{ig})などを直接計算してくれる関数は用意されていないので,多少頑張る必要があります。

8.10 項目情報量・テスト情報量

IRTにかぎらず,統計的推測では推測の精度が問われます。 例えば正規分布の母平均の推測の場合,標準誤差は母分散σ2\sigma^2およびサンプルサイズnnを用いてσ2n\sqrt{\frac{\sigma^2}{n}}と表されました。 この式は,データの数が多いほど推測の精度が高くなるという自然な考えを表しています。

テストや心理尺度では,標準誤差を真値と誤差の関係から考えます。 Section 4.6 で紹介しましたが,潜在特性の測定(古典的テスト理論)では,観測得点xxと真値ttと誤差eeの関係を x=t+e(8.50) x = t + e \qquad(8.50) と表し,信頼性係数を ρ=1σe2σx2(8.51) \rho = 1-\frac{\sigma_e^2}{\sigma_x^2} \qquad(8.51) と定義していました。 標準誤差は「真値がttの人が測定を繰り返した場合に,(平均的に)どの程度観測得点xxがばらつくか」を表しています。 したがって,xN(t,σe2)x\sim N(t,\sigma_e^2)とすると,心理測定における標準誤差はσe\sigma_eにあたります。 信頼性係数ρ\rhoに関しては,上の式を変形させると σe=σx21ρ2(8.52) \sigma_e = \sigma_x^2\sqrt{1-\rho^2} \qquad(8.52) となることから,これが間接的に測定の精度を表す指標だということが言えます。 なお実際に標準誤差を計算する場合には,ρ\rhoの推定値としてα\alpha係数などを代入します。

さて,信頼性係数に基づく上記の標準誤差の式には真値ttは含まれていません。 つまり真値の値が何であれ標準誤差は同じということを意味しています。 ですがこの考え方は割と不自然なものです。 特性値の真値が平均くらいの人では,尺度の項目はほとんどすべてがその人の特性値を測定するのに意味のある(正解であるほど真値ttは高いだろうという判断につながる)項目です。 一方で,特性値の真値が非常に高い人では,実はたいていの項目はその人の特性値を測定するのにほぼ無意味です。 というのも,ある程度難易度の高い項目に正解する人は,難易度の低い項目にはほぼ確実に正解するだろうと言えるためです。 つまり,「ある程度難易度の高い項目」の情報さえあれば「難易度の低い項目」の情報は無くても特に問題ない(真値の推定には影響しない)と言えます。 先ほど,標準誤差はサンプルサイズが多くなるほど小さくなる,ということを説明しました。 本来,特性値の真値によって「その人の特性値を測定するのに意味のある」項目の数は変わるはずだと考えると,標準誤差も個人(の特性値の値)ごとに変わると考えるのが自然な気がしてきます。

IRTでは,この問題を克服するために項目情報量 (item information)およびテスト情報量 (test information)というものを用います。 この項目情報量・テスト情報量は,後述するテストの設計などにおいても大きな力を発揮してくれます。

8.10.1 情報量とは

項目情報量の話に入る前に,少しだけ「情報量」とは何かについてごくごく簡単に説明しておきましょう。 我々は,一般名詞として「情報」という言葉を使うことが多々あります。 この時の「情報」は,何かしらの事柄について「知らなかった」状態から「知っている」状態に遷移させるものという見方ができます。

この「情報」を確率の世界に持ち込むため,以下の例え話における「情報」について考えてみましょう。

いま,Aさんは52枚のトランプから1枚のカードを引きました。 Aさんの向かいに座っているBさんには,もちろんAさんのカードが何かは全くわかりません。 つまりこの段階でBさんがAさんのカードを言い当てる確率は1/521/52です。 ここで,Aさんから「スペードのカードを持っている」という情報が得られたとします。 すると,Bさんの中にあったカードの候補は52枚から13枚に減ります。 この時点では,BさんがAさんのカードを言い当てる確率は一気に4倍になり1/131/13になっています。 「スペードのカードを持っている」という情報はBさんの予測の精度を一気に4倍に引き上げたわけです。

ここで,Bさんは次に「黒のカードを持っている」という情報を得たとします。 この情報では,Bさんの中にあったカードの候補は13枚のまま変わりません。 したがってこの情報は,今のBさんにとっては何の価値も持っていません。 一方で,Bさんが「エースのカードを持っている」という情報を得た場合には,前の情報と合わせて「スペードのエース」であることが確定します。 この場合には,BさんがAさんのカードを当てる確率はさらに一気に13倍となるわけです。

この例え話からは,情報量に関する以下の性質がわかります。

  • 情報は予測の精度を高める: 情報量が多くなるほど,正確に予測ができるようになります。
  • より不確実性を下げる情報のほうが情報量が多い: マークは4通りである一方で,数字は13通りあるため,数字の情報のほうが不確実性をより大きく下げている=情報量が多いと解釈できます。
  • 独立な情報の情報量は加法的である: マークと数字は互いに独立なので,先に「エースを持っている」という情報を得たとしてもカードを当てる確率はやはり13倍(1/521/52から1/41/4)になっていました。

ここまでの話をIRTに置き換えてみます。 θp\theta_pの予測に際して,項目がもつ情報量を考えてみましょう。 図 8.32 には困難度がそれぞれ(0,3)(0,3)の2つの項目のICCが描かれています。

図 8.32: IRTにおける情報の考え方

ある回答者ppは,それまでの項目への解答状況からθ̂p=0\hat{\theta}_p=0だと予測されているとします。 すると,この回答者ppにとって,緑の項目(βi=0\beta_i=0)は正解するかどうかがかなり不確実(50%)です。 そしてこの項目への反応によってθ̂p=0\hat{\theta}_p=0の予測は変化し,多少なりともより正確なものになるでしょう。 そう考えると,緑の項目はθ̂p=0\hat{\theta}_p=0の人にとってそれなりの情報量を持っていると考えられます。

一方で赤い項目(βi=3\beta_i=3)を見ると,θ̂p=0\hat{\theta}_p=0の人はほぼ確実に正解しないだろうという予想がつきます。 つまり,出題時点で赤い項目に対する不確実性はほぼゼロ,ということです。 実際にこの回答者ppが赤い項目に不正解だったとしても,「そうでしょうね」という話なのでθp\theta_pの予測に対して特に新規の情報をもたらすものとは言えません24。 同様に,緑の項目でもθ̂p=3\hat{\theta}_p=3の人にとってはほぼ情報量のない項目,ということになります。 このように,ある項目は,その項目の困難度が回答者のθ\thetaに対してちょうどよいほど,θ\thetaの推定に対して多くの情報量を持っているのです。 ということである項目が持つ項目情報量は,実際にはθ\thetaの値によって変動する関数になっています。 この関数のことを項目情報関数 (item information function [IIF])と呼びます25。 この節のはじめにお話したように,IRTでは真値によって標準誤差が変わるようにすることを考えます。 この後具体的な関数の形を示しますが,項目情報関数をθ\thetaの関数として表すことによってその目的を達成しようというわけです。

8.10.2 項目情報関数

IRTにおける情報量の考え方のイメージを確認したので,具体的な項目情報関数の式を見てみましょう。 二値型モデルの場合,項目情報関数Ii(θ)I_i(\theta)Ii(θ)=P(yi=1)2P(yi=1)(1P(yi=1))(8.53) I_i(\theta)=\frac{P'(y_{i}=1)^2}{P(y_{i}=1)(1-P(y_{i}=1))} \qquad(8.53) という式になります。ここでP(yi=1)P'(y_{i}=1)は,P(yi=1)P(y_{i}=1)θ\thetaで微分した導関数です。 …といってもよくわからないと思うので,実際に見てみましょう。 図 8.33 は, 図 8.32 に示したICCを持つ2項目のロジスティックモデルにおける項目情報関数です。 どちらの線も先程説明したように,正解するかが最も不確実なθp=βi\theta_p=\beta_iのところで情報量が最大になっています。

図 8.33: 項目情報関数

IIFの式 Ii(θ)=P(yi=1)2P(yi=1)(1P(yi=1))(8.54) I_i(\theta)=\frac{P'(y_{i}=1)^2}{P(y_{i}=1)(1-P(y_{i}=1))} \qquad(8.54) を求めるためには,P(yi=1)P(y_{i}=1)の導関数を求める必要があります。 正規累積モデルは項目反応関数の中に積分が含まれているためこれを求めるのはかなり大変なのですが,それと比べるとロジスティックモデルでは解析的にIIFの式を求めることができます。 例えば2パラメータロジスティックモデルの場合,IIFは Ii(θ)=αi2P(yi=1)(1P(yi=1))(8.55) I_i(\theta)=\alpha_i^2P(y_{i}=1)(1-P(y_{i}=1)) \qquad(8.55) となります。 この式からも,IIFはP(yi=1)=0.5P(y_{i}=1)=0.5となるθ=βi\theta=\beta_iの点において最大値をとることがわかります。 同時に,識別力αi\alpha_iが大きいほど項目情報量も多くなることがわかります。

図 8.34 に,識別力が異なる2つの項目(βi=0\beta_i=0; ICCは 図 8.8 ) のロジスティックモデルでのIIFを示しました。 識別力が高い項目ほど,項目特性関数の傾きが大きくなっていました。 そのためθp\theta_pが高い人は正解しθp\theta_pが低い人は不正解する,というメリハリがついていました。 θp\theta_pの推定という視点から言えば,αi\alpha_iが大きい項目ほどその項目への正誤が特性値の推定により確かな情報を与えるという意味で,α\alphaと項目情報量には密接な関係があるわけです。

図 8.34: 異なる識別力を持つ項目のIIF

しかし改めて 図 8.34 をよく見ると,θp=0\theta_p=0付近では識別力の高い緑の項目のほうが高い情報量を持っている一方で,θp\theta_pが0から離れたところではむしろ識別力の低い青い項目のほうが高い情報量を持っています。 そして対応する 図 8.8 を見ると,識別力αi\alpha_iが高くなるほど,θp\theta_pβi\beta_iから離れると急激にP(yi=1)P(y_{i}=1)が0.5から離れていることが分かります。 つまり「正解するかが不確実」な(≒情報量を持ちうる)θp\theta_pの範囲は識別力の高さと反比例しており,その結果θp\theta_pβi\beta_iから離れたところでの項目情報量も少なくなってしまうのです。

この事実は「識別力は高ければ高いほど良いわけではない」ことを教えてくれています。 極端な例として「θp<βi\theta_p<\beta_iの人は確実に間違えるが,θpβi\theta_p\geq\beta_iの人は確実に正解する」項目を考えてみましょう。 その項目は識別力αi\alpha_i\rightarrow\inftyであり, 図 8.35 のようなICCをもつ項目です。

図 8.35: 完全な識別力を持つ項目があったら

この項目は「θp\theta_pβi\beta_iより上か下か」に関しては完全な情報を教えてくれます(Ii(θp=βi)=I_i(\theta_p=\beta_i)=\infty)。 その一方でθp>βi\theta_p>\beta_iの範囲では,θp\theta_pβi\beta_iよりどの程度高かろうとP(yi=1)P(y_{i}=1)は100%(逆もまた然り)なので,「θp\theta_pβi\beta_iよりどの程度上か下か」に関しては全くわかりません(Ii(θpβi)0I_i(\theta_p\neq\beta_i)\simeq0)。 このように,識別力が高すぎる項目は有効なθp\theta_pの範囲が狭すぎるため,多くの回答者にとって意味のない項目になってしまう可能性があるわけです。

多値型モデルでの項目情報量

(いつか書きたいという気持ちだけは常に持っている)

8.10.3 テスト情報関数

ここまで,項目情報関数の定義およびその性質を見てきました。 この「情報」は,θp\theta_pを推定するための「情報」を表しているわけですが,θp\theta_pの推定は尺度・テスト内のすべての項目への回答を総合して行われるものです。 そのため,次は項目情報関数を拡張して尺度全体での情報関数を考えてみます。

情報量には「独立な情報の情報量は加法的である」という性質がありました。 IRTでもこの性質は健在です。 もしテスト全体が局所独立の仮定を満たしているならば,テスト情報関数(test information function [TIF]) は単純に I(θ)=i=1IIi(θ)(8.56) I(\theta) = \sum_{i=1}^{I}I_{i}(\theta) \qquad(8.56) と表せます。 そしてテスト情報量の重要な性質として,真の特性値がθp\theta_pの人における最尤推定量θ̂p\hat{\theta}_pの誤差分散は漸近的に σ(θ̂p|θp)2=1I(θp)(8.57) \sigma^2_{(\hat{\theta}_p|\theta_p)}=\frac{1}{I(\theta_p)} \qquad(8.57) となることが知られています。 ということで,標準誤差は σ(θ̂p|θp)=1I(θp)(8.58) \sigma_{(\hat{\theta}_p|\theta_p)}=\frac{1}{\sqrt{I(\theta_p)}} \qquad(8.58) となります。 これは古典的テスト理論の表記に合わせると,最尤推定値θ̂p=θp+εp\hat{\theta}_p=\theta_p+\varepsilon_pとした場合のσε\sigma_{\varepsilon}に相当するものです。 このように,IRTでは測定の精度を表す標準誤差を,θp\theta_pの関数として求めることができるわけです。

8.10.4 Rで情報量を確認する

項目・テスト情報関数の考え方がわかったところで,実際の項目での関数を確認してみましょう。 これまでICCを出すために使っていたplot()関数によって,IIFも出すことができます。 項目情報関数を出す場合にはtype = "infotrace"とします。

コード 8.36: 項目情報関数
plot(result_2plm, type = "infotrace", facet_items = FALSE)
図 8.36: 項目情報関数

先ほど説明した通り,低識別力の項目1は他の項目と比べて圧倒的に情報量が少ないことがわかります。 また、θp>2\theta_p>2のエリアでは、Q1_10以外のほぼどの項目も情報量が無いため、この10項目では特性値が高い人についてはまともな推定はできなさそうだということがわかります。 この項目情報関数を見れば,例えばどうしても1項目削除しないといけないとしたら,基本的にはQ1_1を削除したら良いが,もしもθp=6\theta_p=-6付近がメインターゲットだとしたらむしろQ1_1は残したほうが良いとか、θp>2\theta_p>2がメインターゲットならば現状の項目セットではどうしようもないので新規項目を追加する必要がある、というように状況に応じた決定をサポートしてくれるでしょう。

具体的に特定のθp\theta_pの値における項目情報量を出したい場合にはiteminfo()関数を使います。

コード 8.37: 項目情報量
# はじめに項目パラメータの情報だけ抽出する必要がある
# 1番目の項目を抽出するなら
item1 <- extract.item(result_grm, item = 1)
# 引数Thetaには,どの値での特性値を出したいかを指定する
iteminfo(item1, Theta = c(-4, -2, 0, 2, 4))
[1] 0.10354648 0.10516020 0.10114364 0.08337954 0.04723855

続いてテスト情報関数です。 こちらは引数をtype = "info"とすることで表示できます。

コード 8.38: テスト情報関数
# 引数typeは同じ
# type = "infotrace"とすると前述の項目情報関数になる
plot(result_2plm, type = "info")
図 8.37: テスト情報関数

尺度(10項目)全体を見ても,項目困難度がまんべんなく存在している-3から0付近の情報量が高い一方で,すべての項目の困難度が0以上なのでθp\theta_pが正のエリアでは情報量がさほど高くないことがわかります。

プロットする際に引数type"infoSE"にすると,標準誤差を合わせて出してくれるようになります。

コード 8.39: テスト情報関数に標準誤差を添えて
# 引数typeは同じ
plot(result_2plm, type = "infoSE")
図 8.38: テスト情報関数に標準誤差を添えて
うまくいかない場合

もしかしたら,以下のようなエラーが出てプロットが表示されないかもしれません。

Error in eval(expr, envir, enclos): latticeExtra package is not available. Please install.

その場合は,指示に従ってinstall.packages("latticeExtra")をしてから再度実行してください。

最も情報量が高いθp=1.4\theta_p=-1.4付近では,標準誤差はだいたい12.820.60\frac{1}{\sqrt{2.82}}\simeq0.60くらいになっています。 仮にθp\theta_pの真値が1.4-1.4の人が大量にいた場合,その人達の推定値は平均で0.6くらいは上下するだろう,ということです。 裏を返せばその程度の精度でしか推定できていないならば,例えばθ̂p=1.3\hat{\theta}_p=-1.3の人とθ̂p=1.5\hat{\theta}_p=-1.5の人がいたとしても,前者のほうが真の特性値が高いとはなかなか言い切れないわけですね。

具体的に特定のθp\theta_pの値におけるテスト情報量を出したい場合にはtestinfo()関数を使います。

コード 8.40: テスト情報量
# 引数Thetaには,どの値での特性値を出したいかを指定する
testinfo(result_grm, Theta = c(-4, -2, 0, 2, 4))
[1] 2.8917688 3.9684301 3.8450407 2.5423962 0.5038195
コード 8.41: 標準誤差
# 標準誤差を出すならばこの逆数のルートを取れば良い
1 / sqrt(testinfo(result_2plm, Theta = c(-4, -2, 0, 2, 4)))
[1] 1.2460379 0.6335903 0.7433266 1.8358362 5.2840351

θp=4\theta_p=4の人では,テスト情報量が0.5程度しかなく,標準誤差が5.28と非常に大きくなっています。 つまりこの10項目で推定を行っても,θ\thetaの真値が4の人に対する推定値は平均で±5\pm 5くらいには変動するわけです。 そう考えると,この10項目ではどうあがいてもθp\theta_pが高い層においてまともな推定値は得られないということが分かります。

8.11 (おまけ)多母集団モデルと特異項目機能

多母集団同時分析のところ Section 7.9 では,多母集団同時分析という枠組みを紹介しました。 当然ながらIRTにおいても,多母集団同時分析を実行することができます。 そもそもSEMやIRTにおいて多母集団モデルを行う理由が何だったのかを思い出してみると,その一つには「集団ごとに異なる回答傾向の原因を探る」という点がありました。 SEMにおける多母集団モデルでは,例えばある(心理)尺度への回答において,特定のグループ(e.g., 男性・公立高校生・日本人)と別のグループ(e.g., 女性・私立高校生・米国人)で平均点に差があったとしたら,

  1. まず測定不変性(因子負荷および切片が同じ)を確認した上で
  2. 因子得点の平均値を集団ごとに自由推定することで差を見る

という手続きを取るのが一般的です。 ここで重要になってくるのが,そもそも回答傾向に違いが生じる原因は2種類あるという点です。 先程の手続きの2つのステップに対応する形で

  1. そもそも項目・尺度が測っている構成概念の意味がグループ間で異なる
  2. 構成概念の意味は同じだが,その特性の強さがグループ間で異なる

という2つの可能性が考えられます。 SEMの多母集団モデルでは,多くの場合2番目に強い関心があることに加えて,ひとつひとつの項目よりは尺度全体としての測定不変性に関心があるためか,1番目についてはあくまでも前提条件扱いであまり重要視されていない気がします。

これに対してIRTの場合,観測されるデータは回答者と項目の交互作用として考えられる側面が強いため,「具体的にどの項目がグループ間で異なる挙動をしているか」というような項目側の要因にも関心があることが多いです。 この「異なるグループ間で項目が異なる挙動をする」ことを,一般的には特異項目機能 (Differential Item Functioning [DIF])と呼びます。 そんなわけで,IRTにおける多母集団モデルは主に以下の2種類の用途で用いられることがあります。

  1. 特性値θ\thetaの分布はグループ間で共通として,グループごとに推定される項目パラメータα,β\alpha, \betaの差異(=DIF)を見る
  2. 項目パラメータα,β\alpha, \betaはグループ間で共通だとして,特性値θ\thetaの分布のグループごとの差異を見る(SEMで言うところの強測定不変モデルと同じ)

mirtパッケージで多母集団IRTモデルを実行するためには,mirt()関数の代わりにmultipleGroup()関数を使用します。 multipleGroup()関数では,mirt()関数に引数に加えてgroupという引数が登場します。 lavaanパッケージのときにも多母集団モデルでは引数groupを指定しましたが,その時とは指定の仕方が異なるのでご注意ください。

コード 8.42: 多母集団IRTモデル(用途1)
# modelに1を与えておけば普通の一次元モデルになる
# 引数groupにはcharacter型かfactor型のベクトルを与える必要がある
result_group <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender))

結果を見るための関数は概ねmirt()のときと同じように使えます。 とりあえず推定された項目パラメータを見てみましょう。

コード 8.43: 項目パラメータのチェック
coef(result_group, simplify = TRUE, IRTpars = TRUE)
$`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はそれぞれ01となっています。 multipleGroup()関数は,デフォルトではこのように「各グループのθ\thetaの母集団分布がそれぞれ標準正規分布である」という設定のもとで項目パラメータを推定します。

特性値θ\thetaの分布のグループごとの差異を見たい場合には,SEMの多母集団モデルと同じ要領で等値制約を表す引数を与える必要があります。 lavaanでの多母集団モデルを実行したときには,因子負荷の等値制約などを引数group.equalで指定しました。 multipleGroup()では,invarianceという引数が用意されています。 この引数には, 表 8.4 に示す値を入れることが出来ます。 free_meansおよびfree_varsという値があるということは,これらを設定しない限り,全てのグループのθ\thetaは標準正規分布に固定されてしまうということなのでご注意ください。

表 8.4: 引数invarianceに指定できるもの
指定 制約されるもの
free_means グループ1以外のθ\thetaの平均値を自由推定する
free_vars グループ1以外のθ\thetaの分散を自由推定する
slopes 項目の識別力パラメータをグループ間で同じとする
intercepts 項目の切片パラメータをグループ間で同じとする
コード 8.44: 多母集団IRTモデル(用途2)
# 引数invarianceで等値制約を指定する
# "free_means" と "free_vars"を入れないとN(0,1)のままなので注意
result_group2 <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender), invariance = c("free_means", "free_vars", "slopes", "intercepts"))

改めて推定されたパラメータを見てみましょう。

コード 8.45: 項目パラメータのチェック
coef(result_group2, simplify = TRUE, IRTpars = TRUE)
$`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×22\times2クロス表に対するχ2\chi^2検定を拡張したMantel-Haenszel検定 (マンテル・ヘンツェル検定: Mantel & Haenszel, 1959に基づく方法です。 手法自体は非常に古いものですが,その使いやすさから近年でもDIF検出の方法として非常に多く用いられているようです Berrío et al., 2020

MH検定では,まず各集団を何らかの変数の値によってKK層に層化します。 DIF検出の場合には,IRTで推定された特性値θ\thetaや合計点などによって層化するのが一般的です。 表 8.5 は,そうして層化されたうちの第kk層について,ある二値項目に対する回答をグループごとに集計したものです。 MH検定では, 表 8.5 のような2×22\times2クロス表がKK個あるときに,それらを全部ひっくるめて独立性の検定を行う手法と言えます。

表 8.5: 第kk(k=1,2,,K)(k=1,2,\cdots,K)におけるクロス表
項目jj
項目ii 当てはまらない 当てはまる
グループA xA0(k)x_{A0}^{(k)} xA1(k)x_{A1}^{(k)} nA(k)n_A^{(k)}
グループB xB0(k)x_{B0}^{(k)} xB1(k)x_{B1}^{(k)} nB(k)n_B^{(k)}
n0(k)n_0^{(k)} n1(k)n_1^{(k)} NN

まずは一つの層から考えることで,具体的な方法を紹介していきましょう。 例えばθ\thetaが最も低い人達が集められた第1層において,グループAの人たちのほうが「当てはまる」と回答した割合が高かったとします。 このように,特定の層に絞ったときに回答傾向に違いが見られるというのも立派なDIFです。 そしてDIFが発生している場合には,2つのグループにおける回答の割合つまりオッズ比 OR(k)=xA1(k)/nA(k)xB1(k)/nB(k)(8.59) \mathrm{OR}^{(k)} = \frac{x_{A1}^{(k)}/n_A^{(k)}}{x_{B1}^{(k)}/n_B^{(k)}} \qquad(8.59) が1ではない,という見方ができます。 MH検定では,このオッズ比を全ての層について統合したもの OR(MH)=k=1KxA1(k)/nA(k)k=1KxB1(k)/nB(k)(8.60) \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)}} \qquad(8.60) が1であるかどうかを検定します。

MH検定の結果が有意であれば,データ全体として一方のグループのほうが「当てはまる」と回答した割合が高かったと言えます。 ICC的に言えば 図 8.39 のような状況であると言えるわけです。

図 8.39: MH検定で検出可能な(均一)DIF

一方で,MH検定では検出できないDIFも存在します。 図 8.40 に示されているDIFでは,θp\theta_pの低い層ではグループAのほうが「当てはまる」割合が高い一方で,θp\theta_pの高い層では逆転が起こっています。 OR(k)\mathrm{OR}^{(k)}を一つ一つ見れば1ではないわけですが,これを統合したOR(MH)\mathrm{OR}^{(\mathrm{MH})}は1になってしまうわけです。 MH検定による方法は, 図 8.39 に示した均一DIFは検出可能な一方で, 図 8.40 のような不均一DIFは検出できないかもしれない方法であることに気をつけましょう。

図 8.40: 不均一DIF

項目パラメータの推定値を比べる方法

IRTに基づくパラメトリックな方法として,まずは項目パラメータの推定値の差を直接検定する方法を紹介します。 Wald検定は,パラメータの推定値に対して用いられる検定の一種です。 一般的に用いられる一変量Wald検定では,パラメータξ\xiの推定値ξ̂\hat{\xi}に対して W=(ξ̂ξ0σξ̂)2(8.61) W = \left(\frac{\hat{\xi}-\xi_0}{\sigma_{\hat{\xi}}}\right)^2 \qquad(8.61) という検定統計量を考えます。ただしξ0\xi_0は帰無仮説における値,σξ̂\sigma_{\hat{\xi}}ξ̂\hat{\xi}の標準誤差です。 検定統計量WWは自由度1のχ2\chi^2分布に従うことが知られているため,これを用いて「ξ̂\hat{\xi}ξ0\xi_0である」という帰無仮説を検定します。

Lord (1980 はこのWald検定によってDIFを検出する方法を提案しました。 DIFが発生している場合には,グループごとに推定された項目パラメータに差があると考えられるため,(8.61)式の中身を「2つのグループでの推定値」に置き換えることにより,例えば項目iiの困難度βi\beta_iであれば W=(βi(A)βi(B)σβi(A)βi(B))2=(βi(A)βi(B)σβi(A)+σβi(B))2(8.62) \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} \qquad(8.62) という形で,同じようにWald検定を実行することができます。

Wald検定で用いるWald統計量WWは「推定値の差」を「標準誤差」で割った形をしているため,パラメータの数が複数であっても同じように定義することができます。 例えば項目iiの識別力と困難度を並べたベクトル𝛏i=[αiβi]\symbf{\xi}_i = \begin{bmatrix}\alpha_i & \beta_i\end{bmatrix}に対して W=𝛏i𝚺1𝛏i(8.63) W = \symbf{\xi}_i\symbf{\Sigma}^{-1}\symbf{\xi}_i^\top \qquad(8.63) という統計量(𝚺\symbf{\Sigma}は標準誤差行列)を考えてやれば,これが自由度2(=パラメータの数)のχ2\chi^2分布に従うため,同じようにして検定ができるのです。

項目特性曲線を比べる方法

IRTでは,推定された項目パラメータをもとに項目特性曲線(ICC)を描くことが出来ました。 ICCは,いわば「得点の期待値」を表したグラフです。 したがって「グループ間でICCのずれが大きい」場合にはDIFが発生していると考えることが出来そうです Raju, 1988。 ICCのずれとは, 図 8.41 でオレンジ色に塗りつぶされた部分の面積のことです。 これが大きい場合には,確かにDIFが発生していると言えそうです。

図 8.41: グループ間でのICCのズレ

グループAにおけるICC(正答確率)は,2パラメータロジスティックモデルならば Pi(A)(θ)=P(ypi=1|θ)=11+exp[αi(A)(θβi(A))](8.64) P_i^{(A)}(\theta) = P(y_{pi}=1|\theta) = \frac{1}{1+\exp[\alpha_i^{(A)}(\theta-\beta_i^{(A)})]} \qquad(8.64) と表せるため,ICCのズレの面積に基づくDIFの指標として,例えば DIFi=|Pi(A)(θ)Pi(B)(θ)|f(θ)dθ(8.65) \mathrm{DIF}_i = \int_{-\infty}^{\infty}\left\vert P_i^{(A)}(\theta) - P_i^{(B)}(\theta)\right\vert f(\theta)d\theta \qquad(8.65) というものを考えることができるわけです(e.g., Chalmers, 2018; Raju et al., 1995。ただしf(θ)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"としたら良いわけです。

コード 8.46: DIFの検出(1)
DIF(result_group, which.par = c("a1", "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の多母集団モデルのときにも登場したAICBICなどの情報量規準およびX2χ2\chi^2)統計量とこれに関連した自由度・pp値が表示されています。 例えばQ1_1行のAICの値-24.578というのは,「Q1_1のみグループごとに項目パラメータが同じ=DIFが無いという制約をかけることで,AICが24.578悪化した」ということを意味します。 したがって,AICのみで判断するならばQ1_1にはDIFが生じている,と言えるわけです。 同様に,χ2\chi^2検定の結果も「その項目にはDIFが発生していないとする(グループ間で等値制約を置く)ことでモデル適合度がどの程度悪化するか」という観点に基づく尤度比検定が行われています27

実際にDIFがありそうと判断された項目についてICCを見るには,引数plotdif = TRUEを与えてください。

コード 8.47: DIFがありそうな項目のICC
# デフォルトではSABICがマイナスになっている項目のプロットを表示する
DIF(result_group, which.par = c("a1", "d"), 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
図 8.42: DIFがありそうな項目のICC

確かにICCにはそこそこのずれがありそうです。

「項目1つ1つについてDIFの有無をチェックする」という意味では,全く逆方向からのアプローチを取ることも可能です。 パラメータ推定時にinvariance = c("free_means","free_vars","slopes","intercepts")を与えると,全ての項目パラメータが同じ,すなわち全ての項目にDIFが無いことを仮定した上での推定を行うことが出来ました。 このときの結果(result_group2)をベースラインとして考えると,DIF()関数には特定の1項目のみから等値制約を外した(=その項目だけDIFがあると仮定した)ときの変化を計算してもらいたいところです。 これを実現するためには,DIF()関数の引数schemeを指定する必要があります。

コード 8.48: DIFの検出(2)
# 等値制約を"drop"して検証してほしい
DIF(result_group2, which.par = c("a1", "d"), scheme = "drop")
      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<.05p < .05)項目は等値制約を外したほうが適合度が良くなっている=DIFの疑いあり,と判断することができます。

続いてはLordのWald検定の結果を出してみます。 ただ,Wald検定を行うためには「推定値の標準誤差行列」が必要になるため,まずは標準誤差を推定するように引数SE = TRUEを付け加えてやり直す必要があります。

コード 8.49: 標準誤差を推定してほしい
# SE = TRUEにするだけ
result_group <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender), SE = TRUE)

推定された標準誤差行列は一応result_group@vcovで見ることができますが,見たところでよくわからないと思うのでそっとしておいてあげてください。

無事標準誤差が推定できたので,これを用いてWald検定に望みます。 Wald検定は,DIF()関数で引数Wald = TRUEを与えることで実行可能です。

コード 8.50: DIFの検出(3):Wald検定
DIF(result_group, which.par = c("a1", "d"), Wald = TRUE)
      groups      W df     p
Q1_1     1,2 28.652  2     0
Q1_2     1,2 36.341  2     0
Q1_3     1,2 22.179  2     0
Q1_4     1,2 16.473  2     0
Q1_5     1,2 17.345  2     0
Q1_6     1,2  1.786  2 0.409
Q1_7     1,2  8.373  2 0.015
Q1_8     1,2  5.129  2 0.077
Q1_9     1,2 19.703  2     0
Q1_10    1,2 13.783  2 0.001

結果の見方はこれまでと基本的に同じです。 つまり検定結果が有意になった(p<.05p < .05)項目はグループ間でパラメータの推定値が有意に異なる=DIFの疑いあり,ということになります。

最後はICCのズレに基づく評価方法です。 こちらは Chalmers (2018 による方法がDRF()という関数に実装されています。 DRF()関数は,デフォルトでは項目特性曲線(ICC)の代わりにテスト特性曲線(TCC)の差について評価を行います28。 項目ごとのDIFを見たい場合には,引数DIF = TRUEを与える必要があるのでご注意ください。 また引数plot = TRUEを与えると,「ICCの差のプロット」を出してくれるようになります。

コード 8.51: DIFの検出(4):ICCのズレをチェック
DRF(result_group, DIF = TRUE, plot = TRUE)
図 8.43: 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に適当な数を指定してください。

コード 8.52: ブートストラップを使って検定
# drawsの値は大きい方が正確になるが,その分時間がかかる
DRF(result_group, DIF = TRUE, draws = 500)
$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
コード 8.53: DIFの検定のプロット
# plot = TRUEをつけると検定はせずに図だけ表示する
DRF(result_group, DIF = TRUE, draws = 500, plot = TRUE)
図 8.44: DIFの検定のプロット

結果の見方はやはり同じです。 つまり検定結果が有意になった(p<.05p < .05)項目はIRFがグループ間でズレている=DIFの疑いあり,ということになります。

DIF検出の実際

ここまでDIFを検出するいくつかの方法を紹介しました。 色々な観点からの評価がある,という意味では,やはり複数の指標を組み合わせて総合的に判断するのが良いのだろうと思います。

それはそれとして,DIFの有無を評価しようとする場合,そもそも関心のあるグループの間で特性値θ\thetaの分布に差が無いということはあまり考えられません(無かったら非常にラクなのですが…)。 したがって,そもそも「θ\thetaの分布が同じだとしてDIFを検出する」というスキームは現実的に結構厳しいのです。 この問題を克服するためには,「グループ間で挙動が変わらない=DIFが無いことが事前にわかっている項目」であるアンカー項目 (anchor items) を加えておくことが考えられます。 アンカー項目の項目パラメータを基準とすることで,θ\thetaの分布がグループで異なる場合であっても特定の項目にDIFが生じているかを評価することが出来るわけです。 具体的にアンカー項目をどのように決定したら良いかはまた難しい話(レビューとして Kopf et al., 2015なので,ここではとりあえずアンカー項目がある場合のやり方だけ紹介しておきます。

multipleGroup()の引数invarianceには,実は項目の名前(データの列名)を与えることも出来ます。

コード 8.54: アンカー項目を用いた多母集団モデル
# Q1_1からQ1_3がアンカー項目だったとします
result_anchor <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender), invariance = c(colnames(dat_bin)[1:3], "free_means", "free_vars"))

あとの流れはほぼ同じなので,紹介はここまでとします。

コード 8.55: 狙った等値制約が置かれていることをチェック
# Q1_1からQ1_3のパラメータだけグループ間で同じ
coef(result_anchor, simplify = TRUE, IRTpars = TRUE)
$`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
コード 8.56: DIFの検出
# もともと等値制約が置かれているQ1_1からQ1_3については結果が出ない
DIF(result_anchor, which.par = c("a1", "d"))
      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モデルとして, Section 8.6 ではGRMやGPCMといったモデルを紹介してきました。 多値反応に対するIRTモデルには他にも様々なものが提案されており,近年注目を集めている(気がする)のがIRTree Böckenholt, 2017 と呼ばれるモデル群です。 これらのモデルでは,回答決定のプロセスとして名前通り「ツリー構造」を仮定します。 例えば「とても嫌い(1)」「嫌い(2)」「どちらでもない(3)」「好き(4)」「とても好き(5)」という5件法の項目があったとしましょう。 この項目への回答決定プロセスの考え方の一例としては, 図 8.45 のように

  1. まず「どちらでもない(3)」を選択(態度保留)するかを二択で考える
  2. 態度を表明する場合には「好き(4-5)」か「嫌い(1-2)」かを二択で考える
  3. 「とても(1,5)」の選択肢を選ぶほど好き or 嫌いかを二択で考える

という感じで3つの二択を順番に行う,というものがあるでしょう。

図 8.45: IRTreeの一例

するとこの多値項目の背後には,実は3つの二値項目(query)が存在していると考えることができ,それぞれが

  1. 「どちらでもない」を選ぶ傾向(中心傾向):θp(M)\theta_p^{(M)}が強い人かどうか
  2. 実際にその対象物が好きかどうか:θp(A)\theta_p^{(A)}
  3. 極端な選択肢を選ぶ傾向(極端傾向):θp(E)\theta_p^{(E)}が強い人かどうか

という異なる個人特性によって動いていると仮定することができそうです(もちろんプロセス単位で完全に独立しているとは思いませんが,ある程度は)30。 各queryに対する反応確率は,普通の正規累積・ロジスティックモデルなどで表現できます。 例えば中心傾向についてのqueryだけで表現できる「どちらでもない」を選ぶ確率は P(ypi=3|θp(M))=αi(M)(θp(M)βi(M))12πexp(z22)dz=𝚽[αi(M)(θp(M)βi(M))](8.66) \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} \qquad(8.66) などとなるわけです。IRTreeモデルでは,項目パラメータにも上付き添字がついています。 したがって,項目ごとに「どちらでもない」および極端な選択肢の選ばれやすさが異なる,という設定がされています。

同様に,「好き(4)」を選ぶ確率は,「『どちらでもない』を選ばない」かつ「好き」かつ「『とても』ではない」の積なので P(ypi=4|θp(M),θp(A)θp(E))=(1𝚽[αi(M)(θp(M)βi(M))])𝚽[αi(A)(θp(A)βi(A))](1𝚽[αi(E)(θp(E)βi(E))])(8.67) \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} \qquad(8.67) という形で表現できるわけです。

IRTreeの目的あるいは利点は,「測定したい回答者の特性値」と「回答傾向」を分離して考えられる点です。 一般的に,リッカート尺度への回答には個人のクセのようなものが混ざっていることが知られています。 そのせいで,「どちらでもない(3)」を選んだ人が「好きと嫌いのちょうど中間」なのか「本当は『どちらかといえば好き』くらいなんだけど,その程度では『好き』は選びづらい」のかは区別がつきにくいものです。 IRTreeモデルならば,上の例で言えば中心傾向と極端傾向の強さを,個人特性θp\theta_pの一種として考えることができます。 そのため,残ったθp(A)\theta_p^{(A)}より純粋な「測定したい回答者の特性値」として解釈できるかもしれないというわけです。 当然ながら,IRTreeモデルでは 図 8.45 に示したような回答決定プロセスの仮定がかなり重要になってきます。 逆に言えば,複数のツリー構造間でモデル比較をすることで,どのプロセスが最も当てはまりそうか,ということを評価したりも出来るかもしれません…。

8.12.2 強制選択形式に対するIRTモデル

IRTreeモデルはいわば,得られた回答データから回答傾向を統計的に分離しようと頑張っているモデルの一種です。 そもそもリッカート尺度では回答傾向の影響を逃れることができないわけですが,その対抗策としては他にも「そもそも中間選択肢をなくそう」という方向を考えることができます。 また,入社試験などの重要な局面では「ウソをついて本当の自分より良い人に見せよう」という思惑が働いて回答が歪められる可能性があったりします (フェイキング: e.g., Jackson et al., 2000。 これを防ぐためには,例えば「友だちが多い」「真面目である」といった複数の短文を同時に提示し,その中から一つだけを選択させたり,ランキングを付けさせるという回答方法が考えられます。 これならば,前述の「中間選択肢を選びがち」といった問題も防げますし,フェイキングの問題もかなり抑制できるでしょう。

複数の選択肢から一つの選択肢を選ばせる回答方法は,強制選択 (Forced Choice) 型項目などと呼ばれます。 また,特に選択肢数が2個の場合は一対比較 (Paired Comparison) と呼ばれることもあります。 そしてこれらの出題形式に対しては,やはり普通のIRTとは少し異なるモデルが考えられています。 最も代表的なモデルが,Thurstonian IRT (TIRT) モデル Brown & Maydeu-Olivares, 2011でしょう。 TIRTモデルでは,各選択肢がそれぞれ独立に「効用」を持っていると考えます。 この「効用」とは,経済学などで出てくるものと同じようなものだと考えてOKです。 あるいは単純に「各選択肢が自分に当てはまっていると感じる程度」というイメージでも構いません。 とにかく強制選択型項目への回答では,回答者は効用を最大化する選択肢を選ぶのだ,と考えているわけです。 ここでは,回答者ppが項目iiの2つの選択肢A,BA, Bに対してそれぞれ効用upi(A),upi(B)u_{pi}^{(A)},u_{pi}^{(B)}を持っているとしましょう。 すると,回答者ppは効用を最大化する選択肢を選ぶため,観測されるデータとしては ypi={Aifupi(A)upi(B)0Bifupi(A)upi(B)<0(8.68)\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} \qquad(8.68) と表わせます。

TIRTモデルでは,選択肢の効用upi(A)u_{pi}^{(A)}を,よくある因子分析モデルの形 upi(A)=τi(A)+fp(A)bi(A)+εpi(A),εpi(A)N(0,σεpi(A))(8.69) 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)}}) \qquad(8.69) と定義します。つまり,平均的な効用τi(A)\tau_{i}^{(A)}と,その選択肢が反映している因子の得点fi(A)f_{i}^{(A)}と因子負荷bp(A)b_p^{(A)}の和によって表される効用の期待値に,そこに正規分布に従う誤差(独自因子得点)εpi(A)\varepsilon_{pi}^{(A)}が加わって効用の実現値が確率的に決定すると考えているわけです。 すると,2つの選択肢A,BA, Bの中からAAを選ぶ確率は最終的に P(ypi=A)=P(upi(A)upi(B)0)=𝚽[αi(AB)+(fi(A)βp(A)fi(B)βp(B))](8.70) \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} \qquad(8.70) と表されます31。ただし表記を簡略化するため αi(AB)=τi(A)τi(A)σεpi(A)+σεpi(B),βi(A)=bi(A)σεpi(A)+σεpi(B),βi(B)=bi(B)σεpi(B)+σεpi(B)(8.71) \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} \qquad(8.71) という形で別の記号を用いて表現しています。 つまりTIRTモデルは,一対比較の場合には,各項目が2つの因子のみに負荷量を持っていると設定したCFAと実質的には同じことをしています。 そして一つの項目に3つ以上の選択肢がある場合は,対ごとに異なる項目を構成しているとみなして分析を行います。 例えば項目iiA,B,C,DA, B, C, Dの4択だとすると,この項目は「AとB」「AとC」「AとD」「BとC」「BとD」「CとD」という6つの二値項目があると考えます。 そして回答者ppが選択肢AAを選んだとすると,「AとB」「AとC」「AとD」については(いずれもAを選んだという)結果が観測され,残りの「BとC」「BとD」「CとD」については欠測扱いにします。 そして観測された3項目についてはそれぞれ P(ypi=A|A,B)=𝚽[αi(AB)+(fi(A)βp(A)fi(B)βp(B))]P(ypi=A|A,C)=𝚽[αi(AC)+(fi(A)βp(A)fi(C)βp(C))]P(ypi=A|A,D)=𝚽[αi(AD)+(fi(A)βp(A)fi(D)βp(D))](8.72) \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} \qquad(8.72) という形でそれぞれ選択肢Aに関する部分が共通しているという設定のもとで定式化されるわけです。

8.12.3 回答時間を利用するIRTモデル

コンピュータによる測定 (Computer Based Testing [CBT]) では,項目反応以外の情報を得ることが容易になりました。 もっとも代表的な情報として,回答時間 (Response Time [RT]) が挙げられるでしょう。 Section 3.2.2 では,回答時間を「適当な回答をしている人を検出する」といった目的で使用しましたが,その一方で回答者の特性値推定に利用するという方向性も相当に研究されています。 そういうわけで,回答時間を利用するIRTモデルを挙げていくとキリがないのですが,ここではシンプルなモデルとして,「当てはまる」「当てはまらない」の二択の項目に対するモデルFerrando & Lorenzo-Seva, 2007を紹介します。

回答時間を利用するIRTモデルでは,そもそも回答者の特性値と回答時間にはどのような関係があるかを考えるところから始まります。 性格特性を尋ねる項目の場合には,「特性の強さが中途半端な人ほど回答時間が長い」という逆U字の関係になることが経験的に知られています。 例えば「友だちが多い」という項目に対して,友だちが多い/少ない自信がある人,つまり「外向性」といった特性が極端な人は,きっとすぐに選択肢を選ぶことでしょう。 これに対して,『友だちはそこそこ居るけど』という人は,きっと「当てはまる」「当てはまらない」のどちらを選ぶかに迷いが生じるはずです。 したがって,同じ「当てはまる」を選んだ人同士であっても,素早く「当てはまる」を選んだ人のほうがより強い特性を持っていると考えるのが自然です。

ということで,この関係を尤度関数に組み込んでIRTモデルを構築していきます。 と言っても考え方は簡単で,回答時間tpit_{pi}についても項目反応ypiy_{pi}の(8.1)式と同じように tpi=g(ωp,γi)(8.73) t_{pi} = g(\omega_p, \gamma_i) \qquad(8.73) という形で,回答者の要因ωp\omega_pと項目の要因γi\gamma_iの交互作用として考えます。 Ferrando & Lorenzo-Seva (2007 では,これを具体的に ln(tpi)=μ+ωp+γi+bδpi+εpi,εpiN(0,σε2)δpi=αi2(θpβi)2(8.74) \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} \qquad(8.74) と置きました。 ここでμ\muは全体的な切片項,ωp\omega_pは回答者ppの回答スピード,γi\gamma_iは項目iiの所要時間を表現したパラメータです。 また,δpi\delta_{pi}の中身には(θpβi)2(\theta_p - \beta_i)^2というものが入っています。この項は,θp=βi\theta_p = \beta_iのとき,つまり反応確率が五分五分のときに最小値を取ります。 係数bbが負の値であれば,回答者の特性値θp\theta_pが項目困難度βi\beta_iに近いほど回答時間の期待値が長くなることで,逆U字の関係を表現しているわけです。

回答時間などの追加情報を用いて特性値を推定するIRTモデルでは,うまく行けばより精度の高い,あるいは妥当性の高い推定を行うことが出来ると期待されています。 ちなみに,IRTに限らず社会調査などの分野では,このように回答内容以外の情報(位置情報や回答デバイスの情報などを含む)をパラデータ Kreuter, 2013 と呼んでいます。 パラデータは,やはり適当な回答の検出に使われるだけでなく,よりノイズの少ないデータを収集するための調査改善など様々な用途で利用されていたりします。

8.12.4 展開型IRTモデル

これまでに紹介したIRTモデルはすべて,θ\thetaの値が高ければ高いほど高い点数をとる(i.e., 「当てはまる」寄りの回答をする,正解する)という仮定がされていました。 項目特性曲線にしても,(αi>0\alpha_i>0ならば)右上がり(単調増加)でした。 ですが質問のタイプによっては「過ぎたるはなお及ばざるが如し」的なものが存在しています。 よくあるのは「態度」に関する項目です。 例えば「消費税を5%にするべきだ」という項目は,「消費税は完全に撤廃してほしい」人も「消費税はこのままで良い」人も「そう思わない」を選択するでしょう。 同様に,心理特性に関しても「カフェで友人と静かに会話をするのは楽しい(因子:外向性)」Stark et al., 2006という項目に対しては「そもそも他人と話したくない人」も「もっとワイワイ騒ぎたい人」も「当てはまらない」を選ぶと思われます。 つまりこうした項目における項目特性曲線は, 図 8.46 のように,あるところを境に単調減少に転じると考えられるわけです。

図 8.46: 中間にピークが来るICC

ロジスティック・正規累積モデルでは,「項目(βi\beta_i)」と「回答者(θp\theta_p)」は共通の軸の上でマッピングされていると考えることが出来ます。 すごく簡単にいえば「困難度がβi\beta_iの項目はθp=βi\theta_p=\beta_iの回答者にとって(易しすぎず難しすぎず)ちょうどよい難易度(50%)の項目である」ということです。 これを先程の項目に置き換えて考えると,βi\beta_iθp\theta_pの距離が近いほど,その人によってちょうどよい程度の(理想の)意見・態度を表している」ということができそうです。 このような「理想点(ideal point)」の考え方に基づくIRTモデルを展開型 (unfolding) IRTモデルと呼びます32

展開型IRTモデルでは,反応確率P(ypi=1)P(y_{pi}=1)が「βi\beta_iθp\theta_pの距離」によって規定されます。 したがって,最も基本的な一次元・二値モデルの場合は P(ypi=1)=exp[12(βiαiθp)2](8.75) P(y_{pi}=1) = \exp\left[-\frac{1}{2}(\beta_i-\alpha_i\theta_p)^2\right] \qquad(8.75) と表されますMaydeu-olivares et al., 2006。 ここでも2つの項目パラメータの意味は概ね通常の場合と同じです。 βi\beta_iはその項目の「位置」を表しており,θp\theta_pがこの値に近い人ほど「当てはまる」と回答する確率が高くなることがわかります。 αi\alpha_iは項目識別力を表していますが,展開型モデルの場合は,正規分布の分散パラメータのように山の狭さを規定します( 図 8.47 )。 ですが結局「θp\theta_pが理想点βi\beta_i付近か否かを区別する強さ」という意味では,これまでに出てきた項目識別力と同じような役割を果たしているわけです。

図 8.47: 識別力の役割(展開型)

mirtでは,itemtype="ideal"と指定することで展開型モデルを実行可能です。

コード 8.57: mirtで展開型モデル
result_ideal <- mirt(dat_bin, model = 1, itemtype = "ideal")

通常のIRTモデルと同じようにcoef()でパラメータの推定値を見ることもできるのですが,どうやらIRTpars = TRUEには対応していないようなのでご注意ください。

コード 8.58: 展開型モデルの推定値のチェック
# IRTpars = TRUEとするとエラーが出る
coef(result_ideal, simplify = TRUE)
$items
         a1      d
Q1_1  0.118 -0.716
Q1_2  0.214 -0.468
Q1_3  0.250 -0.571
Q1_4  0.287 -0.609
Q1_5  0.255 -0.600
Q1_6  0.618 -0.312
Q1_7  0.651 -0.455
Q1_8  0.621 -0.489
Q1_9  0.571 -0.653
Q1_10 0.667 -1.187

$means
F1 
 0 

$cov
   F1
F1  1

また,ICCのプロットも同じように実行可能です。

コード 8.59: 項目特性曲線(展開型)
plot(result_ideal, type = "trace", facet_items = FALSE)
図 8.48: 項目特性曲線(展開型)

例えば1番目の(一番山が緩い)項目は,(αi,βi)=(0.118,0.716)(\alpha_i,\beta_i)=(0.118,-0.716)ということで,θ=0.716/0.1186.068\theta=0.716/0.118\simeq 6.068のあたりで最大値P(ypi=1)=1P(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 等化

Section 8.1 では,IRTを項目依存性と集団依存性をクリアするための理論だと紹介しました。 ですが実際には,IRTモデル(2パラメータのGRM)はカテゴリカル因子分析と同じなのでθp\theta_pはデータをとった集団ごとに標準化されています。 つまり最初に例示したような,異なる集団に対して実施した異なる尺度の比較はこのままではできません。 もちろん反対に,異なる尺度を用いて実施した異なる集団同士の比較もこのままではできません。 こうした比較を可能にするために,IRTでは等化 (equating) またはリンキング (linking) という手続き Kolen & Brennan, 2014を行います34。 等化の技術は,TOEICのような大規模テストでも実際に利用されています。 そのおかげで,毎回異なるはずの難易度によらず「TOEIC ***点」というように点数を履歴書に書けるわけです。

(正直言うと,等化=異なる集団や尺度の比較を考えない限りは,特に心理尺度においてわざわざ因子分析やSEMではなくIRTを選ぶ理由はほぼ無い気すらします)

等化の基本的な考え方を説明しましょう。 ここからはIRTで推定されるパラメータが集団によらない項目パラメータ項目によらない個人特性だという仮定のもとで話を進めます。

共通回答者による等化

まず,個人特性値が「項目によらない」のであれば,どんな尺度で推定した場合でもθp\theta_pの値は同じになるはずです。 いま,ppさんが同じ構成概念を測定する尺度11と尺度22による2回の調査の両方に回答したとします。 この2回のデータをもとにパラメータを推定する場合には,いままでの手順にのっとってそれぞれθpN(0,1)\theta_p\sim N(0,1)と仮定して推定を行います。 その結果

  • 尺度1での推定値θ̂p(1)=0.5\hat{\theta}_p^{(1)}=0.5
  • 尺度2での推定値θ̂p(2)=0.5\hat{\theta}_p^{(2)}=-0.5

だったとしましょう。 仮に尺度1による推定が「項目によらない個人特性」の値とすると,尺度2での推定値は本来の値より1小さい値になっているようです。 ppさんの特性値が尺度1回答時と尺度2回答時で変わっていないとしたら,この1のズレは「尺度2は尺度1よりも特性値を小さめに推定する尺度である」ことが原因と考えられます。 そこで,θ̂p(2)=0.5\hat{\theta}_p^{(2)}=-0.5がどうにかしてθ̂p(1)\hat{\theta}_p^{(1)}と同じ0.50.5になるように調整する(ゲタを履かせる)必要がありそうです。

もう少し規模を広げて,同じようにこの尺度1と尺度2による2回の調査の両方に回答した人が100人いたとします。 その結果,2回の推定値の散布図が 図 8.49 のようになったとします。 そしてこの2回の推定値の回帰直線を求めたら θ̂p(2)=Aθ̂p(1)+B(8.76) \hat{\theta}_p^{(2)} = A\hat{\theta}_p^{(1)} + B \qquad(8.76) であったとします。 この場合,全員のθ̂p(2)\hat{\theta}_p^{(2)}をなるべくθ̂p(1)\hat{\theta}_p^{(1)}に近づくように調整した値θ̂p(21)\hat{\theta}_p^{(2-1)}は, θ̂p(21)=θ̂p(2)BA(8.77) \hat{\theta}_p^{(2-1)} = \frac{\hat{\theta}_p^{(2)}-B}{A} \qquad(8.77) という形で求められそうだとわかります。

図 8.49: 2回の測定における推定値

2パラメータロジスティックモデルでは,項目反応関数は P(ypi=1)=11+exp[αi(θpβi)](8.78) \begin{aligned} P(y_{pi}=1)&= \frac{1}{1+\exp\left[-\alpha_{i}(\theta_p - \beta_i)\right]} \end{aligned} \qquad(8.78) でした。ここで, θp*=θpBAβi*=βiBAαi*=Aαi(8.79) \begin{aligned} \theta_p^{*} &= \frac{\theta_p - B}{A} \\ \beta_i^{*} &= \frac{\beta_i - B}{A} \\ \alpha_i^{*} &= A\alpha_i \end{aligned} \qquad(8.79) とおくと, P(ypi=1)=11+exp[αi*(θp*βi*)]=11+exp[Aαi([θpBA][βiBA])]=11+exp[αi(θpβi)](8.80) \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} \qquad(8.80) となり,もとの項目反応関数と同じになります。 探索的因子分析では,潜在変数のスケールの不定性の問題から,全ての観測変数を標準化していました。 等化では,このスケールの不定性をある種逆手に取ることでα,β\alpha, \betaθ\thetaを適当なスケールに変換しているわけです。 具体的には,2つの尺度に両方とも回答した全員のθ̂p(2)\hat{\theta}_p^{(2)}がなるべくθ̂p(1)\hat{\theta}_p^{(1)}に近づくようにスケールを変換し,これに応じて項目パラメータにも変換をかけます。 こうして得られた変換後の項目パラメータならば,尺度Aと尺度Bでも問題なく比較ができるようになります。

共通項目による等化

先程説明した等化の流れは,共通の項目がある場合でも同じように考えられます。 つまり,2つの集団に対して実施する際に,それぞれの尺度に共通の項目をいくつか入れておけば, βi(21)=βi(2)BAαi(21)=Aαi(2)(8.81) \begin{aligned} \beta_i^{(2-1)} &= \frac{\beta_i^{(2)} - B}{A} \\ \alpha_i^{(2-1)} &= A\alpha_i^{(2)} \end{aligned} \qquad(8.81) となるような等化係数A,BA,Bを求めることができ,それに応じた変換が可能となります。 等化係数A,BA,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) では,コンピュータの力を使って,例えば「θ=(4,2,0,2,4)\theta=(-4,-2,0,2,4)の各点における標準誤差が0.2以下=テスト情報量が25以上になる最小項目数のセット」などの条件にあった最適な項目セットを探し出してくれます35

ただこれだけのテストの設計を行うためにはかなりの事前準備が必要となります。 というのも項目情報関数は項目パラメータが既知でないとわかりません。 また,目的に応じて最適な項目セットを作成するために,なるべく候補となる項目の数は多いほうが良いです。 もちろんすべての項目のパラメータを一気に推定するわけにはいかないので,多くの場合は複数回に分けて集めたデータをもとにパラメータを等化したりして項目プールを作成する必要があるでしょう。 ということで心理尺度の場合にはなかなかそこまで大量の項目を用意すること自体が難しいので,ここまで理論的な尺度構成を行っている先行事例はたぶんほとんど存在しません。 一応用途としては,100-200項目くらい用意した候補の中から「最高の20項目」を選んで心理尺度を作成する,といった方法はありそうな気がします36

8.13.3 適応型テスト

テストや尺度の中には通常様々な困難度の項目が混ざっているため,ほぼ確実に「その人にはほぼ意味のない項目」が含まれることになります。 この問題は,先ほど紹介したテストの設計の考え方を持ってしても避けられません。 そこで自動テスト構成をもっと局所的に行うと「個人に応じて最適な項目セットを選び出す」という考えになります。 つまり視力検査のようにθp\theta_pが高そうな人にはより困難度の高い項目を,θp\theta_pが低そうな人にはより困難度の低い項目を出していけば,効率的にテスト情報量を稼ぐことができそうです。 コンピュータ上で行われるコンピュータ適応型テスト (computerized adaptive testing [CAT]) 37では,この考えに基づいて個人ごとに次の項目を出し分けていくことができます。

ここではIRTに基づく最もシンプルな項目選択法を紹介しましょう。 推定の標準誤差は σ(θ̂p|θp)=1I(θp) \sigma_{(\hat{\theta}_p|\theta_p)}=\frac{1}{\sqrt{I(\theta_p)}} で表されていました。 したがって推定精度をなるべく上げるためには,I(θp)I(\theta_p)を大きくする=項目情報量Ii(θp)I_i(\theta_p)が大きい項目を選んでいけばよいわけです。 ただしIi(θp)I_i(\theta_p)は「真値θp\theta_pにおける項目情報量」のため,真値θp\theta_pがわからない状況ではIi(θp)I_i(\theta_p)の値もわかりません。 そこで,1問解答するたびに推定を行い,その時点での暫定的な推定値θ̂p*\hat{\theta}_p^*における項目情報量Ii(θ̂p*)I_i(\hat{\theta}_p^*)が最大となる項目をとりあえず選び続けたらなんだかうまく行きそうです38

CATはうまく行けば従来式と同程度の推定精度を半分以下の項目数で達成できたりもする,かなり夢のある技術です。 ただすべての回答者にとって最適な項目を用意できないと効果が薄れてしまうため,単純に自動テスト構成するときよりももっと大量の項目を用意しておく必要があります。 また,テストを実施しリアルタイムで次の項目を選択するためのシステムの構築など,実装に向けたハードルはかなり高いものです。 とりあえずは「そういう技術もIRTによって実現できるんだなぁ」くらいに思っておいていただければ良いと思います39


  1. Choi & Asilkalkan (2019 によれば,IRTができるRパッケージは45個もあるらしいです(今はもっと多い?)。細かく目的や用途が異なると思うので,特定の場面ではもっと使い勝手の良いパッケージがあったりするかもしれません。あとは計算の精度や速度の違いなど?↩︎

  2. 「項目反応理論」というと,ほとんどの場合この後紹介する代表的なモデルを指すのですが,モデル式(項目反応関数)がそれだけしか無いのではなく,実際には様々なモデル式による形が考えられるというのは重要なポイントです。「項目反応理論」という言葉はあくまでも項目依存性や集団依存性を解消したパラメータを得るための理論的な方法を考える土台としての理論である,という認識です。↩︎

  3. 頻度論的に言えば,ここでの「確率」は「あるθp\theta_pの一個人が何回も回答した際の当てはまるの割合」や「母集団内でのθp\theta_pの人たちにおける当てはまるの割合」と解釈できます。これはどちらでもOKです。↩︎

  4. もとが学力テストなので「困難度」という呼び方をしています。正解・不正解がない心理尺度でも特に気にせず「困難度」と呼びます。↩︎

  5. もっと近づけるためには1.704倍のほうがよかったり,ズレを極限まで減らすためには0.07056x3+1.5976x0.07056x^3+1.5976xにすると良いBowling et al., 2009ということは分かっていますが,そこまでの精度は必要ない,ということで,わかりやすい1.7が用いられています。↩︎

  6. 正規累積モデルはプロビットモデル,ロジスティックモデルはロジットモデルに対応します。…と考えると「本当は理論的にはどちらを使ったほうが良いか」が明確なケースもあるのかもしれません。実用上はどちらでも良いとされていますが,本当にそのまま1.7倍の変換をするだけで交換可能なのか,というところについてはやや議論があるようです Cho, 2023↩︎

  7. IRTの説明に合わせて誤差のεpi\varepsilon_{pi}の符号を変えています。↩︎

  8. 厳密には,「独立」と「無相関」は別の概念です。理想は独立なのですが,実用上は無相関であれば良いだろう,という感じになっています。厳密な独立とは異なるという意味で,弱局所独立性と呼ぶこともあります。↩︎

  9. θp\theta_pで条件づけた」ということが明確になるようにここだけP(.|θp)P(.|\theta_p)としていますが,これまでのP(.)P(.)と意味は同じです。↩︎

  10. 後半では少しだけ多次元モデルの紹介をしますが,これは一次元モデルの拡張として展開されるものであり,IRTの基本は一次元だといえます。↩︎

  11. もちろん,多因子モデルの場合には,一次元とは限りません。ただ同様に,局所独立性が満たされているならばnn次元性が満たされている,ということはできると思います。↩︎

  12. 問1の答えを使わないと問2が解けないような場合,局所独立性が満たされなくなってしまいます。そこで「問1が解けたら1点,問2も解けたら2点」というように得点化してこれを一つの大きな問題として考えることで局所独立性を確保しようという考えに至るわけです。↩︎

  13. P(ypi=1)=exp(πpi1)P(ypi=0)P(y_{pi}=1)=\exp(\pi_{pi1})P(y_{pi}=0)と表せますが,カテゴリk=0k=0は存在しない,つまりP(ypi=0)=0P(y_{pi}=0)=0です。ということはP(ypi*=1)0,1=P(y^*_{pi}=1)_{0,1}=\inftyとなるため,exp(πpi1)=\exp(\pi_{pi1})=\inftyとなってしまい計算ができないのです。↩︎

  14. 完全には非補償ではない(ある程度なら別の次元の特性値でリカバーできる)という意味で「部分補償モデル」のほうがより正しい名前らしいですが,そこまでこだわるほどの違いはありません。↩︎

  15. あるかはわかりませんが,例えば「2つの特性を両方とも持っている人だけ当てはまると回答する尺度・項目」のようなものを考えた場合には,普通の因子分析よりも非補償型のモデルを適用したほうが良い(回答メカニズムに則している)可能性があります。↩︎

  16. 非補償型モデルは,二値ならばitemtype='PC2PL'またはitemtype='PC3PL'で実行可能です(PCはpartially compensatoryの頭文字)。ただパラメータの推定は結構安定しないので利用する際は気をつけてください。↩︎

  17. psych::fa()のように固有値分解など使えばいいじゃないと思われるかもしれませんが,mirt()関数として他の様々なIRTモデルと同じ文法で実装するためにはこうするのが最も自然なのです。↩︎

  18. NOHARMという名のソフトウェアで出力される推定結果をもとに計算する方法なので,「NOHARMに基づく方法」などと言ったほうが正確かもしれません。↩︎

  19. これは厳密には弱局所独立性 (weak liocl independence)と呼ばれるものです。一方で,ある反応ベクトル𝐲p\symbf{y}_pが得られる確率が各項目の反応確率ypiy_{pi}の積に等しい,という形でも局所独立性は表されるのですが,こちらは強局所独立性 (strong liocl independence)と呼ばれています。本来はこの強局所独立性が必要なのですが,現実世界の多くの場面では,弱局所独立性が満たされている場合には強局所独立性も満たされていることが多いということで,ここでは弱局所独立性の検証を行っていきます。↩︎

  20. DIMTESTについては関数が見つけられていません。一応DIM-Packというフリーソフトを使えば出来そうですが…↩︎

  21. 正確に言えば,「𝐲p\symbf{y}_pという回答をした人の特性値としてθ̂p\hat{\theta}_pという値がどの程度もっともらしいか」という感じですが,あえて曲解させています。↩︎

  22. outfitは重み付けない二乗和,infitは情報量(後述)で重み付けた二乗和のようです。↩︎

  23. もちろんモデル全体での適合度を考えることも可能です。例えば1パラメータと2パラメータの両モデルをlavaanで実行して,compareFit()することでどちらのモデルが良さそうかを判断したりというプロセスが考えられます。↩︎

  24. もちろんこの回答者ppが赤い項目に正解した場合には,θp\theta_pの予測を大きく動かすため,情報量が多くなります。ですがそもそも「回答者ppが赤い項目に正解する」という事象の発生確率がほぼゼロなので,「この項目への回答がθp\theta_pの予測にもたらす情報量の期待値」という観点ではやはりほぼゼロなのです。↩︎

  25. 項目情報「量」という場合には特定のθ\thetaにおける値を指す一方で,項目情報「関数」という場合にはθ\theta全域を対象とすることが多い気がするので,使い分けに気をつけてください。↩︎

  26. MH検定はmirtパッケージには用意されていないようです。RではdifR::difMH()などの関数によってMH検定を行うことができます。↩︎

  27. 項目の数だけ検定を繰り返しているので,実はここで表示されている結果には多重検定の問題が生じています。これを解消するためには,引数p.adjustを指定すると良いです。例えばp.adjust = "holm"とするとHolmの方法によって調整されたpp値を表示してくれるようになります。↩︎

  28. この「テスト全体でのグループ間の挙動の違い」はDIFと同じように特異テスト機能 (Differentian Test Functioning [DTF])と呼ばれたりします。↩︎

  29. ここで紹介するモデルを始め,多くの(特殊な)IRTモデルはmirtパッケージに限らずRのパッケージとしては実装されていなかったりします。そういったモデルについては,自分で尤度関数をプログラムしてoptim()関数による数値最適化を利用する,あるいはstanなどを利用してマルコフ連鎖モンテカルロ法(MCMC)によって推定する,といった方法を取る必要があります。stanを使うためにはまずベイズ統計学を理解する必要があるのがちょっと大変ですが,慣れるとほぼなんでもできるようになります。↩︎

  30. θp\theta_pの上付き添字のM, A, EはそれぞれMiddle/Mid-point, Attitude/Agree, Extremeなどの頭文字です。↩︎

  31. 2つの選択肢の比較のときに,その背後にある効用を比較していて,しかも効用が確率的に変動する,という考え方をサーストンの比較判断の原則 (law of comparative judgment: Thurstone, 1927と呼ぶため,このモデルはThurstonian IRTと名付けられています。↩︎

  32. これに対して通常のIRTモデルはdominance modelと呼びます(日本語訳がわからない…)。↩︎

  33. Rでの実践例までやると膨大になってしまうので,概念の紹介だけにとどめます。気になる人は自分で調べたり直接質問してください。↩︎

  34. RではequateIRTplinkといったパッケージによって等化のメソッドが提供されています。↩︎

  35. RではlpSolveglpkAPIRsymphonyなどのパッケージを使うと自動テスト構成ができるようになります。↩︎

  36. 100-200項目すべてに真面目に回答してくれる人たちを数百人用意する必要があるのですが,それさえクリアできれば…↩︎

  37. RではcatRmirtCATといったパッケージが役に立つかもしれません。↩︎

  38. もちろん他にも項目の選び方は色々あります。ざっくりいうと「推定精度をなるべく高める」方向性と「なるべく人によって異なる項目が出るようにする(項目情報が外部に漏れるのを防ぐため)」という方向性があります。↩︎

  39. IRTによらないCATもあることにはあります。ただIRTに基づくCATのほうが圧倒的マジョリティです。↩︎