多変量解析_08_項目反応理論

Chapter 8 項目反応理論 Chapter 8 項目反応理論 本資料は，神戸大学経営学研究科で 2022 年度より担当している「統計的方法論特殊研究（多 変量解析）」の講義資料です。CC BY-NC 4.0 ライセンスの下に提供されています。 作成者連絡先 神戸大学大学院経営学研究科 分寺杏介（ぶんじ・きょうすけ） mail: bunji@bear.kobe-u.ac.jp website: https://www2.kobe-u.ac.jp/~bunji/ 最後は項目反応理論 (Item Response Theory [IRT]; Lord and Novick, 1968) のお話です。 Contents 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 8.10 8.11 8.12 項目反応理論とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 基本的なモデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 因子分析との関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . R で IRT . . . IRT の前提条件 多値型モデル . . 多次元モデル . . 適合度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 項目情報量・テスト情報量 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . （おまけ）多母集団モデルと特異項目機能 . いろいろな IRT モデル . . . . . . . . . . IRT の使い道 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 3 10 11 16 19 28 35 49 57 72 79 後ほど説明しますが，実は IRT の最もメジャーなモデルの一つは，因子分析と数学的に同じ ものです。そのため lavaan でも IRT 分析はできるのですが，IRT 専用のパッケージも色々 あります*1 。この授業では mirt というパッケージを使用します。 *1 1 Choi and Asilkalkan (2019) によれば，IRT の R パッケージは 45 個もあるらしいです。細かく目的や用 途が異なると思うので，特定の場面ではもっと使い勝手の良いパッケージがあったりするかもしれません。あ とは計算の精度や速度の違いなど？

• https://creativecommons.org/licenses/by-nc/4.0/
• https://www2.kobe-u.ac.jp/~bunji/

8.1 項目反応理論とは 事前準備 1 # install.packages("mirt") 2 library(mirt) 3 dat <- readRDS("chapter04.rds") 4 # 毎回列名を指定するのが面倒なので先に作っておく 5 cols <- paste0("Q1_", 1:10) 8.1 項目反応理論とは 第 1 回では，バーンアウトを例に心理尺度のお話をしました。その時には，現在では 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 人の標準化得点がそ れぞれ 0.2 点と-0.25 点だったとしましょう。2 人はそれぞれ平均点以上・以下だったので， もちろん標準化得点でも A さんのほうが高い得点になっています。ただ，実はまだこの情報 だけでは A さんのほうがよりバーンアウトに近いとは言い切れません。それは，2人の所属 するグループの平均値が同じとは限らないためです。 標準化得点は，当然ながらそのグループの得点分布に依存します。もしかしたら，UWES の 回答者平均点のほうが MBI-GS よりも高い理由は，A さんの所属するグループよりも B さ んの所属するグループの方が平均的にバーンアウト傾向が強いためかもしれません。このよ うに，心理尺度による評価は項目だけでなく集団にも依存します（集団依存性）。 IRT では，項目依存性と集団依存性をクリアした得点を算出するために集団によらない項目 パラメータと項目によらない個人特性値の存在を考えています。仮にそのような項目パラメ 2

8.2 基本的なモデル ータが 𝛽 ，個人特性値が 𝜃 で表せるとしたら，項目 𝑖 に対する回答者 𝑝 の回答は 𝑦𝑝𝑖 = 𝑓(𝜃𝑝 , 𝛽𝑖 ), (8.1) つまり 𝛽𝑖 と 𝜃𝑝 による何らかの関数（項目反応関数; item response function [IRF]）にな るはずです。後ほど IRT の具体的なモデルを紹介しますが，IRT の本質はこのように，回答 を (集団によらない) 項目パラメータと (項目によらない) 個人特性値に分離して分析するこ とにあります*2 。 8.2 基本的なモデル 最もシンプルな IRT モデルを考えるため，一旦ここからは項目への反応を二値（当てはまら ない・あてはまる）に限定して考えてみます。というのも，もともと IRT は「正解・不正解」 の二値で考えるような学力テストなどを想定して作られているもののためです。 8.2.1 正規累積モデル 以前カテゴリカルデータの相関（ポリコリック・ポリシリアル相関）の計算時には，背後に ある連続量を考えるということを説明しました。ここでも，二値反応の場合 𝑦𝑝𝑖 の背後に連 続量 𝜃𝑝 が存在し，𝜃𝑝 が閾値を超えていたら 1，超えていなければ 0 と回答すると考えます。 そして項目の閾値を 𝛽𝑖 に設定してみると， 𝑦𝑝𝑖 = { 0 𝜃𝑝 < 𝛽𝑖 1 𝜃𝑝 ≥ 𝛽𝑖 (8.2) 。𝛽𝑖 の値が大きいほど，𝑦𝑝𝑖 = 1 になる 𝜃𝑝 の範囲が狭くなるため，多く となります（図8.1） の人がその項目には「当てはまらない」と回答するだろう，ということです。 q p ~ N(0, 1) ypi = 0 ypi = 1 bi 図8.1 二値反応と潜在的な連続量の関係 古典的テスト理論でも回帰分析でも因子分析でも，観測値は必ず何らかの誤差を伴って発生 するものです。そこで先程の(8.2)式に誤差 𝜀𝑝𝑖 を追加して 𝑦𝑝𝑖 = { 0 (𝜃𝑝 − 𝜀𝑝𝑖 ) < 𝛽𝑖 1 (𝜃𝑝 − 𝜀𝑝𝑖 ) ≥ 𝛽𝑖 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀2 ) (8.3) とします。誤差を追加することによって，図8.2のように「特性値が 𝜃𝑝 の人が項目 𝑗 に『当 てはまる』と回答する確率*3 」を考えることが出来ます。 *2 3 IRT モデルで分析をしたら一発 OK，というわけではありません。項目依存性や集団依存性を解消したパラ

8.2 基本的なモデル N(q p, s2e ) P(ypi = 0) P(ypi = 1) qp bi 図8.2 誤差を追加すると さらに 𝛽𝑖 と 𝜃𝑝 を一つの関数 𝑓(𝜃𝑝 , 𝛽𝑖 ) として見るために 𝛽𝑖 を移項して， 𝑦𝑝𝑖 = { 0 (𝜃𝑝 − 𝛽𝑖 − 𝜀𝑝𝑖 ) < 0 1 (𝜃𝑝 − 𝛽𝑖 − 𝜀𝑝𝑖 ) ≥ 0 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀2 ) (8.4) とします（図8.3）。当然ながら，図8.2と図8.3は，閾値と分布が平行移動しただけなので 𝑃 (𝑦𝑝𝑖 = 1) は変わりません。 N(q p - b i, s2e ) P(ypi = 1) P(ypi = 0) 0 q p - bi 図8.3 𝛽𝑖 を移項すると 図8.3の色つきの確率分布は，正規分布 𝑁 (𝜃𝑝 − 𝛽𝑖 , 𝜎𝜀 ) と表すことができます。そして式から， 𝑦𝑝𝑖 = 1 となるのは誤差 𝜀𝑝𝑖 が (𝜃𝑝 − 𝛽𝑖 ) 以下のときだとわかります。したがって，𝜀𝑝𝑖 が標 準正規分布に従う (𝜎𝜀2 = 1) とすると，「特性値が 𝜃𝑝 の人が項目 𝑗 に当てはまると回答する 確率」（図8.4の青い部分）は 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 (𝜀𝑝𝑖 ≤ 𝜃𝑝 − 𝛽𝑖 ) (𝜃𝑝 −𝛽𝑖 ) =∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.5) と表すことができます。教科書などで見たことがある形かもしれませんが，この積分の中身 は標準正規分布の確率密度関数です。このように，𝑃 (𝑦𝑝𝑖 = 1) を標準正規分布の累積確率で 表すモデルのことを正規累積モデル (normal ogive model) と呼びます。この後項目パラメ ータを 1 つ増やすので，それに対してこの正規累積モデルは，特に 1 パラメータ正規累積モ デルと呼ばれます。 *3 4 メータを得るための理論的な方法を考える土台として「項目反応理論」というものがある，という感じです。 頻度論的に言えば，ここでの「確率」は「ある 𝜃𝑝 の一個人が何回も回答した際の当てはまるの割合」や「母 集団内での 𝜃𝑝 の人たちにおける当てはまるの割合」と解釈できます。これはどちらでも OK です。

8.2 基本的なモデル e pi ~ N(0, 1) P(ypi = 1) P(ypi = 0) 0 図8.4 q p - bi 最終的な反応確率の分布 先程は 𝜎𝜀 = 1 の場合を見ましたが，もう少し一般化して 𝜎𝜀 が項目ごとに異なる (𝜎𝜀 = 1/𝛼𝑖 ) とします。このとき，誤差の確率分布は 𝜀𝑝𝑖 ∼ 𝑁 (0, 1 𝛼𝑖 ) となるため，両辺を 𝛼𝑖 倍して 𝛼𝑖 𝜀𝑝𝑖 ∼ 𝑁 (0, 1) と表すことができます。𝛼𝑖 𝜀𝑝𝑖 が標準正規分布に従うならば，(8.5)式の正規 累積モデルは 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 [𝜀𝑝𝑖 ≤ (𝜃𝑝 − 𝛽𝑖 )] = 𝑃 [𝛼𝑖 𝜀𝑝𝑖 ≤ 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] 𝛼𝑖 (𝜃𝑝 −𝛽𝑖 ) =∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.6) という形で書き換えることができます。また正規累積モデルが出てきました。(8.5)式との違 いは積分の上限が 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 ) に変わった点です。いま，突拍子もなく 𝛼𝑖 というパラメータ を導入しましたが，これによって項目のパラメータは (𝛽𝑖 , 𝛼𝑖 ) の 2 つになりました。という ことでこの正規累積モデルは，2 パラメータ正規累積モデルと呼ばれます。 8.2.2 項目特性曲線 ここで，2 つの項目パラメータ (𝛽𝑖 , 𝛼𝑖 ) が何を表しているのかを明らかにしていくため， 図8.3を別の形で表してみます。 図8.5の下半分は，横軸に 𝜃𝑝 − 𝛽𝑖 を，縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとったグラフです。図の上半分 にあるように，𝜃𝑝 − 𝛽𝑖 の値が大きくなるほど分布の青い部分が増える，つまり 𝑃 (𝑦𝑝𝑖 = 1) の値が大きくなることを表しています。 続いてパラメータ 𝛽𝑖 の挙動を見るために，横軸に 𝜃𝑝 を，縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとったグラ フを，𝛽𝑖 の値を変えながらいくつか重ねてみます（図8.6）。誤差は平均 0 の正規分布に従う ため，𝑃 (𝑦𝑝𝑖 = 1) = 0.5 になるのは，𝜃𝑝 = 𝛽𝑖 のときです。したがって 𝛽𝑖 の値が大きい項目 ほど，𝜃𝑝 の値が大きくないと「当てはまる」と回答する確率が高くなりません。そのような 意味で，𝛽𝑖 は項目困難度 (item difficulty) と呼ばれます*4 。 次は 𝛼𝑖 の役割を見るために，図8.5よりも誤差分散 (𝛼𝑖 ) が小さい場合や大きい場合を見てみ ましょう。図8.7の左側に示されているように，誤差分散が小さいほど分布の広がりがなくな るため，𝜃𝑝 − 𝛽𝑖 の値の変化に対する 𝑃 (𝑦𝑝𝑖 = 1) の変化が急激になっています。 これを踏まえて，𝛽𝑖 の時と同じように横軸に 𝜃𝑝 をとり，縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとったグラ *4 5 もとが学力テストなので「困難度」という呼び方をしています。正解・不正解がない心理尺度でも特に気にせ ず「困難度」と呼びます。

8.2 基本的なモデル P(ypi = 1) q p - bi 0 図8.5 パラメータと正答率の関係 P(ypi = 1) b i = - 1.5 bi = 0 bi = 2 qp - 1.5 0 2 図8.6 𝛽𝑖 の役割 P(ypi = 1) P(ypi = 1) q p - bi 0 q p - bi 0 図8.7 左：誤差分散が小さい場合，右：誤差分散が大きい場合 6

8.2 基本的なモデル フを，今度は 𝛼𝑖 の値を変えながらいくつかプロットしたものが図8.8です（𝛽𝑖 = 0）。誤差 分散の逆数である 𝛼𝑖 の値が大きい項目ほど，𝜃𝑝 が低い人と高い人の間で 𝑃 (𝑦𝑝𝑖 = 1) が大 きく変化しています。言い換えると 𝛼𝑖 の値が大きい項目ほど，その項目への回答によって 𝜃𝑝 の高低を識別する能力が高いといえます。ということで 𝛼𝑖 のことを項目識別力 (item discrimination) と呼びます。 P(ypi = 1) ai = 0.5 ai = 1 ai = 2 qp 図8.8 𝛼𝑖 の役割 図8.6や図8.8のように，横軸に 𝜃𝑝 をとり，縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとったグラフは，いわば回 答者の特性値に対する項目の特性を表しています。ということで，これらのグラフのことを 項目特性曲線 (item characteristic curve [ICC]) と呼びます。 8.2.3 ロジスティックモデル 2 パラメータ正規累積モデルの式を再掲しておきましょう。 𝛼𝑖 (𝜃𝑝 −𝛽𝑖 ) 𝑃 (𝑦𝑝𝑖 = 1) = ∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.7) 正規累積モデルには，積分計算が含まれています。ある程度予想がつくかもしれませんが， コンピュータは積分計算が比較的得意ではありません。ということで，積分計算がいらない 別のモデルを考えてみましょう。 正規累積モデルがやっていることは，プロビット回帰と同じです。そんな正規累積モデルの 代わりになるものということは，正規分布のように • 累積分布関数が S 字になっていて • 確率密度関数は左右対称の山型 な分布があれば良いわけですが…そのような確率分布として，ロジスティック分布というも のがありました。ロジスティック分布の確率密度および累積分布関数は 𝑓(𝑥) = exp(−𝑥) 2 [1 + exp(−𝑥)] , 𝐹 (𝑥) = 1 1 + exp(−𝑥) (8.8) と表され，確率密度関数を正規分布と重ねて見ると図8.9のようになります。こうしてみると， 裾の重さは違いますが確かにロジスティック分布は正規分布のような左右対称形になってい るようです。 7

8.2 基本的なモデル 0.4 正規分布 ロジスティック分布 f(x) 0.3 0.2 0.1 0.0 -4 -2 図8.9 0 X 2 4 正規分布とロジスティック分布 また，ロジスティック分布の 𝑥 をおよそ 1.7 倍すると，累積分布関数が正規分布とかなり近 くなることが知られています (図8.10)*5 。 1.00 F(x) 0.75 正規分布 1.7倍ロジスティック分布 0.50 0.25 0.00 -4 -2 図8.10 0 X 2 4 1.7 倍ロジスティック分布 ということで，正規分布の代わりにロジスティック分布を使った (2 パラメータ) ロジスティ ックモデルは 𝑃 (𝑦𝑝𝑖 = 1) = 1 1 + exp [−1.7𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.9) と表されます。いま説明したように，この 1.7 はあくまでも「正規累積モデルに近づけるた めの定数」なので，特に正規累積モデルとの比較をする予定が無いならば 1.7 は無くても問 題ありません。ただ，結果を報告する際には「正規累積モデルとロジスティックモデルのど ちらを使用したのか」「ロジスティックモデルならば 1.7 倍した値なのか」がハッキリわか るようにしておきましょう。論文などでは，モデル式を書いておけば一発です。 以後の説明では，特に正規累積モデルと比較することもないので1.7は省略していきます。 また，1 パラメータロジスティックモデルはラッシュモデル (Rasch model) と呼ばれるこ ともあります。これは，IRT とは独立に Rasch という人が理論的に「1 パラメータロジステ ィックモデルと同じ関数」を導出した (Rasch, 1980) ことに由来します。ということで「ラ ッシュモデル」にはやや高尚な理念があったりしますが，私もよく理解できていないのでそ この説明はしません。とりあえず「ラッシュモデル」という名前が出てきたら「困難度だけ *5 8 もっと近づけるためには 1.704 倍のほうがよかったり， ズレを極限まで減らすためには 0.07056𝑥3 +1.5976𝑥 にすると良い (Bowling et al., 2009) ということは分かっていますが，そこまでの精度は必要ない，という ことで，わかりやすい 1.7 が用いられています。

8.2 基本的なモデル のロジスティックモデルなんだな」と理解できていれば良いでしょう。厳密な使い分けとし ては，ラッシュモデルは「正規累積モデルの近似」を目的としていないので，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 パラメータモデルに加えて「当て推量」の パラメータ 𝑐𝑖 が追加されます。 𝑃 (𝑦𝑝𝑖 = 1) = 𝑐𝑖 + *6 1 − 𝑐𝑖 1 − exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.10) 正規累積モデルはプロビットモデル，ロジスティックモデルはロジットモデルに対応します。…と考えると 「本当は理論的にはどちらを使ったほうが良いか」が明確なケースもあるのかもしれません。実用上はどちら でも良いと思います。 9

8.3 因子分析との関係 𝑐𝑖 は項目特性曲線の下限を操作します。𝜃𝑝 がどんなに低い人でも，必ず 𝑃 (𝑦𝑝𝑖 = 1) は 𝑐𝑖 以上の値になります。例えば 4 択問題であれば，どんなに能 力が低くても，勘で答えれば正答率は 1/4 になります。これが 𝑐𝑖 です。 4 パラメータ 4 パラメータモデルでは，3 パラメータモデルに加えて「上限」のパラ メータ 𝑑𝑖 が追加されます。 𝑃 (𝑦𝑝𝑖 = 1) = 𝑐𝑖 + 𝑑𝑖 − 𝑐 𝑖 1 − exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.11) 4 パラメータモデルでは，𝜃𝑝 がどんなに高い人でも，必ず 𝑃 (𝑦𝑝𝑖 = 1) は 𝑑𝑖 以 下の値になります。どんなに能力がある人でも 100% にはならないもの，例え ばバスケのフリースローなどをイメージすると良いかもしれません。 5 パラメータ 5 パラメータモデルでは，4 パラメータモデルに加えて「非対称性」の パラメータ 𝑒𝑖 が追加されます。 𝑃 (𝑦𝑝𝑖 = 1) = 𝑐𝑖 + 𝑑𝑖 − 𝑐 𝑖 (1 − exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )]) 𝑒𝑖 (8.12) 通常の IRT モデルでは，𝑃 (𝑦𝑝𝑖 = 1) の動き方は 0.5 を軸に対称になっていま す。これに対して，例えば「最初は急激に成功率が上がるが，上に行くほどき つくなる」的なことを表現したい場合には非対称性を考慮する必要があるでし ょう。 ただ，これらのモデルは安定した推定に必要なサンプルサイズが莫大であったり，そ もそも推定量に一致性がないなどの問題を抱えていることから，実務場面で使われる ことはほとんどありません（強いて言えば，当て推量パラメータを「選択肢数分の 1」 に固定した 3 パラメータモデルなどはありえる） 。ので， 「そういうのもあるんだなぁ」 くらいに思っておけば良いと思います。 8.3 因子分析との関係 Chapter 7: SEM の多母集団同時分析のところで説明したように，標準化していない（＝切 片が 0 でない）因子分析モデル（1 因子モデル）は 𝑦𝑝𝑖 = 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 − 𝜀𝑝𝑖 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀 ) (8.13) という形になっています*7 。そしてカテゴリカル因子分析では，観測変数 𝑦𝑝𝑖 はその背後に ある連続量によって決まる，という考え方をしました。心理尺度の項目が 2 件法だとしたら， 観測変数 𝑦𝑝𝑖 とその背後の連続量の関係は図8.3と同じようにして図8.11のように表すことが できます。 したがって，標準正規分布に従う誤差 𝜀𝑝𝑖 が 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 より小さいときに 𝑦𝑝𝑖 = 1 となるので， その確率は 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 (𝜀𝑝𝑖 ≤ 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 ) (𝜏𝑖 +𝑓𝑝 𝑏𝑖 ) =∫ −∞ *7 10 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 IRT の説明に合わせて誤差の 𝜀𝑝𝑖 の符号を変えています。 (8.14)

8.4 R で IRT N(t i + fpbi, se ) P(ypi = 0) 0 図8.11 P(ypi = 1) t i + f pb i 因子分析モデルの背後に となります。これを 2 パラメータ正規累積モデル（8.6式）に対応させて見ると， 𝜏𝑖 ) 𝑏𝑖 = 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 ) 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 = 𝑏𝑖 (𝑓𝑝 + (8.15) ということで，(𝑓𝑝 , 𝑏𝑖 , 𝜏𝑖 ) = (𝜃𝑝 , 𝛼𝑖 , 𝛼𝑖 𝛽𝑖 ) としてみればこれら 2 つのモデルは完全に同一だ ということがわかります。 IRT では，(8.15)式に基づく（因子分析的な）パラメータ化を行うことがあります。つまり 2 パラメータロジスティックモデルで言えば 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] 1 = 1 + exp [−(𝛼𝑖 𝜃𝑝 − 𝜏𝑖 )] 𝑃 (𝑦𝑝𝑖 = 1) = (8.16) という形で，困難度パラメータの代わりに切片パラメータを求めるということです。これも どちらの設定を使用しても良いのですが，ソフトウェアによってはこの切片パラメータを使 用する方法がデフォルトになっていることもあるので，出力を見る際は気をつけましょう。 なお IRT の枠組みでは，どちらかというと困難度パラメータ 𝛽𝑖 を用いる解釈のほうがよく 使われます。推定結果で切片パラメータが出てきた場合には，自分で 𝛽𝑖 = 𝜏𝑖 𝛼𝑖 と変換する必 要があるかもしれません。 8.4 R で IRT IRT の基本を学んだところで，いよいよ実際に R で IRT を実行してみます。これまでの説 明に合わせて，まずは二値で試してみるため，それ用のデータを作成します。 二値データの作成 1 # 使用する部分だけ取り出す 2 # 細かいことは考えずに10項目をそのまま使います 3 dat_bin <- dat[, cols] 4 # 1から3は0，4から6は1に変換して二値化 5 dat_bin <- (dat_bin >= 4) 6 # as.numericで数字に変換 7 dat_bin <- apply(dat_bin, 2, as.numeric) 推定には，mirt パッケージの中にある mirt() 関数を使用します。基本的な引数は以下の 3 つです。 11

8.4 R で IRT data 推定に用いるデータです。lavaan の関数と異なり，推定に使う変数のみが入ったデ ータを用意する必要があります。 model lavaan のように各変数と因子の関係を記述したモデル式です。先程説明したよう に，IRT は探索的因子分析と同じものなので，モデル式を記述することで多次元の IRT モデルも実行可能になる，というわけです（詳細は後ほど）。変数名が連番になっ ている場合，後述のようにハイフンでまとめて指定できるのでおすすめです。 itemtype モデルの指定です。たぶん正規累積モデルは選べず，すべてロジスティック関数 ベースのモデルです。これまでに紹介したモデルでは 'Rasch' ラッシュモデル，つまり 1 パラメータモデルです。後述のものに合わせ て'1PL' と指定することはできないので気をつけてください。 '2PL' 2 パラメータモデルです。同じように'3PL', '4PL' も指定可能です。 となります。他にも色々ありますが，後ほどもう少し紹介します。 R で項目反応理論 1 # モデル式の指定 2 # lavaanとは違い，ふつうにイコールで結びます 3 model <- "f_1 = Q1_1-Q1_10" 4 # ハイフンの指定は以下の書き方の簡略化 5 # model <- 'f_1 = Q1_1 + Q1_2 + Q1_3 + (中略) + Q1_10' 1 result_2plm <- mirt(dat_bin, model = model, itemtype = "2PL") 得られた推定値を見る場合には，coef() 関数を使います。 推定値を見る 1 coef(result_2plm) 1 $Q1_1 2 3 a1 d g u par 0.416 1.262 0 1 4 5 $Q1_2 6 7 a1 d g u par 1.023 2.359 0 1 8 9 $Q1_3 10 11 a1 d g u par 1.027 1.911 0 1 12 13 14 12 $Q1_4 a1 d g u

8.4 R で IRT 15 par 1.033 1.737 0 1 16 17 $Q1_5 18 19 a1 d g u par 1.01 1.791 0 1 20 21 $Q1_6 22 23 a1 d g u par 1.203 1.911 0 1 24 25 $Q1_7 26 27 a1 d g u par 1.367 1.667 0 1 28 29 $Q1_8 30 31 a1 d g u par 1.225 1.582 0 1 32 33 $Q1_9 34 35 a1 d g u par 1.339 1.365 0 1 36 37 $Q1_10 38 39 a1 d g u par 1.27 0.031 0 1 40 41 $GroupPars 42 43 MEAN_1 COV_11 par 0 1 デフォルトでは，1 項目ごとの list 型で返すため少し見にくいです。引数 simplify = TRUE とすることで，データフレームの形で表示してくれるようになります。 推定値を見る 1 coef(result_2plm, simplify = TRUE) 1 $items 2 13 a1 d g u 3 Q1_1 0.416 1.262 0 1 4 Q1_2 1.023 2.359 0 1 5 Q1_3 1.027 1.911 0 1 6 Q1_4 1.033 1.737 0 1 7 Q1_5 1.010 1.791 0 1 8 Q1_6 1.203 1.911 0 1 9 Q1_7 1.367 1.667 0 1

8.4 R で IRT 10 Q1_8 1.225 1.582 0 1 11 Q1_9 1.339 1.365 0 1 12 Q1_10 1.270 0.031 0 1 13 14 $means 15 f_1 16 0 17 18 $cov 19 20 f_1 f_1 1 最初の $item 部分に表示されているものが項目パラメータです。左から a1 識別力パラメータの推定値 d 切片パラメータの推定値 g 当て推量パラメータ (guessing) の推定値（2 パラメータモデルなら全て 0） u 上限パラメータ (upper bound) の推定値（2 パラメータモデルなら全て 1） となっています。このように，mirt はデフォルトでは切片パラメータを出力します。困難度 パラメータを見たい場合には d 列を a1 列で割るか，引数 IRTpars = TRUE をつけて出力し ましょう。 困難度パラメータを見る 1 coef(result_2plm, simplify = TRUE, IRTpars = TRUE) 1 $items 2 a Q1_1 0.416 -3.031 0 1 4 Q1_2 1.023 -2.307 0 1 5 Q1_3 1.027 -1.861 0 1 6 Q1_4 1.033 -1.681 0 1 7 Q1_5 1.010 -1.773 0 1 8 Q1_6 1.203 -1.588 0 1 9 Q1_7 1.367 -1.220 0 1 10 Q1_8 1.225 -1.292 0 1 11 Q1_9 1.339 -1.019 0 1 12 Q1_10 1.270 -0.025 0 1 13 14 $means 15 f_1 16 0 17 14 b g u 3

8.4 R で IRT 18 $cov 19 20 f_1 f_1 1 これで d 列の代わりに b（困難度パラメータ）列が表示されるようになりました。出力を見 ると，例えば Q1_1 では，1.262/0.416 ≃ 3.031 ということで，確かに d 列を a1 列で割った 値（の符号を入れ替えたもの）が表示されています。 $item の下に表示されている $means と $cov はそれぞれ特性値（因子得点）𝜃 の平均値と 分散です。この時点では一般的な因子分析の制約（𝜃 ∼ 𝑁 (0, 1)）に沿った値のみが表示され ているため特に注目する必要はありません。後ほど紹介する多次元モデルや多母集団モデル になると，ここに表示される値も重要な意味を持つことになります。 plot() 関数を使うと，引数 type に応じて，推定結果に基づく様々なプロットを全項目ま とめて出すことができます。ICC を出したい場合には type = "trace" としてあげてくだ さい。 項目特性曲線 1 plot(result_2plm, type = "trace") Item Probability Functions -6 -4 -2 0 2 4 6 Q1_9 Q1_10 Q1_5 Q1_6 1.0 0.8 0.6 0.4 0.2 0.0 Q1_7 Q1_8 P (q ) 1.0 0.8 0.6 0.4 0.2 0.0 Q1_1 Q1_2 Q1_3 Q1_4 1.0 0.8 0.6 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 q 図8.12 項目特性曲線 デフォルトではすべての項目の ICC を縦横に並べてくれるのですが，項目数が多くなると見 るのが大変そうです。この場合，引数 facet_items = FALSE とすると，すべての項目のプ ロットを重ねて表示してくれます。また引数 which.items を指定すると，特定の項目だけ の ICC を出すことも出来ます。 15

8.5 IRT の前提条件 1 # Q1_1のICCを出すならitem = 1 2 plot(result_2plm, type = "trace", which.items = 1:5, facet_items = FALSE) Item Probability Functions P1 1.0 P (q ) 0.8 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 0.6 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 q 図8.13 最初の 5 項目の ICC ということで図8.12や8.13を見ると，識別力がダントツで低い項目 1 では，ICC の傾きが緩 やかになっていることがわかります。 項目パラメータを固定した上での特性値 𝜃𝑝 の推定値を出す場合は，fscores() 関数を使い ます。このあたりは因子分析や CFA での「項目パラメータを推定」→「項目パラメータを 固定して個人特性値を推定」と全く同じです。 theta を出す 1 head(fscores(result_2plm)) 1 f_1 2 [1,] -1.537468218 3 [2,] -0.007789283 4 [3,] 5 [4,] -0.643698451 6 [5,] 0.085994618 7 [6,] 0.788036256 0.226340769 8.5 IRT の前提条件 IRT，特に 2 パラメータモデルは基本的には因子分析と同じようなものなのですが，理論的 背景の違いなどから留意点も少し異なってきます。ここでは，IRT を適用するにあたって重 要となる 2 つの前提条件について見ていきます。 16

8.5 IRT の前提条件 8.5.1 局所独立性 局所独立性 (local independence) の仮定は，因子分析のときと概ね同じです。探索的因子 分析では，独自因子の間が完全に無相関であるという仮定をおき，その中でデータの相関行 列を最もよく復元できる因子負荷行列と独自因子の分散 (𝐁⊤ 𝐁 + 𝚿) を求めていました。独 自因子の間に相関があるということは，その当該項目の間に何か別の共通因子があるという ことです。 局所独立性をもう少し固い言い方で説明すると， 「𝜃𝑝 で条件づけたとき，項目𝑖と項目𝑗 (𝑖 ≠ 𝑗) への回答は完全に独立」となります。項目への回答は二値なので，2 × 2 クロス表をイメージ してみましょう。もし局所独立性が満たされている場合，𝜃𝑝 が特定の値の人を 1000 人ほど 集めて（つまり 𝜃𝑝 で条件づけて）表8.1のように 2 つの項目 𝑖, 𝑗 の回答のクロス表を作ると， 連関がゼロになっています。項目 𝑖 に「当てはまる」と回答した人を見ても， 「当てはまらな い」と回答した人を見ても，項目 𝑗 に「当てはまる」と回答した人の割合は全体の割合と同 じ 6:4 になっています。 表8.1 局所独立な 2 項目 項目 𝑗 項目 𝑖 当てはまらない 当てはまる 計 当てはまらない 180 420 600 当てはまる 120 280 400 計 300 700 1000 一方で，𝜃𝑝 以外の別の共通因子がある場合，クロス表に連関が現れます。例えば「他者への 助けを求める傾向」の尺度において，宗教観の影響を受けると思われる項目（e.g., 困ったと きは神に祈る）どうしのクロス表を見てみると，表8.2のように • 項目 𝑖 に「当てはまる」と答えた人では 340/400 = 85% の人が項目 𝑗 でも「当ては まる」と回答 • 項目 𝑖 に「当てはまらない」と答えた人では 360/600 = 60% の人が項目 𝑗 でも「当 てはまる」と回答 しているため，項目 𝑖 の回答と項目 𝑗 の回答が関連を持っていることになります。項目 𝑖 に 「当てはまる」と答えた人のほうが項目 𝑗 でも「当てはまる」と回答する確率が高い，という ことです。 表8.2 局所独立ではない 2 項目 項目 𝑗 当てはまらない 項目 𝑖 当てはまる 計 17 当てはまらない 当てはまる 計 240 360 600 60 340 400 300 700 1000

8.5 IRT の前提条件 局所独立の仮定とは，このようにクロス表を作成したときに，どの項目ペア間でも，どの特 性値の値でも連関がない，ということを指しています*8 。 局所独立性がなぜ重要なのでしょうか。その理由の一つは計算上の問題です。IRT では基本 的に最尤法によってパラメータを推定します。もしも局所独立性が満たされていれば，𝜃𝑝 で 条件づけたときの各項目に対する反応確率は独立なので，回答者 𝑝 に対する尤度関数は単純 にすべての項目をかけ合わせた 𝐼 𝑃 (𝐲𝑝 |𝜃𝑝 ) = ∏ 𝑃 (𝑦𝑝𝑖 = 1|𝜃𝑝 )𝑦𝑝𝑖 𝑃 (𝑦𝑝𝑖 = 0|𝜃𝑝 )1−𝑦𝑝𝑖 (8.17) 𝑖=1 と表せます*9 。これに対してもしも局所独立が満たされていない場合（表8.2） ，𝑃 (𝑦𝑝𝑗 = 1|𝜃𝑝 ) の値が 𝑦𝑝𝑖 によって変わることになります。そのため単純に掛け算ができなくなってしまい， もっと複雑な計算が必要になってしまうのです。 ということで，IRT モデルは「局所独立性が満たされている」という仮定にもとづいて定式 化されています。したがって，データを収集する時点で局所独立性が満たされていることを しっかりと確認しておく必要があります。なお後半では，局所独立性が満たされているかを 事後的に評価する指標を紹介します。 8.5.2 一次元性 これまで紹介してきた IRT モデルはいわば一因子の因子分析モデルです。因子分析では「項 目をグルーピングする」ことに重きを置きがちだったので，因子数は「いい感じにグループ ができる数」にすることが一般的でした。これに対して IRT はもともと単一のテストに対 する単一の能力を測定することが主たる目的と言えるため，基本的には因子数は 1 となりま す*10 。したがって，IRT の重要な前提条件としてすべての項目が同じ能力・特性のみを反映 していること（一次元性）が必要となります。 一次元性を検証するためのフレームワークは，多くの場合因子分析と共通です。つまり， Chapter 6: 因子分析のところで紹介したスクリープロットや平行分析，MAP などの基準に よって「因子数 1」が提案されたら良い，ということです。ですが，因子分析的な次元性の確 認方法では，項目数が増えるとなかなか「因子数 1 が最適」とはならないケースが出てきま す。したがって IRT では，「一次元とみなしても問題ないか」といった別の視点から一次元 性の評価を行うことがあります。これは，例えば一因子の CFA を適用した際の適合度指標 （RMSEA や CFA, AGFI など）が許容範囲にあるかによって確認することが出来ます。 また IRT では，局所独立性の観点から一次元性を確認する指標もいくつか提案されていま す。一次元性が満たされているということは，回答行動に影響する共通因子が一つ（𝜃𝑝 ）だ けということです。したがって，Chapter 6（図 6.6）でも少し説明したように，一次元性が *8 厳密には，「独立」と「無相関」は別の概念です。理想は独立なのですが，実用上は無相関であれば良いだろ う，という感じになっています。厳密な独立とは異なるという意味で，弱局所独立性と呼ぶこともあります。 *9 「𝜃 で条件づけた」ということが明確になるようにここだけ 𝑃 (.|𝜃 ) としていますが，これまでの 𝑃 (.) と 𝑝 𝑝 意味は同じです。 *10 次回には少しだけ多次元モデルの紹介をしますが，これは一次元モデルの拡張として展開されるものであり， IRT の基本は一次元だといえます。 18

8.6 多値型モデル 満たされていない場合には，共通因子で説明出来ない部分（独自因子）の間に共分散が発生 することになります。そのように考えると，一次元性は局所独立と同じようなものだとみな すことができるわけです。厳密には局所独立性が満たされている場合，一次元性も満たされ ていると言って良いだろう，ということになります*11 。ということで，より洗練された一 次元性確認の方法として，例えば DIMTEST(Stout et al., 2001) や DETECT(Zhang and Stout, 1999) という名前の方法が提案されていたりします。 局所独立と一次元性 一次元性が満たされているだけでは局所独立性が満たされるとは限りません。例えば • 変数 𝑥 の平均値を求めよ • 変数 𝑥 の分散を求めよ という 2 つの問題があった場合，この 2 つはいずれも「統計の知識」を問うている という意味では同じ潜在特性に関する項目であり，一次元性を満たしているといえま す。ですが，分散を計算する過程では必ず平均値を計算する必要があるため，項目1 に間違えた人はほぼ確実に項目2に間違えます。したがって，この 2 問は一次元性は 満たされているが局所独立性は満たされていない状態だといえます。このように，あ る項目への回答が他の項目への回答に直接的に影響を与える場合を実験的独立性の欠 如と呼びます。心理尺度でいえば， 前の項目で『どちらかといえば当てはまる』以上を選んだ方にお聞きします。前 の項目で『どちらかといえば当てはまらない』以下を選んだ方は，この項目では 『全く当てはまらない』を選んでください。 といった指示を出した場合には，実験的独立性が欠如していると言えそうです。他に は，回答時間が足りない場合の後ろの方の項目も，「前の項目に回答出来ていない人 は次の項目にも回答できていない」ため，実験的独立性が欠如しているとみなすこと が出来ます。 ということで局所独立性を正確に表すと，一次元性に加えてこの実験的独立性が満た されていることが必要といえます (南風原, 2000; Lord and Novick, 1968)。 8.6 多値型モデル ここまで，二値の項目反応データに対する基本的なモデルを紹介し，R での実行方法を示し ました。しかし実際の場面では，リッカート式尺度ならばたいてい 5 件法から 7 件法の多値 項目です。そしてテストの場面でも，部分点や組問*12 を考えると多値型のモデルが必要とな るでしょう。ということで，代表的な 2 つの多値型モデルを紹介したいと思います。どちら を使っても本質的には違いは無いので，説明を読んでしっくり来たほうを使えば良いと思い ます。 *11 もちろん，多因子モデルの場合には，一次元とは限りません。ただ同様に，局所独立性が満たされているなら ば 𝑛 次元性が満たされている，ということはできると思います。 *12 問 1 の答えを使わないと問 2 が解けないような場合，局所独立性が満たされなくなってしまいます。そこで 「問 1 が解けたら 1 点，問 2 も解けたら 2 点」というように得点化してこれを一つの大きな問題として考える ことで局所独立性を確保しようという考えに至るわけです。 19

8.6 多値型モデル 8.6.1 段階反応モデル 段階反応モデル (Graded response model [GRM]: Samejima, 1969) では，複数の二値 IRT モデルを組み合わせて多値反応を表現します。項目 𝑖 に対する回答者 𝑝 の回答 𝑦𝑝𝑖 = 𝑘 (𝑘 = 1, 2, ⋯ , 𝐾) について「𝑘 以上のカテゴリを選ぶ確率」を考えると，これはまだ二値（𝑘 以上 or 𝑘 より小さい）なので，2 パラメータロジスティックモデルならば 𝑃 (𝑦𝑝𝑖 ≥ 𝑘) = 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] (8.18) と表せます。なお，困難度パラメータは「項目 𝑖 のカテゴリ 𝑘」ごとに用意されるため 𝛽𝑖𝑘 と しています。 異なるカテゴリを選ぶ確率は当然排反な事象なので，二値モデルを組み合わせると「ちょう ど 𝑘 番目のカテゴリを選ぶ確率」は 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 ≥ 𝑘) − 𝑃 (𝑦𝑝𝑖 ≥ 𝑘 + 1) (8.19) と表すことができます。ただし，端のカテゴリに関する確率はそれぞれ 𝑃 (𝑦𝑝𝑖 ≥ 1) = 1 𝑃 (𝑦𝑝𝑖 ≥ 𝐾 + 1) = 0 (8.20) とします。 「1 以上のカテゴリを選ぶ確率」は当然 1 になるため，これに対する困難度パラメ ータ 𝛽𝑖1 = −∞ となります。同様に， 「𝐾 + 1 以上のカテゴリを選ぶ確率」は 0 になります。 𝑃 (𝑦𝑝𝑖 = 𝐾|𝜃𝑝 ) を考えるときだけ仮想的な「さらに上のカテゴリ」を想定する，ということ です。このように GRM は二値型モデルの素直な拡張のため，数学的にはカテゴリカル因子 分析と同じものだとみなすことができます。 図8.14に，4 件法の項目に対する段階反応モデル（(𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4 ]) = (1, [−1.5, 0, 2])）の イメージを示しました。左右で色の塗り方や分布の位置は変えていません。点線の位置だけ です。 左は，𝜃𝑝 の位置に応じて各カテゴリの反応確率が変化する状況を表しています。反応確率は 全カテゴリ合計で 1 となっているため，点線の位置での各カテゴリの長さの比率がそのまま カテゴリ選択確率となっています。 右は，𝛽𝑖𝑘 の位置に点線を引いています。カテゴリの困難度 𝛽𝑖𝑘 は二値モデルの時と同じよ うに 𝑃 (𝑦𝑝𝑖 ≥ 𝑘) = 0.5 になる 𝜃𝑝 の値を表しています。そのため，上の確率密度分布を見る と，色が変わるポイントはすべての分布で同じになっていることがわかります。したがって， 𝜃𝑝 の値が大きい（正規分布が右に動く）ほど，また 𝛽𝑖𝑘 の値が小さい（閾値の点線が左に動 く）ほど，上位カテゴリを選択する確率が高くなっていくわけです。 続いて，図8.15には，いくつかの項目パラメータのセットにおける項目特性曲線 (ICC) を示 しました。図8.14の下半分では 𝑃 (𝑦𝑝𝑖 ≥ 𝑘) のプロットを示しましたが，通常「項目特性曲線」 というと 𝑃 (𝑦𝑝𝑖 = 𝑘) のプロットです。二値項目の場合はどちらでも同じですが，多値項目 の場合は異なるため改めて ICC を示しています。なお，多値型モデルの場合は図8.15の一本 一本の線のことを項目反応カテゴリ特性曲線 (item response category characteristic curve [IRCCC]) と呼ぶこともあります。が，そこまで細かい名前を覚える必要はないでしょう。 20

8.6 多値型モデル P(ypi ³ k) P(ypi ³ k) 0.5 qp qp - 1.5 0 図8.14 P(ypi = k) 3 2 段階反応モデル P(ypi = k) (ai, b i2, b i3, b i4) (1, - 1.5, 0, 2) 1 0 (ai, b i2, b i3, b i4) (1, - 0.5, 1, 3) 4 1 4 3 2 2 qp - 0.75 P(ypi = k) 1 qp 0.25 P(ypi = k) (ai, b i2, b i3, b i4) (0.5, - 1.5, 0, 2) 4 2 (ai, b i2, b i3, b i4) (2, - 1.5, 0, 2) 1 3 4 1 2 3 2 qp - 0.75 qp - 0.75 1 図8.15 1 GRM の項目特性曲線 左上のプロットは，図8.14と同じ項目パラメータ（(𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4 ]) = (1, [−1.5, 0, 2])）で の IRCCC です。基本的には，このように 𝜃𝑝 の値に応じて左から順に各カテゴリの選択確率 のピークが訪れます。点線は中間カテゴリの選択確率が最も高くなる箇所ですが，これは具 体的にカテゴリ困難度 𝛽𝑖𝑘 と 𝛽𝑖,𝑘+1 のちょうど中間になります。 右 上 の プ ロ ッ ト は， 左 上 か ら す べ て の カ テ ゴ リ 困 難 度 の 値 を 1 ず つ 大 き く し ま し た （(𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4 ]) = (1, [−0.5, 1, 3])）。二値モデルの時と同じように，すべてのカテゴ リ困難度が同じ値だけ変化するときには IRCCC は平行移動します。 21

8.6 多値型モデル 左 下 の プ ロ ッ ト は， 左 上 か ら 識 別 力 を 半 分 に し ま し た （(𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4 ]) = (0.5, [−0.5, 1, 3])） 。 ま た 右 下 の プ ロ ッ ト は， 左 上 か ら 識 別 力 を 2 倍 に し て い ま す （(𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4] ) = (2, [−0.5, 1, 3])）。識別力は各 IRCCC の山のきつさを調整していま す。識別力が低い場合，𝜃𝑝 の大小に対して各カテゴリの選択確率があまり変化しなくなるた め，反対に「あるカテゴリを選択した」という情報が 𝜃𝑝 の推定に及ぼす影響も小さくなって しまいます。 また，左下のプロットのように，どの 𝜃𝑝 の値においてもあるカテゴリの選択確率が最大とな らないこともありえます。識別力が低い場合や，カテゴリ困難度が近すぎる場合，あるいは そもそもカテゴリ数が多い場合にはそのような現象が起こりえます。ただ何も問題ではない ので，そのようなカテゴリが出現した場合には「選ばれにくいんだなぁ」くらいに思ってお いてください。 8.6.2 一般化部分採点モデル 有名な多値型モデルのもう一つが，一般化部分採点モデル (generalized partial credit model [GPCM]: Muraki, 1992) です。GPCM では，隣のカテゴリとの関係を別の視点から考えま す。ということで，まずは「隣のカテゴリとの関係」を二値（2 パラメータロジスティック） モデルで考えてみましょう。二値モデルでは， 「カテゴリ 0 を選ぶ確率」と「カテゴリ 1 を選 ぶ確率」の和は当然 1 になります。そして，「カテゴリ 1 を選ぶ（e.g., 『当てはまる』を選 ぶ，正解する）」確率は 𝑃 (𝑦𝑝𝑖 = 1) 𝑃 (𝑦𝑝𝑖 = 1) + 𝑃 (𝑦𝑝𝑖 = 0) 1 = 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] 𝑃 (𝑦𝑝𝑖 = 1) = = (8.21) exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] 1 + exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] となります。つまりこの式は，隣のカテゴリとの二択において，当該カテゴリを選ぶ確率と いう見方ができるわけです。これを多値 (𝑘 = 1, 2, ⋯ , 𝐾) に置き換えると，「カテゴリ 𝑘 − 1 と 𝑘 の二択だったら 𝑘 の方を選ぶ確率」に拡張することができそうです。ということで， ∗ GPCM ではこの確率 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 を 𝑃 (𝑦𝑝𝑖 = 𝑘) 𝑃 (𝑦𝑝𝑖 = 𝑘) + 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) 1 = 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] ∗ 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 = = = exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] (8.22) 1 + exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] exp [𝜋𝑝𝑖𝑘 ] 1 + exp [𝜋𝑝𝑖𝑘 ] と表します。なお，この後の式展開を考えて 𝜋𝑝𝑖𝑘 = 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 ) と置き換えています。 ここから，全カテゴリの中でカテゴリ 𝑘 を選ぶ確率を導出しましょう。(8.22) 式を 𝑃 (𝑦𝑝𝑖 = 22

8.6 多値型モデル 𝑘) = の形に変形すると， 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) となります。最後の右辺の 𝑃 1−𝑃 ∗ 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 ∗ 1 − 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 (8.23) の部分は 2 つのカテゴリ間のオッズを表しているため， ∗ 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 = exp(𝜋𝑝𝑖𝑘 ) ∗ 1 − 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 (8.24) 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) exp(𝜋𝑝𝑖𝑘 ) (8.25) を用いて という漸化式が得られ，特定のカテゴリ 𝑘 を選ぶ確率が表現できるようになります。例えば 4 カテゴリの場合には， 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 (𝑦𝑝𝑖 = 1) 𝑃 (𝑦𝑝𝑖 = 2) = 𝑃 (𝑦𝑝𝑖 = 1) exp(𝜋𝑝𝑖2 ) 𝑃 (𝑦𝑝𝑖 = 3) = 𝑃 (𝑦𝑝𝑖 = 1) exp(𝜋𝑝𝑖2 ) exp(𝜋𝑝𝑖3 ) (8.26) 𝑃 (𝑦𝑝𝑖 = 4) = 𝑃 (𝑦𝑝𝑖 = 1) exp(𝜋𝑝𝑖2 ) exp(𝜋𝑝𝑖3 ) exp(𝜋𝑝𝑖4 ) という形で 𝑃 (𝑦𝑝𝑖 = 1) を起点に「一個上のカテゴリに上がるときのオッズ」の積によって すべてのカテゴリの反応確率が表されるわけです。ということで，この式を一般化して 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) exp(𝜋𝑝𝑖𝑘 ) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 2) exp(𝜋𝑝𝑖,𝑘−1 ) exp(𝜋𝑝𝑖𝑘 ) ⋮ 𝑘 = 𝑃 (𝑦𝑝𝑖 = 1) ∏ exp(𝜋𝑝𝑖𝑐 ) (8.27) 𝑐=2 𝑘 = 𝑃 (𝑦𝑝𝑖 = 1) exp (∑ 𝜋𝑝𝑖𝑐 ) 𝑐=2 と書き下します。ここで 𝑃 (𝑦𝑝𝑖 = 1) はこのまま残ってしまいますが，GPCM の考え方 ではこれ以上計算できません*13 。ですが先程の式を見ると，すべてのカテゴリの反応確 率は「𝑃 (𝑦𝑝𝑖 = 1) の何倍」という形になっています。ので一旦 𝜋𝑝𝑖1 = 0 とすることで 𝑃 (𝑦𝑝𝑖 = 1) = exp(𝜋𝑝𝑖1 ) = 1（基準）と仮置きして， 𝑘 𝑃 (𝑦𝑝𝑖 = 𝑘) = exp (∑ 𝜋𝑝𝑖𝑐 ) (8.28) 𝑐=1 とまとめてしまいます。つまり，各カテゴリの反応確率を「一番下のカテゴリの何倍か」という 形で表してみます。こによって，ひとまず 𝑃 (𝑦𝑝𝑖 = 𝑘) が定式化されましたが，𝑃 (𝑦𝑝𝑖 = 1) = 1 と設定している以上 𝑃 (𝑦𝑝𝑖 = 𝑘) の全カテゴリの総和が 1 になるかわかりません。 *13 𝑃 (𝑦𝑝𝑖 = 1) = exp(𝜋𝑝𝑖1 )𝑃 (𝑦𝑝𝑖 = 0) と表せますが，カテゴリ 𝑘 = 0 は存在しない，つまり 𝑃 (𝑦𝑝𝑖 = 0) = 0 ∗ = 1) です。ということは 𝑃 (𝑦𝑝𝑖 0,1 = ∞ となるため，exp(𝜋𝑝𝑖1 ) = ∞ となってしまい計算ができないので す。 23

8.6 多値型モデル ということで最後に，𝑃 (𝑦𝑝𝑖 = 1) の値が何であっても全カテゴリでの総和が 1 になるように するため，総和で割ることを考えます。 𝑃 (𝑦𝑝𝑖 = 𝑘) 𝐾 ∑𝑙=1 𝑃 (𝑦𝑝𝑖 = 𝑙) = 𝑃 (𝑦𝑝𝑖 𝑃 (𝑦𝑝𝑖 = 𝑘) = 1) + 𝑃 (𝑦𝑝𝑖 = 2) + ⋯ + 𝑃 (𝑦𝑝𝑖 = 𝐾) 𝑘 exp (∑𝑐=1 𝜋𝑝𝑖𝑐 ) = 1 2 𝐾 exp (∑𝑙=1 𝜋𝑝𝑖𝑙 ) + exp (∑𝑙=1 𝜋𝑝𝑖𝑙 ) + ⋯ + exp (∑𝑙=1 𝜋𝑝𝑖𝑙 ) (8.29) 𝑘 = exp (∑𝑐=1 𝜋𝑝𝑖𝑐 ) 𝐾 𝑙 ∑𝑙=1 exp (∑𝑐=1 𝜋𝑝𝑖𝑙 ) となります。すると分子も分母もすべての項に共通の exp(𝜋𝑝𝑖1 ) がかかっているため，この 項の値が何であっても全体の割合に占める各カテゴリの反応確率の割合は変わらずに済みま す。ということで，ようやく GPCM の式が導出できました。図8.16に，ここまでの考え方を まとめてみました。特定のカテゴリの反応確率 𝑃 (𝑦𝑝𝑖 = 𝑘) が直接計算できないため，いった ん一番下のカテゴリの反応確率 𝑃 (𝑦𝑝𝑖 = 1) = 1 とおいて，他のカテゴリは「前のカテゴリの 何倍」という形で表します。最後に全体の総和が 1 になるように調整をしたものが GPCM における反応確率というわけです。 図8.16 GPCM のイメージ 図8.17には，GPCM での項目特性曲線をプロットしました。GPCM では，カテゴリ困難度 𝛽𝑖𝑘 は「カテゴリ 𝑘 − 1 と 𝑘 の二択」における選択確率がちょうど 50% になるポイントを表 しています。したがって，隣接するカテゴリの反応確率はちょうど 𝛽𝑖𝑘 のところで大小関係 が入れ替わります。𝛽𝑖𝑘 はそのように解釈できるわけです。 24

8.6 多値型モデル P(ypi = k) (ai, b i2, b i3, b i4) (1, - 1.5, 0, 2) 1 4 3 2 qp - 1.5 0 2 図8.17 GPCM の項目特性曲線 GRM では，カテゴリ困難度は必ず単調増加になっている必要がありますが，GPCM では単 調増加で無くても問題ありません。図8.18にカテゴリ困難度が単調増加にならないケースを プロットしてみました。この場合， 「カテゴリ 2 がカテゴリ 1 を上回る」よりも先に「カテゴ リ 3 がカテゴリ 2 を上回る」ため，結果的にカテゴリ 2 が最大になることがありません。 P(ypi = k) 1 (ai, b i2, b i3, b i4) (1, 0, - 1.5, 2) 3 4 2 qp - 1.5 図8.18 0 2 カテゴリ困難度が単調増加にならないケース 8.6.3 R で多値型モデル それでは，mirt で 2 つの多値型モデルを実行してみましょう。…といっても引数 itemtype に graded(GRM) か gpcm を指定したら良いだけなので難しいことはありません。 25 1 # 多値型なのでdatをそのまま突っ込む 2 result_grm <- mirt(dat[, cols], model = "f_1 = Q1_1-Q1_10", itemtype = "g raded")

8.6 多値型モデル 1 coef(result_grm, IRTpars = TRUE, simplify = TRUE) 1 $items 2 a b1 b2 b3 b4 b5 3 Q1_1 0.572 -6.382 -3.936 -2.323 -1.040 1.269 4 Q1_2 1.181 -4.051 -2.765 -2.087 -0.840 0.821 5 Q1_3 1.238 -3.245 -2.243 -1.643 -0.564 1.019 6 Q1_4 1.074 -3.293 -2.174 -1.648 -0.706 0.415 7 Q1_5 1.116 -3.951 -2.472 -1.664 -0.460 1.216 8 Q1_6 1.062 -3.988 -2.698 -1.726 -0.440 1.429 9 Q1_7 1.177 -3.454 -2.084 -1.314 -0.178 1.485 10 Q1_8 1.097 -3.671 -2.186 -1.386 -0.023 1.744 11 Q1_9 1.365 -3.306 -2.016 -1.039 -0.312 0.921 12 Q1_10 1.200 -2.226 -1.048 -0.056 0.461 1.586 13 14 $means 15 f_1 16 0 17 18 $cov 19 20 f_1 f_1 1 dat の各項目は 6 件法なので，カテゴリ困難度 b が 5 つ出力されました。二値型モデルの時 と同じように，項目特性曲線を出力することももちろん可能です。 項目特性曲線 1 # カテゴリ数のぶんだけ線を引いてくれる 2 plot(result_grm, type = "trace") また，itemplot(., type = "threshold") 関数を使うと，図8.14の下半分のような「カテ ゴリ 𝑘 以上の反応をする確率」のプロットを出すこともできます。 項目特性曲線 1 # itemplot()は1項目しかプロットできない 2 # 何番目の項目をプロットするかを引数itemで指定 3 itemplot(result_grm, item = 6, type = "threshold") 同じように GPCM でもやってみましょう。やり方も結果の見方も同じなので詳しい説明は 省略します。 26 1 # 多値型なのでdatをそのまま突っ込む 2 result_gpcm <- mirt(dat[, cols], model = "f_1 = Q1_1-Q1_10", itemtype = " gpcm")

8.6 多値型モデル Item Probability Functions -6 -4 -2 0 2 4 6 Q1_9 Q1_10 Q1_5 Q1_6 1.0 0.8 0.6 0.4 0.2 0.0 Q1_7 Q1_8 P (q ) 1.0 0.8 0.6 0.4 0.2 0.0 Q1_1 Q1_2 Q1_3 Q1_4 1.0 0.8 0.6 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 q 図8.19 GRM で項目特性曲線 Probability Thresholds for Item 6 1.0 0.8 P(x > 1) P(x > 2) P(x > 3) P(x > 4) P(x > 5) P (q ) 0.6 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 q 図8.20 1 coef(result_grm, IRTpars = TRUE, simplify = TRUE) 1 $items 2 27 GRM で項目特性曲線 2 a b1 b2 b3 b4 b5 3 Q1_1 0.572 -6.382 -3.936 -2.323 -1.040 1.269 4 Q1_2 1.181 -4.051 -2.765 -2.087 -0.840 0.821 5 Q1_3 1.238 -3.245 -2.243 -1.643 -0.564 1.019 6 Q1_4 1.074 -3.293 -2.174 -1.648 -0.706 0.415 7 Q1_5 1.116 -3.951 -2.472 -1.664 -0.460 1.216 P1 P2 P3 P4 P5 P6

8.7 多次元モデル 8 Q1_6 1.062 -3.988 -2.698 -1.726 -0.440 1.429 9 Q1_7 1.177 -3.454 -2.084 -1.314 -0.178 1.485 10 Q1_8 1.097 -3.671 -2.186 -1.386 -0.023 1.744 11 Q1_9 1.365 -3.306 -2.016 -1.039 -0.312 0.921 12 Q1_10 1.200 -2.226 -1.048 -0.056 0.461 1.586 13 14 $means 15 f_1 16 0 17 18 $cov 19 20 f_1 f_1 1 項目特性曲線 1 # カテゴリ数のぶんだけ線を引いてくれる 2 plot(result_gpcm, type = "trace") Item Probability Functions -6 -4 -2 0 2 4 6 Q1_9 Q1_10 Q1_5 Q1_6 1.0 0.8 0.6 0.4 0.2 0.0 Q1_7 Q1_8 P (q ) 1.0 0.8 0.6 0.4 0.2 0.0 Q1_1 Q1_2 Q1_3 P1 P2 P3 P4 P5 P6 Q1_4 1.0 0.8 0.6 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 q 図8.21 GPCM でも項目特性曲線 8.7 多次元モデル これまでの IRT モデルでは一次元性の仮定に基づく一次元モデルのみを扱ってきました。で すが多くの心理尺度やテストでは，複数の次元の特性を同時に測定したいというニーズがあ ります。先程紹介したようにIRTはカテゴリカル因子分析と数学的には同値なので，IRT で も多次元モデルを考えていきましょう。 28

8.7 多次元モデル 多次元因子分析モデルでは，項目反応を 𝑇 𝑦𝑝𝑖 = 𝜏𝑖 + ∑ 𝑓𝑝𝑡 𝑏𝑡𝑖 − 𝜀𝑝𝑖 (8.30) 𝑡=1 と表しました。これと同じようにして，多次元 IRT モデル（2 パラメータロジスティックモ デル）での項目反応関数は 𝑃 (𝑦𝑝𝑖 = 1) = 1 1+ 𝑇 exp [− ∑𝑡=1 𝛼𝑡𝑖 𝜃𝑝𝑡 − 𝜏𝑖 ] (8.31) となります。項目パラメータは，識別力が次元の数（𝑇 個）と，切片が一つです。また，項 目の全体的な識別力を評価する際には，次元ごとの識別力を一つにまとめた多次元識別力 (multidimensional discrimination: 豊田, 2013) √ √ 𝑇 MDISC𝑖 = √∑ 𝛼2𝑡𝑖 ⎷ 𝑡=1 (8.32) を用いることもあります。これを用いると，切片パラメータ 𝜏𝑖 を困難度に変換することもで きるようになります： 𝛽𝑖 = −𝜏𝑖 . MDISC𝑖 (8.33) 二次元までならば，ICC と同じように 𝛉𝑝 と 𝑃 (𝑦𝑝𝑖 = 1) の関係をプロットすることも可能で す。図8.22の左は，とある二次元項目における項目反応曲面 (item response surface [IRS]) です。X 軸，Y 軸がそれぞれの 𝜃𝑡 の値を表しており，Z 軸の値が反応確率 𝑃 (𝑦𝑝𝑖 = 1) を表 しています。この項目では，([𝛼𝑖1 , 𝛼𝑖2 ], 𝜏𝑖 ) = ([1.5, 0.5], 0.7) となっていることから，𝜃𝑝2 の 値はさほど反応確率に影響を及ぼしていないことが図からもわかります。なお図8.22の右は， 2 つの 𝜃𝑝𝑡 の値に対して 𝑃 (𝑦𝑝𝑖 = 1) が等しくなるラインを表しています。 図8.22 29 項目反応曲面 (Reckase, 2019, Fig. 12.1)

8.7 多次元モデル 8.7.1 補償モデルと非補償モデル ここまで紹介した多次元モデルでは，反応確率に対して 𝜃𝑝𝑡 が及ぼす影響は加法的といえま す。例えば ([𝛼𝑖1 , 𝛼𝑖2 ], 𝜏𝑖 ) = ([1.5, 1], 0) の項目の場合， 𝛉𝑝 = [−1 1] の人と 𝛉𝑝 = [0 −0.5] の人では反応確率は同じになります。また，ある次元の特性値がとても低い場合でも，別の次 元の特性値がとても高い場合には反応確率が高くなる可能性があります。このような次元間 の関係を反映した先程の多次元 IRT モデルは補償型 (compensatory) モデルと呼ばれます。 一方で，複数の次元の要素がすべて揃って初めて反応確率が高くなるようなケースも考えら れます。例えば数学のテストが英語で実施される場合， 「英語」と「数学」の二次元構造ではあ りますが，補償型モデルのように「英語が苦手でも数学力が高ければ解答できる」とは考えに くいです。英語で書かれた問題文を理解するだけの英語力と，計算を実行できるだけの数学 力の両方が揃って初めて正答できるでしょう。このような状況では，一方の特性値の低さを もう一方の特性値の高さでカバーすることはできません。こうした次元間の関係を反映した モデルを非補償 (noncompensatory) モデルあるいは部分補償 (partially compensatory) モデルと呼びます*14 。 最もシンプルな非補償モデルの項目反応関数は 𝑇 𝑃 (𝑦𝑝𝑖 = 1) = ∏ 𝑡=1 1 1 + exp [−𝛼𝑖𝑡 (𝜃𝑝𝑡 − 𝛽𝑖𝑡 )] (8.34) と表されます。各次元に関して独立に算出された「その次元の閾値を超える確率」を計算し， 最後にすべての積をとっています。そのため，一つでも特性値の低い次元があれば反応確率 は一気に小さくなってしまう，というわけです。非補償モデルの項目反応曲面は図8.23の左 のようになります。図8.22と比べると，両方の次元の 𝜃𝑝 が高くないと 𝑃 (𝑦𝑝𝑖 = 1) が高くな らないことが見て取れます。 図8.23 *14 30 非補償モデルの項目反応曲面 (Reckase, 2019, Fig. 12.2) 完全には非補償ではない（ある程度なら別の次元の特性値でリカバーできる）という意味で「部分補償モデ ル」のほうがより正しい名前らしいですが，そこまでこだわるほどの違いはありません。

8.7 多次元モデル 補償モデルと非補償モデルを比べると，補償モデルのほうがよく用いられているような印象 です。その理由の一つとして，非補償モデルのほうがパラメータの推定が難しいという点が 挙げられると思います。上の式を見るとわかるように，非補償モデルでは困難度パラメータ 𝛽𝑖𝑡 が次元の数だけあるためそもそものパラメータ数が補償型よりも多くなっています。ま た，非補償モデルは項目反応関数の全体の形が積になっています。コンピュータは足し算よ りも掛け算のほうが苦手なので，そういった意味でも計算の難易度が高いのです。また，心 理尺度的には補償型のほうが従来の因子分析モデルと同じ形になっているため，扱いやすい というのもあると思います。とはいえ，多次元の考え方として非補償型もあるんだというこ とを頭の片隅に置いておくと，いつかどこかで何かの参考になるかもしれません。 R でやってみる ここからは補償型の多次元 IRT モデルを mirt で実行する方法を紹介します*15 。補償型モ デルはカテゴリカル因子分析と同値なので，因子分析と同じように「探索的モデル」と「検 証的モデル」の 2 種類があります。ただ探索的モデルでは因子の回転に応じて 𝜃𝑝𝑡 の意味が 変化してしまうため実用上はなかなか使いにくいような気もします。 【探索的補償型モデル】 探索的モデルを実行する場合には，引数 model に「次元数」を与え るだけです。itemtype はこれまでと同じように自由に選んでください。二値だったら 2PL とか，多値だったら graded とか gpcm とか。 探索的補償型モデル 1 # GRMでやってみましょう 2 # このあたりまで来ると推定にもそこそこ時間がかかるかもしれません。 3 result_expl_comp <- mirt(dat[, cols], model = 2, itemtype = "graded") 1 # 推定値を見る 2 # 多次元の場合引数IRTparsは使えないので注意 3 coef(result_expl_comp, simplify = T) 1 $items 2 a2 d1 d2 d3 d4 d5 Q1_1 -0.251 0.876 3.879 2.433 1.457 0.665 -0.785 4 Q1_2 -0.817 1.699 5.685 3.977 3.034 1.234 -1.208 5 Q1_3 -0.991 2.402 5.788 4.099 3.032 1.048 -1.891 6 Q1_4 -0.775 0.903 3.646 2.414 1.831 0.781 -0.465 7 Q1_5 -0.757 1.523 5.105 3.271 2.221 0.612 -1.628 8 Q1_6 -1.449 -0.152 4.668 3.195 2.067 0.537 -1.707 9 Q1_7 -1.645 -0.184 4.602 2.822 1.798 0.252 -2.011 10 Q1_8 -1.329 -0.047 4.263 2.559 1.625 0.028 -2.047 11 Q1_9 -1.918 -0.131 5.226 3.255 1.695 0.508 -1.501 12 Q1_10 -1.436 *15 31 a1 3 0.000 2.864 1.361 0.075 -0.595 -2.040 非補償型モデルは，二値ならば itemtype='PC2PL'または itemtype='PC3PL'で実行可能です（PC は partially compensatory の頭文字）。ただパラメータの推定は結構安定しないので利用する際は気をつけて ください。

8.7 多次元モデル 13 14 $means 15 F1 F2 16 0 0 17 18 $cov 19 F1 F2 20 F1 1 0 21 F2 0 1 結果を見ると，識別力パラメータが次元の数 (a1, a2) と，切片パラメータが（カテゴリ数-1） 個 (d1-d5) 推定されています。その下の $mean と $cov には，それぞれ 𝜃𝑝𝑡 の平均と共分散 行列が出ています。これを見る限り，まだ因子間相関が 0 なので回転前の因子負荷（識別力） が出ているようです。さらに詳しく見ると，最後の項目の a2 が 0.000 となっています。実 は mirt での探索的モデルの推定の際には，このように一部の因子負荷を 0 に固定した状態 で推定を行っています*16 。「探索的」と紹介していましたが，実際には「モデルが識別でき る最小限の制約を置いて推定する」ということを行っているのです。IRT に限らず CFA も そうなのですが，複数の因子がある場合に識別性を持たせるには，最低でも 𝑇 (𝑇 −1) 2 個の因 子負荷を固定する必要があることが知られています。そのため，2 因子の場合は第 2 因子の 因子負荷を一つ 0 にしているわけです。同様に，3 因子の場合にはさらに第 3 因子の因子負 荷を 2 つ 0 にする必要があり，4 因子の場合にはさらにさらに第 4 因子の因子負荷を 3 つ 0 にする必要があり…という要領です。 というわけで，いくら探索的とは言えこのままでは各次元の解釈が出来ないので，回転しま しょう。回転させる場合は，psych::fa() の時と同じように引数 rotate を指定してあげ ます。 psych::fa() のデフォルトである oblimin 回転をしてみる 1 coef(result_expl_comp, rotate = "oblimin", simplify = TRUE) 1 2 Rotation: 1 $items 2 32 oblimin a1 a2 d1 d2 d3 d4 d5 3 Q1_1 -0.146 0.953 3.879 2.433 1.457 0.665 -0.785 4 Q1_2 0.037 1.872 5.685 3.977 3.034 1.234 -1.208 5 Q1_3 -0.107 2.634 5.788 4.099 3.032 1.048 -1.891 6 Q1_4 0.351 1.021 3.646 2.414 1.831 0.781 -0.465 7 Q1_5 0.057 1.680 5.105 3.271 2.221 0.612 -1.628 *16 psych::fa() のように固有値分解を使えばいいじゃないと思われるかもしれませんが，mirt() 関数として 他の様々な IRT モデルと同じ方法でパラメータを推定するためにはこうするのが最も自然なのです。

8.7 多次元モデル 8 Q1_6 1.474 -0.050 4.668 3.195 2.067 0.537 -1.707 9 Q1_7 1.679 -0.070 4.602 2.822 1.798 0.252 -2.011 10 Q1_8 1.311 0.051 4.263 2.559 1.625 0.028 -2.047 11 Q1_9 1.920 0.007 5.226 3.255 1.695 0.508 -1.501 12 Q1_10 1.394 0.110 2.864 1.361 0.075 -0.595 -2.040 13 14 $means 15 F1 F2 16 0 0 17 18 $cov 19 F1 20 F1 1.00 0.35 21 F2 0.35 1.00 F2 回転後の識別力は，いい感じに最初の 5 項目と最後の 5 項目に分かれてくれました。 【検証的補償型モデル】 検証的モデルを行う際には，引数 model に各次元と項目の関係を 記述してあげましょう。 1 # model式の右辺は列番号だけでもOKです 2 model <- " 3 f_1 = 1-5 4 f_2 = 6-10 5 " 6 7 result_conf_comp <- mirt(dat[, cols], model = model, itemtype = 8 coef(result_2dim, simplify = TRUE) 1 $items "graded") 2 a1 d1 d2 d3 d4 d5 3 Q1_1 0.895 0.000 3.865 2.425 1.453 0.665 -0.779 4 Q1_2 1.882 0.000 5.693 3.979 3.032 1.228 -1.212 5 Q1_3 2.590 0.000 5.789 4.095 3.026 1.042 -1.888 6 Q1_4 1.109 0.000 3.575 2.353 1.780 0.750 -0.465 7 Q1_5 1.701 0.000 5.118 3.275 2.219 0.605 -1.632 8 Q1_6 0.000 1.463 4.670 3.197 2.070 0.537 -1.712 9 Q1_7 0.000 1.606 4.538 2.779 1.770 0.249 -1.980 10 Q1_8 0.000 1.330 4.263 2.557 1.621 0.025 -2.047 11 Q1_9 0.000 1.979 5.311 3.310 1.719 0.509 -1.532 12 Q1_10 0.000 1.421 2.853 1.352 0.071 -0.596 -2.032 13 14 $means 15 f_1 f_2 16 33 a2 0 0

8.7 多次元モデル 17 18 $cov 19 f_1 f_2 20 f_1 1 0 21 f_2 0 1 きちんと項目 1-5 の次元 2 の識別力（a2）と項目 6-10 の次元 1 の識別力（a1）が 0 に固定 されています。先程の探索的モデルと見比べても，因子の順序こそ違えど概ね同じような推 定値が得られていますね。 ちなみに，探索的モデルでも検証的モデルでも，ICC は一次元モデルの時と同じように itemplot() によって出すことが出来ます。 多次元でも plot() 1 # 多値型項目だとすべてのカテゴリのIRCCCが重なるのでめっちゃ見にくい 2 # 全項目やると重いので，which.itemで少しずつやるのがおすすめ 3 plot(result_conf_comp, type = "trace", which.item = 4) Q1_4 1.0 0.8 1.0 0.8 0.6 0.6 P (q ) 0.4 0.2 0.4 0.0 6 4 2 q2 0 -2 -4 -6 -6 -4 -2 0 2 4 6 0.2 q1 0.0 図8.24 項目反応カテゴリ特性曲線 多次元でも itemplot() 1 # 引数rotを使うとグラフを回すこともできる 2 # 項目4はa2=0なので，a1方向に見るとIRCCCが見やすい 3 plot(result_conf_comp, type = "trace", which.item = 4, rot = list(xaxis = -90, yaxis = 0, zaxis = 0)) 34

8.8 適合度 Q1_4 1.0 1.0 0.8 0.8 0.6 (q ) 0.6 0.4 0.4 0.2 0.2 0.0 -6 -4 -2 0 2 4 64 20 -2 -4 6 -6 q2 q1 図8.25 0.0 グラフを回してみる 8.8 適合度 IRT でも SEM と同じように適合度を考えることができます。特に「個人」と「項目」の両 方に関心があることの多い IRT では，それぞれに対しての適合度を考える必要があります。 もちろんこれから紹介する適合度の考え方は，そのまま因子分析や SEM にも応用しようと 思えばできるので，そういった意味でも参考にしてみてください。 8.8.1 局所独立性の確認 まずはじめに，データとモデルの適合度の一つとして，データが局所独立の仮定に適合して いるかをチェックします。8.5.1項で説明したように，局所独立性が満たされている場合には， 𝜃𝑝 で条件づけたときの 2 項目の回答の間に相関が無いはずです。 𝑋 2 統計量 𝑋 2 統計量 (Chen and Thissen, 1997) では，カテゴリ変数のクロス表に対して用いられる 𝜒2 検定の枠組みを利用して局所独立性を評価します。ある二値項目 𝑖 の項目パラメータ (𝛼𝑖 , 𝛽𝑖 ) が分かっている場合，その項目の母集団全体での期待正答率は ∞ 𝐸[𝑃 (𝑦𝑝𝑖 = 1)] = ∫ 𝑃 (𝑦𝑝𝑖 = 1|𝜃)𝑓(𝜃)𝑑𝜃 (8.35) −∞ という形で求めることができます。ここで 𝑓(𝜃) は特性値の確率分布（ふつうは標準正規分 布）における 𝑥 = 𝜃 での確率密度を表しています。同様に，その項目の誤答率の期待値を求 めるならば 𝑃 (𝑦𝑝𝑖 = 1|𝜃) の部分を 1 − 𝑃 (𝑦𝑝𝑖 = 1|𝜃) にしたら良いだけです。 この考え方を 2 項目に広げると，母集団全体において項目 (𝑖, 𝑗) に両方とも正解する確率は もし2項目が独立ならば ∞ 𝐸[𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 1)] = ∫ −∞ 35 𝑃 (𝑦𝑝𝑖 = 1|𝜃)𝑃 (𝑦𝑝𝑗 = 1|𝜃)𝑓(𝜃)𝑑𝜃 (8.36)

8.8 適合度 と表せます。ということは，サンプルサイズが 𝑁 の場合，項目 (𝑖, 𝑗) に両方とも正解する人 数の期待値 𝐸11 は 𝑁 × 𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 1) です。同様にして，一方の項目にだけ正解する 人数および両方の項目に誤答する人数の期待値も求めると，表8.3のようになります。 表8.3 期待度数（独立な 2 項目の分布） 項目 𝑗 当てはまらない 当てはまらない 項目 𝑖 当てはまる 当てはまる 計 𝐸00 = 𝐸01 = 𝑁 × 𝑃 (𝑦𝑝𝑖 = 0 ∧ 𝑦𝑝𝑗 = 0) 𝑁 × 𝑃 (𝑦𝑝𝑖 = 0 ∧ 𝑦𝑝𝑗 = 1) 𝐸10 = 𝐸11 = 𝑁 × 𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 0) 𝑁 × 𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 1) - - 計 1000 この期待度数分布はもし2項目が局所独立ならばこうなるだろうという状態を表しており， 表8.1に相当するものを 𝜃 について周辺化して作成したものといえます。つまり局所独立では なくなるほど，実際のデータでのクロス表は（表8.2のように）この期待度数から離れるだろ う，ということです。 クロス表での連関の検定は（学部の統計でやっているかもしれない）𝜒2 検定です。その検定 統計量は 1 1 𝜒2 = ∑ ∑ 𝑠=0 𝑡=0 (𝑂𝑠𝑡 − 𝐸𝑠𝑡 )2 𝐸𝑠𝑡 (8.37) で計算されます。ここで 𝑂𝑠𝑡 は，実際のデータで 𝑦𝑝𝑖 = 𝑠 ∧ 𝑦𝑝𝑗 = 𝑡 だった人数を表してい ます。 こうして計算された 𝜒2 統計量は，自由が (行数-1) × (列数-1) の 𝜒2 分布に従うことが知ら れています。二値データの場合は 1 × 1 = 1 ということです。そしてここまでのプロセスは そのまま多値項目に対しても拡張可能です。例えば dat の項目のように 6 件法どうしの場合 は，自由度は 5 × 5 = 25 となるわけです。ということで，項目のペア毎に 𝜒2 統計量を計算 して統計的仮説検定を行い，有意だったペアの間は局所独立ではない可能性が疑われる，と いうことになります。 𝑄3 統計量 𝑄3 統計量 (Yen, 1984) では，回帰分析でいうところの偏相関係数にあたるものを利用しま す。偏相関係数は「変数 𝑧 の影響を取り除いた時の 𝑥 と 𝑦 の相関係数」で，疑似相関の影響 を検討する際などに利用されるものです。これを IRT の文脈に当てはめると「変数 𝜃𝑝 の影 響を取り除いた時の 𝑦𝑝𝑖 と 𝑦𝑝𝑗 の相関係数」ということで，これが 0 であれば局所独立だろ うと言えそうです。 回帰分析では，𝑥 と 𝑦 をそれぞれ 𝑧 で回帰した時の予測値 𝑥,̂ 𝑦 ̂ に対して，残差の相関 𝑟(𝑥−𝑥,𝑦− ̂ 𝑦)̂ のことを偏相関係数と呼んでいました。IRT（ロジスティックモデル）もその中身 はロジスティック回帰なので，同じようにして 𝑦𝑝𝑖 と 𝑦𝑝𝑗 をそれぞれ 𝜃𝑝 で予測した時の予測 36

8.8 適合度 値（期待正答率）を 𝑃 ̂ (𝑦𝑝𝑖 = 1) = 1 (8.38) 1 + exp [−𝛼𝑖̂ (𝜃𝑝̂ − 𝛽𝑖̂ )] というように，𝛼𝑖 , 𝛽𝑖 , 𝜃𝑝 にそれぞれ推定値 𝛼𝑖̂ , 𝛽𝑖̂ , 𝜃𝑝̂ を入れることで計算が出来ます。する と，残差を 𝑑𝑝𝑖 = 𝑦𝑝𝑖 − 𝑃 ̂ (𝑦𝑝𝑖 = 1) と求めることができます。図8.26に 𝑑𝑝𝑖 を表してみまし た。この項目には A さんも B さんも正解しているため，ふたりとも 𝑦 軸の値は 1 です。です が他の全項目から推定された 2 人の 𝜃𝑝̂ はそれぞれ −2, 2 だったとします。この場合，A さん は他の項目の出来から考えるとこの問題にはほぼ間違えるはずなのに正解しています。とい うことで 𝑑𝑝𝑖 の値が大きくなっています。このようになる理由としては，やはり 𝜃𝑝 以外の別 の要因によって正解できた，と考えるのが妥当でしょう。これを全体に広げてみた時に，𝑑𝑝𝑖 が高い人ほど別の項目でも同様に 𝑑𝑝𝑗 の値が大きくなっているとしたら，この 2 つの項目に は何か正解に寄与する 𝜃𝑝 以外の別の要因がある，と考えることができるわけです。 1.00 A B ypi 0.75 dpi 0.50 0.25 0.00 -4 -2 0 q^p 図8.26 2 4 残差 ということで 𝑄3 統計量は 𝑄3 = 𝑟𝐝𝑖 ,𝐝𝑗 (8.39) という形になります。𝑄3 統計量は相関係数なので，サンプルサイズに依存しない効果量の 指標として見ることができ，絶対値が 0.2 を超えてくると怪しいと判断できるようです。 R でやってみる mirt では residuals という関数でいま紹介した指標を出してくれます。引数 type で，出 してほしい統計量を指定できます。ちょっと厄介なのですが，𝜒2 統計量を出してほしい場合 には type = "LD" と指定してください。LD は local dependence の頭文字ですが，そうい ったら 𝑄3 だって LD だろ，と思われるかもしれません。これは Chen and Thissen (1997) の書き方に則っているのだと思います。我慢してください。 X2 統計量を出す 1 # 引数df.p = TRUEとすると，自由度とp値を出してくれる 2 residuals(result_2plm, type = "LD", df.p = TRUE) 1 Degrees of freedom (lower triangle) and p-values: 2 37

8.8 適合度 3 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 4 Q1_1 NA 0 0 0.019 0 0.000 0.000 0.093 0.824 0.004 5 Q1_2 1 NA 0 0.000 0 0.172 0.000 0.208 0.000 0.000 6 Q1_3 1 1 NA 0.000 0 0.006 0.003 0.000 0.000 0.000 7 Q1_4 1 1 1 NA 0 0.001 0.193 0.000 0.001 0.485 8 Q1_5 1 1 1 1.000 NA 0.003 0.000 0.004 0.000 0.019 9 Q1_6 1 1 1 1.000 1 10 Q1_7 1 1 1 1.000 1 1.000 11 Q1_8 1 1 1 1.000 1 1.000 1.000 12 Q1_9 1 1 1 1.000 1 1.000 1.000 1.000 13 Q1_10 1 1 1 1.000 1 1.000 1.000 1.000 1.000 NA 0.000 0.002 0.099 0.477 NA 0.032 0.019 0.550 NA 0.045 0.121 NA 0.000 NA 14 15 LD matrix (lower triangle) and standardized values. 16 17 18 19 Upper triangle summary: Min. 1st Qu. -0.109 -0.067 Median -0.026 Mean 3rd Qu. 0.002 Max. 0.063 0.219 20 21 Q1_1 Q1_2 Q1_3 Q1_4 Q1_8 Q1_9 22 Q1_1 NA 0.144 0.111 0.048 0.082 -0.098 -0.090 -0.034 0.005 23 Q1_2 50.489 NA 0.171 0.077 0.130 -0.028 -0.071 -0.026 -0.100 24 Q1_3 29.999 71.332 NA 0.102 0.219 -0.055 -0.061 -0.078 -0.109 25 Q1_4 5.504 14.503 25.105 NA 0.094 -0.067 -0.026 -0.079 -0.065 26 Q1_5 16.380 41.085 116.151 21.292 NA -0.059 -0.080 -0.058 -0.104 27 Q1_6 23.577 28 Q1_7 19.698 12.128 29 Q1_8 2.820 8.112 9.698 30 Q1_9 0.049 24.351 29.035 10.195 26.099 2.723 31 Q1_10 8.358 25.794 16.877 0.507 32 1.863 1.583 7.479 10.768 8.934 Q1_5 Q1_7 NA 0.101 0.063 0.033 1.695 15.588 24.903 NA 0.044 0.048 4.621 NA 0.041 5.519 4.012 NA 0.357 2.411 37.068 14.878 15.199 0.487 8.539 Q1_6 5.500 Q1_10 33 Q1_1 -0.059 34 Q1_2 -0.103 35 Q1_3 -0.083 36 Q1_4 0.014 37 Q1_5 -0.048 38 Q1_6 -0.014 39 Q1_7 0.012 40 Q1_8 0.031 41 Q1_9 0.123 42 Q1_10 NA 上の行列の上三角が 𝑝 値です。したがってこれが 0.05 を下回っている場合は「局所独立で はない」という判断になるのですが，見た感じかなり多くのペアで有意になっています。こ れは，統計的仮説検定なのでサンプルサイズの影響を受けているためです。また上の行列の 下三角は自由度です。今回は二値に変換したデータに対して行っているので自由度は 1 にな っています。もしも 6 件法のまま (GRM や GPCM で) 推定した結果に対して出した場合に 38

8.8 適合度 は，自由度は先程説明した通り 5 × 5 = 25 となります。 そして下の行列は実際に計算された 𝜒2 統計量（下三角）および標準化した値（上三角；た ぶんクラメールの連関係数）です。ということで，この行列の上三角について大きいペアが， 局所独立からより強く離れているといえます。 続いて 𝑄3 統計量です。 Q3 統計量を出す 1 # 仮説検定検定は無いので引数df.pは書いても無視される 2 residuals(result_2plm, type = "Q3") 1 Q3 summary statistics: 2 3 Min. 1st Qu. -0.235 -0.166 Median -0.101 Mean 3rd Qu. -0.067 0.033 Max. 0.244 4 5 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 6 Q1_1 1.000 0.150 0.109 0.033 0.074 -0.149 -0.148 -0.072 -0.027 7 Q1_2 0.150 1.000 0.188 0.056 0.131 -0.101 -0.173 -0.097 -0.209 8 Q1_3 0.109 0.188 1.000 0.082 0.244 -0.149 -0.170 -0.182 -0.235 9 Q1_4 0.033 0.056 0.082 1.000 0.071 -0.166 -0.122 -0.186 -0.175 10 Q1_5 0.074 0.131 0.244 0.071 1.000 -0.152 -0.197 -0.152 -0.226 11 Q1_6 -0.149 -0.101 -0.149 -0.166 -0.152 1.000 0.054 12 Q1_7 -0.148 -0.173 -0.170 -0.122 -0.197 0.054 1.000 -0.041 -0.047 13 Q1_8 -0.072 -0.097 -0.182 -0.186 -0.152 0.005 -0.041 14 Q1_9 -0.027 -0.209 -0.235 -0.175 -0.226 -0.048 -0.047 -0.045 1.000 15 Q1_10 -0.115 -0.193 -0.187 -0.068 -0.145 -0.114 -0.105 -0.069 0.040 16 0.005 -0.048 1.000 -0.045 Q1_10 17 Q1_1 -0.115 18 Q1_2 -0.193 19 Q1_3 -0.187 20 Q1_4 -0.068 21 Q1_5 -0.145 22 Q1_6 -0.114 23 Q1_7 -0.105 24 Q1_8 -0.069 25 Q1_9 0.040 26 Q1_10 1.000 ところどころ高めの相関が出ており，完全な局所独立ではないといえそうです。というのも 今回のデータは，本来 2 因子に分かれるはずの 10 項目を無理やり 1 因子とみなして分析し ています。そもそもデータが一次元性を持っていないはずなので，局所独立性も保たれては いない，ということになります。 （本来こういう時には一旦立ち止まり，項目の内容を吟味し たりモデルを変更したりといった手段を考えるべきですが，ここでは気にせず進めます。） ちなみに，引数 suppress を与えると，絶対値が suppress 以下の値が NA に置き換わるため， 39

8.8 適合度 項目がたくさんある場合でも「閾値を超えているペア」を探しやすくなるかもしれません。 1 residuals(result_2plm, type = "Q3", suppress = 0.2) 1 Q3 summary statistics: 2 3 Min. 1st Qu. -0.235 -0.166 Median -0.101 Mean 3rd Qu. -0.067 Max. 0.033 0.244 4 5 Q1_1 6 Q1_1 7 Q1_2 8 Q1_3 9 Q1_4 10 Q1_5 11 Q1_6 12 Q1_7 13 Q1_8 14 Q1_9 15 Q1_10 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 1 1.000 -0.209 1.000 0.244 -0.235 1.000 -0.226 1 0.244 1 1 1 -0.209 -0.235 -0.226 1.000 1 8.8.2 個人適合度 続いて個人の適合度です。IRT における個人適合度の基本的な考え方としては「𝜃𝑝 が高いの に困難度が低い項目に間違えるのはおかしい」「𝜃𝑝 が低いのに困難度が高い項目に正解する のはおかしい」というものです。するとこれは，図8.26的な考え方になります。A さんのほ うが予測値と実測値の間の差 𝑑𝑝𝑖 が大きいことは，前述の通り「𝜃𝑝 以外に回答に影響を与え る要因がある（のでモデルが間違っている） 」という可能性の他に，「（もしモデルは正しい と仮定したら）回答者の方に何か原因があるのでは」と考える，というわけです。 回答者の方の原因としては，心理尺度でいえば以下のようなものが考えられます (Meijer et al., 2016)。 • 回答者の言語能力が不十分なために項目の内容を十分に理解できていない可能性。低 学齢への実施や第一言語以外での実施で起こる可能性がある。 • 心理尺度が測定しようとしている心理特性に対する認識がない状態（Traitedness）。 例えば「そんなこと考えたこともない」といった場合，項目の意図を他の人とは異な るニュアンスで捉えてしまい，結果的に「大多数が当てはまると回答する項目に当て はまらないと回答する」といった反応につながる。 • 回答のモチベーションが低い可能性。やる気が無いので適当に回答している？ • 特定の心理特性を持っている人は平均的な人より一貫性の低い回答をする可能性があ る。例えば精神的な問題を抱えていたり，希死念慮を持ち合わせている人はそうでな い人と比べて適合度が低いという先行研究もあるらしい。 ということで，個人適合度が低いイコール適当な回答者（だから除外してもよい）とも言い 切れないのですが，特徴的な回答者をあぶり出すなどの目的でも個人適合度は使い道があり 40

8.8 適合度 そうです。 実際に個人適合度の指標として提案されているものは相当数ある (Meijer and Sijtsma, 2001; Karabatsos, 2003) のですが，わかりやすい＆統計的仮説検定のフレームワークが整 備されていると言った理由から，𝑧ℎ （またの名を 𝑙𝑧 ）という統計量 (Drasgow et al., 1984) が用いられることが多いようです (Meijer et al., 2016)。 𝑧ℎ (𝑙𝑧 ) 統計量 𝑧ℎ (𝑙𝑧 ) 統計量では，ある回答者の実際の反応パタン 𝐲𝑝 = [𝑦𝑝1 𝑦𝑝2 ⋯ 𝑦𝑝𝐼 ] と「考えられ る全反応パタン」での尤度を比較することにより，その回答者の反応パタンの「普通さ」を評 価します。まず，回答者 𝑝 の特性値の最尤推定値 𝜃𝑝̂ に関して，その時の尤度を計算します。 局所独立性が満たされているならば，尤度は単純に全ての項目をかけ合わせたら良いので， 𝐼 𝑙0 = ∏ 𝑃 (𝑦𝑝𝑖 = 1)𝑦𝑝𝑖 𝑃 (𝑦𝑝𝑖 = 0)1−𝑦𝑝𝑖 (8.40) 𝑖=1 となります。ただ掛け算は（コンピュータ的に）大変なので，対数尤度 𝐼 𝑙𝑙0 = ∑ 𝑦𝑝𝑖 ln 𝑃 (𝑦𝑝𝑖 = 1) + (1 − 𝑦𝑝𝑖 ) ln 𝑃 (𝑦𝑝𝑖 = 0) (8.41) 𝑖=1 を使います。この尤度 𝑙𝑙0 は，「特性値が 𝜃𝑝̂ の人が 𝐲𝑝 という回答をする確率」のようなも の*17 です。ここで 𝐲𝑝 以外の回答パタンに目を向けてみると，もしかしたら「𝐲𝑝 よりも 𝜃𝑝̂ と整合的な回答パタン」があるかもしれません。もちろん反対に「𝜃𝑝̂ と全く合わない回答パ タン」というものもあるかもしれません。 そこで，考えられる全回答パタン（2 件法 10 項目ならば 210 通り）の中で 𝐲𝑝 という回答パ タンが 𝜃𝑝̂ とどの程度整合的かを相対的にチェックすることを考えます。全回答パタン内で 𝐲𝑝 という回答パタンと 𝜃𝑝̂ の整合性（尤度）が平均くらいであれば，その回答パタンはまあ 起こりうるものだろうとみなせる一方で，もし整合性が低い場合には，𝐲𝑝 という回答パタン は起こりにくい（ズレが大きい）だろうと判断できる，ということです。 実際に「全回答パタンでの対数尤度の期待値」は（各項目への回答は離散変数なので） 𝐼 𝐸ℎ = ∑ 𝑃 (𝑦𝑝𝑖 = 0) ln 𝑃 (𝑦𝑝𝑖 = 0) + 𝑃 (𝑦𝑝𝑖 = 1) ln 𝑃 (𝑦𝑝𝑖 = 1) (8.42) 𝑖=1 として求めることができます。つまり 𝑙𝑙0 − 𝐸ℎ が「平均的な回答パタンでの尤度と実際の回 答パタン 𝐲𝑝 での尤度の差」であり，これがマイナス方向に大きいほどモデルとデータの適 合度が悪いということを意味します。 ここで全回答パタンにおける対数尤度の分散を 𝐼 𝜎ℎ2 = ∑ {𝑃 (𝑦𝑝𝑖 = 1)𝑃 (𝑦𝑝𝑖 = 0) [ln 𝑖=1 *17 41 2 𝑃 (𝑦𝑝𝑖 = 1) ] } 𝑃 (𝑦𝑝𝑖 = 0) (8.43) 正確に言えば， 「𝐲𝑝 という回答をした人の特性値として 𝜃𝑝̂ という値がどの程度もっともらしいか」という感 じですが，あえて曲解させています。

8.8 適合度 として求めることができるのですが，これを用いると，𝑙𝑙0 を全回答パタン内で標準化した 値 𝑧ℎ = 𝑙𝑙0 −𝐸ℎ 𝜎ℎ が経験的に標準正規分布に従うということが知られています。ということで， |𝑧ℎ | > 1.96 となった回答者はどうやら「特性値が 𝜃𝑝̂ のもとでは平均的ではない回答パタン である」と判断することができるわけです。 図8.27に，𝑧ℎ 統計量の考え方を表しました。ヒストグラムは母集団からサンプリングされた 様々な 𝜃𝑝̂ の人で計算した 𝑧ℎ をまとめたものです。これが近似的に正規分布に従っていると いうことで，この外側 5% を棄却域とみなして検定を行うことができます。B さんの 𝑧ℎ は， この分布の中で割と平均的なところにあるため，B さんの回答は「 （B さんの推定値 𝜃B̂ にお ける）ふつう」に近いとみなせます。一方 A さんの 𝑧ℎ は，この分布の中でもかなり小さい ̂ からすると当ては 方にあります。ということは，A さんの回答パタンは A さんの推定値 𝜃A まりが悪く（つまり困難度の低い項目に「当てはまらない」と回答し，困難度の高い項目で 「当てはまる」と回答している可能性が高く），なんなら母集団全体の中でも相当悪いほうだ と判断でき，何らかの問題があるのではないかと疑いをかけられるわけです。 図8.27 個人適合度の考え方 ちなみに，正規分布的には上側の外れ値も「ふつうではない」ということになりますが，こ ちらはどう考えたら良いかよくわかりません。理論的には「特性値 𝜃𝑝̂ に対してあまりにも適 合しすぎている（例えば困難度が 𝜃 ̂ 以下の項目では全て「当てはまる」と回答し，𝜃 ̂ 以上の 項目ではすべて「当てはまらない」と回答しているなど） 」という感じなのですが，だからと いって即座に問題とは言えなさそうです。ということで，基本的には 𝑧ℎ < −1.96 の回答者 に着目しておけば良いと思います。 R でやってみる ※ 個人適合度と回答パタンの関係がわかりやすいので，ここだけ GRM で推定した結果を使 用します。 mirt では，personfit() という関数によって 𝑧ℎ 統計量を出すことが出来ます。 個人適合度を出す 1 42 personfit(result_grm)

8.8 適合度 1 outfit z.outfit infit z.infit Zh 2 1 0.3043410 -2.51418387 0.2901436 -2.630507534 1.11739286 3 2 0.4319025 -1.53744653 0.4991522 -1.313386687 1.54458648 4 3 0.5417319 -1.17103185 0.6082584 -0.964496321 1.30478354 5 4 0.6961604 -0.71567303 0.7021298 -0.707746942 0.38159507 6 5 0.5979698 -0.95977207 0.5796844 -1.045595986 0.06132109 7 6 1.6611227 8 7 0.2629479 -2.06106332 0.2970130 -1.974928887 9 8 1.2820317 10 9 0.4296002 -0.79327083 0.4483505 -0.806286109 0.88382209 11 10 0.6018210 -0.83208882 0.6646485 -0.682369844 0.51748259 12 11 0.7508775 -0.39482320 0.8851335 -0.095554151 0.66704225 13 12 0.9285353 0.04354098 1.1524578 0.453682964 0.30674141 14 13 1.9537886 1.83531480 1.9490471 1.865766458 -1.36629763 15 14 1.0385754 0.23709522 1.0081517 0.168761287 -0.60996496 16 15 0.8936027 -0.09915012 0.9407077 0.005928737 0.05273411 17 16 0.6633546 -0.62999115 0.7489752 -0.432539405 0.73371137 18 17 0.5941524 -1.07615247 0.6027536 -1.058618799 0.98045292 19 18 1.7516559 1.78554390 1.7275073 1.788745281 -1.75844255 20 19 1.4892405 1.27332749 1.5781245 1.468740378 -1.18728452 21 20 1.3794907 0.81105690 1.5238636 1.050442554 -0.07950879 22 0.99612493 2.2816175 0.82619048 1.2682515 1.585552919 -0.57702558 1.15924498 0.803602305 -1.12899488 [ reached 'max' / getOption("max.print") -- omitted 2412 rows ] outfit や infit というものは，𝑧ℎ よりもシンプルに予測値と実測値の間の差 𝑑𝑝𝑖 の二乗和 をもとにした適合度指標 (Smith et al., 1995) です*18 。ということでこれらについても大き いほど適合度が悪いと判断できます。実際にこれらの指標の相関を見てみると，いずれもか なり高い相関になっていることがわかります。 指標間の相関 1 # z_hだけは方向が違う（小さいほど悪い）ので他の指標と負の相関 2 cor(personfit(result_grm)) 1 outfit z.outfit infit z.infit Zh 2 outfit 1.0000000 0.9029153 0.9767846 0.8895739 -0.8593484 3 z.outfit 0.9029153 1.0000000 0.8962045 0.9898224 -0.8803762 4 infit 0.9767846 0.8962045 1.0000000 0.9083434 -0.8420095 5 z.infit 0.8895739 0.9898224 0.9083434 1.0000000 -0.8617816 6 Zh -0.8593484 -0.8803762 -0.8420095 -0.8617816 1.0000000 では実際に，当てはまりが悪い人の回答パタンを見てみましょう。 *18 43 outfit は重み付けない二乗和，infit は情報量（後述）で重み付けた二乗和のようです。

8.8 適合度 1 # 引数stats.only = FALSEにすると，統計量に加えて回答パタンも残してくれる 2 pf <- personfit(result_grm, stats.only = FALSE) 3 # Zhの値が最小の人を見てみる 4 pf[which.min(pf[, "Zh"]), ] 1 2 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 581 3 4 1 1 infit 1 1 1 z.infit 6 6 6 6 outfit z.outfit 6 5.804625 5.127515 Zh 581 5.552154 5.122718 -7.275079 いま扱っているモデルは 10 項目 1 因子のモデルです。そして 10 項目の識別力は，いずれも 正の値になっています。つまり，基本的にはすべての項目間には正の相関があるため，この 回答者のように「ある項目では 1，別の項目では 6」という回答の付け方には違和感があり ます。 1 # Zhの値が-1.96より小さい人 2 subset(pf, Zh < -1.96) 1 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 outfit 2 45 6 5 3 2 3 3 6 3 6 4 1.362101 3 75 6 6 3 3 1 6 6 5 6 1 2.777980 4 96 6 6 4 1 1 6 4 5 5 6 2.782816 5 98 1 1 1 1 2 3 6 3 4 3 2.006426 6 122 6 6 6 6 6 6 6 3 1 2 3.089210 7 123 6 6 6 5 6 6 1 3 3 2 2.065644 8 126 4 6 6 4 6 2 3 2 2 2 1.459869 9 172 5 2 1 6 3 6 6 6 6 6 4.695997 10 186 6 5 2 2 2 2 2 6 5 4 1.817030 11 226 1 1 1 4 6 5 6 5 6 6 4.784881 12 z.outfit infit z.infit Zh 13 45 0.8901956 1.401457 0.9775542 -2.193202 14 75 2.4882677 2.751967 2.5542871 -3.279082 15 96 2.5103934 2.855421 2.6749604 -2.107685 16 98 2.2242279 2.099394 2.4577322 -3.172528 17 122 2.3332034 3.034305 2.3832158 -2.030530 18 123 1.7906200 1.955780 1.7068100 -1.971380 19 126 1.1442225 1.314690 0.8576219 -1.980124 20 172 3.3697866 3.651095 2.8513846 -3.188176 21 186 1.8103247 1.791534 1.7810377 -2.201609 22 226 3.9067830 4.276645 3.7176248 -4.057956 23 [ reached 'max' / getOption("max.print") -- omitted 121 rows ] 同様にして，𝑧ℎ < −1.96 の人たちを見てみると，項目ごとにかなりバラバラな回答をしてい る人たちがあぶり出されています。ということで，𝑧ℎ 統計量によって怪しい回答者をあぶり 出すことができそうだということがわかりました。 44

8.8 適合度 ちなみに反対に 𝑧ℎ > 1.96 の人たちを見たところですべての項目で同じような値をつけてお り，やはり「適合度が悪い」という判断はなかなか難しそうです（ストレートライン気味な 感じはありますが）。 1 # Zhの値が1.96より大きい人 2 subset(pf, Zh > 1.96) 1 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 outfit 2 119 5 4 4 4 4 4 4 4 3 2 0.1059937 3 136 6 5 5 6 5 5 5 5 5 5 0.1749942 4 221 2 2 1 1 2 2 1 2 2 1 0.3892825 5 674 5 5 5 4 5 5 4 4 5 2 0.2410515 6 684 4 4 5 5 4 4 4 4 3 2 0.1893306 7 694 5 5 5 4 5 5 4 4 3 2 0.2993084 8 739 5 5 5 6 4 4 4 4 3 3 0.2791157 9 804 5 5 5 5 4 5 5 5 5 3 0.1527831 10 826 5 5 4 4 4 4 4 3 3 1 0.2959759 11 834 5 5 4 5 5 4 4 4 3 3 0.1765065 12 855 5 5 5 5 5 4 4 4 4 3 0.1063125 13 1055 5 5 4 4 5 5 4 4 4 2 0.1968235 14 1115 5 4 5 5 4 4 4 4 4 2 0.1644982 15 z.outfit infit z.infit Zh 16 119 -3.933211 0.1122620 -3.899088 2.556666 17 136 -2.032175 0.2246921 -1.890473 2.065611 18 221 -1.326634 0.4313157 -1.407809 2.117805 19 674 -2.270639 0.2776608 -2.150842 1.970796 20 684 -3.039619 0.2021914 -2.982029 2.001422 21 694 -2.123527 0.3096588 -2.123376 2.017472 22 739 -2.191891 0.3059683 -2.111215 1.971686 23 804 -2.557928 0.1559106 -2.653563 1.963533 24 826 -2.472784 0.3014073 -2.463963 2.015447 25 834 -2.919233 0.1741386 -3.000314 2.105461 26 855 -3.292733 0.1044904 -3.412590 2.040042 27 1055 -2.751655 0.2080324 -2.738016 1.975189 28 1115 -3.111734 0.1829595 -3.020623 2.079906 29 [ reached 'max' / getOption("max.print") -- omitted 6 rows ] ここで紹介した 𝑧ℎ 統計量のような個人適合度は，モデルベースで適当回答者を除外するの に使うこともできます。もちろん適合度の値だけで機械的に除外するのは危ないですが，こ れによって一貫しない回答をしている人などをある程度除外することができれば，項目パラ メータの推定もより安定するようになるでしょう。 45

8.8 適合度 8.8.3 項目適合度 IRT では項目一つ一つの性能（識別力や困難度）に関心があることが多いため，適合度も項 目レベルで考えます*19 。基本的な考え方としては，推定されたパラメータに基づく ICC が， データとどの程度一致するかを確認していきます。適合度が悪い項目が見つかった場合には， 項目の内容を吟味して削除するかを検討したり，場合によってはモデル自体を修正（因子数 や因子構造の修正など）する必要があるかもしれません。項目適合度は，具体的には以下の ような手順によって算出していきます (Orlando and Thissen, 2000) 。 1. 項目パラメータと 𝜃𝑝 を推定する 2. 回答者を 𝜃𝑝 の値で並び替える 3. 回答者を 𝜃𝑝 の値によっていくつかのグループに分ける 4. 各グループでの反応確率を計算する 5. 4. で求めた反応確率と「モデル上の期待反応確率」を 𝜒2 的な統計量や視覚的方法を 用いて比較する 視覚的なチェック ということで，まずは視覚的な比較を見てみましょう。mirt では，項目適合度を見るための 関数として itemfit() が用意されています。これに引数 empirical.plot を与えることで プロットを見ることができます。 項目適合度を視覚的にチェック 1 # 引数empirical.plotに項目番号か列名を指定します 2 itemfit(result_2plm, empirical.plot = 8) 図8.28を見ると，各カテゴリごとに ICC が表示されています。これに重なっている小さいマ ルの Y 座標は，手順 3 で分割した各グループ（ここでは 10 グループ）における 𝑃 (𝑦𝑝𝑖 = 1) を表しています（多値型の場合はカテゴリ選択確率が表示される） 。また，X 座標はそのグル ープにおける 𝜃𝑝 の平均値を表しています。ということで，各カテゴリにおける実線に小さい マルが近いほどデータとモデルの当てはまりが良い，と判断することができるわけです。 ただカテゴリ数や項目数が多くなってくるとすべてを目視で確認するのはかなり大変です。 また視覚的なチェックだけではどの項目が最も乖離しているのかを判断するのも大変でしょ う。ということで，「𝜒2 的な統計量」にもとづいて怪しい項目にアタリをつけるのが良さそ うです。 カイ二乗的な統計量 𝜒2 統計量は先ほど局所独立性の検証で登場したものです。つまり何らかの方法によってグル ープに分けられたクロス表に対する期待度数 𝐸 と観測度数 𝑂 について， 𝜒2 = ∑ *19 46 (𝑂 − 𝐸)2 𝐸 (8.44) もちろんモデル全体での適合度を考えることも可能です。例えば 1 パラメータと 2 パラメータの両モデルを lavaan で実行して，compareFit() することでどちらのモデルが良さそうかを判断したりというプロセス が考えられます。

8.8 適合度 Empirical plot for item 8 1.0 0.8 P (q ) 0.6 0.4 0.2 0.0 -4 -2 0 2 4 q 図8.28 ICC とデータの反応確率 と計算したものが，帰無仮説：モデルがデータに適合しているならば小さい値におさまるは ずだという考え方です。 Bock (1972) は，𝜒2 統計量を計算するために「頻出回答パターン」によるグルーピングを行 いました。つまり項目 𝑖 = 1 の適合度を見るときには「項目 1 以外の各項目の回答パターン」 ×「項目 1 の正誤」でクロス表を作り，十分に観測度数 𝑂 が計算できるだけの数が揃ってい るカテゴリについてズレの総和をとったのです。 Yen (1984) はもう少しこれをシンプルに考えて，𝜃𝑝 の値によってグループ分けすることに しました。さらに Yen は，グループの数も 10 に固定して算出した指標を 𝑄1 と名付けまし た。グループ 𝑔 (𝑔 = 1, 2, ⋯ , 10) における特性値 𝜃𝑝̂ の平均値を 𝜃𝑔̄ とおくと，期待正答率は （2 パラメータロジスティックモデルならば） 𝐸𝑖𝑔 = 1 (8.45) 1 + exp [−𝛼𝑖̂ (𝜃𝑔̄ − 𝛽𝑖̂ )] と求めることができます。すると，二値項目における 𝑄1 統計量は 2 10 [(1 − 𝑂𝑖𝑔 ) − (1 − 𝐸𝑖𝑔 )] (𝑂𝑖𝑔 − 𝐸𝑖𝑔 )2 ⎞ ⎟ + 𝐸𝑖𝑔 1 − 𝐸𝑖𝑔 ⎠ 10 2 2 (1 − 𝐸𝑖𝑔 )(𝑂𝑖𝑔 − 𝐸𝑖𝑔 ) + 𝐸𝑖𝑔 (𝑂𝑖𝑔 − 𝐸𝑖𝑔 ) ) = ∑ 𝑁𝑔 ( 𝐸𝑖𝑔 (1 − 𝐸𝑖𝑔 ) 𝑔=1 ⎜ 𝑄1 = ∑ 𝑁𝑔 ⎛ 𝑔=1 ⎝ 10 = ∑ 𝑁𝑔 ( 𝑔=1 (8.46) (𝑂𝑖𝑔 − 𝐸𝑖𝑔 )2 ) 𝐸𝑖𝑔 (1 − 𝐸𝑖𝑔 ) となります。なお 1 行目の右辺は，第 1 項が正答セルにおける乖離度を表しています。同様 に，第 2 項は (1 − 𝑂𝑖𝑔 ), (1 − 𝐸𝑖𝑔 ) に置き換わっていることから，誤答セルにおける乖離度を 表していることがわかります。したがって，𝑄1 統計量は「グループ」×「項目 𝑖 の正誤」と 47

8.8 適合度 いう 10 × 2 クロス表における 𝜒2 統計量を表しているのです。そしてこの 𝑄1 統計量は自由 度が 10− パラメータ数の 𝜒2 分布に従うため，これを用いて検定ができるわけです。 同様に，差ではなく比を用いた統計量として 𝐺2 統計量というものもあります。 10 𝐺2 = ∑ 𝑁𝑔 [𝑂𝑖𝑔 ln 𝑔=1 𝑂𝑖𝑔 1 − 𝑂𝑖𝑔 + (1 − 𝑂𝑖𝑔 ) ln ] 𝐸𝑖𝑔 1 − 𝐸𝑖𝑔 (8.47) ほかにもいくつか適合度指標は提案されているのですが，基本的には前述の通り「期待反応 確率と実際の反応確率の差」を何らかの形で統計量に落とし込んでいます。細かな設定と しては 1. グループの数 2. 𝜃𝑝 の推定方法（基本的には最尤法，ただベイズ推定などでも可） 3. 各グループの特性値の代表値の決め方（平均値 or 中央値） 4. 期待反応確率の計算方法（個人の 𝜃 ̂ における期待反応確率 𝐸 𝑝 𝑖𝑝 のグループ平均とい う手もあり） といったところで様々なオプションが考えられるようです (加藤他, 2014)。 R でやってみる 項目適合度をを出す場合は itemfit() 関数を使います。引数 fit_stats に，どの統計量を 出してほしいかを指定します。なお，X2 を指定すると内部では 𝑄1 統計量（i.e., グループ数 10）を出してきます。 項目適合度をチェック 1 # まとめて複数の統計量を出すこともできます 2 itemfit(result_2plm, fit_stats = c("X2", "G2")) 1 item X2 df.X2 RMSEA.X2 p.X2 G2 df.G2 RMSEA.G2 p.G2 2 1 Q1_1 375.251 8 0.137 0 458.736 6 0.176 0 3 2 Q1_2 73.800 8 0.058 0 99.817 5 0.088 0 4 3 Q1_3 92.664 8 0.066 0 135.273 5 0.104 0 5 4 Q1_4 102.569 8 0.070 0 153.430 5 0.111 0 6 5 Q1_5 126.263 8 0.078 0 169.111 5 0.116 0 7 6 Q1_6 157.177 8 0.088 0 191.087 5 0.124 0 8 7 Q1_7 151.358 8 0.086 0 219.612 4 0.149 0 9 8 Q1_8 125.809 8 0.078 0 154.268 6 0.101 0 10 9 Q1_9 151.495 8 0.086 0 207.812 5 0.129 0 11 10 Q1_10 572.082 8 0.170 0 754.816 5 0.248 0 G2 の自由度 (df.G2) がいろいろな値になっていますがあまり気にしなくて OK です。今回 も例によって検定の結果 (p.*) は全て有意となりました。このように，サンプルサイズが大 きい場合これらの統計量はあまり参考にならない可能性が高いです。 …じゃあどうするんだ，という話になるわけですが，一つの方法としては効果量的な考え方 48

8.9 項目情報量・テスト情報量 をとるために「実際の差 (𝑂𝑖𝑔 − 𝐸𝑖𝑔 ) や比 えば 1 10 10 ∑𝑔=1 𝑂𝑖𝑔 𝐸𝑖𝑔 を計算する」という方法が考えられます。例 |𝑂𝑖𝑔 − 𝐸𝑖𝑔 | が大きい項目は，基本的に視覚的にプロットを見ても差があるで しょう。ただ (𝑂𝑖𝑔 − 𝐸𝑖𝑔 ) などを直接計算してくれる関数は用意されていないので，多少頑 張る必要があります。 8.9 項目情報量・テスト情報量 IRT にかぎらず，統計的推測では推測の精度が問われます。例えば正規分布の母平均の推測 2 の場合，標準誤差は母分散 𝜎2 およびサンプルサイズ 𝑛 を用いて √ 𝜎𝑛 と表されました。こ の式は，データの数が多いほど推測の精度が高くなるという自然な考えを表しています。 テストや心理尺度では，標準誤差を真値と誤差の関係から考えます。Chapter 4 で紹介しま したが，潜在特性の測定（古典的テスト理論）では，観測得点 𝑥 と真値 𝑡 と誤差 𝑒 の関係を 𝑥=𝑡+𝑒 (8.48) と表し，信頼性係数を 𝜌 =1− 𝜎𝑒2 𝜎𝑥2 (8.49) と定義していました。標準誤差は「真値が 𝑡 の人が測定を繰り返した場合に， （平均的に）ど の程度観測得点 𝑥 がばらつくか」を表しています。したがって，𝑥 ∼ 𝑁 (𝑡, 𝜎𝑒2 ) とすると，心理 測定における標準誤差は 𝜎𝑒 にあたります。信頼性係数 𝜌 に関しては，上の式を変形させると 𝜎𝑒 = 𝜎𝑥2 √1 − 𝜌2 (8.50) となることから，これが間接的に測定の精度を表す指標だということが言えます。なお実際 に標準誤差を計算する場合には，𝜌 の推定値として 𝛼 係数などを代入します。 さて，信頼性係数に基づく上記の標準誤差の式には真値 𝑡 は含まれていません。つまり真値 の値が何であれ標準誤差は同じということを意味しています。ですがこの考え方は割と不自 然なものです。特性値の真値が平均くらいの人では，尺度の項目はほとんどすべてがその人 の特性値を測定するのに意味のある（正解であるほど真値 𝑡 は高いだろうという判断につな がる）項目です。一方で，特性値の真値が非常に高い人では，実はたいていの項目はその人 の特性値を測定するのにほぼ無意味です。というのも，ある程度難易度の高い項目に正解す る人は，難易度の低い項目にはほぼ確実に正解するだろうと言えるためです。つまり，「あ る程度難易度の高い項目」の情報さえあれば「難易度の低い項目」の情報は無くても特に問 題ない（真値の推定には影響しない）と言えます。先程，標準誤差はサンプルサイズが多く なるほど小さくなる，ということを説明しました。本来，特性値の真値によって「その人の 特性値を測定するのに意味のある」項目の数は変わるはずであり，それに伴って標準誤差も 変わるはずなのです。 IRT では，この問題を克服するために項目情報量 (item information) およびテスト情報量 (test information) というものを用います。この項目情報量・テスト情報量は，後述するテ ストの設計などにおいても大きな力を発揮してくれます。 49

8.9 項目情報量・テスト情報量 8.9.1 情報量とは 項目情報量の話に入る前に，少しだけ「情報量」とは何かについてごくごく簡単に説明して おきましょう。我々は，一般名詞として「情報」という言葉を使うことが多々あります。こ の時の「情報」は，何かしらの事柄について「知らなかった」状態から「知っている」状態 に遷移させるものという見方ができます。 この「情報」を確率の世界に持ち込んで考えてみます。いま，A さんは 52 枚のトランプから 1 枚のカードを引きました。A さんの向かいに座っている B さんには，もちろん A さんのカ ードが何かは全くわかりません。つまりこの段階で B さんが A さんのカードを言い当てる 確率は 1/52 です。ここで，A さんから「スペードのカードを持っている」という情報が得 られたとします。すると，B さんの中にあったカードの候補は 52 枚から 13 枚に減ります。 この時点では，B さんが A さんのカードを言い当てる確率は一気に 4 倍になり 1/13 になっ ています。「スペードのカードを持っている」という情報は B さんの予測の精度を一気に 4 倍に引き上げたわけです。 ここで，B さんは次に「黒のカードを持っている」という情報を得たとします。この情報で は，B さんの中にあったカードの候補は 13 枚のまま変わりません。したがってこの情報は， 今の B さんにとっては何の価値も持っていません。一方で，B さんが「エースのカードを持 っている」という情報を得た場合には，前の情報と合わせて「スペードのエース」であるこ とが確定します。この場合には，B さんが A さんのカードを当てる確率は一気に 13 倍とな るわけです。 以上の例え話からは，情報量に関する以下の性質がわかります。 情報は予測の精度を高める 情報量が多くなるほど，正確に予測ができるようになります。 より不確実性を下げる情報のほうが情報量が多い マークは 4 通りである一方で，数字は 13 通りあるため，数字の情報のほうが不確実性をより大きく下げている＝情報量が多い と解釈できます。 独立な情報の情報量は加法的である マークと数字は互いに独立なので，先に「エースを持 っている」という情報を得たとしてもカードを当てる確率はやはり 13 倍になってい ました。 ここまでの話を IRT に置き換えてみます。𝜃𝑝 の予測に際して，項目がもつ情報量を考えて みましょう。図8.29には困難度がそれぞれ (0, 3) の項目が描かれています。 ある回答者 𝑝 は，それまでの項目への解答状況から 𝜃𝑝̂ = 0 だと予測されているとします。す ると，この回答者 𝑝 にとって，緑の項目は正解するかどうかがかなり不確実（50%）です。そ してこの項目への反応によって 𝜃𝑝̂ = 0 の予測は変化し，多少なりともより正確なものにな るでしょう。そう考えると，緑の項目は 𝜃𝑝̂ = 0 の人にとってそれなりの情報量を持っている と考えられます。 一方で赤い項目を見ると，𝜃𝑝̂ = 0 の人はほぼ確実に正解しないだろうという予想がつきま す。つまり，出題時点で赤い項目に対する不確実性はほぼゼロ，ということです。実際にこ 50

8.9 項目情報量・テスト情報量 P(ypi = 1) bi = 0 bi = 3 qp 0 図8.29 IRT における情報の考え方 の回答者 𝑝 が赤い項目に不正解だったとしても，「そうでしょうね」という話なので 𝜃𝑝 の予 測に対して特に新規の情報をもたらすものとは言えません*20 。同様に，緑の項目でも 𝜃𝑝̂ = 3 の人にとってはほぼ情報量のない項目，ということになります。このように，ある項目は，そ の項目の困難度が回答者の 𝜃 に対してちょうどよいほど，𝜃 の推定に対して多くの情報量を 持っているのです。ということである項目が持つ項目情報量は，実際には 𝜃 の値によって変 動する関数になっています。この関数のことを項目情報関数 (item information function [IIF]) と呼びます*21 。この節のはじめにお話したように，IRT では真値によって標準誤差が 変わるようにすることを考えます。この後具体的な関数の形を示しますが，項目情報関数を 𝜃 の関数として表すことによってその目的を達成しようというわけです。 8.9.2 項目情報関数 IRT における情報量の考え方のイメージを確認したので，具体的な項目情報関数の式を見て みましょう。二値型モデルの場合，項目情報関数 𝐼𝑖 (𝜃) は 𝐼𝑖 (𝜃) = 𝑃 ′ (𝑦𝑖 = 1)2 𝑃 (𝑦𝑖 = 1)(1 − 𝑃 (𝑦𝑖 = 1)) (8.51) という式になります。ここで 𝑃 ′ (𝑦𝑖 = 1) は，𝑃 (𝑦𝑖 = 1) を 𝜃 で微分した導関数です。…とい ってもよくわからないと思うので，実際に見てみましょう。図8.30は，図8.29に示した ICC を持つ 2 項目のロジスティックモデルにおける項目情報関数です。どちらの線も先程説明し たように，正解するかが最も不確実な 𝜃𝑝 = 𝛽𝑖 のところで情報量が最大になっています。 IIF の式 𝐼𝑖 (𝜃) = 𝑃 ′ (𝑦𝑖 = 1)2 𝑃 (𝑦𝑖 = 1)(1 − 𝑃 (𝑦𝑖 = 1)) (8.52) を求めるためには，𝑃 (𝑦𝑖 = 1) の導関数を求める必要があります。正規累積モデルは項目反 応関数の中に積分が含まれているためこれを求めるのはかなり大変なのですが，それと比べ *20 もちろんこの回答者 𝑝 が赤い項目に正解した場合には，𝜃𝑝 の予測を大きく動かすため，情報量が多くなりま す。ですがそもそも「回答者 𝑝 が赤い項目に正解する」という事象の発生確率がほぼゼロなので， 「この項目 への回答が 𝜃𝑝 の予測にもたらす情報量の期待値」という観点ではやはりほぼゼロなのです。 *21 項目情報「量」という場合には特定の 𝜃 における値を指す一方で，項目情報「関数」という場合には 𝜃 全域 を対象とすることが多い気がするので，使い分けに気をつけてください。 51

8.9 項目情報量・テスト情報量 Ij(q ) bi = 0 bi = 3 qp 3 0 図8.30 項目情報関数 るとロジスティックモデルでは解析的に IIF の式を求めることができます。例えば 2 パラメ ータロジスティックモデルの場合，IIF は 𝐼𝑖 (𝜃) = 𝛼2𝑖 𝑃 (𝑦𝑖 = 1)(1 − 𝑃 (𝑦𝑖 = 1)) (8.53) となります。この式からも，IIF は𝑃 (𝑦𝑖 = 1) = 0.5 となる𝜃 = 𝛽𝑖 の点において最大値をとる ことがわかります。同時に，識別力𝛼𝑖 が大きいほど項目情報量も多くなることがわかります。 図8.31に，識別力が異なる 2 つの項目 (𝛽𝑖 = 0; ICC は図8.8) のロジスティックモデルでの IIF を示しました。識別力が高い項目ほど，項目特性関数の傾きが大きくなっていました。そ のため 𝜃𝑝 が高い人は正解し 𝜃𝑝 が低い人は不正解する，というメリハリがついていました。 𝜃𝑝 の推定という視点から言えば，𝛼𝑖 が大きい項目ほどその項目への正誤が特性値の推定に より確かな情報を与えるという意味で，𝛼 と項目情報量には密接な関係があるわけです。 Ij(q ) ai = 1 ai = 0.5 qp 0 図8.31 異なる識別力を持つ項目の IIF 図8.31をよく見ると，𝜃𝑝 = 0 付近では識別力の高い緑の項目のほうが高い情報量を持ってい る一方で，𝜃𝑝 が 0 から離れたところではむしろ識別力の低い青い項目のほうが高い情報量 を持っています。図8.8を見ると，識別力 𝛼𝑖 が高くなるほど，𝜃𝑝 が 𝛽𝑖 から離れると急激に 𝑃 (𝑦𝑖 = 1) が 0.5 から離れていきます。つまり「正解するかが不確実」な 𝜃𝑝 の範囲は識別力 の高さと反比例しており，その結果 𝜃𝑝 が 𝛽𝑖 から離れたところでの項目情報量も少なくなっ てしまうのです。 52

8.9 項目情報量・テスト情報量 この事実は「識別力は高ければ高いほど良いわけではない」ことを教えてくれています。例 えば「𝜃𝑝 < 𝛽𝑖 の人は確実に間違えるが，𝜃𝑝 ≥ 𝛽𝑖 の人は確実に正解する」項目があったとす ると，その項目は図8.32のように識別力 𝛼𝑖 → ∞ になる項目です。 P(ypi = 1) ai ®¥ qp bi 図8.32 完全な識別力を持つ項目があったら この項目は「𝜃𝑝 が 𝛽𝑖 より上か下か」に関しては完全な情報を教えてくれます (𝐼𝑖 (𝜃𝑝 = 𝛽𝑖 ) = ∞)。その一方で，𝜃𝑝 が 𝛽𝑖 よりどの程度高かろうと 𝑃 (𝑦𝑖 = 1) は 100% なので， 「𝜃𝑝 が 𝛽𝑖 よ りどの程度上か下か」に関しては全くわかりません (𝐼𝑖 (𝜃𝑝 ≠ 𝛽𝑖 ) ≃ 0)。このように，識別力 が高すぎる項目は有効な 𝜃𝑝 の範囲が狭すぎるため，多くの回答者にとって意味のない項目 になってしまう可能性があるわけです。 8.9.3 テスト情報関数 ここまで，項目情報関数の定義およびその性質を見てきました。この「情報」は，𝜃𝑝 を推定 するための「情報」を表しているわけですが，𝜃𝑝 の推定は尺度・テスト内のすべての項目へ の回答を総合して行われるものです。そのため，次は項目情報関数を拡張して尺度全体での 情報関数を考えてみます。 情報量には「独立な情報の情報量は加法的である」という性質がありました。IRT でもこの 性質は健在です。もしテスト全体が局所独立の仮定を満たしているならば，テスト情報関数 (test information function [TIF]) は単純に 𝐼 𝐼(𝜃) = ∑ 𝐼𝑖 (𝜃) (8.54) 𝑖=1 と表せます。そしてテスト情報量の重要な性質として，真の特性値が 𝜃𝑝 の人における最尤推 定量 𝜃𝑝̂ の誤差分散は漸近的に 𝜎(2𝜃 ̂ 𝑝 |𝜃𝑝 ) = 1 𝐼(𝜃𝑝 ) (8.55) となることが知られています。ということで，標準誤差は 𝜎 (𝜃 ̂ 𝑝 |𝜃𝑝 ) 53 = 1 √𝐼(𝜃𝑝 ) (8.56)

8.9 項目情報量・テスト情報量 となります。これは古典的テスト理論の表記に合わせると，最尤推定値 𝜃𝑝̂ = 𝜃𝑝 + 𝜀𝑝 とした 場合の 𝜎𝜀 に相当するものです。このように，IRT では測定の精度を表す標準誤差を，𝜃𝑝 の 関数として求めることができるわけです。 R で確認してみる 項目・テスト情報関数の考え方がわかったところで，実際の項目での関数を確認してみまし ょう。これまで ICC を出すために使っていた plot() 関数によって，IIF も出すことができ ます。項目情報関数を出す場合には type = "infotrace" とします。 項目情報関数 1 plot(result_2plm, type = "infotrace", facet_items = FALSE) Item Information 0.4 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 I (q ) 0.3 0.2 0.1 0.0 -6 -4 -2 0 2 4 6 q 図8.33 項目情報関数 先ほど説明した通り，低識別力の項目 1 は他の項目と比べて圧倒的に情報量が少ないことが わかります。また、𝜃𝑝 > 2 のエリアでは、Q1_10 以外のほぼどの項目も情報量が無いため、 この 10 項目では特性値が高い人についてはまともな推定はできなさそうだということがわ かります。この項目情報関数を見れば，例えばどうしても 1 項目削除しないといけないとし たら，基本的には Q1_1 を削除したら良いが，もしも 𝜃𝑝 = −6 付近がメインターゲットだと したらむしろ Q1_1 は残したほうが良いとか、𝜃𝑝 > 2 がメインターゲットならば現状の項目 セットではどうしようもないので新規項目を追加する必要がある、というように状況に応じ た決定をサポートしてくれるでしょう。 具体的に特定の 𝜃𝑝 の値における項目情報量を出したい場合には iteminfo() 関数を使い ます。 54

8.9 項目情報量・テスト情報量 項目情報量 1 # はじめに項目パラメータの情報だけ抽出する必要がある 2 # 1番目の項目を抽出する 3 item1 <- extract.item(result_grm, item = 1) 4 # 引数Thetaには，どの値での特性値を出したいかを指定する 5 iteminfo(item1, Theta = c(-4, -2, 0, 2, 4)) 1 [1] 0.10354648 0.10516020 0.10114364 0.08337954 0.04723855 続いてテスト情報関数です。こちらは引数を type = "info" とすることで表示できます。 テスト情報関数 1 # 引数typeは同じ 2 # type = "infotrace"とすると前述の項目情報関数になるので注意 3 plot(result_2plm, type = "info") Test Information 2.5 I (q ) 2.0 1.5 1.0 0.5 0.0 -6 -4 -2 0 2 4 6 q 図8.34 テスト情報関数 尺度（10 項目）全体を見ても，項目困難度がまんべんなく存在している-3 から 0 付近の情報 量が高い一方で，すべての項目の困難度が 0 以上なので 𝜃𝑝 が正のエリアでは情報量がさほ ど高くないことがわかります。 プロットする際に引数 type を"infoSE" にすると，標準誤差を合わせて出してくれるよう になります。 55

8.9 項目情報量・テスト情報量 テスト情報関数に標準誤差を添えて 1 # 引数typeは同じ 2 plot(result_2plm, type = "infoSE") Test Information and Standard Errors 12 2.5 10 2.0 1.5 6 SE(q ) I (q ) 8 1.0 4 0.5 2 0.0 -6 -4 -2 0 2 4 6 q 図8.35 テスト情報関数 最も情報量が高い 𝜃𝑝 = −1.4 付近では，標準誤差はだいたい √1 2.82 ≃ 0.60 くらいになって います。仮に 𝜃𝑝 の真値が −1.4 の人が大量にいた場合，その人達の推定値は平均で 0.6 くら いは上下するだろう，ということです。裏を返せばその程度の精度でしか推定できていない ならば，例えば 𝜃𝑝̂ = −1.3 の人と 𝜃𝑝̂ = −1.5 の人がいたとしても，前者のほうが真の特性値 が高いとはなかなか言い切れないわけですね。 具体的に特定の 𝜃𝑝 の値におけるテスト情報量を出したい場合には testinfo() 関数を使い ます。 テスト情報量 1 # 引数Thetaには，どの値での特性値を出したいかを指定する 2 testinfo(result_grm, Theta = c(-4, -2, 0, 2, 4)) 1 [1] 2.8917688 3.9684301 3.8450407 2.5423962 0.5038195 1 # 標準誤差を出すならばこの逆数のルートを取れば良い 2 1 / sqrt(testinfo(result_2plm, Theta = c(-4, -2, 0, 2, 4))) 1 [1] 1.2460379 0.6335903 0.7433266 1.8358362 5.2840351 𝜃𝑝 = 4 の人では，テスト情報量が 0.5 程度しかなく，標準誤差が 5.28 と非常に大きくなっ 56

8.10 （おまけ）多母集団モデルと特異項目機能 ています。つまりこの 10 項目で推定を行っても，𝜃 の真値が 4 の人に対する推定値は平均で ±5 くらいには変動するわけです。そう考えると，この 10 項目ではどうあがいても 𝜃𝑝 が高 い層においてまともな推定値は得られないということが分かります。 8.10 （おまけ）多母集団モデルと特異項目機能 Chapter 7: SEM の中で，多母集団同時分析という枠組みを紹介しました。当然ながら 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. 特性値 𝜃 の分布はグループ間で共通として，グループごとに推定される項目パラメー タ 𝛼, 𝛽 の差異（=DIF）を見る 2. 項目パラメータ 𝛼, 𝛽 はグループ間で共通だとして，特性値 𝜃 の分布のグループごとの 差異を見る（SEM で言うところの強測定不変モデルと同じ） mirt パッケージで多母集団 IRT モデルを実行するためには，mirt() 関数の代わりに multipleGroup() 関数を使用します。multipleGroup() 関数では，mirt() 関数に引数に 加えて group という引数が登場します。lavaan パッケージのときにも多母集団モデルでは 引数 group を指定しましたが，その時とは指定の仕方が異なるのでご注意ください。 57

8.10 （おまけ）多母集団モデルと特異項目機能 多母集団 IRT モデル (用途 1) 1 # modelに1を与えておけば普通の一次元モデルになる 2 # 引数groupにはcharacter型かfactor型のベクトルを与える必要がある 3 result_group <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender)) 結果を見るための関数は概ね mirt() のときと同じように使えます。とりあえず推定された 項目パラメータを見てみましょう。 1 coef(result_group, simplify = TRUE, IRTpars = TRUE) 1 $`1` 2 $items 3 a b g u 4 Q1_1 0.447 -2.061 0 1 5 Q1_2 1.057 -1.751 0 1 6 Q1_3 0.971 -1.545 0 1 7 Q1_4 1.169 -1.271 0 1 8 Q1_5 1.073 -1.385 0 1 9 Q1_6 1.249 -1.455 0 1 10 Q1_7 1.281 -1.086 0 1 11 Q1_8 1.273 -1.111 0 1 12 Q1_9 1.045 -0.930 0 1 13 Q1_10 1.156 0.203 0 1 14 15 $means 16 F1 17 0 18 19 $cov 20 21 F1 F1 1 22 23 24 $`2` 25 $items 26 58 a b g u 27 Q1_1 0.294 -4.870 0 1 28 Q1_2 0.853 -3.029 0 1 29 Q1_3 0.962 -2.178 0 1 30 Q1_4 0.858 -2.118 0 1 31 Q1_5 0.871 -2.184 0 1 32 Q1_6 1.234 -1.613 0 1 33 Q1_7 1.449 -1.268 0 1 34 Q1_8 1.228 -1.369 0 1

8.10 （おまけ）多母集団モデルと特異項目機能 35 Q1_9 36 Q1_10 1.322 -0.127 0 1 1.574 -1.041 0 1 37 38 $means 39 F1 40 0 41 42 $cov 43 44 F1 F1 1 結果を見ると，$`1`と $`2`という 2 つに分かれて，それぞれ項目パラメータの推定値が出て います。これは，引数 group で指定したグループの値ごとに項目パラメータをそれぞれ推定 した結果です。 （2 つのグループで 𝜃 の分布が同じという前提のもとで）もしも推定された項 目パラメータが大きく異なる場合には，DIF が発生していると考えられるわけです。また，各 グループにおいて $means および $cov はそれぞれ 0，1 となっています。multipleGroup() 関数は，デフォルトではこのように「各グループの 𝜃 の母集団分布がそれぞれ標準正規分布 である」という設定のもとで項目パラメータを推定します。 特性値 𝜃 の分布のグループごとの差異を見たい場合には，等値制約を表す引数を与える必要 があります。lavaan での多母集団モデルを実行したときには，因子負荷の等値制約などを 引数 group.equal で指定しました。multipleGroup() では，invariance という引数が用 意されています。この引数には，以下に示す値を入れることが出来ます。free_means およ び free_vars という値があるということは，これらを設定しない限り，全てのグループの𝜃 は標準正規分布に固定されてしまうということなのでご注意ください。 指定 制約されるもの free_means グループ 1 以外の 𝜃 の平均値を自由推定する free_vars グループ 1 以外の 𝜃 の分散を自由推定する slopes 項目の識別力パラメータをグループ間で同じとする intercepts 項目の切片パラメータをグループ間で同じとする 多母集団 IRT モデル (用途 2) 1 # 引数invarianceで等値制約を指定する 2 # "free_means" と "free_vars"を入れないとN(0,1)のままなので注意 3 result_group2 <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender), invariance = c("free_means", "free_vars", "slopes", "intercepts")) 改めて推定されたパラメータを見てみましょう。 1 59 coef(result_group2, simplify = TRUE, IRTpars = TRUE)

8.10 （おまけ）多母集団モデルと特異項目機能 1 $`1` 2 $items 3 a b g u 4 Q1_1 0.445 -2.584 0 1 5 Q1_2 1.072 -1.978 0 1 6 Q1_3 1.067 -1.555 0 1 7 Q1_4 1.041 -1.414 0 1 8 Q1_5 1.038 -1.482 0 1 9 Q1_6 1.133 -1.396 0 1 10 Q1_7 1.294 -1.003 0 1 11 Q1_8 1.162 -1.078 0 1 12 Q1_9 1.272 -0.792 0 1 13 Q1_10 1.208 0.242 0 1 14 15 $means 16 F1 17 0 18 19 $cov 20 21 F1 F1 1 22 23 24 $`2` 25 $items 26 a b g u 27 Q1_1 0.445 -2.584 0 1 28 Q1_2 1.072 -1.978 0 1 29 Q1_3 1.067 -1.555 0 1 30 Q1_4 1.041 -1.414 0 1 31 Q1_5 1.038 -1.482 0 1 32 Q1_6 1.133 -1.396 0 1 33 Q1_7 1.294 -1.003 0 1 34 Q1_8 1.162 -1.078 0 1 35 Q1_9 1.272 -0.792 0 1 36 Q1_10 1.208 0.242 0 1 37 38 $means 39 F1 40 0.4 41 42 $cov 43 F1 44 F1 1.013 確かに項目パラメータの推定値がグループ間で同じになり，$`2`の $means および $cov の 値が変わりました。この結果は（2 つのグループで DIF が発生している項目が無いという前 60

8.10 （おまけ）多母集団モデルと特異項目機能 提のもとで）グループ 2 のほうが 𝜃 の平均値が高いことを示しているわけです。 8.10.1 DIF の評価 ここからは，DIF が発生しているかを評価するいくつかの方法を紹介します。とりあえず二 値 IRT モデルについての方法を紹介していきますが，多くの方法はすでに多値モデルに対す る拡張も提案されているはずなので，興味があれば探してみてください。 ノンパラメトリックな方法 DIF の検出法は，大きく分けると IRT モデルに基づく，つまり項目パラメータの推定値に よって行われるパラメトリックな方法と，IRT を使用しないノンパラメトリックな方法に分 けられます。ノンパラメトリックな方法として最も有名なのが，2 × 2 クロス表に対する 𝜒2 検定を拡張した Mantel-Haenszel 検定 (マンテル・ヘンツェル検定: Mantel and Haenszel, 1959) に基づく方法です。手法自体は非常に古いものですが，その使いやすさから近年でも DIF 検出の方法として非常に多く用いられているようです (Berrío et al., 2020)。 MH 検定では，まず各集団を何らかの変数の値によって 𝐾 層に層化します。DIF 検出の場合 には，IRT で推定された特性値 𝜃 や合計点などによって層化するのが一般的です。表8.5は， そうして層化されたうちの第 𝑘 層について，ある二値項目に対する回答をグループごとに集 計したものです。MH 検定では，表8.5のような 2 × 2 クロス表が 𝐾 個あるときに，それら を全部ひっくるめて独立性の検定を行う手法と言えます。 表8.5 第 𝑘 層 (𝑘 = 1, 2, ⋯ , 𝐾) におけるクロス表 項目 𝑖 当てはまらない 当てはまる 計 グループ A 𝑥𝐴0 (𝑘) 𝑥𝐴1 (𝑘) 𝑛𝐴 (𝑘) グループ B 𝑥𝐵0 (𝑘) 𝑥𝐵1 (𝑘) 𝑛𝐵 計 𝑛0 (𝑘) 𝑛1 (𝑘) 𝑁 (𝑘) (𝑘) まずは一つの層から考えることで，具体的な方法を紹介していきましょう。例えば 𝜃 が最も 低い人達が集められた第 1 層において，グループ A の人たちのほうが「当てはまる」と回答 した割合が高かったとします。このように，特定の層に絞ったときに回答傾向に違いが見ら れるというのも立派な DIF です。そして DIF が発生している場合には，2 つのグループに おける回答の割合つまりオッズ比 (𝑘) OR = (𝑘) (𝑘) (𝑘) (𝑘) 𝑥𝐴1 /𝑛𝐴 (8.57) 𝑥𝐵1 /𝑛𝐵 が 1 ではない，という見方ができます。MH 検定では，このオッズ比を全ての層について統 合したもの OR (MH) が 1 であるかどうかを検定します。 61 = 𝐾 (𝑘) (𝑘) 𝐾 (𝑘) (𝑘) ∑𝑘=1 𝑥𝐴1 /𝑛𝐴 ∑𝑘=1 𝑥𝐵1 /𝑛𝐵 (8.58)

8.10 （おまけ）多母集団モデルと特異項目機能 MH 検定の結果が有意であれば，データ全体として一方のグループのほうが「当てはまる」 と回答した割合が高かったと言えます。ICC 的に言えば図8.36のような状況であると言える わけです。 P(ypi = 1) Group A Group B qp 図8.36 MH 検定で検出可能な（均一）DIF 一方で，MH 検定では検出できない DIF も存在します。図8.37に示されている DIF では，𝜃𝑝 の低い層ではグループ A のほうが「当てはまる」割合が高い一方で，𝜃𝑝 の高い層では逆転が 起こっています。OR (𝑘) (MH) を一つ一つ見れば 1 ではないわけですが，これを統合した OR は 1 になってしまうわけです。MH 検定による方法は，図8.36に示した均一 DIF は検出可能 な一方で，図8.37のような不均一 DIF は検出できない方法であることに気をつけましょう。 P(ypi = 1) Group A Group B qp 図8.37 不均一 DIF 項目パラメータの推定値を比べる方法 IRT に基づくパラメトリックな方法として，まずは項目パラメータの推定値の差を直接検定 する方法を紹介します。Wald 検定は，パラメータの推定値に対して用いられる検定の一種 です。一般的に用いられる一変量 Wald 検定では，パラメータ 𝜉 の推定値 𝜉 ̂ に対して 2 𝜉 ̂ − 𝜉0 ) 𝑊 =( 𝜎𝜉 ̂ (8.59) という検定統計量を考えます。ただし 𝜉0 は帰無仮説における値，𝜎𝜉 ̂ は 𝜉 ̂ の標準誤差です。 検定統計量 𝑊 は自由度 1 の 𝜒2 分布に従うことが知られているため，これを用いて「𝜉 ̂ が 𝜉0 である」という帰無仮説を検定します。 Lord (1980) はこの Wald 検定によって DIF を検出する方法を提案しました。DIF が発生 している場合には，グループごとに推定された項目パラメータに差があると考えられるため， 62

8.10 （おまけ）多母集団モデルと特異項目機能 (8.59) 式の中身を「2 つのグループでの推定値」に置き換えることにより，例えば項目 𝑖 の 困難度 𝛽𝑖 であれば (𝐴) (𝐵) 2 𝛽 − 𝛽𝑖 𝑊 =( 𝑖 ) 𝜎𝛽(𝐴) −𝛽(𝐵) 𝑖 =( (𝐴) 𝛽𝑖 𝜎𝛽(𝐴) 𝑖 𝑖 (𝐵) (8.60) 2 − 𝛽𝑖 ) + 𝜎𝛽(𝐵) 𝑖 という形で，同じように Wald 検定を実行することができます。 Wald 検定で用いる Wald 統計量 𝑊 は「推定値の差」を「標準誤差」で割った形をしている ため，パラメータの数が複数であっても同じように定義することができます。例えば項目 𝑖 の識別力と困難度を並べたベクトル 𝛏𝑖 = [𝛼𝑖 𝛽𝑖 ] に対して 𝑊 = 𝛏𝑖 𝚺−1 𝛏⊤ 𝑖 (8.61) という統計量（𝚺 は標準誤差行列）を考えてやれば，これが自由度 2（＝パラメータの数）の 𝜒2 分布に従うため，同じようにして検定ができるのです。 項目特性曲線を比べる方法 IRT では，推定された項目パラメータをもとに項目特性曲線（ICC）を描くことが出来まし た。ICC は，いわば「得点の期待値」を表したグラフです。したがって「グループ間で ICC のずれが大きい」場合には DIF が発生していると考えることが出来そうです (Raju, 1988)。 ICC のずれとは，図8.38でオレンジ色に塗りつぶされた部分の面積のことです。これが大き い場合には，確かに DIF が発生していると言えそうです。 P(ypi = 1) Group A Group B qp 図8.38 グループ間での ICC のズレ グループ A における ICC（正答確率）は，2 パラメータロジスティックモデルならば (𝐴) 𝑃𝑖 (𝜃) = 𝑃 (𝑦𝑝𝑖 = 1|𝜃) = 1 1+ (𝐴) exp[𝛼𝑖 (𝜃 (𝐴) − 𝛽𝑖 )] (8.62) と表せるため，ICC のズレの面積に基づく DIF の指標として，例えば ∞ DIF𝑖 = ∫ (𝐴) ∣𝑃𝑖 (𝐵) (𝜃) − 𝑃𝑖 (𝜃)∣ 𝑓(𝜃)𝑑𝜃 (8.63) −∞ というものを考えることができるわけです (e.g., Chalmers, 2018; Raju et al., 1995)。ただ し 𝑓(𝜃) は 𝜃 の母集団分布（普通は標準正規分布）を指しています。 63

8.10 （おまけ）多母集団モデルと特異項目機能 R でやってみる DIF の検出法は相当たくさんあります（レビューとしてKleinman and Teresi, 2016; Berrío et al., 2020など）。ということで検出法の紹介はこのあたりでおしまいにして，mirt パッケ ージに実装されている DIF 検出法を試してみましょう*22 。これまでに紹介した方法と異な るのですが，DIF() 関数を使うと，モデル比較および尤度比検定の観点から DIF の有無を評 価してくれます。DIF() 関数のデフォルトの挙動では，第一引数で与えたモデルをベースラ インに特定の1項目のみに等値制約をかけたときの変化を計算します。ここで等値制約がか かるのは引数 which.par で指定したパラメータです。したがって，2 パラメータモデルの場 合は which.par = c("a1", "d") と指定することで「その項目の全ての項目パラメータの DIF を丸ごと評価する」ことを意味します。（あまり無いケースですが）例えば識別力には 興味がなく，困難度（切片）の DIF だけに関心がある場合には which.par = "d" としたら 良いわけです。 DIF の検出 1 DIF(result_group, which.par = c("a1", "d")) 1 converged AIC SABIC HQ BIC X2 df p 2 Q1_1 TRUE -24.578 -19.339 -20.363 -12.985 28.578 2 0 3 Q1_2 TRUE -32.846 -27.608 -28.632 -21.253 36.846 2 0 4 Q1_3 TRUE -18.162 -12.924 -13.948 -6.569 22.162 2 0 5 Q1_4 TRUE -13.108 -7.870 -8.894 -1.515 17.108 2 0 6 Q1_5 TRUE -13.712 -8.473 -9.497 -2.119 17.712 2 0 7 Q1_6 TRUE 2.217 7.456 6.432 13.810 1.783 2 0.41 8 Q1_7 TRUE -4.187 1.052 0.028 7.406 8.187 2 0.017 -1.149 4.090 3.066 10.444 5.149 9 Q1_8 TRUE 10 Q1_9 TRUE -15.366 -10.128 -11.151 11 Q1_10 TRUE -9.558 -4.320 -5.344 -3.773 19.366 2.035 13.558 2 0.076 2 0 2 0.001 結果を見てみると，SEM の多母集団モデルのときにも登場した AIC や BIC などの情報量規 準および X2（𝜒2 ）統計量とこれに関連した自由度・𝑝 値が表示されています。例えば Q1_1 行 の AIC の値-24.578 というのは，「Q1_1 のみグループごとに項目パラメータが同じ＝ DIF が無いという制約をかけることで，AIC が 24.578 悪化した」ということを意味します。し たがって，AICのみで判断するならばQ1_1 には DIF が生じている，と言えるわけです。同 様に，𝜒2 検定の結果も「その項目には DIF が発生していないとすることでモデル適合度が どの程度悪化するか」という観点に基づく尤度比検定が行われています*23 。 実際に DIF がありそうと判断された項目について ICC を見るには，引数 plotdif = TRUE を与えてください。 64 *22 MH 検定は mirt パッケージには用意されていないようです。R では difR::difMH() などの関数によって MH 検定を行うことができます。 *23 項目の数だけ検定を繰り返しているので，実はここで表示されている結果には多重検定の問題が生じていま す。これを解消するためには，引数 p.adjust を指定すると良いです。例えば p.adjust = "holm" とする と Holm の方法によって調整された 𝑝 値を表示してくれるようになります。

8.10 （おまけ）多母集団モデルと特異項目機能 DIF がありそうな項目の ICC 1 # デフォルトではSABICがマイナスになっている項目のプロットを表示する 2 DIF(result_group, which.par = c("a1", "d"), plotdif = TRUE) Item Probability Functions -6 -4 -2 0 2 4 6 Q1_5 Q1_9 Q1_10 1.0 0.8 0.6 0.4 P (q ) 0.2 0.0 Q1_1 Q1_2 Q1_3 cat1:1 cat1:2 Q1_4 1.0 0.8 0.6 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 q 1 converged AIC SABIC HQ BIC X2 df p 2 Q1_1 TRUE -24.578 -19.339 -20.363 -12.985 28.578 2 0 3 Q1_2 TRUE -32.846 -27.608 -28.632 -21.253 36.846 2 0 4 Q1_3 TRUE -18.162 -12.924 -13.948 -6.569 22.162 2 0 5 Q1_4 TRUE -13.108 -7.870 -8.894 -1.515 17.108 2 0 6 Q1_5 TRUE -13.712 -8.473 -9.497 -2.119 17.712 2 0 7 Q1_6 TRUE 2.217 7.456 6.432 13.810 1.783 2 0.41 8 Q1_7 TRUE -4.187 1.052 0.028 7.406 8.187 2 0.017 9 Q1_8 TRUE -1.149 4.090 3.066 10.444 5.149 2 0.076 10 Q1_9 TRUE -15.366 -10.128 -11.151 11 Q1_10 TRUE -9.558 -4.320 -5.344 -3.773 19.366 2.035 13.558 2 0 2 0.001 確かに ICC にはそこそこのずれがありそうです。 「項 目 1 つ 1 つ に つ い て DIF の 有 無 を チ ェ ッ ク す る」 と い う 意 味 で は， 全 く 逆 方 向 か ら の ア プ ロ ー チ を 取 る こ と も 可 能 で す。 パ ラ メ ー タ 推 定 時 に invariance = c("free_means","free_vars","slopes","intercepts") を与えると，全ての項目パラ メータが同じ，すなわち全ての項目に DIF が無いことを仮定した上での推定を行うことが 出来ました。このときの結果（result_group2）をベースラインとして考えると，DIF() 関 数には特定の1項目のみから等値制約を外したときの変化を計算してもらいたいところです。 これを実現するためには，DIF() 関数の引数 scheme を指定する必要があります。 65

8.10 （おまけ）多母集団モデルと特異項目機能 DIF の検出 (2) 1 # 等値制約を"drop"して検証してほしい 2 DIF(result_group2, which.par = c("a1", "d"), scheme = "drop") 1 converged AIC SABIC HQ BIC X2 df p 2 Q1_1 TRUE -10.691 -5.452 -6.476 0.902 14.691 2 0.001 3 Q1_2 TRUE -8.581 -3.343 -4.367 3.012 12.581 2 0.002 4 Q1_3 TRUE -0.172 5.066 4.042 11.420 4.172 2 0.124 5 Q1_4 TRUE 1.756 6.995 5.971 13.349 2.244 2 0.326 6 Q1_5 TRUE 2.772 8.010 6.986 14.365 1.228 2 0.541 7 Q1_6 TRUE -2.267 2.972 1.948 9.326 6.267 2 0.044 8 Q1_7 TRUE 1.468 6.706 5.683 13.061 2.532 2 0.282 9 Q1_8 TRUE 0.590 5.829 4.805 12.183 3.41 2 0.182 10 Q1_9 TRUE -1.618 3.620 2.596 9.975 5.618 2 11 Q1_10 TRUE 3.055 8.293 7.270 14.648 0.945 2 0.623 0.06 結果の見方は先程と同じです。情報量規準の値がマイナスになっている・検定が有意になっ ている（𝑝 < .05）項目は DIF の疑いあり，と判断することができます。 続いては Lord の Wald 検定の結果を出してみます。ただ，Wald 検定を行うためには「推定 値の標準誤差行列」が必要になるため，まずは標準誤差を推定するように引数 SE = TRUE を 付け加えてやり直す必要があります。 標準誤差を推定してほしい 1 # SE = TRUEにするだけ 2 result_group <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender), SE = TRUE) 推定された標準誤差行列は一応 result_group@vcov で見ることができますが，見たところ でよくわからないと思うのでそっとしておいてあげてください。 無事標準誤差が推定できたので，これを用いて Wald 検定に望みます。Wald 検定は，DIF() 関数で引数 Wald = TRUE を与えることで実行可能です。 Wald 検定 1 DIF(result_group, which.par = c("a1", "d"), Wald = TRUE) 1 66 W df p 2 Q1_1 28.652 2 0 3 Q1_2 36.341 2 0 4 Q1_3 22.179 2 0 5 Q1_4 16.473 2 0 6 Q1_5 17.345 2 0

8.10 （おまけ）多母集団モデルと特異項目機能 7 Q1_6 1.786 2 0.409 8 Q1_7 8.373 2 0.015 9 Q1_8 5.129 2 0.077 10 Q1_9 19.703 11 Q1_10 13.783 2 0 2 0.001 結果の見方はこれまでと基本的に同じです。つまり検定結果が有意になった（𝑝 < .05）項目 は DIF の疑いあり，ということになります。 最後は ICC のズレに基づく評価方法です。こちらはChalmers (2018) による方法が DRF() という関数に実装されています。DRF() 関数は，デフォルトでは項目特性曲線（ICC）の代 わりにテスト特性曲線（TCC）の差について評価を行います*24 。項目ごとの DIF を見た い場合には，引数 DIF = TRUE を与える必要があるのでご注意ください。また引数 plot = TRUE を与えると，「ICC の差のプロット」を出してくれるようになります。 ICC のズレをチェック 1 DRF(result_group, DIF = TRUE, plot = TRUE) Signed DIF -6 -4 -2 0 2 4 6 0.3 Q1_9 Q1_10 Q1_5 Q1_6 0.2 0.1 0.0 Q1_7 Q1_8 0.3 sDIFq 0.2 0.1 0.0 0.3 Q1_1 Q1_2 Q1_3 Q1_4 0.2 0.1 0.0 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 q 結果を見ると，3 種類の DIF が表示されています。 sDIF 符合つき（signed）のズレの面積（単純に差を積分したもの）を表しています。図8.38の ようにズレが上側と下側にある場合，これらが打ち消し合って sDIF の値は小さくな ります。つまり sDIF は「総合的に見てどちらのグループのほうが平均点が高いか」を *24 67 この「テスト全体でのグループ間の挙動の違い」は DIF と同じように特異テスト機能 (Differentian Test Functioning [DTF]) と呼ばれたりします。

8.10 （おまけ）多母集団モデルと特異項目機能 表していると言えるでしょう。 uDIF 符号なし（unsigned）のズレの面積（差の絶対値を積分したもの: 8.63式）を表して います。したがって，sDIF よりも大きな値が表示されているはずです。 dDIF (8.63) 式とは別のやり方で符号なしのズレを求めるために，絶対値の代わりに差を二 乗することを考えた DIF です。二乗した分を戻すため，最後に全体にルートをかけま す。 「二乗平均のルートをとる」という考え方は，いわゆる逸脱度（deviance）的な考 え方なので dDIF と名付けられているのだと思います。多分。uDIF と比べると，ズ レの面積が同じだとしても，大きくズレている箇所があるほど，二乗されている分だ け dDIF のほうが大きな値になると思われます。 いずれにしても，この (s|u|d)DIF の値が大きい項目は DIF の疑いあり，と言えるわけで す。ICC のズレに基づく検定を行うためには「ICC のズレの標準誤差」のようなものが必 要になりますが，これはパラメータの標準誤差だけでは対応できないようです。そのため， DRF() 関数では代わりにブートストラップ標準誤差を用いて検定を行います。これを行うた めには，引数 draws に適当な数を指定してください。 ブートストラップを使って検定 1 # drawsの値は大きい方が正確になるが，その分時間がかかる 2 DRF(result_group, DIF = TRUE, draws = 500) 1 $sDIF 2 sDIF CI_2.5 CI_97.5 X2 df p 3 1 0.096 0.057 0.133 27.274 1 0 4 2 0.086 0.055 0.117 31.302 1 0 5 3 0.077 0.047 0.109 23.465 1 0 6 4 0.066 0.031 0.103 13.69 1 0 7 5 0.069 0.035 0.102 16.354 1 0 8 6 0.022 -0.010 0.058 1.677 1 0.195 9 7 0.050 0.009 0.086 7.019 1 0.008 10 8 0.041 0.004 0.076 5.075 1 0.024 11 9 0.070 0.029 0.110 12.506 1 0 12 10 0.078 0.038 0.119 1 0 13.62 13 14 $uDIF 15 68 uDIF CI_2.5 CI_97.5 X2 df p 16 1 0.096 0.059 0.133 27.701 2 0 17 2 0.086 0.055 0.117 26.172 2 0 18 3 0.077 0.048 0.110 18.848 1 0 19 4 0.067 0.036 0.105 14.664 2 0.001 20 5 0.069 0.038 0.102 13.637 2 0.001 21 6 0.022 0.004 0.060 1.972 1 22 7 0.050 0.020 0.089 7.249 2 0.027 23 8 0.041 0.011 0.077 5.324 2 0.16 0.07

8.10 （おまけ）多母集団モデルと特異項目機能 24 9 0.079 0.048 0.121 20.382 2 25 10 0.078 0.044 0.122 14.286 2 0.001 0 26 27 $dDIF 28 dDIF CI_2.5 CI_97.5 29 1 0.105 0.067 0.149 30 2 0.112 0.069 0.162 31 3 0.087 0.054 0.130 32 4 0.094 0.048 0.146 33 5 0.089 0.045 0.139 34 6 0.027 0.005 0.080 35 7 0.054 0.021 0.101 36 8 0.048 0.013 0.096 37 9 0.084 0.050 0.130 38 10 0.084 0.050 0.135 1 # plot = TRUEをつけると検定はせずに図だけ表示する 2 DRF(result_group, DIF = TRUE, draws = 500, plot = TRUE) Signed DIF -6 -4 -2 0 2 4 6 Q1_9 Q1_10 Q1_5 Q1_6 0.4 0.2 0.0 Q1_7 Q1_8 sDIFq 0.4 0.2 0.0 Q1_1 Q1_2 Q1_3 Q1_4 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 q 結果の見方はやはり同じです。つまり検定結果が有意になった（𝑝 < .05）項目は DIF の疑 いあり，ということになります。 DIF 検出の実際 ここまで DIF を検出するいくつかの方法を紹介しました。色々な観点からの評価がある，と いう意味では，やはり複数の指標を組み合わせて総合的に判断するのが良いのだろうと思い ます。 69

8.10 （おまけ）多母集団モデルと特異項目機能 それはそれとして，DIF の有無を評価しようとする場合，そもそも関心のあるグループの間 で特性値 𝜃 の分布に差が無いということはあまり考えられません（無かったら非常にラクな のですが…） 。したがって，そもそも「𝜃 の分布が同じだとして DIF を検出する」というスキ ームは現実的に結構厳しいのです。この問題を克服するためには， 「グループ間で挙動が変わ らない＝ DIF が無いことが事前にわかっている項目」であるアンカー項目 (anchor items) を加えておくことが考えられます。アンカー項目の項目パラメータを基準とすることで，𝜃 の分布がグループで異なる場合であっても特定の項目に DIF が生じているかを評価するこ とが出来るわけです。具体的にアンカー項目をどのように決定したら良いかはまた難しい話 (レビューとして Kopf et al., 2015) なので，ここではとりあえずアンカー項目がある場合の やり方だけ紹介しておきます。 multipleGroup() の引数 invariance には，実は項目の名前（データの列名）を与えるこ とも出来ます。 アンカー項目を用いた多母集団モデル 1 # Q1_1からQ1_3がアンカー項目だったとします 2 result_anchor <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender), invariance = c("free_means", "free_vars", colnames(dat_bin)[1:3])) あとの流れはほぼ同じなので，紹介はここまでとします。 1 # Q1_1からQ1_3のパラメータだけグループ間で同じ 2 coef(result_anchor, simplify = TRUE, IRTpars = TRUE) 1 $`1` 2 $items 3 a 4 Q1_1 0.444 -2.359 0 1 5 Q1_2 1.014 -1.831 0 1 6 Q1_3 0.962 -1.456 0 1 7 Q1_4 1.165 -1.274 0 1 8 Q1_5 1.061 -1.396 0 1 9 Q1_6 1.255 -1.450 0 1 10 Q1_7 1.291 -1.082 0 1 11 Q1_8 1.279 -1.108 0 1 12 Q1_9 1.056 -0.924 0 1 13 Q1_10 1.167 14 15 $means 16 F1 17 0 18 19 70 b g u $cov 0.202 0 1

8.10 （おまけ）多母集団モデルと特異項目機能 20 21 F1 F1 1 22 23 24 $`2` 25 $items 26 a b g u 27 Q1_1 0.444 -2.359 0 1 28 Q1_2 1.014 -1.831 0 1 29 Q1_3 0.962 -1.456 0 1 30 Q1_4 0.949 -1.176 0 1 31 Q1_5 0.957 -1.246 0 1 32 Q1_6 1.345 -0.734 0 1 33 Q1_7 1.579 -0.418 0 1 34 Q1_8 1.353 -0.502 0 1 35 Q1_9 1.736 -0.204 0 1 36 Q1_10 1.441 0.625 0 1 37 38 $means 39 F1 40 0.741 41 42 $cov 43 F1 44 F1 0.827 1 # もともと等値制約が置かれているQ1_1からQ1_3については結果が出ない 2 DIF(result_anchor, which.par = c("a1", "d")) 1 converged AIC SABIC HQ BIC 2 Q1_1 TRUE 0.000 0.000 0.000 0.000 X2 df 0 0 NaN p 3 Q1_2 TRUE 0.000 0.000 0.000 0.000 0 0 NaN 4 Q1_3 TRUE 0.000 0.000 0.000 0.000 0 0 NaN 5 Q1_4 TRUE -1.403 3.835 2.812 10.190 5.403 2 0.067 6 Q1_5 TRUE 0.571 5.809 4.785 12.164 3.429 2 0.18 7 Q1_6 TRUE -18.971 -13.733 -14.757 -7.378 22.971 2 0 8 Q1_7 TRUE -15.059 -9.821 -10.844 -3.466 19.059 2 0 9 Q1_8 TRUE -14.729 -9.490 -10.514 -3.136 18.729 2 0 10 Q1_9 TRUE -14.360 -9.122 -10.146 -2.767 2 0 11 Q1_10 TRUE -10.458 -5.219 -6.243 18.36 1.135 14.458 2 0.001 結果を見て分かるように，アンカー項目を適切に設定しないと，その他の項目の DIF の評価 も適切に行うことができません。そういった意味では，DIF の検出というのは結構職人技だ ったりするのです。 71

8.11 いろいろな IRT モデル 8.11 いろいろな IRT モデル IRT の本質は，(8.1)式に集約されていると個人的には思っています。それは項目への反応 （データ）は回答者と項目の交互作用によって確率的に決定するという点です。そのような挙 動をとるモデルの形の一つとして，因子分析と同じようなモデル（正規累積・ロジスティッ ク）や多値反応に対するモデル（GRM・GPCM など）を紹介してきました。しかし世の中 にはもっと様々な IRT モデルが存在しています。全てを紹介するのは不可能ですが，ここで はいくつかのモデルのさわりだけを（個人的な関心に基づいて）紹介していきたいと思いま す*25 。 8.11.1 多値反応に対するその他の IRT モデル 多値反応に対する IRT モデルとして，8.6節では GRM や GPCM といったモデルを紹介し てきました。多値反応に対する IRT モデルには他にも様々なものが提案されており，近年注 目を集めている（気がする）のが IRTree (Böckenholt, 2017) と呼ばれるモデル群です。こ れらのモデルでは，回答決定のプロセスとして名前通り「ツリー構造」を仮定します。例え ば「とても嫌い (1)」 「嫌い (2)」 「どちらでもない (3)」 「好き (4)」 「とても好き (5)」という 5 件法の項目があったとしましょう。この項目への回答決定プロセスの考え方の一例としては， 図8.39のように 1. まず「どちらでもない (3)」を選択（態度保留）するかを二択で考える 2. 態度を表明する場合には「好き (4-5)」か「嫌い (1-2)」かを二択で考える 3.「とても (1,5)」の選択肢を選ぶほど好き or 嫌いかを二択で考える という感じで３つの二択を順番に行う，というものがあるでしょう。 図8.39 IRTree の一例 するとこの多値項目の背後には，実は３つの二値項目（query）が存在していると考えるこ *25 ここで紹介するモデルを始め，多くの（特殊な）IRT モデルは mirt パッケージに限らず R のパッケージ としては実装されていなかったりします。そういったモデルについては，自分で尤度関数をプログラムして optim() 関数による数値最適化を利用する，あるいは stan などを利用してマルコフ連鎖モンテカルロ法 （MCMC）によって推定する，といった方法を取る必要があります。stan を使うためにはまずベイズ統計学 を理解する必要があるのがちょっと大変ですが，慣れるとほぼなんでもできるようになります。 72

8.11 いろいろな IRT モデル とができ，それぞれが (𝑀) 1.「どちらでもない」を選ぶ傾向（中心傾向）：𝜃𝑝 2. が強い人かどうか (𝐴) 実際にその対象物が好きかどうか：𝜃𝑝 (𝐸) 3. 極端な選択肢を選ぶ傾向（極端傾向）：𝜃𝑝 が強い人かどうか という異なる個人特性によって動いていると仮定することができそうです（もちろんプロセ ス単位で完全に独立しているとは思いませんが，ある程度は）*26 。各 query に対する反応確 率は，普通の正規累積・ロジスティックモデルなどで表現できます。例えば中心傾向につい ての query だけで表現できる「どちらでもない」を選ぶ確率は (𝑀) 𝑃 (𝑦𝑝𝑖 = (𝑀) 3|𝜃𝑝 ) 𝛼𝑖 (𝑀) (𝜃𝑝 (𝑀) −𝛽𝑖 ) =∫ −∞ = (𝑀) 𝚽 [𝛼𝑖 (𝑀) (𝜃𝑝 − 𝑧2 1 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.64) (𝑀) 𝛽𝑖 )] などとなるわけです。IRTree モデルでは，項目パラメータにも上付き添字がついています。 したがって，項目ごとに「どちらでもない」および極端な選択肢の選ばれやすさが異なる，と いう設定がされています。 同様に， 「好き (4)」を選ぶ確率は， 「 『どちらでもない』を選ばない」かつ「好き」かつ「 『と ても』ではない」の積なので (𝑀) 𝑃 (𝑦𝑝𝑖 = 4|𝜃𝑝 (𝑀) (1 − 𝚽 [𝛼𝑖 (𝐴) (𝐸) , 𝜃𝑝 𝜃𝑝 ) = (𝑀) (𝜃𝑝 (𝑀) − 𝛽𝑖 (𝐴) )]) 𝚽 [𝛼𝑖 (𝐴) (𝜃𝑝 (𝐴) (𝐸) − 𝛽𝑖 )] (1 − 𝚽 [𝛼𝑖 (𝐸) (𝜃𝑝 (𝐸) − 𝛽𝑖 )]) (8.65) という形で表現できるわけです。 IRTree の目的あるいは利点は，「測定したい回答者の特性値」と「回答傾向」を分離して考 えられる点です。一般的に，リッカート尺度への回答には個人のクセのようなものが混ざっ ていることが知られています。そのせいで，「どちらでもない (3)」を選んだ人が「好きと嫌 いのちょうど中間」なのか「本当は『どちらかといえば好き』くらいなんだけど，その程度 では『好き』は選びづらい」のかは区別がつきにくいものです。IRTree モデルならば，上の 例で言えば中心傾向と極端傾向の強さを，個人特性 𝜃𝑝 の一種として考えることができます。 (𝐴) そのため，残った 𝜃𝑝 はより純粋な「測定したい回答者の特性値」として解釈できるかもし れないというわけです。当然ながら，IRTree モデルでは図8.39に示したような回答決定プロ セスの仮定がかなり重要になってきます。逆に言えば，複数のツリー構造間でモデル比較を することで，どのプロセスが最も当てはまりそうか，ということを評価したりも出来るかも しれません…。 8.11.2 強制選択形式に対する IRT モデル IRTree モデルはいわば，得られた回答データから回答傾向を統計的に分離しようと頑張って いるモデルの一種です。そもそもリッカート尺度では回答傾向の影響を逃れることができな *26 𝜃𝑝 の上付き添字の M, A, E はそれぞれ Middle/Mid-point, Attitude/Agree, Extreme などの頭文字で す。 73

8.11 いろいろな IRT モデル いわけですが，その対抗策としては他にも「そもそも中間選択肢をなくそう」という方向を 考えることができます。また，入社試験などの重要な局面では「ウソをついて本当の自分よ り良い人に見せよう」という思惑が働いて回答が歪められる可能性があったりします (フェ イキング：e.g., Jackson et al., 2000)。これを防ぐためには，例えば「友だちが多い」 「真面 目である」といった複数の短文を同時に提示し，その中から一つだけを選択させたり，ラン キングを付けさせるという回答方法が考えられます。これならば，前述の「中間選択肢を選 びがち」といった問題も防げますし，フェイキングの問題もかなり抑制できるでしょう。 複数の選択肢から一つの選択肢を選ばせる回答方法は，強制選択 (Forced Choice) 型項目な どと呼ばれます。また，特に選択肢数が２個の場合は一対比較 (Paired Comparison) と呼 ばれることもあります。そしてこれらの出題形式に対しては，やはり普通の IRT とは少し異 なるモデルが考えられています。最も代表的なモデルが，Thurstonian IRT (TIRT) モデル (Brown and Maydeu-Olivares, 2011) でしょう。TIRT モデルでは，各選択肢がそれぞれ独 立に「効用」を持っていると考えます。この「効用」とは，経済学などで出てくるものと同 じようなものだと思って OK です。あるいは単純に「各選択肢が自分に当てはまっていると 感じる程度」というイメージでも構いません。とにかく強制選択型項目への回答では，回答 者は効用を最大化する選択肢を選ぶのだ，と考えているわけです。ここでは，回答者 𝑝 が項 (𝐴) (𝐵) 目 𝑖 の２つの選択肢 𝐴, 𝐵 に対してそれぞれ効用 𝑢𝑝𝑖 , 𝑢𝑝𝑖 を持っているとしましょう。する と，回答者 𝑝 は効用を最大化する選択肢を選ぶため，観測されるデータとしては 𝑦𝑝𝑖 = { 𝐴 𝐵 (𝐴) (𝐵) if 𝑢𝑝𝑖 − 𝑢𝑝𝑖 ≥ 0 (𝐴) (𝐵) if 𝑢𝑝𝑖 − 𝑢𝑝𝑖 < 0 (8.66) と表わせます。 (𝐴) TIRT モデルでは，選択肢の効用 𝑢𝑝𝑖 を，よくある因子分析モデルの形 (𝐴) (𝐴) 𝑢𝑝𝑖 = 𝜏𝑖 (𝐴) (𝐴) + 𝑓 𝑝 𝑏𝑖 (𝐴) (𝐴) (8.67) 𝑝𝑖 (𝐴) と定義します。つまり，平均的な効用 𝜏𝑖 (𝐴) + 𝜀𝑝𝑖 , 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀(𝐴) ) (𝐴) と，その選択肢が反映している因子の得点 𝑓𝑖 と因子負荷 𝑏𝑝 の和によって表される効用の期待値に，そこに正規分布に従う誤差（独自 (𝐴) 因子得点）𝜀𝑝𝑖 が加わって効用の実現値が確率的に決定すると考えているわけです。すると， ２つの選択肢 𝐴, 𝐵 の中から 𝐴 を選ぶ確率は最終的に (𝐴) (𝐵) 𝑃 (𝑦𝑝𝑖 = 𝐴) = 𝑃 (𝑢𝑝𝑖 − 𝑢𝑝𝑖 ≥ 0) (𝐴𝐵) = 𝚽 [𝛼𝑖 (𝐴) (𝐴) 𝛽𝑝 + (𝑓𝑖 (8.68) (𝐵) (𝐵) 𝛽𝑝 )] − 𝑓𝑖 と表されます*27 。ただし表記を簡略化するため (𝐴𝐵) 𝛼𝑖 (𝐴) = (𝐴) 𝜏𝑖 − 𝜏 𝑖 , √𝜎𝜀(𝐴) + 𝜎𝜀(𝐵) 𝑝𝑖 𝑝𝑖 (𝐴) 𝛽𝑖 (𝐴) = 𝑏𝑖 + 𝜎𝜀(𝐵) √𝜎𝜀(𝐴) 𝑝𝑖 𝑝𝑖 , (𝐵) 𝛽𝑖 (𝐵) = 𝑏𝑖 + 𝜎𝜀(𝐵) √𝜎𝜀(𝐵) 𝑝𝑖 𝑝𝑖 という形で別の記号を用いて表現しています。つまり TIRT モデルは，一対比較の場合には， 各項目が２つの因子のみに負荷量を持っていると設定した CFA と実質的には同じことをし *27 74 ２つの選択肢の比較のときに，その背後にある効用を比較していて，しかも効用が確率的に変動する，とい う考え方をサーストンの比較判断の原則 (law of comparative judgment: Thurstone, 1927) と呼ぶため， このモデルは Thurstonian IRT と名付けられています。

8.11 いろいろな IRT モデル ています。そして一つの項目に３つ以上の選択肢がある場合は，対ごとに異なる項目を構成 しているとみなして分析を行います。例えば項目 𝑖 が 𝐴, 𝐵, 𝐶, 𝐷 の４択だとすると，この項 目は「A と B」 「A と C」 「A と D」 「B と C」 「B と D」 「C と D」という６つの二値項目が あると考えます。そして回答者 𝑝 が選択肢 𝐴 を選んだとすると， 「A と B」 「A と C」 「A と D」については（いずれも A を選んだという）結果が観測され，残りの「B と C」「B と D」 「C と D」については欠測扱いにします。そして観測された３項目についてはそれぞれ (𝐴𝐵) + (𝑓𝑖 (𝐴𝐶) + (𝑓𝑖 (𝐴𝐷) + (𝑓𝑖 𝑃 (𝑦𝑝𝑖 = 𝐴|𝐴, 𝐵) = 𝚽 [𝛼𝑖 𝑃 (𝑦𝑝𝑖 = 𝐴|𝐴, 𝐶) = 𝚽 [𝛼𝑖 𝑃 (𝑦𝑝𝑖 = 𝐴|𝐴, 𝐷) = 𝚽 [𝛼𝑖 (𝐴) (𝐴) 𝛽𝑝 − 𝑓𝑖 (𝐵) (𝐵) 𝛽𝑝 )] (𝐴) (𝐴) 𝛽𝑝 − 𝑓𝑖 (𝐴) (𝐴) 𝛽𝑝 − 𝑓𝑖 (𝐶) (𝐶) 𝛽𝑝 )] (8.69) (𝐷) (𝐷) 𝛽𝑝 )] という形でそれぞれ選択肢 A に関する部分が共通しているという設定のもとで定式化される わけです。 8.11.3 回答時間を利用する IRT モデル コンピュータによる測定 (Computer Based Testing [CBT]) では，項目反応以外の情報を得 ることが容易になりました。もっとも代表的な情報として，回答時間 (Response Time [RT]) が挙げられるでしょう。Chapter 4 では，回答時間を「適当な回答をしている人を検出する」 といった目的で使用しましたが，その一方で回答者の特性値推定に利用するという方向性も 相当に研究されています。そういうわけで，回答時間を利用する IRT モデルを挙げていくと キリがないのですが，ここではシンプルなモデルとして， 「当てはまる」 「当てはまらない」の 二択の項目に対するモデル (Ferrando and Lorenzo-Seva, 2007) を紹介します。 回答時間を利用する IRT モデルでは，そもそも回答者の特性値と回答時間にはどのような関 係があるかを考えるところから始まります。性格特性を尋ねる項目の場合には， 「特性の強さ が中途半端な人ほど回答時間が長い」という逆U字の関係になることが経験的に知られてい ます。例えば「友だちが多い」という項目に対して，友だちが多い／少ない自信がある人，つ まり「外向性」といった特性が極端な人は，きっとすぐに選択肢を選ぶことでしょう。これ に対して，『友だちはそこそこ居るけど』という人は，きっと「当てはまる」「当てはまらな い」のどちらを選ぶかに迷いが生じるはずです。したがって，同じ「当てはまる」を選んだ 人同士であっても，素早く「当てはまる」を選んだ人のほうがより強い特性を持っていると 考えるのが自然です。 ということで，この関係を尤度関数に組み込んで IRT モデルを構築していきます。と言って も考え方は簡単で，回答時間 𝑡𝑝𝑖 についても項目反応 𝑦𝑝𝑖 の(8.1)式と同じように 𝑡𝑝𝑖 = 𝑔(𝜔𝑝 , 𝛾𝑖 ) (8.70) と， 回答者の要因 𝜔𝑝 と項目の要因 𝛾𝑖 の交互作用として考えます。 Ferrando and Lorenzo-Seva (2007) では，これを具体的に ln(𝑡𝑝𝑖 ) = 𝜇 + 𝜔𝑝 + 𝛾𝑖 + 𝑏𝛿𝑝𝑖 + 𝜀𝑝𝑖 , 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀2 ) 𝛿𝑝𝑖 = √𝛼2𝑖 (𝜃𝑝 − 𝛽𝑖 )2 (8.71) と置きました。ここで 𝜇 は全体的な切片項，𝜔𝑝 は回答者 𝑝 の回答スピード，𝛾𝑖 は項目 𝑖 の 所要時間を表現したパラメータです。また，𝛿𝑝𝑖 の中身には (𝜃𝑝 − 𝛽𝑖 )2 というものが入って 75

8.11 いろいろな IRT モデル います。この項は，𝜃𝑝 = 𝛽𝑖 のとき，つまり反応確率が五分五分のときに最小値を取ります。 係数 𝑏 が負の値であれば，回答者の特性値 𝜃𝑝 が項目困難度 𝛽𝑖 に近いほど回答時間の期待値 が長くなることで，逆 U 字の関係を表現しているわけです。 回答時間などの追加情報を用いて特性値を推定する IRT モデルでは，うまく行けばより精度 の高い，あるいは妥当性の高い推定を行うことが出来ると期待されています。ちなみに，IRT に限らず社会調査などの分野では，このように回答内容以外の情報（位置情報や回答デバイ スの情報などを含む）をパラデータ (Kreuter, 2013) と呼んでいます。パラデータは，やは り適当な回答の検出に使われるだけでなく，よりノイズの少ないデータを収集するための調 査改善など様々な用途で利用されていたりします。 8.11.4 展開型 IRT モデル 前半で紹介した IRT モデルはすべて，𝜃 の値が高ければ高いほど高い点数をとる（i.e., 「当 てはまる」寄りの回答をする，正解する）という仮定がされていました。項目特性曲線にし ても， （𝛼𝑖 > 0 ならば）右上がり（単調増加）でした。ですが質問のタイプによっては「過ぎ たるはなお及ばざるが如し」的なものが存在しています。よくあるのは「態度」に関する項 目です。例えば「消費税を5%にするべきだ」という項目は，「消費税は完全に撤廃してほし い」人も「消費税はこのままで良い」人も「そう思わない」を選択するでしょう。また，心 理特性に関しても「カフェで友人と静かに会話をするのは楽しい（因子：外向性）」(Stark et al., 2006) という項目に対しては「そもそも他人と話したくない人」も「もっとワイワイ 騒ぎたい人」も「当てはまらない」を選ぶと思われます。つまりこうした項目における項目 特性曲線は，図8.40のように，あるところを境に単調減少に転じると考えられるわけです。 図8.40 個人適合度の考え方 ロジスティック・正規累積モデルでは，「項目 (𝛽𝑖 )」と「回答者 (𝜃𝑝 )」は共通の軸の上でマ ッピングされていると考えることが出来ます。すごく簡単にいえば「困難度が 𝛽𝑖 の項目は 𝜃𝑝 = 𝛽𝑖 の回答者にとって（易しすぎず難しすぎず）ちょうどよい難易度（50%）の項目で ある」ということです。これを先程の項目に置き換えて考えると， 「𝛽𝑖 と 𝜃𝑝 の距離が近いほ ど，その人によってちょうどよい程度の（理想の）意見・態度を表している」ということが できそうです。このような「理想点 (ideal point)」の考え方に基づく IRT モデルを展開型 76

8.11 いろいろな IRT モデル (unfolding) IRT モデルと呼びます*28 。 展開型 IRT モデルでは，反応確率 𝑃 (𝑦𝑝𝑖 = 1) が「𝛽𝑖 と 𝜃𝑝 の距離」によって規定されます。 したがって，最も基本的な一次元・二値モデルの場合は 1 𝑃 (𝑦𝑝𝑖 = 1) = exp [− (𝛽𝑖 − 𝛼𝑖 𝜃𝑝 )2 ] 2 (8.72) と表されます (Maydeu-olivares et al., 2006)。ここでも 2 つの項目パラメータの意味は概 ね通常の場合と同じです。𝛽𝑖 はその項目の「位置」を表しており，𝜃𝑝 がこの値に近い人ほ ど「当てはまる」と回答する確率が高くなることがわかります。𝛼𝑖 は項目識別力を表してい ますが，展開型モデルの場合は，正規分布の分散パラメータのように山の狭さを規定します （図8.41） 。ですが結局「𝜃𝑝 が理想点 𝛽𝑖 付近か否かを区別する強さ」という意味では，これま でに出てきた項目識別力と同じような役割を果たしているわけです。 P(ypi = 1) ai = 0.5 ai = 1 ai = 2 qp q p = bi 図8.41 識別力の役割（展開型） mirt では，itemtype="ideal" と指定することで展開型モデルを実行可能です。 1 result_ideal <- mirt(dat_bin, model = 1, itemtype = "ideal") 通常の IRT モデルと同じように coef() でパラメータの推定値を見ることもできるのです が，どうやら IRTpars = TRUE には対応していないようなのでご注意ください。 1 # IRTpars = TRUEとするとエラーが出る 2 coef(result_ideal, simplify = TRUE) 1 $items 2 77 a1 d 3 Q1_1 0.118 -0.716 4 Q1_2 0.214 -0.468 5 Q1_3 0.250 -0.571 6 Q1_4 0.287 -0.609 7 Q1_5 0.255 -0.600 8 Q1_6 0.618 -0.312 *28 これに対して通常の IRT モデルは dominance model と呼びます（日本語訳がわからない…）。

8.11 いろいろな IRT モデル 9 Q1_7 0.651 -0.455 10 Q1_8 0.621 -0.489 11 Q1_9 0.571 -0.653 12 Q1_10 0.667 -1.187 13 14 $means 15 F1 16 0 17 18 $cov 19 20 F1 F1 1 また，ICC のプロットも同じように実行可能です。 1 # 10番目の項目 2 plot(result_ideal, type = "trace", facet_items = FALSE) Item Probability Functions P1 1.0 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 P (q ) 0.8 0.6 0.4 0.2 0.0 -6 -4 -2 0 2 4 6 q 図8.42 項目特性曲線（展開型） 例えば 1 番目の（一番山が緩い）項目は，(𝛼𝑖 , 𝛽𝑖 ) = (0.118, −0.716) ということで，𝜃 = 0.716/0.118 ≃ 6.068 のあたりで最大値 𝑃 (𝑦𝑝𝑖 = 1) = 1 となります。そのため，𝜃 が 6 まで しか表示されていない図8.42では右上がりのように見えています。 多値の展開型 IRT モデル 多値の展開型 IRT モデルの代表的なものは，一般化段階展開型モデル (Generalized graded unfolding model［GGUM］: Roberts, 2019) です。このモデルは，GPCM をベースに展開型 IRT モデルを構成しています。 例えば 4 カテゴリの項目がある場合には，「全く当てはまらない（下）」「あまり当て 78

8.12 IRT の使い道 はまらない（下） 」 〜 「あまり当てはまらない（上） 」 「全く当てはまらない（上） 」とい う要領で 8 個のカテゴリがあると考えます。こうすることで，隣同士のカテゴリの比 較は今まで通りできるので GPCM 的な方法が適用できる，というわけです。そして 「全く当てはまらない（下）」と「全く当てはまらない（上）」の確率を合わせたもの を「全く当てはまらない」の選択確率として表します。このままだとパラメータ数が 多すぎる（4 カテゴリに対してカテゴリ困難度が 7 個必要になってしまう）ので，理 想点を軸にして各カテゴリ閾値間の距離は左右で等しいという仮定をおいた上で定式 化しています。 なお mirt では，itemtype = "ggum" と指定することで実行可能です。 8.12 IRT の使い道 最後に，IRT の考え方によって出来るようになることをいくつか紹介します*29 。 8.12.1 等化 本 Chapter の最初では，IRT を項目依存性と集団依存性をクリアするための理論だと紹介 しました。ですが実際には，IRT モデル（2 パラメータの GRM）はカテゴリカル因子分析と 同じなので 𝜃𝑝 はデータをとった集団ごとに標準化されています。つまり前回の最初に示し たような，異なる集団に対して実施した異なる尺度の比較はこのままではできません。もち ろん反対に，異なる尺度を用いて実施した異なる集団同士の比較もこのままではできません。 こうした比較を可能にするために，IRT では等化 (equating) またはリンキング (linking) という手続き (Kolen and Brennan, 2014) を行います*30 。等化の技術は，TOEIC のような 大規模テストでも実際に利用されています。そのおかげで，毎回異なるはずの難易度によら ず「TOEIC *** 点」というように点数を履歴書に書けるわけです。 等化の基本的な考え方を説明しましょう。ここからは IRT で推定されるパラメータが集団に よらない項目パラメータと項目によらない個人特性値だという仮定のもとで話を進めます。 【共通受験者による等化】 まず，個人特性値が「項目によらない」のであれば，どんな尺度 で推定した場合でも 𝜃𝑝 の値は同じになるはずです。いま，𝑝 さんは同じ構成概念を測定する 尺度 1 と尺度 2 による 2 回の調査の両方に回答しました。この 2 回のデータをもとにパラメ ータを推定する場合には，いままでの手順にのっとってそれぞれ 𝜃𝑝 ∼ 𝑁 (0, 1) と仮定して推 定を行います。その結果 (1) • 尺度 A での推定値 𝜃𝑝̂ = 0.5 (2) • 尺度 B での推定値 𝜃𝑝̂ = −0.5 だったとします。尺度 A を基準に考えると，尺度 B での推定値は本来の値より 1 小さい値 (2) になっているようです。ですが個人特性値は「項目によらない」であるべきなので，𝜃𝑝̂ も (1) 𝜃𝑝̂ と同じ 0.5 になるように調整する必要がありそうです。 79 *29 R での実践例までやると膨大になってしまうので，概念の紹介だけにとどめます。気になる人は自分で調べ *30 たり直接質問してください。 R では equateIRT や plink といったパッケージによって等化のメソッドが提供されています。

8.12 IRT の使い道 もう少し規模を広げて，同じようにこの尺度 1 と尺度 2 による 2 回の調査の両方に回答した 人が 100 人いたとします。その結果，2 回の推定値の散布図が図8.43のようになったとしま す。そしてこの 2 回の推定値の回帰直線を求めたら (2) (1) 𝜃𝑝̂ = 𝐴𝜃𝑝̂ + 𝐵 (2) (8.73) (1) であったとします。この場合，全員の 𝜃𝑝̂ をなるべく 𝜃𝑝̂ に近づけるためには， (2−1) 𝜃𝑝̂ (2) = 𝜃𝑝̂ − 𝐵 𝐴 (8.74) と変換してあげるのがよさそうだとわかります。 4 3 qp (2 ) 2 1 0 -1 -2 -1 0 (1 ) 1 2 3 qp 図8.43 2 回の測定における推定値 2 パラメータロジスティックモデルでは，項目反応関数は 𝑃 (𝑦𝑝𝑖 = 1) = 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.75) でした。ここで， 𝜃𝑝 − 𝐵 𝐴 𝛽 −𝐵 𝛽𝑖∗ = 𝑖 𝐴 𝛼∗𝑖 = 𝐴𝛼𝑖 𝜃𝑝∗ = (8.76) とおくと， 𝑃 (𝑦𝑝𝑖 = 1) = = 1 1 + exp [−𝛼∗𝑖 (𝜃𝑝∗ − 𝛽𝑖∗ )] 1 1 + exp [−𝐴𝛼𝑖 ([ = 𝜃𝑝 −𝐵 𝐴 ] − [ 𝛽𝑖𝐴−𝐵 ])] (8.77) 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] となり，もとの項目反応関数と同じになります。因子分析の Chapter では，潜在変数のスケ ールの不定性の問題から，全ての観測変数を標準化していました。等化では，このスケール 80

8.12 IRT の使い道 の不定性をある種逆手に取ることで 𝛼, 𝛽 や 𝜃 を適当なスケールに変換しているわけです。具 (2) (1) 体的には，2 つの尺度に両方とも回答した全員の 𝜃𝑝̂ がなるべく 𝜃𝑝̂ に近づくようにスケー ルを変換し，これに応じて項目パラメータにも変換をかけます。こうして得られた変換後の 項目パラメータならば，尺度 A と尺度 B でも問題なく比較ができるようになります。 【共通項目による等化】 先程説明した等化の流れは，共通の項目がある場合でも同じように 考えられます。つまり，2 つの集団に対して実施する際に，それぞれの尺度に共通の項目を いくつか入れておけば， (2−1) 𝛽𝑖 (2−1) 𝛼𝑖 (2) 𝛽𝑖 − 𝐵 𝐴 (2) = 𝐴𝛼𝑖 = (8.78) となるような等化係数 𝐴, 𝐵 を求めることができ，それに応じた変換が可能となります。等 化係数 𝐴, 𝐵 を求める方法としては，推定値の要約統計量を用いる方法 (mean-mean 法, mean-sigma 法) や， 「変換後の ICC がなるべく重なるようにする」方法 (Stocking-Lord 法: Stocking and Lord, 1983, Haebara 法: Haebara, 1980) が一般的です。 等化は IRT の可能性を大きく広げてくれる重要な技術ですが，その適用には様々なハードル があります。何より 𝜃 が等質であるという大前提が重要です。もしも尺度 1 と尺度 2 がそも そも異なる構成概念を測定しているならば，パラメータを調整したところで同じモノサシの 上で比較することはナンセンスです。本来は全部まとめてデータを取った上で一次元性の検 証ができれば良いのですが，等化をしたい場合にはそんなことができないので，理論的側面 などから測定内容の等質性を担保してあげる必要があります。また，等化係数はパラメータ の推定値をもとに計算されます。つまり，等化の係数には少なからず誤差がある，というこ とです。厄介なのは，尺度 2 での項目パラメータの推定値それ自体にも誤差があるという点 です。つまり等化後のパラメータには「等化前の誤差＋等化係数の誤差」という 2 つの誤差 が重なっているので，等化をしない場合と比べるとどうしても推定値は不安定になりやすい です。理論的には，等化は繰り返し行うことも可能です（尺度 4→ 尺度 3→ 尺度 2→ 尺度 1） が，実際に等化を繰り返すとなかなか厳しいことになるので注意が必要です。 8.12.2 テストの設計 テストの一番の目的は，回答者の特性値 𝜃 を精度良く推定することです。ですがテストの目 的によって「どの範囲の 𝜃 をどの程度の精度で推定できたら良いか」は変わってきます。例 えばふつうの心理尺度では，比較的広い範囲の特性値をそれなりの精度で推定したいと思う でしょう。この場合には，困難度が高い項目から低い項目までをまんべんなく用意するのが 良さそうです。一方で，心療内科が扱うような「診断」を目的とした尺度や合格／不合格を 決めるテストなどの場合，ある特性の値が具体的にどの程度かよりも「基準値以上か以下か」 に強い関心があったりします。このような場合には，困難度が基準値周辺の，識別力ができ るだけ高い項目ばかりを用意するほうが良さそうです。このように最適な項目の特性はテス トの目的に応じて異なります。 IRT では項目・テスト情報関数が 𝜃 の関数になっているおかげで，様々な 𝜃 の値に合わせて オーダーメイドで尺度を構成することができます。幸いなことにテスト情報関数は局所独立 の仮定が満たされていれば項目情報関数の和なので，他の項目との関係を気にすることなく， 81

8.12 IRT の使い道 一つ一つの項目に対して「その項目を追加したらテスト情報量がどの程度増加するか」のみ を考えたらよいわけです。特に自動テスト構成 (automated test assembly) では，コンピュ ータの力を使って，例えば「𝜃 = (−4, −2, 0, 2, 4) の各点における標準誤差が 0.2 以下＝テス ト情報量が 25 以上になる最小項目数のセット」などの条件にあった最適な項目セットを探 し出してくれます*31 。 ただこれだけのテストの設計を行うためにはかなりの事前準備が必要となります。というの も項目情報関数は項目パラメータが既知でないとわかりません。また，目的に応じて最適な 項目セットを作成するために，なるべく候補となる項目の数は多いほうが良いです。もちろ んすべての項目のパラメータを一気に推定するわけにはいかないので，多くの場合は複数回 に分けて集めたデータをもとにパラメータを等化したりして項目プールを作成する必要があ るでしょう。ということで心理尺度の場合にはなかなかそこまで大量の項目を用意すること 自体が難しいので，ここまで理論的な尺度構成を行っている先行事例はたぶんほとんど存在 しません。一応用途としては，100-200 項目くらい用意した候補の中から「最高の 20 項目」 を選んで心理尺度を作成する，といった方法はありそうな気がします*32 。 8.12.3 適応型テスト テストや尺度の中には通常様々な困難度の項目が混ざっているため，ほぼ確実に「その人に はほぼ意味のない項目」が含まれることになります。この問題は，先ほど紹介したテストの 設計の考え方を持ってしても避けられません。そこで自動テスト構成をもっと局所的に行う と「個人に応じて最適な項目セットを選び出す」という考えになります。つまり視力検査の ように 𝜃𝑝 が高そうな人にはより困難度の高い項目を，𝜃𝑝 が低そうな人にはより困難度の低 い項目を出していけば，効率的にテスト情報量を稼ぐことができそうです。コンピュータ上 で行われるコンピュータ適応型テスト (computerized adaptive testing [CAT]) *33 では， この考えに基づいて個人ごとに次の項目を出し分けていくことができます。 ここでは IRT に基づく最もシンプルな項目選択法を紹介しましょう。推定の標準誤差は 𝜎 (𝜃 ̂ 𝑝 |𝜃𝑝 ) = 1 (8.79) √𝐼(𝜃𝑝 ) で表されていました。したがって推定精度をなるべく上げるためには，𝐼(𝜃𝑝 ) を大きくする ＝項目情報量 𝐼𝑖 (𝜃𝑝 ) が大きい項目を選んでいけばよいわけです。ただし 𝐼𝑖 (𝜃𝑝 ) は「真値 𝜃𝑝 における項目情報量」のため，真値 𝜃𝑝 がわからない状況では 𝐼𝑖 (𝜃𝑝 ) の値もわかりません。そ こで，1 問解答するたびに推定を行い，その時点での暫定的な推定値 𝜃𝑝∗̂ における項目情報量 𝐼𝑖 (𝜃𝑝∗̂ ) が最大となる項目をとりあえず選び続けたらなんだかうまく行きそうです*34 。 *31 R では lpSolve や glpkAPI，Rsymphony などのパッケージを使うと自動テスト構成ができるようになりま *32 100-200 項目すべてに真面目に回答してくれる人たちを数百人用意する必要があるのですが，それさえクリ アできれば… R では catR や mirtCAT といったパッケージが役に立つかもしれません。 す。 *33 *34 82 もちろん他にも項目の選び方は色々あります。ざっくりいうと「推定精度をなるべく高める」方向性と「なる べく人によって異なる項目が出るようにする（項目情報が外部に漏れるのを防ぐため）」という方向性があり ます。

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

参考文献 参考文献 Berrío, Á. I., Gómez-Benito, J., and Arias-Patiño, E. M. (2020). Developments and trends in research on methods of detecting differential item functioning. Educational Research Review, 31, 100340. https://doi.org/10.1016/j.edurev.2020.100340 Bock, R. D. (1972). Estimating item parameters and latent ability when responses are scored in two or more nominal categories. Psychometrika, 37 (1), 29–51. https://doi. org/10.1007/BF02291411 Böckenholt, U. (2017). Measuring response styles in Likert items. Psychological Methods, 22(1), 69–83. https://doi.org/10.1037/met0000106 Bowling, S. R., Khasawneh, M. T., Kaewkuekool, S., and Cho, B. R. (2009). A logistic approximation to the cumulative normal distribution. Journal of Industrial Engineering and Management, 2(1), 114–127. https://doi.org/10.3926/jiem.2009.v2n1.p114127 Brown, A. and Maydeu-Olivares, A. (2011). Item response modeling of forced-choice questionnaires. Educational and Psychological Measurement, 71(3), 460–502. https: //doi.org/10.1177/0013164410375112 Chalmers, R. P. (2018). Model-based measures for detecting and quantifying response bias. Psychometrika, 83(3), 696–732. https://doi.org/10.1007/s11336-018-9626-9 Chen, W.-H. and Thissen, D. (1997). Local dependence indexes for item pairs using item response theory. Journal of Educational and Behavioral Statistics, 22(3), 265–289. https://doi.org/10.2307/1165285 Choi, Y.-J. and Asilkalkan, A. (2019). R packages for item response theory analysis: Descriptions and features. Measurement: Interdisciplinary Research and Perspectives, 17 (3), 168–175. https://doi.org/10.1080/15366367.2019.1586404 Drasgow, F., Levine, M. V., and Williams, E. A. (1984). Appropriateness measurement with polychotomous item response models and standardized indices.Technical Report ADA141365, Model Based Measurement Laboratory, University of Illinois Ferrando, P. J. and Lorenzo-Seva, U. (2007). An item response theory model for incorporating response time data in binary personality items. Applied Psychological Measurement, 31(6), 525–543. https://doi.org/10.1177/0146621606295197 84

• https://doi.org/10.1016/j.edurev.2020.100340
• https://doi.org/10.1007/BF02291411
• https://doi.org/10.1037/met0000106
• https://doi.org/10.3926/jiem.2009.v2n1.p114-127
• https://doi.org/10.1177/0013164410375112
• https://doi.org/10.1007/s11336-018-9626-9
• https://doi.org/10.2307/1165285
• https://doi.org/10.1080/15366367.2019.1586404
• https://doi.org/10.1177/0146621606295197

参考文献 Haebara, T. (1980). Equating logistic ability scales by a weighted least squares method. Japanese Psychological Research, 22(3), 144–149. https://doi.org/10.4992/psycholr es1954.22.144 南風原朝和 (2000). 局所独立性，https://www.p.u-tokyo.ac.jp/~haebara/local_ind/ Jackson, D. N., Wroblewski, V. R., and Ashton, M. C. (2000). The impact of faking on employment tests: Does forced choice offer a solution? Human Performance, 13(4), 371–388. https://doi.org/10.1207/S15327043HUP1304_3 Karabatsos, G. (2003). Comparing the aberrant response detection performance of thirty-six person-fit statistics. Applied Measurement in Education, 16(4), 277–298. https://doi.org/10.1207/S15324818AME1604_2 加藤健太郎・山田剛史・川端一光 (2014). R による項目反応理論 オーム社 Kleinman, M. and Teresi, J. A. (2016). Differential item functioning magnitude and impact measures from item response theory models. Psychological test and assessment modeling, 58(1), 79–98. Kolen, M. J. and Brennan, R. L. (2014). Test equating, scaling, and linking (3rd ed.). Springer. Kopf, J., Zeileis, A., and Strobl, C. (2015). Anchor selection strategies for DIF analysis. Educational and Psychological Measurement, 75(1), 22–56. https://doi.org/10.1177/ 0013164414529792 Kreuter, F. (Ed.) (2013). Improving surveys with paradata: Analytic use of process information, Hoboken, New Jersey. John Wiley & Sons. Lord, F. M. (1980). Applications of item response theory to practical testing problems, Hillsdale, N.J. L. Erlbaum Associates. Lord, F. M. and Novick, M. R. (1968). Statistical theories of mental test scores. AddisonWesley. Mantel, N. and Haenszel, W. (1959). Statistical aspects of the analysis of data from retrospective studies of disease. Journal of the National Cancer Institute, https: //doi.org/10.1093/jnci/22.4.719 Maydeu-olivares, A., Hernández, A., and Mcdonald, R. P. (2006). A multidimensional ideal point item response theory model for binary data. Multivariate Behavioral Research, 41(4), 445–472. https://doi.org/10.1207/s15327906mbr4104_2 Meijer, R. R. and Sijtsma, K. (2001). Methodology review: Evaluating person fit. Applied Psychological Measurement, 25(2), 107–135. https://doi.org/10.1177/014662101220 31957 85

• https://doi.org/10.4992/psycholres1954.22.144
• https://www.p.u-tokyo.ac.jp/~haebara/local_ind/
• https://doi.org/10.1207/S15327043HUP1304_3
• https://doi.org/10.1207/S15324818AME1604_2
• https://doi.org/10.1177/0013164414529792
• https://doi.org/10.1093/jnci/22.4.719
• https://doi.org/10.1207/s15327906mbr4104_2
• https://doi.org/10.1177/01466210122031957

参考文献 Meijer, R. R., Niessen, A. S. M., and Tendeiro, J. N. (2016). A practical guide to check the consistency of item response patterns in clinical research through personfit statistics: Examples and a computer program. Assessment, 23(1), 52–62. https: //doi.org/10.1177/1073191115577800 Muraki, E. (1992). A generalized partial credit model: Application of an EM algorithm. Applied Psychological Measurement, 16(2), 159–176. https://doi.org/10.1177/014662 169201600206 Orlando, M. and Thissen, D. (2000). Likelihood-based item-fit indices for dichotomous item response theory models. Applied Psychological Measurement, 24(1), 50–64. https: //doi.org/10.1177/01466216000241003 Raju, N. S. (1988). The area between two item characteristic curves. Psychometrika, 53(4), 495–502. https://doi.org/10.1007/BF02294403 Raju, N. S., van der Linden, W. J., and Fleer, P. F. (1995). IRT-based internal measures of differential functioning of items and tests. Applied Psychological Measurement, 19(4), 353–368. https://doi.org/10.1177/014662169501900405 Rasch, G. (1980). Probabilistic models for some intelligence and attainment tests. University of Chicago Press. (Original work published 1960) Reckase, M. D. (2019). Logistic multidimensional models. In van der Linden, W. J. (Ed.) Handbook of item response theory volume one: Models (pp. 189–209.). CRC Press. Roberts, J. S. (2019). Generalized graded unfolding model. In van der Linden, W. J. (Ed.) Handbook of item response theory volume one: Models (pp. 369–392.). CRC Press. Samejima, F. (1969). Estimation of latent ability using a response pattern of graded scores. Psychometrika, 34(S1), 1–97. https://doi.org/10.1007/BF03372160 Smith, R. M., Schumacker, R. E., and Busch, M. J. (1995). Using item mean squares to evaluate fit to the rasch model. A paper presented at the Annual Meeting of the American Educational Research Association Stark, S., Chernyshenko, O. S., Drasgow, F., and Williams, B. A. (2006). Examining assumptions about item responding in personality assessment: Should ideal point methods be considered for scale development and scoring? Journal of Applied Psychology, 91(1), 25. https://doi.org/10.1037/0021-9010.91.1.25 Stocking, M. L. and Lord, F. M. (1983). Developing a common metric in item response theory. Applied Psychological Measurement, 7 (2), 201–210. https://doi.org/10.1177/ 014662168300700208 Stout, W., Froelich, A. G., and Gao, F. (2001). Using resampling methods to produce an 86

• https://doi.org/10.1177/1073191115577800
• https://doi.org/10.1177/014662169201600206
• https://doi.org/10.1177/01466216000241003
• https://doi.org/10.1007/BF02294403
• https://doi.org/10.1177/014662169501900405
• https://doi.org/10.1007/BF03372160
• https://doi.org/10.1037/0021-9010.91.1.25
• https://doi.org/10.1177/014662168300700208

参考文献 improved DIMTEST prodecure. In Boomsma, A., van Duijn, M. A. J., and Snijders, T. A. B. (Eds.) Essays on item response theory. Springer. Thurstone, L. L. (1927). A law of comparative judgment. Psychological Review, 34(4), 273–286. https://doi.org/10.1037/h0070288 豊田秀樹（編）(2013). 項目反応理論［中級編］ 朝倉書店 Yen, W. M. (1984). Effects of local item dependence on the fit and equating performance of the three-parameter logistic model. Applied Psychological Measurement, 8(2), 125– 145. https://doi.org/10.1177/014662168400800201 Zhang, J. and Stout, W. (1999). The theoretical DETECT index of dimensionality and its application to approximate simple structure. Psychometrika, 64(2), 213–249. https://doi.org/10.1007/BF02294536 87

• https://doi.org/10.1037/h0070288
• https://doi.org/10.1177/014662168400800201
• https://doi.org/10.1007/BF02294536