2.6K Views
July 25, 23
スライド概要
神戸大学経営学研究科で2022年度より担当している「統計的方法論特殊研究(多変量解析)」の講義資料「08_項目反応理論」です。
2025/2/5: 全体的に内容に手を入れました。RmarkdownからQuartoに変更しました。
神戸大学経営学研究科准教授 分寺杏介(ぶんじ・きょうすけ)です。 主に心理学的な測定・教育測定に関する研究を行っています。 講義資料や学会発表のスライドを公開していきます。 ※スライドに誤りを見つけた方は,炎上させずにこっそりお伝えいただけると幸いです。
Chapter 8 項目反応理論 Chapter 8 項目反応理論 データから「回答者の特性値」と「項目の特性値」を同時に推定する分析法である項目反応理論 について,前半では理論的な考え方と mirt パッケージでのやり方を説明し,後半では実践上の いくつかのトピック(多値,多次元モデル・テスト情報量・様々な拡張など)を解説しています。 本資料は,神戸大学経営学研究科で 2022 年度より担当している「統計的方法論特殊研究(多 変量解析) 」の講義資料です。CC BY-NC 4.0 ライセンスの下に提供されています。 作成者連絡先 神戸大学大学院経営学研究科 分寺杏介(ぶんじ・きょうすけ) mail: [email protected] website: https://www2.kobe-u.ac.jp/~bunji/ この Chapter では項目反応理論(Item Response Theory [IRT]: Lord & Novick, 1968)のお話をしていきます。 後ほど説明しますが,実は IRT の最もメジャーなモデルの一つは,因子分析と数学的に同じ ものです。そのため lavaan でも基本的な IRT 分析はできるのですが,IRT 専用のパッケージ も色々あります*1 。この講義では mirt というパッケージを使用します。 事前準備 1 2 3 4 # install.packages("mirt") library(mirt) dat <- readRDS("chapter04.rds") # 毎回列名を指定するのが面倒なので先に作っておく *1 Choi & Asilkalkan(2019)によれば,IRT ができる R パッケージは 45 個もあるらしいです(今はもっと多 い?)。細かく目的や用途が異なると思うので,特定の場面ではもっと使い勝手の良いパッケージがあったりす るかもしれません。あとは計算の精度や速度の違いなど? 1
8.1 項目反応理論とは 5 2 cols <- paste0("Q1_", 1:10) 8.1 項目反応理論とは 項目反応理論は,古典的テスト理論と同じく,テスト理論と呼ばれるカテゴリに属するものの 一つです。テスト理論の目的をとてもシンプルに言うと,以下の 2 つです。 • 回答者の潜在的な特性(心理特性の強さ・能力など)を測定する • 測定に使用した項目(心理尺度の各項目・テストの問題など)の性能を評価する 1 点目は直感的にもすぐ分かると思います。2 点目について少し補足しておくと,「項目の性 能」とは,単純にはテストの問題の難易度などを指しています。もしも難しすぎてほとんど誰も 正解していない問題があったとしたら,その問題には,いわば床効果が出ているわけなのであま り良い問題とは言えません。また,因子分析的に見たときに因子負荷が低い項目があったとした ら,その項目への回答は,測定したい「回答者の潜在的な特性」を反映していないので,注意が 必要かもしれません。テスト理論では,そのような形で項目の良し悪しを評価することで,「回 答者の潜在的な特性」をより良く測定することを目指します。 そんなテスト理論の中で,項目反応理論が目指しているところを理解するために,仮想的な状 況を考えてみましょう。Chapter 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 人の標準化得点を計算してみたところ,そ れぞれ A さんが 0.2 点,B さんが-0.25 点だったとしましょう。2 人はそれぞれ平均点以上・以
8.1 項目反応理論とは 3 下だったので,もちろん標準化得点でも A さんのほうが高い得点になっています。ただ,実はま だこの情報だけでは A さんのほうがよりバーンアウトに近いとは言い切れません。それは,2 人 の所属するグループの平均値が同じとは限らないためです。 標準化得点は,当然ながらそのグループの得点分布に依存します。もしかしたら,UWES の 回答者平均点のほうが MBI-GS よりも高い理由は,UWES のほうが高得点になりやすい項目な のではなく,単に A さんの所属するグループよりも B さんの所属するグループの方が平均的に バーンアウト傾向が強いためかもしれません。このように,心理尺度得点に基づく評価は項目だ けでなく集団にも依存します(集団依存性) 。 IRT では,項目依存性と集団依存性をクリアした得点を算出するために,まず「集団によらな い項目パラメータ」と「項目によらない個人特性値」の存在を考えています。もちろんこれらの パラメータは,単純な和得点・平均値からはわかりません。が,仮にそのようなパラメータがあ るとして,項目パラメータが 𝛽,個人特性値が 𝜃 で表せるとしたら,項目 𝑖 に対する回答者 𝑝 の 反応(回答)は 𝑦𝑝𝑖 = 𝑓(𝜃𝑝 , 𝛽𝑖 ), (8.1) つまり 𝛽𝑖 と 𝜃𝑝 による何らかの関数(項目反応関数: item response function [IRF])になる はずです。項目反応関数の形はいろいろと考えられますが,例えば多くの心理尺度に対しては, 𝜃𝑝 の値が大きいほど「当てはまる」寄りの回答をする(= 𝑦𝑝𝑖 の値が大きくなる)と考えられま す。そのようなデータを IRT で分析する場合には,項目反応関数として「𝜃𝑝 に対して単調増加 の関数」を持ってきたら良いわけです。…という感じで IRT では,観測されたデータにおいて想 定される 𝑦𝑝𝑖 , 𝜃𝑝 , 𝛽𝑖 の関係性をうまく表現できるような項目反応関数が色々と提案されてきまし た。後ほど IRT の代表的なモデルを紹介しますが,IRT の本質はこのように,回答を (集団によ らない) 項目パラメータと (項目によらない) 個人特性値に分離して分析することにあります*2 。 ということで,一般的な IRT の利用場面では, 1. 特定の関数形(IRF)を仮定する 2. そのモデルのもとで項目パラメータを推定する 3. 推定した項目パラメータをもとに,項目の性能を評価する 4. 合わせて個人特性値を推定する 5. 推定した個人特性値を使ってあとはご自由に といったことが行われます。それでは,まずは「とりあえず IRT で分析」ができるようになる ために必要な最低限の知識を学んでいきましょう。 *2 「項目反応理論」というと,ほとんどの場合この後紹介する代表的なモデルを指すのですが,モデル式(項目反 応関数)がそれだけしか無いのではなく,実際には様々なモデル式による形が考えられるというのは重要なポイ ントです。「項目反応理論」という言葉はあくまでも項目依存性や集団依存性を解消したパラメータを得るため の理論的な方法を考える土台としての理論である,という認識です。
8.2 基本的なモデル 4 8.2 基本的なモデル 最もシンプルな IRT モデルを考えるため,一旦ここからは項目への反応を二値(当てはまら ない・あてはまる)に限定して考えてみます。というのも,もともと IRT は「正解・不正解」の 二値で考えるような学力テストなどを想定して作られているもののためです。 8.2.1 正規累積モデル 以前カテゴリカルデータの相関(ポリコリック・ポリシリアル相関)の計算時には,背後にあ る連続量を考えるという説明をしました。ここでも,二値反応 𝑦𝑝𝑖 の背後に連続量 𝜃𝑝 が存在し, 𝜃𝑝 が閾値を超えていたら 1,超えていなければ 0 と回答すると考えます。そして項目の閾値を 𝛽𝑖 と表すと, 𝑦𝑝𝑖 = { 0 𝜃𝑝 < 𝛽𝑖 1 𝜃𝑝 ≥ 𝛽𝑖 (8.2) となります( 図 8.1 ) 。𝛽𝑖 の値が大きい項目ほど,𝑦𝑝𝑖 = 1 になる 𝜃𝑝 の範囲が狭くなるため, 「当 てはまらない」と回答する人の割合が高くなるだろう,ということです。 q p ~ N(0, 1) ypi = 0 ypi = 1 bi 図 8.1: 二値反応と潜在的な連続量の関係 ただし,このままでは能力 𝜃𝑝 を持つ 𝑝 さんは難易度が 𝛽𝑖 以下の問題には,どんなジャンル・ 内容の問題でも 100% 正解することになってしまいます(𝜃𝑝 ≥ 𝛽𝑖 なので) 。しかし,もちろん客 観的に見た難易度が全く同じであっても,内容によって 𝑝 さんが答えられる・答えられない問題 があるはずです。さらに言えば,𝑝 さんの体調や時の運(たまたま昨日一夜漬けした内容に入っ ていた,など)によっても答えられるかどうかは変わるはずです。というように,心理尺度やテ ストへの回答では,常に「その他の要因」による変動を考えてきました。言い換えると,古典的 テスト理論でも回帰分析でも因子分析でも,観測値は必ず何らかの誤差を伴って発生するものだ と考えているわけです。そこで先程の (8.2) 式に誤差項 𝜀𝑝𝑖 を追加して 𝑦𝑝𝑖 = { 0 (𝜃𝑝 − 𝜀𝑝𝑖 ) < 𝛽𝑖 1 (𝜃𝑝 − 𝜀𝑝𝑖 ) ≥ 𝛽𝑖 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀2 ) (8.3) とします。なお誤差項は,回帰分析などと同じく,最も一般的な「平均 0 の正規分布に従う」と 仮定してあります。
8.2 基本的なモデル 5 誤差項を追加することによって,「特性値が 𝜃𝑝 の人が項目 𝑗 に『当てはまる』と回答する確 率*3 」を考えることが出来ます。図 8.1 では「全ての人の 𝜃𝑝 」の分布を示していた一方で,図 8.2 では「𝜃𝑝 が特定の値の人」における分布になっている点に注意してください。 q p - e pi ~ N(q p, s2e ) P(ypi = 0) P(ypi = 1) qp bi 図 8.2: 誤差を追加すると さらに 𝛽𝑖 と 𝜃𝑝 を一つの関数 𝑓(𝜃𝑝 , 𝛽𝑖 ) として見るために 𝛽𝑖 を移項して, 0 𝑦𝑝𝑖 = { 1 (𝜃𝑝 − 𝛽𝑖 − 𝜀𝑝𝑖 ) < 0 (𝜃𝑝 − 𝛽𝑖 − 𝜀𝑝𝑖 ) ≥ 0 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀2 ) (8.4) とします( 図 8.3 ) 。当然ながら,図 8.2 と 図 8.3 は,閾値と分布が平行移動しただけなので 𝑃 (𝑦𝑝𝑖 = 1) は変わりません。 q p - b i - e pi ~ N(q p - b i, s2e ) P(ypi = 1) P(ypi = 0) 0 q p - bi 図 8.3: 𝛽𝑖 を移項すると 図 8.3 の色つきの確率分布は,正規分布 𝑁 (𝜃𝑝 − 𝛽𝑖 , 𝜎𝜀2 ) と表すことができます。そして式か ら,𝑦𝑝𝑖 = 1 となるのは誤差 𝜀𝑝𝑖 が (𝜃𝑝 − 𝛽𝑖 ) 以下のときだとわかります。したがって,仮に 𝜀𝑝𝑖 が標準正規分布に従う(𝜎𝜀2 = 1)としてみると, 「特性値が 𝜃𝑝 の人が項目 𝑗 に当てはまると回答 *3 頻度論的に言えば,ここでの「確率」は「ある 𝜃 𝑝 の一個人が何回も回答した際の当てはまるの割合」や「母集 団内での 𝜃𝑝 の人たちにおける当てはまるの割合」と解釈できます。これはどちらでも OK です。
8.2 基本的なモデル 6 する確率」 (図 8.4 の青い部分)は 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 (𝜀𝑝𝑖 ≤ 𝜃𝑝 − 𝛽𝑖 ) (𝜃𝑝 −𝛽𝑖 ) =∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.5) と表すことができます。教科書などで見たことがある形かもしれませんが,この積分の中身は標 準正規分布の確率密度関数です。このように,𝑃 (𝑦𝑝𝑖 = 1) を標準正規分布の累積確率で表すモデ ルのことを正規累積モデル (normal ogive model) と呼びます。この後項目パラメータを 1 つ 増やすので,それに対してこの正規累積モデルは,特に 1 パラメータ正規累積モデルと呼ばれ ます。 e pi ~ N(0, 1) P(ypi = 1) P(ypi = 0) 0 q p - bi 図 8.4: 最終的な反応確率の分布 ここまでは,途中で話を簡単にするために勝手に 𝜎𝜀 = 1 として進めていましたが,ここから はもう少し一般化して 𝜎𝜀2 が項目ごとに異なる (𝜎𝜀2 = 1/𝛼𝑖 ) とします。このとき,誤差の確率 分布は 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝛼1 ) となるため,両辺を 𝛼𝑖 倍して 𝛼𝑖 𝜀𝑝𝑖 ∼ 𝑁 (0, 1) と表すことができます。 𝑖 𝛼𝑖 𝜀𝑝𝑖 が標準正規分布に従うならば,(8.5) 式の正規累積モデルは 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 [𝜀𝑝𝑖 ≤ (𝜃𝑝 − 𝛽𝑖 )] = 𝑃 [𝛼𝑖 𝜀𝑝𝑖 ≤ 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] 𝛼𝑖 (𝜃𝑝 −𝛽𝑖 ) =∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.6) という形で書き換えることができます。また正規累積モデルが出てきました。(8.5) 式との違い は積分の上限が 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 ) に変わった点だけです。いま,突拍子もなく 𝛼𝑖 というパラメータを 導入しましたが,とりあえずこれによって項目 𝑖 が持つパラメータは (𝛽𝑖 , 𝛼𝑖 ) の 2 つになりまし た。ということで,(8.6) 式の正規累積モデルは 2 パラメータ正規累積モデルと呼ばれます。 8.2.2 項目特性曲線 ここで,2 つの項目パラメータ (𝛽𝑖 , 𝛼𝑖 ) が何を表しているのかを明らかにしていくため,図 8.3 を別の形で表してみます。 図 8.5 の下半分は,横軸に 𝜃𝑝 − 𝛽𝑖 を,縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとったグラフです。図の上半分に
8.2 基本的なモデル 7 P(ypi = 1) q p - bi 0 図 8.5: パラメータと正答率の関係 あるように,𝜃𝑝 − 𝛽𝑖 の値が大きくなるほど分布の青い部分が増える,つまり 𝑃 (𝑦𝑝𝑖 = 1) の値が 大きくなることを表しています。 これを踏まえて,パラメータ 𝛽𝑖 の挙動を見るために,横軸に 𝜃𝑝 を,縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとっ たグラフを,𝛽𝑖 の値を変えながらいくつか重ねてみます(図 8.6 ) 。誤差は平均 0 の正規分布に 従うため,𝑃 (𝑦𝑝𝑖 = 1) = 0.5 になるのは,𝜃𝑝 = 𝛽𝑖 のときです。したがって 𝛽𝑖 の値が大きい項目 ほど,𝜃𝑝 の値が大きくないと「当てはまる」と回答する確率が高くなりません。そのような意味 で,𝛽𝑖 は項目困難度 (item difficulty) と呼ばれます*4 。 次は 𝛼𝑖 の役割を見るために,図 8.6 よりも誤差分散 (𝛼𝑖 ) が小さい場合や大きい場合を見てみ ましょう。図 8.7 の左側に示されているように,誤差分散が小さいほど分布の広がりがなくなる ため,𝜃𝑝 − 𝛽𝑖 の値の変化に対する 𝑃 (𝑦𝑝𝑖 = 1) の変化が急激になっています。 これを踏まえて,𝛽𝑖 の時と同じように横軸に 𝜃𝑝 をとり,縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとったグラフ を,今度は 𝛼𝑖 の値を変えながらいくつかプロットしたものが図 8.8 です(𝛽𝑖 = 0) 。誤差分散の 逆数である 𝛼𝑖 の値が大きい項目ほど,𝜃𝑝 が低い人と高い人の間で 𝑃 (𝑦𝑝𝑖 = 1) が大きく変化し ています。言い換えると 𝛼𝑖 の値が大きい項目ほど,その項目への回答によって 𝜃𝑝 の高低を識別 する能力が高いといえます。ということで 𝛼𝑖 のことを項目識別力 (item discrimination) と 呼びます。 図 8.6 や 図 8.8 のように,横軸に 𝜃𝑝 をとり,縦軸に 𝑃 (𝑦𝑝𝑖 = 1) をとったグラフは,いわば 回答者の特性値に対する項目の特性を表しています。ということで,これらのグラフのことを項 *4 もとが学力テストなので「困難度」という呼び方をしています。正解・不正解がない心理尺度でも特に気にせず 「困難度」と呼びます。
8.2 基本的なモデル 8 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: 左:誤差分散が小さい場合,右:誤差分散が大きい場合 P(ypi = 1) ai = 0.5 ai = 1 ai = 2 qp 図 8.8: 𝛼𝑖 の役割
8.2 基本的なモデル 9 目特性曲線 (item characteristic curve [ICC]) と呼びます。 8.2.3 ロジスティックモデル 2 パラメータ正規累積モデルの式を再掲しておきましょう。 𝛼𝑖 (𝜃𝑝 −𝛽𝑖 ) 𝑃 (𝑦𝑝𝑖 = 1) = ∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 正規累積モデルには,積分計算が含まれています。ある程度想像がつくかもしれませんが,コン ピュータは積分計算が比較的得意ではありません。ということで,積分計算がいらない別のモデ ルを考えてみましょう。 正規累積モデルがやっていることは,プロビット回帰と同じです。そんな正規累積モデルの代 わりになるものということは,正規分布のように • 累積分布関数が S 字になっていて • 確率密度関数は左右対称の山型 な分布があれば良いわけですが…そのような確率分布として,ロジスティック分布というもの がありました。ロジスティック分布の確率密度および累積分布関数は 𝑓(𝑥) = exp(−𝑥) [1 + exp(−𝑥)] 2 , 𝐹 (𝑥) = 1 1 + exp(−𝑥) (8.7) と表され,確率密度関数を正規分布と重ねて見ると図 8.9 のようになります。こうしてみると, 裾の重さは違いますが確かにロジスティック分布は正規分布のような左右対称形になっているよ うです。 0.4 正規分布 ロジスティック分布 f(x) 0.3 0.2 0.1 0.0 -4 -2 0 X 2 4 図 8.9: 正規分布とロジスティック分布 また,ロジスティック分布の 𝑥 をおよそ 1.7 倍すると,累積分布関数が正規分布とかなり近く なることが知られています ( 図 8.10 )*5 。 *5 もっと近づけるためには 1.704 倍のほうがよかったり,ズレを極限まで減らすためには 0.07056𝑥3 + 1.5976𝑥 にすると良い(Bowling et al., 2009)ということは分かっていますが,そこまでの精度は必要ない,というこ
8.2 基本的なモデル 10 1.00 F(x) 0.75 正規分布 1.7倍ロジスティック分布 0.50 0.25 0.00 -4 -2 0 X 2 4 図 8.10: 1.7 倍ロジスティック分布 ということで,正規分布の代わりにロジスティック分布を使った (2 パラメータ) ロジスティッ クモデルは 𝑃 (𝑦𝑝𝑖 = 1) = 1 1 + exp [−1.7𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.8) と表されます。いま説明したように,この 1.7 はあくまでも「正規累積モデルに近づけるための 定数」なので,特に正規累積モデルとの比較をする予定が無いならば 1.7 は無くても問題ありま せん。ただ,結果を報告する際には「正規累積モデルとロジスティックモデルのどちらを使用し たのか」「ロジスティックモデルならば 1.7 倍した値なのか」がハッキリわかるようにしておき ましょう。論文などでは,モデル式を書いておけば一発です。 (以後の説明では,特に正規累積モデルと比較することもないので 1.7 は省略していきます) また,1 パラメータロジスティックモデルはラッシュモデル (Rasch model) と呼ばれること もあります。これは,IRT とは独立に Rasch という人が理論的に「1 パラメータロジスティック モデルと同じ関数」を導出した(Rasch, 1960/1980)ことに由来します。ということで「ラッ シュモデル」にはやや高尚な理念があったりしますが,私もよく理解できていないのでそこの説 明はしません。とりあえず「ラッシュモデル」という名前が出てきたら「困難度だけのロジス ティックモデルなんだな」と理解できていれば良いでしょう。厳密な使い分けとしては,ラッ シュモデルは「正規累積モデルの近似」を目的としていないので,1.7 倍しない 1 パラメータモ デルを指す,という場合もあるようです。 8.2.4 モデルの使い分け 正規累積モデルとロジスティックモデルのどちらを使用するかは,どちらでも問題ありませ ん(ほぼ同じ関数なので)*6 。コンピュータが今ほど発達していなかった昔では,計算にかかる とで,わかりやすい 1.7 が用いられています。 *6 正規累積モデルはプロビットモデル,ロジスティックモデルはロジットモデルに対応します。…と考えると「本 当は理論的にはどちらを使ったほうが良いか」が明確なケースもあるのかもしれません。実用上はどちらでも良 いとされていますが,本当にそのまま 1.7 倍の変換をするだけで交換可能なのか,というところについてはやや 議論があるようです(Cho, 2023)。
8.2 基本的なモデル 11 時間がクリティカルな問題であったなどの理由でロジスティックモデルのほうが良かったかもし れませんが,現代のコンピュータの性能であればその差は無視してよいレベルでしょう。とはい え,今でも計算の容易さからロジスティックモデルを使用するケースのほうが多いような気がし ます。 パラメータ数に関しては,考え方の違いによるところが大きい気がします。ざっくりと 1 パラ メータモデルと 2 パラメータモデルを比べると, • 1 パラメータモデルのほうがサンプルサイズが少なくても推定が安定する • 2 パラメータモデルのほうが推定の正確さ(モデル適合度)が高くなりやすい といえます。また,項目特性曲線的に考えてみると,識別力の異なる 2 つの項目間に関して, 「特性値が低い人では項目 A のほうが当てはまりやすい一方で,特性値が高い人では項目 B のほ うが当てはまりやすい」といった逆転現象のようなものが起こります(図 8.8) 。困難度パラメー タは「その項目の難しさ」を表すと言ったものの,識別力が異なる 2 パラメータモデルではシン プルな「難しさ」として考えるのは案外難しいのかもしれません。 それに対して 1 パラメータモデルでは,項目特性曲線どうしが交差することは決してない (図 8.6 )ので,どの項目間でも当てはまりやすさの大小関係が変わることはありません。した がって,困難度パラメータを文字通りの「難しさ」として捉えることができます。 …という感じでいろいろ御託を並べてみましたが,基本的には迷ったら 2 パラメータ(ロジス ティック)モデルを選んでおけばとりあえずは安泰です。後述しますが,2 パラメータモデルは 因子分析モデルと数学的に等価である一方で,1 パラメータモデルは「因子負荷をすべて 1 に固 定した因子分析モデル」に相当します。そう考えると「2 パラメータモデルでいっか」という気 持ちになりますね。 INFO パラメータの多いモデル 実は世の中にはもっとパラメータ数の多いモデルが存在します。 【3 パラメータ】 3 パラメータモデルでは,2 パラメータモデルに加えて「当て推量」のパ ラメータ 𝑐𝑖 が追加されます。 𝑃 (𝑦𝑝𝑖 = 1) = 𝑐𝑖 + 1 − 𝑐𝑖 1 − exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.9) 𝑐𝑖 は項目特性曲線の下限を操作します。𝜃𝑝 がどんなに低い人でも,必ず 𝑃 (𝑦𝑝𝑖 = 1) は 𝑐𝑖 以上の値になります。例えば 4 択問題であれば,どんなに能力が低くても,勘 で答えれば正答率は 1/4 になります。これが 𝑐𝑖 です。 【4 パラメータ】 4 パラメータモデルでは,3 パラメータモデルに加えて「上限」のパラ メータ 𝑑𝑖 が追加されます。 𝑃 (𝑦𝑝𝑖 = 1) = 𝑐𝑖 + 𝑑𝑖 − 𝑐𝑖 1 − exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.10)
8.3 因子分析との関係 12 4 パラメータモデルでは,𝜃𝑝 がどんなに高い人でも,必ず 𝑃 (𝑦𝑝𝑖 = 1) は 𝑑𝑖 以下の 値になります。どんなに能力がある人でも 100% にはならないもの,例えばバスケ のフリースローなどをイメージすると良いかもしれません。 【5 パラメータ】 5 パラメータモデルでは,4 パラメータモデルに加えて「非対称性」のパ ラメータ 𝑒𝑖 が追加されます。 𝑃 (𝑦𝑝𝑖 = 1) = 𝑐𝑖 + 𝑑𝑖 − 𝑐𝑖 𝑒 (1 − exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )]) 𝑖 (8.11) 通常の IRT モデルでは,𝑃 (𝑦𝑝𝑖 = 1) の動き方は 0.5 を軸に対称になっています。こ れに対して,例えば「最初は急激に成功率が上がるが,上に行くほどきつくなる」的 なことを表現したい場合には非対称性を考慮する必要があるでしょう。 ただ,これらのモデルは安定した推定に必要なサンプルサイズが莫大であったり,そもそ も推定量に一致性がないなどの問題を抱えていることから,実務場面で使われることはほ とんどありません(強いて言えば,当て推量パラメータを「選択肢数分の 1」に固定した 3 パラメータモデルなどはありえる)。ので,「そういうのもあるんだなぁ」くらいに思って おけば良いと思います。 8.3 因子分析との関係 IRT の基本的なモデルが分かったところで,IRT モデルと因子分析の特定のモデルが数学的に は等価であること(Takane & de Leeuw, 1987)を紹介しておきたいと思います。このことは, 心理尺度の分析において因子分析だけでなくときには IRT も使用できることを意味しています。 それぞれの分析手法は出自が異なる以上,分析の焦点や技法にも違いがあるのですが,適切な設 定のもとでは,因子分析と IRT のいいとこ取り(あるいは併用)ができるということは理解し ておくと良いかもしれません。 Chapter 7 の 多母集団同時分析のところで説明したように,標準化していない(=切片が 0 で ない)因子分析モデル(1 因子モデル)は 𝑦𝑝𝑖 = 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 − 𝜀𝑝𝑖 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀 ) (8.12) という形になっています*7 。そしてカテゴリカル因子分析では,観測変数 𝑦𝑝𝑖 はその背後にある 連続量によって決まる,という考え方をしていました。心理尺度の項目が 2 件法だとしたら,観 測変数 𝑦𝑝𝑖 とその背後の連続量の関係は図 8.3 と同じようにして 図 8.11 のように表すことがで きます。 したがって,標準正規分布に従う誤差 𝜀𝑝𝑖 が 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 より小さいときに 𝑦𝑝𝑖 = 1 となるので, *7 IRT の説明に合わせて誤差の 𝜀 𝑝𝑖 の符号を変えています。
8.4 R で IRT 13 N(t i + fpbi, se ) P(ypi = 0) 0 P(ypi = 1) t i + f pb i 図 8.11: 因子分析モデルの背後に その確率は 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 (𝜀𝑝𝑖 ≤ 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 ) (𝜏𝑖 +𝑓𝑝 𝑏𝑖 ) =∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.13) となります。これを 2 パラメータ正規累積モデル(8.6 式)に対応させて見ると, 𝜏𝑖 ) 𝑏𝑖 = 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 ) 𝜏𝑖 + 𝑓𝑝 𝑏𝑖 = 𝑏𝑖 (𝑓𝑝 + (8.14) ということで,(𝑓𝑝 , 𝑏𝑖 , 𝜏𝑖 ) = (𝜃𝑝 , 𝛼𝑖 , 𝛼𝑖 𝛽𝑖 ) としてみればこれら 2 つのモデルは完全に同一だとい うことがわかります。 IRT では,(8.14) 式に基づく(因子分析的な)パラメータ化を行うことがあります。つまり 2 パラメータロジスティックモデルで言えば 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] 1 = 1 + exp [−(𝛼𝑖 𝜃𝑝 − 𝜏𝑖 )] 𝑃 (𝑦𝑝𝑖 = 1) = (8.15) という形で,困難度パラメータの代わりに切片パラメータを求めるということです。これもどち らの設定を使用しても良いのですが,ソフトウェアによってはこの切片パラメータを使用する方 法がデフォルトになっていることもあるので,出力を見る際は気をつけましょう。なお IRT の 枠組みでは,どちらかというと困難度パラメータ 𝛽𝑖 を用いる解釈のほうがよく使われます。推 定結果で切片パラメータが出てきた場合には,自分で 𝛽𝑖 = 𝛼𝑖 と変換する必要があるかもしれま 𝜏 𝑖 せん。 8.4 R で IRT IRT の基本を学んだところで,いよいよ実際に R で IRT を実行してみます。これまでの説明 に合わせて,まずは二値で試してみるため,それ用のデータを作成します。元々の dat は 6 件法 への回答データでしたが,各選択肢は 1 から 3 が「当てはまらない」寄り,4 から 6 が「当ては まる」寄りでした。そこで,「どちらかと言えば当てはまると回答したかどうか」の二値データ
8.4 R で IRT
に変換したものを使用していきます。
二値データの作成
1
2
3
4
5
6
7
# 使用する部分だけ取り出す
# 細かいことは考えずに10項目をそのまま使います
dat_bin <- dat[, cols]
# 1から3は0,4から6は1に変換して二値化
dat_bin <- (dat_bin >= 4)
# as.numericで数字に変換
dat_bin <- apply(dat_bin, 2, as.numeric)
そして IRT での分析には,mirt パッケージの中にある mirt() 関数を使用します。基本的な
引数は以下の 3 つです。他にも色々ありますが,後ほどもう少し紹介します。
data 推定に用いるデータです。lavaan の関数と異なり,推定に使う変数のみが入ったデータ
を用意する必要があります。
model
lavaan のように各変数と因子の関係を記述したモデル式です。先程説明したように,
IRT は因子分析と同じものなので,モデル式を記述することで多次元の IRT モデルも実
行可能になる,というわけです(詳細は後ほど)。変数名が連番になっている場合,後述
のようにハイフンでまとめて指定できるのでおすすめです。
itemtype モデルの指定です。たぶん正規累積モデルは選べず,すべてロジスティック関数ベー
スのモデルです。これまでに紹介したモデルでは
• 'Rasch':ラッシュモデル,つまり 1 パラメータモデルです。後述のものに合わせて
'1PL' と指定することはできないので気をつけてください。
• '2PL':2 パラメータモデルです。同じように '3PL', '4PL' も指定可能です。
R ではじめての項目反応理論
1
2
3
4
5
6
7
# モデル式の指定
# lavaanとは違い,ふつうにイコールで結びます
model <- "f_1 = Q1_1-Q1_10"
# ハイフンの指定は以下の書き方の簡略化
# model <- 'f_1 = Q1_1 + Q1_2 + Q1_3 + (中略) + Q1_10'
result_2plm <- mirt(dat_bin, model = model, itemtype = "2PL")
得られた推定値を見る場合には,coef() 関数を使います。
推定値を見る
1
coef(result_2plm)
14
8.4 R で IRT 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 15 $Q1_1 a1 d g u par 0.416 1.262 0 1 $Q1_2 a1 d g u par 1.023 2.359 0 1 $Q1_3 a1 d g u par 1.027 1.911 0 1 $Q1_4 a1 d g u par 1.033 1.737 0 1 $Q1_5 a1 d g u par 1.01 1.791 0 1 $Q1_6 a1 d g u par 1.203 1.911 0 1 $Q1_7 a1 d g u par 1.367 1.667 0 1 $Q1_8 a1 d g u par 1.225 1.582 0 1 $Q1_9 a1 d g u par 1.339 1.365 0 1 $Q1_10 a1 d g u par 1.27 0.031 0 1 $GroupPars MEAN_1 COV_11 par 0 1 デフォルトでは,1 項目ごとの list 型で返すため少し見にくく,また今後の取り扱いがやや 面倒になりがちです。引数 simplify = TRUE とすることで,データフレームの形で返してくれ
8.4 R で IRT 16 るようになります。 推定値をシンプルに見る 1 1 2 3 4 5 6 7 8 9 10 11 12 13 coef(result_2plm, simplify = TRUE) $items a1 d g u Q1_1 0.416 1.262 0 1 Q1_2 1.023 2.359 0 1 Q1_3 1.027 1.911 0 1 Q1_4 1.033 1.737 0 1 Q1_5 1.010 1.791 0 1 Q1_6 1.203 1.911 0 1 Q1_7 1.367 1.667 0 1 Q1_8 1.225 1.582 0 1 Q1_9 1.339 1.365 0 1 Q1_10 1.270 0.031 0 1 14 15 16 17 18 19 20 $means f_1 0 $cov f_1 f_1 1 最初の $item 部分に表示されているものが項目パラメータです。左から a1 識別力パラメータの推定値 d 切片パラメータの推定値 g 当て推量パラメータ (guessing) の推定値(2 パラメータモデルなら全て 0) u 上限パラメータ (upper bound) の推定値(2 パラメータモデルなら全て 1) となっています。このように,mirt はデフォルトでは切片パラメータを出力します。困難度 パラメータを見たい場合には d 列を a1 列で割るか,引数 IRTpars = TRUE をつけて出力しま しょう。 困難度パラメータを見る 1 1 2 3 coef(result_2plm, simplify = TRUE, IRTpars = TRUE) $items Q1_1 a b g u 0.416 -3.031 0 1
8.4 R で IRT 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 17 Q1_2 1.023 -2.307 0 1 Q1_3 1.027 -1.861 0 1 Q1_4 1.033 -1.681 0 1 Q1_5 1.010 -1.773 0 1 Q1_6 1.203 -1.588 0 1 Q1_7 1.367 -1.220 0 1 Q1_8 1.225 -1.292 0 1 Q1_9 1.339 -1.019 0 1 Q1_10 1.270 -0.025 0 1 $means f_1 0 $cov f_1 f_1 1 これで d 列の代わりに b(困難度パラメータ)列が表示されるようになりました。出力を見る と,例えば Q1_1 では,1.262/0.416 ≃ 3.031 ということで,確かに d 列を a1 列で割った値(の 符号を入れ替えたもの)が表示されています。 $item の下に表示されている $means と $cov はそれぞれ特性値(因子得点)𝜃 の平均値と分 散です。この時点では一般的な因子分析の制約(平均 0,分散 1)に沿った値のみが表示されて いるだけなので特に注目する必要はありません。後ほど紹介する多次元モデルや多母集団モデル になると,ここに表示される値も重要な意味を持つことになります。 plot() 関数を使うと,引数 type に応じて,推定結果に基づく様々なプロットを全項目まと めて出すことができます。図 8.6 のような ICC を出したい場合には type = "trace"としてあ げてください。 項目特性曲線 1 plot(result_2plm, type = "trace") デフォルトではすべての項目の ICC を縦横に並べてくれるのですが,項目数が多くなると見 るのが大変そうです。この場合,引数 facet_items = FALSE とすると,すべての項目のプロッ トを重ねて表示してくれます。また引数 which.items を指定すると,特定の項目だけの ICC を 出すことも出来ます。 特定の項目の ICC だけを重ねて表示 1 plot(result_2plm, type = "trace", which.items = 1:5, facet_items = FALSE) ということで 図 8.12 や 図 8.13 を見ると,識別力がダントツで低い項目 1 では,ICC の傾き が緩やかになっていることがわかります。
8.4 R で IRT 18 Item Probability Functions -6 -4 -2 0 2 4 6 Q1_10 Q1_5 Q1_6 Q1_7 Q1_8 P (q ) 0.8 0.4 0.0 Q1_9 Q1_1 Q1_2 Q1_3 Q1_4 0.8 0.4 0.0 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 q 図 8.12: 項目特性曲線 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 項目の項目特性曲線 0.8 0.4 0.0
8.5 IRT の前提条件 19 IRT による項目分析では,まずこのように推定されたパラメータをもとに ICC を描くのが重 要です。これは Chapter 4 で出てきたトレースラインの図と同じように解釈することができるの で,例えば • 𝜃 の広い領域で ICC がほぼ上か下に貼り付いている(𝛽𝑖 がとても小さいか大きい) • ICC の傾きが緩やかである(𝛼𝑖 がとても小さい) ような項目は,𝜃 の測定という面からはあまり良くない項目であろう,ということがわかりま す(除外するかは要検討)。この他にも,複数の項目を組み合わせたときにどうなるかという分 析などもあるわけですが,それはおいおい紹介していきます。 話は変わって,項目パラメータを固定した上での特性値 𝜃𝑝 の推定値を出す場合は,fscores() 関数を使います。このあたりは因子分析や CFA での「項目パラメータを推定」→「項目パラメー タを固定して個人特性値を推定」と全く同じです。 theta の推定値を出す 1 head(fscores(result_2plm)) 1 2 3 4 5 6 7 f_1 [1,] -1.537468218 [2,] -0.007789283 [3,] 0.226340769 [4,] -0.643698451 [5,] 0.085994618 [6,] 0.788036256 8.5 IRT の前提条件 IRT,なかでも 2 パラメータモデルは基本的には因子分析と同じようなものなのですが,理論 的背景の違いなどから,因子分析を行う際とは留意すべき点も少し異なってきます。ここでは, IRT を適用するにあたって重要となる 2 つの前提条件について見ていきます。 8.5.1 局所独立性 局所独立性 (local independence) の仮定は,因子分析のときと概ね同じです。探索的因子 分析では,独自因子の間が完全に無相関であるという仮定をおき,その中でデータの相関行列を 最もよく復元できる因子負荷行列と独自因子の分散 (𝐁⊤ 𝐁 + 𝚿) を求めていました。独自因子の 間に相関があるということは,その当該項目の間に何か別の共通因子があるということです。 局所独立性をもう少し固い言い方で説明すると, 「𝜃𝑝 で条件づけたとき,項目 𝑖 と項目 𝑗 (𝑖 ≠ 𝑗) への回答は完全に独立」となります。項目への回答は二値なので,2 × 2 クロス表をイメージして みましょう。𝜃𝑝 が特定の値の人を 1000 人ほど集めて(つまり 𝜃𝑝 で条件づけて)表 8.1 のよう に 2 つの項目 𝑖, 𝑗 の回答のクロス表を作ると,もし局所独立性が満たされているならば連関はゼ
8.5 IRT の前提条件 20 表 8.1: 局所独立な 2 項目 項目 𝑗 項目 𝑖 当てはまらない 当てはまる 計 当てはまらない 180 420 600 当てはまる 120 280 400 計 300 700 1000 表 8.2: 局所独立ではない 2 項目 項目 𝑗 項目 𝑖 当てはまらない 当てはまる 計 当てはまらない 240 60 600 当てはまる 360 340 400 計 300 700 1000 ロになります。項目 𝑖 に「当てはまる」と回答した人を見ても, 「当てはまらない」と回答した人 を見ても,項目 𝑗 に「当てはまる」と回答した人の割合は全体の割合と同じ 6:4 になっています。 一方で,𝜃𝑝 以外の別の共通因子がある場合,クロス表に連関が現れます。例えば「他者への助 けを求める傾向」の尺度において,宗教観の影響を受けると思われる項目(e.g., 困ったときは神 に祈る)どうしのクロス表を見てみると,表 8.2 のように • 項目 𝑖 に「当てはまる」と答えた人では 340/400 = 85% の人が項目 𝑗 でも「当てはまる」 と回答 • 項目 𝑖 に「当てはまらない」と答えた人では 360/600 = 60% の人が項目 𝑗 でも「当ては まる」と回答 しているため,項目 𝑖 の回答と項目 𝑗 の回答が関連を持っていることになります。項目 𝑖 に 「当てはまる」と答えた人のほうが項目 𝑗 でも「当てはまる」と回答する確率が高い,ということ です。 局所独立の仮定とは,このようにクロス表を作成したときに,どの項目ペア間でも,どの特性 値の値でも連関がない,ということを指しています*8 。 局所独立性がなぜ重要なのでしょうか。その理由の一つは計算上の問題です。IRT では基本的 に最尤法によってパラメータを推定します。もしも局所独立性が満たされていれば,𝜃𝑝 で条件づ けたときの各項目に対する反応確率は独立なので,回答者 𝑝 に対する尤度関数は単純にすべての *8 厳密には, 「独立」と「無相関」は別の概念です。理想は独立なのですが,実用上は無相関であれば良いだろう, という感じになっています。厳密な独立とは異なるという意味で,弱局所独立性と呼ぶこともあります。
8.5 IRT の前提条件 21 項目をかけ合わせた 𝐼 𝑃 (𝐲𝑝 |𝜃𝑝 ) = ∏ 𝑃 (𝑦𝑝𝑖 = 1|𝜃𝑝 )𝑦𝑝𝑖 𝑃 (𝑦𝑝𝑖 = 0|𝜃𝑝 )1−𝑦𝑝𝑖 (8.16) 𝑖=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 では,局所独立性の観点から一次元性を確認する指標もいくつか提案されています。 一次元性が満たされているということは,回答行動に影響する共通因子が一つ(𝜃𝑝 )だけという ことです。したがって一次元性が満たされていない場合には,Section 6.3.1(図 6.6)で少し説 明したように,共通因子で説明出来ない部分(独自因子)の間に共分散が発生することになりま す。そのように考えると,一次元性は局所独立と同じようなものだとみなすことができるわけで す。厳密には局所独立性が満たされている場合,一次元性も満たされていると言って良いだろ う,ということになります*11 。ということで,より洗練された一次元性確認の方法として,例え *9 「𝜃 で条件づけた」ということが明確になるようにここだけ 𝑃 (.|𝜃 ) としていますが,これまでの 𝑃 (.) と意 𝑝 𝑝 味は同じです。 *10 後半では少しだけ多次元モデルの紹介をしますが,これは一次元モデルの拡張として展開されるものであり, IRT の基本は一次元だといえます。 *11 もちろん,多因子モデルの場合には,一次元とは限りません。ただ同様に,局所独立性が満たされているならば
8.6 多値型モデル 22 ば DIMTEST(Stout et al., 1996)や DETECT(Zhang & Stout, 1999)という名前の方法が 提案されていたりします。 INFO 局所独立と一次元性 一次元性が満たされているだけでは局所独立性が満たされるとは限りません。例えば • 変数 𝑥 の平均値を求めよ • 変数 𝑥 の分散を求めよ という 2 つの問題があった場合,この 2 つはいずれも「統計の知識」を問うているという 意味では同じ潜在特性に関する項目であり,一次元性を満たしているといえます。ですが, 分散を計算する過程では必ず平均値を計算する必要があるため,項目 1 に間違えた人はほ ぼ確実に項目 2 に間違えます。したがって,この 2 問は一次元性は満たされているが局 所独立性は満たされていない状態だといえます。このように,ある項目への回答が他の項 目への回答に直接的に影響を与える場合を実験的独立性の欠如と呼びます。心理尺度でい えば, • 前の項目で『どちらかといえば当てはまる』以上を選んだ方にお聞きします。 • 前の項目で『どちらかといえば当てはまらない』以下を選んだ方は,この項目 では『全く当てはまらない』を選んでください。 といった指示を出した場合には,実験的独立性が欠如していると言えそうです。他には, 回答時間が足りない場合の後ろの方の項目も, 「前の項目に回答出来ていない人は次の項目 にも回答できていない」ため,実験的独立性が欠如しているとみなすことが出来ます。 ということで局所独立性を正確に表すと,一次元性に加えてこの実験的独立性が必要とい えます(南風原,2000; Lord & Novick, 1968) 。 8.6 多値型モデル ここまで,二値の項目反応データに対する基本的なモデルを紹介し,R での実行方法を示しま した。しかし実際の場面では,リッカート式尺度ならばたいてい 5 件法から 7 件法の多値項目で す。そしてテストの場面でも,部分点や組問*12 を考えると多値型のモデルが必要となるでしょ う。ということで,代表的な 2 つの多値型モデルを紹介したいと思います。どちらを使っても本 質的には違いは無いので,説明を読んでしっくり来たほうを使えば良いと思います。 𝑛 次元性が満たされている,ということはできると思います。 *12 問 1 の答えを使わないと問 2 が解けないような場合,局所独立性が満たされなくなってしまいます。そこで「問 1 が解けたら 1 点,問 2 も解けたら 2 点」というように得点化してこれを一つの大きな問題として考えること で局所独立性を確保しようという考えに至るわけです。
8.6 多値型モデル 23 8.6.1 段階反応モデル 段階反応モデル(Graded response model [GRM]: Samejima, 1969)では,複数の二値 IRT モデルを組み合わせて多値反応を表現します。項目 𝑖 に対する回答者 𝑝 の回答 𝑦𝑝𝑖 = 𝑘 (𝑘 = 1, 2, ⋯ , 𝐾) について「𝑘 以上のカテゴリを選ぶ確率」を考えると,これはまだ二値(𝑘 以上 or 𝑘 より小さい)なので,2 パラメータロジスティックモデルならば 𝑃 (𝑦𝑝𝑖 ≥ 𝑘) = 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] (8.17) と表せます。なお,困難度パラメータは「項目 𝑖 のカテゴリ 𝑘」ごとに用意されるため 𝛽𝑖𝑘 とし ています。 異なるカテゴリを選ぶ確率は当然排反な事象なので,二値モデルを組み合わせると「ちょうど 𝑘 番目のカテゴリを選ぶ確率」は 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 ≥ 𝑘) − 𝑃 (𝑦𝑝𝑖 ≥ 𝑘 + 1) (8.18) と表すことができます。ただし,端のカテゴリに関する確率はそれぞれ 𝑃 (𝑦𝑝𝑖 ≥ 1) = 1 𝑃 (𝑦𝑝𝑖 ≥ 𝐾 + 1) = 0 (8.19) とします。「1 以上のカテゴリを選ぶ確率」は当然 1 になるため,これに対する困難度パラメー タ 𝛽𝑖1 = −∞ となります。同様に,「𝐾 + 1 以上のカテゴリを選ぶ確率」は 0 になります。 𝑃 (𝑦𝑝𝑖 = 𝐾|𝜃𝑝 ) を考えるときだけ仮想的な「さらに上のカテゴリ」を想定する,ということです。 このように GRM は二値型モデルの素直な拡張のため,数学的には正規累積型の GRM はカテゴ リカル因子分析と同じものだとみなすことができます(Takane & de Leeuw, 1987)。 図 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]) と呼ぶこともあります。が,そこまで細かい名前を覚える必要はないでしょう。
8.6 多値型モデル 24 P(ypi ³ k) P(ypi ³ k) 0.5 qp qp - 1.5 0 0 2 図 8.14: 段階反応モデル P(ypi = k) P(ypi = k) (ai, b i2, b i3, b i4) (1, - 1.5, 0, 2) 1 3 (ai, b i2, b i3, b i4) (1, - 0.5, 1, 3) 4 1 2 4 3 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: GRM の項目特性曲線 1
8.6 多値型モデル 25 左上のプロットは,図 8.15 と同じ項目パラメータ((𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4 ]) = (1, [−1.5, 0, 2]))で の IRCCC です。基本的には,このように 𝜃𝑝 の値に応じて左から順に各カテゴリの選択確率の ピークが訪れます。点線は中間カテゴリの選択確率が最も高くなる箇所ですが,これは具体的に カテゴリ困難度 𝛽𝑖𝑘 と 𝛽𝑖,𝑘+1 のちょうど中間になります。 右 上 の プ ロ ッ ト は, 左 上 か ら す べ て の カ テ ゴ リ 困 難 度 の 値 を 1 ず つ 大 き く し ま し た ((𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4 ]) = (1, [−0.5, 1, 3])) 。二値モデルの時と同じように,すべてのカテゴリ困難 度が同じ値だけ変化するときには IRCCC は平行移動します。 左下のプロットは,左上から識別力を半分にしました((𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4 ]) = (0.5, [−1.5, 0, 2])) 。 ま た 右 下 の プ ロ ッ ト は, 左 上 か ら 識 別 力 を 2 倍 に し て い ま す ((𝛼𝑖 , [𝛽𝑖2 , 𝛽𝑖3 , 𝛽𝑖4] ) = (2, [−1.5, 0, 2]))。識別力は各 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.20) exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] 1 + exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] となります。つまりこの式は,隣のカテゴリとの二択において,当該カテゴリを選ぶ確率という 見方ができるわけです。これを多値 (𝑘 = 1, 2, ⋯ , 𝐾) に置き換えると,「カテゴリ 𝑘 − 1 と 𝑘 の 二択だったら 𝑘 の方を選ぶ確率」に拡張することができそうです。ということで,GPCM では
8.6 多値型モデル 26 ∗ この確率 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 を 𝑃 (𝑦𝑝𝑖 = 𝑘) 𝑃 (𝑦𝑝𝑖 = 𝑘) + 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) 1 = 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] ∗ 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 = = = exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] (8.21) 1 + exp [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 )] exp [𝜋𝑝𝑖𝑘 ] 1 + exp [𝜋𝑝𝑖𝑘 ] と表します。なお,この後の式展開を考えて 𝜋𝑝𝑖𝑘 = 𝛼𝑖 (𝜃𝑝 − 𝛽𝑖𝑘 ) と置き換えています。 ここから,全カテゴリの中でカテゴリ 𝑘 を選ぶ確率を導出しましょう。(8.21) 式を 𝑃 (𝑦𝑝𝑖 = 𝑘) = の形に変形すると, 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) ∗ 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 ∗ 1 − 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 (8.22) 𝑃 となります。最後の右辺の 1−𝑃 の部分は 2 つのカテゴリ間のオッズを表しているため, を用いて ∗ 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 = exp(𝜋𝑝𝑖𝑘 ) ∗ 1 − 𝑃 (𝑦𝑝𝑖 = 𝑘)𝑘−1,𝑘 (8.23) 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) exp(𝜋𝑝𝑖𝑘 ) (8.24) という漸化式が得られ,特定のカテゴリ 𝑘 を選ぶ確率が表現できるようになります。例えば 4 カ テゴリの場合には, 𝑃 (𝑦𝑝𝑖 = 1) = 𝑃 (𝑦𝑝𝑖 = 1) 𝑃 (𝑦𝑝𝑖 = 2) = 𝑃 (𝑦𝑝𝑖 = 1) exp(𝜋𝑝𝑖2 ) 𝑃 (𝑦𝑝𝑖 = 3) = 𝑃 (𝑦𝑝𝑖 = 1) exp(𝜋𝑝𝑖2 ) exp(𝜋𝑝𝑖3 ) (8.25) 𝑃 (𝑦𝑝𝑖 = 4) = 𝑃 (𝑦𝑝𝑖 = 1) exp(𝜋𝑝𝑖2 ) exp(𝜋𝑝𝑖3 ) exp(𝜋𝑝𝑖4 ) という形で 𝑃 (𝑦𝑝𝑖 = 1) を起点に「一個上のカテゴリに上がるときのオッズ」の積によってすべ てのカテゴリの反応確率が表されるわけです。ということで,この式を一般化して 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 1) exp(𝜋𝑝𝑖𝑘 ) = 𝑃 (𝑦𝑝𝑖 = 𝑘 − 2) exp(𝜋𝑝𝑖,𝑘−1 ) exp(𝜋𝑝𝑖𝑘 ) ⋮ 𝑘 = 𝑃 (𝑦𝑝𝑖 = 1) ∏ exp(𝜋𝑝𝑖𝑐 ) (8.26) 𝑐=2 𝑘 = 𝑃 (𝑦𝑝𝑖 = 1) exp (∑ 𝜋𝑝𝑖𝑐 ) 𝑐=2 と書き下します。ここではまだ 𝑃 (𝑦𝑝𝑖 = 1) が残っていますが,GPCM の考え方ではこれ以上
8.6 多値型モデル 27 計算できません*13 。改めて (8.26) 式を見ると,すべてのカテゴリの反応確率は「𝑃 (𝑦𝑝𝑖 = 1) の何倍」という形になっていることがわかります。そこで一旦 𝜋𝑝𝑖1 = 0 として得られる 𝑃 (𝑦𝑝𝑖 = 1) = exp(𝜋𝑝𝑖1 ) = 1 を (8.26) 式に代入して, 𝑘 𝑃 (𝑦𝑝𝑖 = 𝑘) = exp (∑ 𝜋𝑝𝑖𝑐 ) (8.27) 𝑐=1 とまとめてしまいます。これは,各カテゴリの反応確率を「一番下のカテゴリが選択される確率 の何倍か」という形で表すことを意味しています。これによって,ひとまず 𝑃 (𝑦𝑝𝑖 = 𝑘) が定式 化されましたが,𝑃 (𝑦𝑝𝑖 = 1) = 1 と設定している以上 𝑃 (𝑦𝑝𝑖 = 𝑘) の全カテゴリの総和が 1 にな るかわかりません。 ということで最後に,𝑃 (𝑦𝑝𝑖 = 1) の値が何であっても全カテゴリでの総和が 1 になるように するため,総和で割ることを考えます。実際に 𝑃 (𝑦𝑝𝑖 = 𝑘) 𝑃 (𝑦𝑝𝑖 = 𝑘) = 𝐾 𝑃 (𝑦𝑝𝑖 = 1) + 𝑃 (𝑦𝑝𝑖 = 2) + ⋯ + 𝑃 (𝑦𝑝𝑖 = 𝐾) ∑𝑙=1 𝑃 (𝑦𝑝𝑖 = 𝑙) 𝑘 exp (∑𝑐=1 𝜋𝑝𝑖𝑐 ) = exp (∑1𝑙=1 𝜋𝑝𝑖𝑙 ) + exp (∑2𝑙=1 𝜋𝑝𝑖𝑙 ) + ⋯ + exp (∑𝐾 𝜋 ) 𝑙=1 𝑝𝑖𝑙 (8.28) exp (∑𝑐=1 𝜋𝑝𝑖𝑐 ) 𝑘 = ∑𝑙=1 exp (∑𝑐=1 𝜋𝑝𝑖𝑙 ) 𝐾 𝑙 とすると,分子も分母もすべての項に共通の exp(𝜋𝑝𝑖1 ) がかかっているため,この項の値が何で あっても全体の割合に占める各カテゴリの反応確率の割合は変わらずに済みます。ということ で,ようやく GPCM の式が導出できました。図 8.16 に,ここまでの考え方をまとめてみまし た。特定のカテゴリの反応確率 𝑃 (𝑦𝑝𝑖 = 𝑘) が直接計算できないため,いったん一番下のカテゴ リの反応確率 𝑃 (𝑦𝑝𝑖 = 1) = 1 とおいて,他のカテゴリは「前のカテゴリの何倍」という形で表 します。最後に全体の総和が 1 になるように調整をしたものが GPCM における反応確率という わけです。 図 8.17 は,GPCM での項目特性曲線をプロットしたものです。GPCM では,カテゴリ困難 度 𝛽𝑖𝑘 は「カテゴリ 𝑘 − 1 と 𝑘 の二択」における選択確率がちょうど 50% になるポイントを表 しています。したがって,隣接するカテゴリの反応確率はちょうど 𝛽𝑖𝑘 のところで大小関係が入 れ替わります。𝛽𝑖𝑘 はそのように解釈できるわけです。 GRM では,カテゴリ困難度は必ず単調増加になっている必要がありますが,GPCM では単 調増加で無くても問題ありません。図 8.18 にカテゴリ困難度が単調増加にならないケースをプ ロットしてみました。この場合,「カテゴリ 2 がカテゴリ 1 を上回る」よりも先に「カテゴリ 3 がカテゴリ 2 を上回る」ため,結果的にカテゴリ 2 が最大になることがありません。 *13 𝑃 (𝑦 = 1) = exp(𝜋 𝑝𝑖 𝑝𝑖1 )𝑃 (𝑦𝑝𝑖 = 0) と表せますが,カテゴリ 𝑘 = 0 は存在しない,つまり 𝑃 (𝑦𝑝𝑖 = 0) = 0 ∗ = 1) です。ということは 𝑃 (𝑦𝑝𝑖 0,1 = ∞ となるため,exp(𝜋𝑝𝑖1 ) = ∞ となってしまい計算ができないのです。
8.6 多値型モデル 28 図 8.16: GPCM のイメージ 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 の項目特性曲線
8.6 多値型モデル 29 P(ypi = k) (ai, b i2, b i3, b i4) (1, 0, - 1.5, 2) 1 3 4 2 qp - 1.5 0 2 図 8.18: カテゴリ困難度が単調増加にならないケース 8.6.3 R で多値型モデル それでは,mirt で 2 つの多値型モデルを実行してみましょう。…といっても引数 itemtype に graded(GRM) か gpcm を指定したら良いだけなので難しいことはありません。 多値型モデル (GRM) 1 # 多値型なのでdatをそのまま突っ込む 2 result_grm <- mirt(dat[, cols], model = "f_1 = Q1_1-Q1_10", itemtype = "grade d") 係数をチェック 1 1 2 3 4 5 6 7 8 9 10 11 coef(result_grm, IRTpars = TRUE, simplify = TRUE) $items Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 a b1 b2 b3 b4 b5 0.572 -6.382 -3.936 -2.323 -1.040 1.269 1.181 -4.051 -2.765 -2.087 -0.840 0.821 1.238 -3.245 -2.243 -1.643 -0.564 1.019 1.074 -3.293 -2.174 -1.648 -0.706 0.415 1.116 -3.951 -2.472 -1.664 -0.460 1.216 1.062 -3.988 -2.698 -1.726 -0.440 1.429 1.177 -3.454 -2.084 -1.314 -0.178 1.485 1.097 -3.671 -2.186 -1.386 -0.023 1.744 1.365 -3.306 -2.016 -1.039 -0.312 0.921
8.6 多値型モデル 12 13 14 15 16 17 18 19 20 30 Q1_10 1.200 -2.226 -1.048 -0.056 0.461 1.586 $means f_1 0 $cov f_1 f_1 1 dat の各項目は 6 件法なので,カテゴリ困難度 b が 5 つ出力されました。 二値型モデルの時と同じように,項目特性曲線を出力することももちろん可能です。 GRM で項目特性曲線 1 2 # カテゴリ数のぶんだけ線を引いてくれる plot(result_grm, type = "trace") Item Probability Functions -6 -2 0 2 4 6 Q1_10 Q1_5 Q1_6 Q1_7 Q1_8 P (q ) 0.8 0.4 0.0 Q1_9 Q1_1 Q1_2 Q1_3 Q1_4 0.8 0.4 0.0 P1 P2 P3 P4 P5 P6 0.8 0.4 0.0 -6 -2 0 2 4 6 -6 -2 0 2 4 6 q 図 8.19: GRM で項目特性曲線 また,itemplot(., type = "threshold") 関数を使うと,図 8.14 の下半分のような「カテ ゴリ 𝑘 以上の反応をする確率」のプロットを出すこともできます。
8.6 多値型モデル
31
GRM で項目特性曲線 2
1
2
3
# itemplot()は1項目しかプロットできない
# 何番目の項目をプロットするかを引数itemで指定
itemplot(result_grm, item = 6, type = "threshold")
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: GRM で項目特性曲線 2
同じように GPCM でもやってみましょう。やり方も結果の見方も同じなので詳しい説明は省
略します。
多値型モデル (GPCM)
1
# 多値型なのでdatをそのまま突っ込む
2
result_gpcm <- mirt(dat[, cols], model = "f_1 = Q1_1-Q1_10", itemtype = "gpc
m")
係数をチェック
1
1
2
3
4
5
6
7
8
9
10
coef(result_gpcm, IRTpars = TRUE, simplify = TRUE)
$items
Q1_1
Q1_2
Q1_3
Q1_4
Q1_5
Q1_6
Q1_7
Q1_8
a
b1
b2
b3
b4
b5
0.222 -5.012 -2.300 -1.027 -3.288 -0.347
0.681 -2.878 -1.254 -2.419 -1.103 0.536
0.672 -2.103 -1.046 -2.050 -0.935 0.791
0.484 -2.141 -0.254 -2.357 -0.891 -0.851
0.582 -3.195 -1.219 -1.986 -0.831 0.978
0.421 -3.016 -1.949 -2.352 -1.175 1.566
0.491 -3.108 -1.031 -1.960 -0.787 1.528
0.467 -3.297 -0.963 -2.341 -0.371 1.841
8.7 多次元モデル 11 12 13 14 15 16 17 18 19 20 32 Q1_9 0.553 -3.276 -1.887 -0.518 -0.992 Q1_10 0.438 -1.844 -0.932 1.298 -1.047 0.491 0.940 $means f_1 0 $cov f_1 f_1 1 GPCM で項目特性曲線 1 2 # カテゴリ数のぶんだけ線を引いてくれる plot(result_gpcm, type = "trace") Item Probability Functions -6 -2 0 2 4 6 Q1_10 Q1_5 Q1_6 Q1_7 Q1_8 P (q ) 0.8 0.4 0.0 Q1_9 Q1_1 Q1_2 Q1_3 Q1_4 0.8 0.4 0.0 P1 P2 P3 P4 P5 P6 0.8 0.4 0.0 -6 -2 0 2 4 6 -6 -2 0 2 4 6 q 図 8.21: GPCM で項目特性曲線 8.7 多次元モデル これまでの IRT モデルでは一次元性の仮定に基づく一次元モデルのみを扱ってきました。で すが多くの心理尺度やテストでは,複数の次元の特性を同時に測定したいというニーズがありま す。先程紹介したように IRT はカテゴリカル因子分析と数学的には同値なので,IRT でも多次 元モデルを考えていきましょう。
8.7 多次元モデル 33 多次元因子分析モデルでは,項目反応を 𝑇 𝑦𝑝𝑖 = 𝜏𝑖 + ∑ 𝑓𝑝𝑡 𝑏𝑡𝑖 − 𝜀𝑝𝑖 (8.29) 𝑡=1 と表しました。これと同じように考えるならば,多次元 IRT モデル(2 パラメータロジスティッ クモデル)での項目反応関数は 𝑃 (𝑦𝑝𝑖 = 1) = 1 1 + exp [− ∑𝑇𝑡=1 𝛼𝑡𝑖 𝜃𝑝𝑡 − 𝜏𝑖 ] (8.30) と表せそうです。 (8.30) 式のモデルにおける項目パラメータは,識別力が次元の数(𝑇 個)と,切片が一つです。 また,項目の全体的な識別力を評価する際には,次元ごとの識別力を一つにまとめた多次元識別 力(multidimensional discrimination: 豊田,2013) 𝑇 MDISC𝑖 = √∑ 𝛼2𝑡𝑖 (8.31) 𝑡=1 を用いることもあります。これを用いると,切片パラメータ 𝜏𝑖 を困難度に変換することもできる ようになります: 𝛽𝑖 = −𝜏𝑖 . MDISC𝑖 (8.32) 二次元までならば,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: 項目反応曲面(Reckase, 2019, fig. 12.1)
8.7 多次元モデル 34 8.7.1 補償モデルと非補償モデル ここまで紹介した多次元モデルでは 図 8.22 からも分かるように,反応確率に対して 𝜃𝑝𝑡 が及 ぼす影響は加法的といえます。例えば ([𝛼𝑖1 , 𝛼𝑖2 ], 𝜏𝑖 ) = ([1.5, 1], 0) の項目の場合,𝜽𝑝 = [−1 の人と 𝜽𝑝 = [0 1] −0.5] の人では反応確率は同じになります。また,ある次元の特性値がとても 低い場合でも,別の次元の特性値がとても高い場合には反応確率が高くなる可能性があります。 このような次元間の関係を反映した先程の多次元 IRT モデルは補償型 (compensatory) モデ ルと呼ばれます。 一方で,複数の次元の要素がすべて揃って初めて反応確率が高くなるようなケースも考えられ ます。例えば数学のテストが日本語ではなく英語で実施される場合,これは「英語」と「数学」 の 2 つの能力が試される二次元構造ではありますが,補償型モデルのように「英語が苦手でも 数学力が高ければ解答できる」とは考えにくいです。英語で書かれた問題文を理解するだけの英 語力と,計算を実行できるだけの数学力の両方が揃って初めて正答できるでしょう。このような 状況では,一方の特性値の低さをもう一方の特性値の高さでカバーすることはできません。こう した次元間の関係を反映したモデルを非補償 (noncompensatory) モデルあるいは部分補償 (partially compensatory) モデルと呼びます*14 。 最もシンプルな非補償モデルの項目反応関数は 𝑇 1 𝑡=1 1 + exp [−𝛼𝑖𝑡 (𝜃𝑝𝑡 − 𝛽𝑖𝑡 )] 𝑃 (𝑦𝑝𝑖 = 1) = ∏ (8.33) と表されます。各次元に関して独立に算出された「その次元の閾値を超える確率」を計算し,最 後にすべての積をとっています。そのため,一つでも特性値の低い次元があれば反応確率は一気 に小さくなってしまう,というわけです。非補償モデルの項目反応曲面は 図 8.23 の左のように なります。図 8.22 と比べると,両方の次元の 𝜃𝑝 が高くないと 𝑃 (𝑦𝑝𝑖 = 1) が高くならないこと が見て取れます。 補償モデルと非補償モデルを比べると,補償モデルのほうがよく用いられているような印象で す。その理由の一つとして,非補償モデルのほうがパラメータの推定が難しいという点が挙げら れると思います。上の式を見るとわかるように,非補償モデルでは困難度パラメータ 𝛽𝑖𝑡 が次元 の数だけあるためそもそものパラメータ数が補償型よりも多くなっています。また,非補償モデ ルは項目反応関数の全体の形が積になっています。コンピュータは足し算よりも掛け算のほうが 苦手なので,そういった意味でも計算の難易度が高いのです。また,心理尺度的には補償型のほ うが従来の因子分析モデルと同じ形になっているため,扱いやすいというのもあると思います。 とはいえ,多次元の考え方として非補償型もあるんだということを頭の片隅に置いておくと,い つかどこかで何かの参考になるかもしれません*15 。 *14 完全には非補償ではない(ある程度なら別の次元の特性値でリカバーできる)という意味で「部分補償モデル」 のほうがより正しい名前らしいですが,そこまでこだわるほどの違いはありません。 *15 あるかはわかりませんが,例えば「2 つの特性を両方とも持っている人だけ当てはまると回答する尺度・項目」 のようなものを考えた場合には,普通の因子分析よりも非補償型のモデルを適用したほうが良い(回答メカニズ ムに則している)可能性があります。
8.7 多次元モデル
図 8.23: 非補償モデルの項目反応曲面(Reckase, 2019, fig. 12.2)
8.7.2 R でやってみる
ここからは補償型の多次元 IRT モデルを mirt で実行する方法を紹介します*16 。補償型モデ
ルはカテゴリカル因子分析と同値なので,因子分析と同じように「探索的モデル」と「検証的モ
デル」の 2 種類があります。ただ探索的モデルでは因子の回転に応じて 𝜃𝑝𝑡 の意味が変化してし
まうため実用上はなかなか使いにくいような気もします。
探索的補償型モデル
探索的モデルを実行する場合には,引数 model に「次元数」を与えるだけです。itemtype は
これまでと同じように自由に選んでください。二値だったら 2PL とか,多値だったら graded と
か gpcm とか。
探索的補償型モデル
1
2
# GRMでやってみましょう
# このあたりまで来ると推定にもそこそこ時間がかかるかもしれません。
3
result_expl_comp <- mirt(dat[, cols], model = 2, itemtype = "graded")
推定値のチェック
1
2
# 多次元の場合引数IRTparsは使えないので注意
coef(result_expl_comp, simplify = T)
*16 非補償型モデルは,二値ならば itemtype='PC2PL' または itemtype='PC3PL' で実行可能です( PC は
partially compensatory の頭文字)。ただパラメータの推定は結構安定しないので利用する際は気をつけてくだ
さい。
35
8.7 多次元モデル 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 36 $items a1 a2 d1 d2 d3 d4 d5 Q1_1 -0.251 -0.877 3.879 2.433 1.457 0.665 -0.785 Q1_2 -0.817 -1.700 5.685 3.977 3.034 1.234 -1.208 Q1_3 -0.990 -2.403 5.789 4.099 3.032 1.048 -1.891 Q1_4 -0.775 -0.903 3.646 2.414 1.831 0.781 -0.465 Q1_5 -0.757 -1.523 5.105 3.271 2.221 0.612 -1.628 Q1_6 -1.449 0.151 4.667 3.195 2.067 0.537 -1.707 Q1_7 -1.645 0.183 4.602 2.822 1.798 0.252 -2.011 Q1_8 -1.329 0.047 4.263 2.559 1.625 0.028 -2.047 Q1_9 -1.918 0.131 5.227 3.255 1.695 0.508 -1.501 Q1_10 -1.436 0.000 2.864 1.361 0.075 -0.595 -2.040 $means F1 F2 0 0 $cov F1 F2 F1 1 0 F2 0 1 結果を見ると,識別力パラメータが次元の数 ( a1, a2) と,切片パラメータが(カテゴリ数-1) 個 ( d1-d5) 推定されています。その下の $mean と $cov には,それぞれ 𝜃𝑝𝑡 の平均と共分散行 列が出ています。これを見る限り,まだ因子間相関が 0 なので回転前の因子負荷(識別力)が出 ているようです。さらに詳しく見ると,最後の項目の a2 が 0.000 となっています。実は mirt での探索的モデルの推定の際には,このように一部の因子負荷を 0 に固定した状態で推定を行っ ています*17 。「探索的」と紹介していましたが,実際には「モデルが識別できる最小限の制約を 置いて推定する」ということを行っているのです。IRT に限らず CFA もそうなのですが,複数 の因子がある場合に識別性を持たせるには,最低でも 𝑇 (𝑇 −1) 個の因子負荷を固定する必要があ 2 ることが知られています。そのため,2 因子の場合は第 2 因子の因子負荷を一つ 0 にしているわ けです。同様に,3 因子の場合にはさらに第 3 因子の因子負荷を 2 つ 0 にする必要があり,4 因 子の場合にはさらにさらに第 4 因子の因子負荷を 3 つ 0 にする必要があり…という要領です。 このままでは,いくら探索的とは言え各次元の解釈が出来ません。mirt() で得られるパラ メータの推定値はいわば「初期解」なので,ここからさらに回転させてあげれば良さそうです。 回転させる場合は,psych::fa() の時と同じように引数 rotate を指定してあげます。 *17 psych::fa() のように固有値分解など使えばいいじゃないと思われるかもしれませんが,mirt() 関数として 他の様々な IRT モデルと同じ文法で実装するためにはこうするのが最も自然なのです。
8.7 多次元モデル 37 識別力の回転 1 # psych::fa()のデフォルトであるoblimin回転をしてみる 2 coef(result_expl_comp, rotate = "oblimin", simplify = TRUE) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Rotation: oblimin $items a1 a2 d1 d2 d3 d4 d5 Q1_1 -0.146 0.953 3.879 2.433 1.457 0.665 -0.785 Q1_2 0.037 1.872 5.685 3.977 3.034 1.234 -1.208 Q1_3 -0.107 2.634 5.789 4.099 3.032 1.048 -1.891 Q1_4 0.351 1.021 3.646 2.414 1.831 0.781 -0.465 Q1_5 0.057 1.680 5.105 3.271 2.221 0.612 -1.628 Q1_6 1.474 -0.050 4.667 3.195 2.067 0.537 -1.707 Q1_7 1.678 -0.070 4.602 2.822 1.798 0.252 -2.011 Q1_8 1.311 0.051 4.263 2.559 1.625 0.028 -2.047 Q1_9 1.920 0.007 5.227 3.255 1.695 0.508 -1.501 Q1_10 1.394 0.110 2.864 1.361 0.075 -0.595 -2.040 $means F1 F2 0 0 $cov F1 F2 F1 1.00 0.35 F2 0.35 1.00 回転後の識別力は,いい感じに最初の 5 項目と最後の 5 項目に分かれてくれました。 検証的補償型モデル 検証的モデルを行う際には,lavaan に似た要領で,引数 model に各次元と項目の関係を記述 してあげる必要があります。 検証的補償型モデル 1 2 3 4 5 6 7 # model式の右辺は列番号だけでもOKです model <- " f_1 = 1-5 f_2 = 6-10 " result_conf_comp <- mirt(dat[, cols], model = model, itemtype = "graded")
8.7 多次元モデル 8 coef(result_conf_comp, simplify = TRUE) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 38 $items a1 a2 d1 d2 d3 d4 d5 Q1_1 0.895 0.000 3.865 2.425 1.453 0.665 -0.779 Q1_2 1.882 0.000 5.693 3.979 3.032 1.228 -1.212 Q1_3 2.590 0.000 5.789 4.095 3.026 1.042 -1.888 Q1_4 1.109 0.000 3.575 2.353 1.780 0.750 -0.465 Q1_5 1.701 0.000 5.118 3.275 2.219 0.605 -1.632 Q1_6 0.000 1.463 4.670 3.197 2.070 0.537 -1.712 Q1_7 0.000 1.606 4.538 2.779 1.770 0.249 -1.980 Q1_8 0.000 1.330 4.263 2.557 1.621 0.025 -2.047 Q1_9 0.000 1.979 5.311 3.310 1.719 0.509 -1.532 Q1_10 0.000 1.421 2.853 1.352 0.071 -0.596 -2.032 $means f_1 f_2 0 0 $cov f_1 f_2 f_1 1 0 f_2 0 1 きちんと項目 1-5 の次元 2 の識別力( a2)と項目 6-10 の次元 1 の識別力( a1)が 0 に固定 されています。先程の探索的モデルと見比べても,因子の順序こそ違えど概ね同じような推定値 が得られていますね。 ち な み に, 多 次 元 モ デ ル に お け る ICC も, 一 次 元 モ デ ル の 時 と 同 じ よ う に plot() や itemplot() によって出力可能です。 多次元でも項目反応カテゴリ特性曲線 1 2 3 4 # もちろん非補償型でも探索的でも可能です # 多値型項目だとすべてのカテゴリのIRCCCが重なるのでめっちゃ見にくい # 全項目やると重いので,which.itemで少しずつやるのがおすすめ plot(result_conf_comp, type = "trace", which.item = 4) また,引数 rot をうまく指定すると,プロットを回すことができます。うまく回せば,多少は 見やすくなるかもしれません。例えば確認的モデルでは,各項目の識別力を一つだけ非ゼロにし ていました。そのため,識別力がゼロである次元の方向は潰すように回転させると良いかもしれ ません(見やすくなるとは言っていない) 。
8.8 次元性の検証 39 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 0 q 2 -2-4 -6 -6 0.2 4 6 0 2 -2 -4 q1 0.0 図 8.24: 項目反応カテゴリ特性曲線 多次元でも itemplot() 1 2 # 項目4はa2=0だったので,a1方向に見るとIRCCCが見やすいかも? plot(result_conf_comp, type = "trace", which.item = 4, rot = list(xaxis = -90, yaxis = 0, zaxis = 0)) Q1_4 1.0 0.8 1.0 0.8 0.6 0.6 (q ) 0.4 0.4 0.2 0.2 0.0 -6 -4 -2 0 2 4 642 0-2 -6 6-4 q1 q2 0.0 図 8.25: グラフを回してみる 8.8 次元性の検証 Section 8.5.2 で少し触れたように,IRT では基本的にそのテストが測定しようとする次元数 は理論的に決まっていることが多く,そのため「何因子が最適化」を探索的に決める方法よりは 「𝑛 因子で問題ないか」を検証する方法が採用されることが多いと思います。ここでは,そのよ
8.8 次元性の検証 40 うな IRT の文脈で用いられることの多い(と思われる)いくつかの次元性検証法を紹介してい きます。 8.8.1 パラメトリックな方法 パラメトリックな方法では, (多次元・一次元の)IRT や因子分析モデルに基づいてパラメータ の推定を行い,その結果を用いて次元性を評価していきます。例えば確認的因子分析によって推 定を行い,SEM 的なモデル適合度をチェックする方法はパラメトリックな検証法の一つと言え るでしょう。その他にも,IRT の文脈で用いられている検証法がいくつか存在します。代表的と 思われる方法の一つが NOHARM(Gessaroli & De Champlain, 1996)です*18 。この方法では, 項目の各ペアでの残差共分散を,推定されたパラメータおよび観測されたデータ(各項目の正答 率)から計算していきます。そして,全ペアに対して計算された残差共分散(を 𝑧 変換したもの) の二乗和(を何倍かしたもの)が 𝜒2 分布に従うことを利用して,「残差共分散が全てゼロであ る」という帰無仮説に対する検定を行います。そして帰無仮説が棄却されなかった場合,「その 次元数で全ての項目ペア間の相関関係を表現できていないとは言えない」ということで,その次 元数で OK,という判断がくだされます。考え方としては Bartlett 検定に近いかもしれません。 8.8.2 ノンパラメトリックな方法 IRT の文脈で提案されたノンパラメトリックな次元性評価法の中でも代表的なものが,ここで 紹介する DIMTEST(Stout et al., 1996)や DETECT(Zhang & Stout, 1999)です。これら の方法では,「その次元(共通因子)数があれば局所独立性が満たされるか」という視点から評 価を行います。Section 8.5.1 では,局所独立性を「どの項目ペア間でも,どの特性値の値でも連 関がない」ことと表現していましたが,これを数式で表すと Cov(𝐲𝑖 , 𝐲𝑗 |𝜃) = 0 for all 𝑖, 𝑗, 𝜃 (8.34) となります*19 。したがって,𝜃 で条件づけた各項目ペアの(条件付き)共分散を計算していけば, 局所独立性を満たすために必要な次元数を考えることができるわけです。 具体的な方法の説明に移る前に,次元性と条件付き共分散のイメージを確認しておきます。 図 8.26 には,2 因子モデルにおいて,なるべく単純構造になるように回転した後の各項目の因 子負荷をベクトルとして表したものです。また 𝚯𝐶1 , 𝚯𝐶2 と表記されているベクトルは,各クラ スタ内の項目の因子負荷の平均を取ったようなものを表しています。同様に 𝚯𝑇 𝑇 と表記されて いるベクトルは,クラスタをまたいだ全項目のベクトルの平均を取ったもので,端的に言えばこ れが「テスト全体で測定しようとしている特性のベクトル」というわけです。よくある例として *18 NOHARM という名のソフトウェアで出力される推定結果をもとに計算する方法なので, 「NOHARM に基づ く方法」などと言ったほうが正確かもしれません。 *19 これは厳密には弱局所独立性 (weak liocl independence) と呼ばれるものです。一方で,ある反応ベクト ル 𝐲𝑝 が得られる確率が各項目の反応確率 𝑦𝑝𝑖 の積に等しい,という形でも局所独立性は表されるのですが,こ ちらは強局所独立性 (strong liocl independence) と呼ばれています。本来はこの強局所独立性が必要なの ですが,現実世界の多くの場面では,弱局所独立性が満たされている場合には強局所独立性も満たされているこ とが多いということで,ここでは弱局所独立性の検証を行っていきます。
8.8 次元性の検証 41 は,𝚯𝑇 𝑇 が「数学の学力」,𝚯𝐶1 , 𝚯𝐶2 がそれぞれ「代数」「幾何」といった感じです。図 8.26 では,𝚯𝐶1 , 𝚯𝐶2 は完全に無相関というわけではないですが,さほど強くもないので,一次元性 が満たされているとは言いづらい感じがします。 図 8.26: 2 次元のテストの因子負荷プロット(Stout et al., 1996, fig. 2) このテストが一次元性を満たしているかを確認するためには,「𝚯𝑇 𝑇 で条件付けられた」共分 散を考えることになります。そして条件付き共分散は,図 8.26 においてはベクトル 𝚯𝑇 𝑇 に向 かって伸ばした足の長さの積のような形で表されます。これを理解するために,図 8.26 に少し 描き足してみました。図 8.27 に追加されたオレンジの線は,𝚯𝑇 𝑇 = 𝚯𝐶1 +𝚯𝐶2 となる(次元間 2 で重みが同じ)場合に 𝚯𝑇 𝑇 の値が同じになる 𝚯𝐶1 , 𝚯𝐶2 の組みあわせを表しています。 図 8.27: 条件付き共分散(Stout et al., 1996, fig. 2 に少し加筆) この図をもとに,項目ペアに対して計算できる条件付き共分散を考えてみましょう。まず,ク ラスター 𝐶1 の中の 2 つの項目の条件付き共分散は正の値になることが期待されます。𝚯𝑇 𝑇 の 値が同じ 2 人(A さんと B さん)について考えてみると,A さんはクラスター 𝐶1 の項目に正解 するのに必要な能力(𝚯2 )の値が高いためクラスター 𝐶1 の項目には正解する確率が高いと予想
8.8 次元性の検証 42 される一方で,B さんは 𝚯2 の値が低いので,正解する確率は低そうです。これをまとめると, 「𝚯𝑇 𝑇 で条件づけたとき,クラスター 𝐶1 内のある項目に正解する人ほど,同じクラスター内の 別の項目にも正解する確率が高い」わけなので,条件付き共分散は正の値になります。 同様に,クラスター 𝐶1 の中のある項目と,別のクラスター 𝐶2 の中のある項目の条件付き共 分散を考えてみると,A さんはクラスター 𝐶1 の項目には正解するがクラスター 𝐶2 の項目には 正解できないと予想されます。そして B さんは反対に,クラスター 𝐶2 の項目にのみ正解できる 可能性が高いでしょう。したがって,「𝚯𝑇 𝑇 で条件づけたとき,クラスター 𝐶1 内のある項目に 正解する人ほど,別のクラスター内のある項目に正解する確率は低い」ことになり,条件付き共 分散は負の値になります。 以上をまとめると,図 8.27 における条件付き共分散の符号は,2 つの項目の因子負荷ベクトル が「ベクトル 𝚯𝑇 𝑇 から見て」同じ方向にあるかどうかに対応するわけです。そして,仮に全て の項目が近い方向のベクトルで表される場合には,当然 𝚯𝑇 𝑇 ベクトルも同じような方向になる ため,条件付き共分散も 0 に近くなっているはずです。という感じで,(8.34) 式が局所独立性, そして一次元性を表しています。 DIMTEST DIMTEST(Stout et al., 1996)では,まずテストを「(次元性の)評価用 (assessment subtest: AT)」と「データの分割用 (partitioning subtest: PT)」の 2 つに分割します。この分割では,理 論的に想定される分かれ方(図 8.26 の 𝚯𝐶1 , 𝚯𝐶2 に相当)でわけたり,探索的因子分析をもと に分けたりします。つまり,一旦「このテストが分かれる(2 次元だ)としたらこうなるだろう」 という形で分けてみるのです。もしこのテストが実際には一次元だとしたら,そのように分けた としても結局全ての項目のベクトルの向きは近いので,片方のテストの得点によって条件付けた としても,もう一方のテストの中の共分散は 0 に近いはずです。一方で,テストが(𝚯𝐶1 , 𝚯𝐶2 のように)分かれる場合には,例えば 𝚯𝐶1 で条件づけた際の 𝚯𝐶2 内の条件付き共分散は全て正 になってしまいます。 以上の考え方に基づいて,DIMTEST では以下の統計量を計算します。 ̂ , 𝐲 |𝑆 = 𝑘). 𝑇 = ∑ ∑ Cov(𝐲 𝑖 𝑗 𝑃𝑇 (8.35) 𝑘 (𝑖,𝑗∈𝐴𝑇 ) ここで 𝑆𝑃 𝑇 は PT における正答数を表しています。したがって,(8.35) 式は「PT における正答 数で(= 𝚯𝑃 𝑇 で)条件づけたときの AT 内の共分散」を求めているわけです。統計量 𝑇 が 0 に 近いほど局所独立が満たされている(そのテストは一次元である)と言えるので,あとはこれを 標準化した値を用いて標準正規分布に基づく 𝑧 検定(Stout, 1987)を行います。 DIMTEST の成功のカギは,AT と PT の分割にあるわけですが,どのように分けたら良いか については Stout et al.(1996)などを参照してください。 DETECT DETECT も条件付き共分散をベースにした方法です。以下では,二値データに対して考案さ れたもの(Zhang & Stout, 1999)を紹介します。多値データの場合にはこれを発展させた方法
8.8 次元性の検証 43 (Zhang, 2007)が提案されています。 DETECT でも,まずはテストをいくつかのクラスター(𝐶1, 𝐶2, ⋯ , 𝐶𝑘)に分割します。そ して以下のようにテスト全体で推定された特性値 𝚯𝑇 𝑇 によって条件づけた共分散の統計量を計 算します。 𝐷= 2 ∑ 𝛿 𝐸[Cov(𝐲𝑖 , 𝐲𝑗 |𝚯𝑇 𝑇 )]. 𝐼(𝐼 − 1) 1≤𝑖≤𝑗≤𝐼 𝑖,𝑗 (8.36) ここで重要なのは指示変数 𝛿𝑖,𝑗 です。これは, • 2 つの項目 𝑖, 𝑗 が同じクラスターにある場合には 1 • 2 つの項目 𝑖, 𝑗 が異なるクラスターにある場合には-1 になります。 DETECT が計算しようとしているものを 図 8.26 をもとに考えてみます。図 8.26 は,あるテ ストの真のクラスターを表しているとします。つまりあるテストはきれいに 2 つのクラスターに 分かれている,ということです。このときに,この真のクラスターと一致するようにテストを分 けた上で DETECT の統計量 𝐷 を計算してみましょう。 • 2 つ の 項 目 𝑖, 𝑗 が 同 じ ク ラ ス タ ー に あ る 場 合 に は, 条 件 付 き 共 分 散 (の 期 待 値) 𝐸[Cov(𝐲𝑖 , 𝐲𝑗 |𝚯𝑇 𝑇 )] も 𝛿𝑖,𝑗 も正の値 • 2 つ の 項 目 𝑖, 𝑗 が 異 な る ク ラ ス タ ー に あ る 場 合 に は, 条 件 付 き 共 分 散 (の 期 待 値) 𝐸[Cov(𝐲𝑖 , 𝐲𝑗 |𝚯𝑇 𝑇 )] も 𝛿𝑖,𝑗 も負の値なので,積は正の値になる ということで,全ての項目ペア間で計算される値(∑ の中身)が全て正の値になります。この ように,統計量 𝐷 は(クラスター数も含めて)最も適切なクラスター基づいて計算したときに 最大値を取る値です。そして,そもそも全ての項目のベクトルの方向が近ければ(次元性が満た されていれば),条件付き共分散(の期待値)は全て 0 に近いので,その平均である統計量 𝐷 も 0 に近い値になっているはずです。 ということで,この統計量を用いて一次元性を評価する場合には,「統計量 𝐷 が最大でいくつ になるか」を見たら良いと言えます。例えば「理論的に分かれるとしたらこう」で分けてみたり, 探索的因子分析の結果をもとに分けてみたり,あるいはとにかく手当たり次第にランダムに様々 な分け方をひたすら試してみたりして,統計量 𝐷 が高々どれくらいの値に収まるかを確認しま しょう。一次元性の検証という意味では,もちろん 𝐷 の値は小さいほどよいのですが,二値デー タの場合例えば Roussos & Ozbek(2006)では最大値が 0.2 以下なら一次元とみなしてよいの ではないか,と言われていたりします。 そ ん な DETECT に つ い て は, R で 計 算 す る た め の 関 数 と し て sirt パ ッ ケ ー ジ に conf.detect() および expl.detect() というものが用意されています*20 。これらは二値・多 値のどちらにも対応していますが,ここでは二値データへの適用例を見てみます。 まずは conf.detect() です。conf は”confirmatory” のことで,その名の通り分析者が指定 *20 DIMTEST については関数が見つけられていません。一応 DIM-Pack というフリーソフトを使えば出来そう ですが…
8.8 次元性の検証 したクラスター構造のもとで DETECT を計算します。そのため,事前に想定されるクラスター がある(数学のテストの中に複数の単元がある,心理尺度の中に複数の下位概念がある)ような 場合に使える方法です。 DETECT で次元性の検証 1 2 3 4 5 6 7 8 9 10 11 # install.packages("sirt") library(sirt) 1 2 3 4 5 6 7 8 9 10 11 12 13 ----------------------------------------------------------Confirmatory DETECT Analysis Conditioning on 1 Scores Bandwidth Scale: 1.1 DETECT Calculation Score 1 ----------------------------------------------------------NScores Mean SD Min Max DETECT Unweighted 1 1.064 NA 1.064 1.064 DETECT Weighted 1 1.064 NA 1.064 1.064 ASSI Unweighted 1 0.511 NA 0.511 0.511 ASSI Weighted 1 0.511 NA 0.511 0.511 RATIO Unweighted 1 0.664 NA 0.664 0.664 RATIO Weighted 1 0.664 NA 0.664 0.664 # テスト全体で測定している特性値 # シンプルに和得点などを使用しても良い fs <- fscores(result_2plm) # 各項目がどのクラスタに属するか # 今回は元のデータの想定に合わせて設定しておきます cluster <- rep(1:2, each = 5) res_detect <- conf.detect(dat_bin, score = fs, itemcluster = cluster) # conf.detect()の返り値には各ペアの残差共分散などが含まれている 結果の中の DETECT Unweighted にかかれている値が (8.36) 式で計算される 𝐷 の値です。こ の値が 0.2 より小さければ,そのテストは一次元であるとみなせそうなのですが…今回使用した 10 項目は元々 2 次元だったため,「一次元ではない」と正しく判定されています(一般に 1 を超 えると「かなり強く多次元である」と言えるようです) 。 続いて,本来一次元であることが想定されるケースとして,項目間の相関が高い(=内的一貫 性の高い)Q1_6 から Q1_10 の 5 項目のみを使って conf.detect() してみましょう。 一次元性が示されそうな例 1 2 3 4 # 今度は和得点でやってみます sum_s <- apply(dat_bin[, 6:10], 1, sum) # 各項目がどのクラスタに属するか # 今回は適当に奇数・偶数で分けてみます 44
8.8 次元性の検証 5 6 cluster <- c(1, 2, 1, 2, 1) res_detect <- conf.detect(dat_bin[, 6:10], score = sum_s, itemcluster = cluster) 1 2 3 4 5 6 7 8 9 10 11 ----------------------------------------------------------Confirmatory DETECT Analysis Conditioning on 1 Score Bandwidth Scale: 1.1 ----------------------------------------------------------unweighted weighted DETECT 0.425 0.425 ASSI 0.200 0.200 RATIO 0.140 0.140 MADCOV100 3.041 3.041 MCOV100 -3.041 -3.041 DETECT Unweighted の値は先程よりは小さくなっていますが,それでも 0.2 よりは大きく なっています…理由はよく分からないのですが,もしかしたら項目数が少なすぎる場合には小さ くなりにくいのかもしれません(要出典) 。 ここまでは,引数 itemcluster によって項目のクラスタリングを与えていました。一次元性 を検証するためには,本当はこれをすべてのパターン(𝐼 項目を 𝐾 個のクラスタに分けるとした らおよそ 𝐾 𝐼 通り)で計算して,その最大値を出す必要があります。ですがこれはどう考えても 非効率な(というかすぐ無理になる)話です。そこで,各クラスタ数における「𝐷 が最大になる 分割」を探索的に求めることを考えましょう。 expl.detect() では,階層的クラスター分析を用いて,指定したクラスタ数の分割の中で 𝐷 が最大になるものを探索的に決定してくれます。考え方はシンプルで,条件付き共分散(の期待 値)𝐸[Cov(𝐲𝑖 , 𝐲𝑗 |𝚯𝑇 𝑇 )] が大きいペアは近くにある=同じクラスターに属すると考えて,距離 の近い項目から順に,最終的に 𝐾 個のクラスターにまとまるまでくっつけていき,そのクラス ター構造の元で計算した 𝐷 を出力してくれます。ということで早速試してみましょう。 探索的 DETECT 1 2 3 4 fs <- fscores(result_2plm) # nclustersは想定される最大クラスター数 # 想定よりやや大きめの値にしておくのが無難 res_detect_expl <- expl.detect(dat_bin, score = fs, nclusters = 5) 1 2 3 4 5 Pairwise Estimation of Conditional Covariances ........................................................... Nonparametric ICC estimation ........................................................... 45
8.8 次元性の検証 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 46 Nonparametric Estimation of conditional covariances 5% 10% 15% 20% 25% 30% 35% 40% 45% 50% 55% 60% 65% 70% 75% 80% 85% 90% 95% Pairwise Estimation of Conditional Covariances ........................................................... Nonparametric ICC estimation ........................................................... Nonparametric Estimation of conditional covariances 5% 10% 15% 20% 25% 30% 35% 40% 45% 50% 55% 60% 65% 70% 75% 80% 85% 90% 95% DETECT (unweighted) Optimal Cluster Size is 4 (Maximum of DETECT Index) N.Cluster N.items N.est N.val size.cluster DETECT.est ASSI.est RATIO.est 2 10 1216 1216 5-5 1.133 0.511 0.674 3 10 1216 1216 5-3-2 1.590 0.778 0.947 4 10 1216 1216 5-2-1-2 1.659 0.867 0.988 5 10 1216 1216 4-1-2-1-2 1.644 0.867 0.979 MADCOV100.est MCOV100.est DETECT.val ASSI.val RATIO.val MADCOV100.val 1 1.68 -1.284 1.136 0.556 0.680 1.669 2 1.68 -1.284 1.543 0.822 0.924 1.669 3 1.68 -1.284 1.612 0.822 0.966 1.669 4 1.68 -1.284 1.583 0.733 0.949 1.669 MCOV100.val 1 -1.283 2 -1.283 3 -1.283 4 -1.283 1 2 3 4 結果を見ると,このデータ( Q1_1 から Q1_10)に対してはクラスター数 4 のときが最大値に なるようです。図 8.28 にその時の分かれ方が示されています。 (やっぱり項目数が少ないとクラ スター数が多めになるのか?) 一点注意として,expl.detect() の内部では,デフォルトではデータを半分に分割して, 1. 前半のデータで条件付き共分散を計算→クラスター構造を決定 2. 後半のデータで 𝐷 を計算 という流れで計算を行っています。これは多分オーバーフィッティングを避けるためだと思い ます。そのため,N.Cluster が 2 のところの DETECT.est(手順 1 で使用したデータで計算され た 𝐷)も DETECTR.val(手順 2 で使用したデータで計算された 𝐷)も,conf.detect() での結 果 1.064 とは異なる値になっています( N.est および N.val がいずれも元の dat の行数の半
8.9 適合度 47 Q1_7 Q1_6 Q1_8 Q1_10 Q1_9 Q1_2 Q1_1 Q1_5 Q1_3 Q1_4 0.00 Height 0.08 Cluster Dendogram with 4 Clusters d stats::hclust (*, "ward.D") 図 8.28: 𝐷 が最大になるクラスターのデンドログラム 分になっている) 。 …ということで,今回のデータではきれいな一次元性を確認することはできなかったのです が,このような方法を用いると,「一次元とみなして分析しても問題ないか」を評価することが できるのです。 8.9 適合度 IRT でも SEM と同じように適合度を考えることができます。特に「個人」と「項目」の両方 に関心があることの多い IRT では,それぞれに対しての適合度を考える必要があります。もち ろんこれから紹介する適合度の考え方は,そのまま因子分析や SEM にも応用しようと思えばで きるので,そういった意味でも参考にしてみてください。 8.9.1 局所独立性の確認 まずはじめに,データとモデルの適合度の一つとして,データが局所独立の仮定に適合してい るかをチェックします。Section 8.5.1 で説明したように,局所独立性が満たされている場合に は,𝜃𝑝 で条件づけたときの 2 項目の回答の間に相関が無いはずです。 𝑋 2 統計量 𝑋 2 統計量(Chen & Thissen, 1997)では,カテゴリ変数のクロス表に対して用いられる 𝜒2 検定の枠組みを利用して局所独立性を評価します。ある二値項目 𝑖 の項目パラメータ (𝛼𝑖 , 𝛽𝑖 ) が 分かっている場合,その項目の母集団全体での期待正答率は ∞ 𝐸[𝑃 (𝑦𝑝𝑖 = 1)] = ∫ −∞ 𝑃 (𝑦𝑝𝑖 = 1|𝜃)𝑓(𝜃)𝑑𝜃 (8.37)
8.9 適合度 48 表 8.3: 期待度数(独立な 2 項目の分布) 項目 𝑗 項目 𝑖 当てはまらない 当てはまる 計 当てはまらない 当てはまる 計 𝐸00 = 𝐸01 = 𝑁 × 𝑃 (𝑦𝑝𝑖 = 0 ∧ 𝑦𝑝𝑗 = 0) 𝐸10 = 𝑁 × 𝑃 (𝑦𝑝𝑖 = 0 ∧ 𝑦𝑝𝑗 = 1) 𝐸11 = 𝑁 × 𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 0) 𝑁 × 𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 1) – – – – 𝑁 という形で求めることができます。ここで 𝑓(𝜃) は特性値の確率分布(ふつうは標準正規分布)に おける 𝑥 = 𝜃 での確率密度を表しています。同様に,その項目の誤答率の期待値を求めるならば 𝑃 (𝑦𝑝𝑖 = 1|𝜃) の部分を 1 − 𝑃 (𝑦𝑝𝑖 = 1|𝜃) にしたら良いだけです。 この考え方を 2 項目に広げると,母集団全体において項目 (𝑖, 𝑗) に両方とも正解する確率はも し 2 項目が独立ならば ∞ 𝐸[𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 1)] = ∫ 𝑃 (𝑦𝑝𝑖 = 1|𝜃)𝑃 (𝑦𝑝𝑗 = 1|𝜃)𝑓(𝜃)𝑑𝜃 (8.38) −∞ と表せます。ということは,サンプルサイズが 𝑁 の場合,項目 (𝑖, 𝑗) に両方とも正解する人数の 期待値 𝐸11 は 𝑁 × 𝑃 (𝑦𝑝𝑖 = 1 ∧ 𝑦𝑝𝑗 = 1) です。同様にして,一方の項目にだけ正解する人数お よび両方の項目に誤答する人数の期待値も求めると,表 8.3 のようになります。 この期待度数分布はもし 2 項目が局所独立ならばこうなるだろうという状態を表しており, 表 8.1 に相当するものを 𝜃 について周辺化して作成したものといえます。つまり局所独立ではな くなるほど,実際のデータでのクロス表は(表 8.2 のように)この期待度数から離れるだろう, ということです。 クロス表での連関の検定は(学部の統計でやっているかもしれない)𝜒2 検定です。その検定 統計量は 1 1 (𝑂𝑠𝑡 − 𝐸𝑠𝑡 )2 𝐸𝑠𝑡 𝑠=0 𝑡=0 𝜒2 = ∑ ∑ (8.39) で計算されます。ここで 𝑂𝑠𝑡 は,実際のデータで 𝑦𝑝𝑖 = 𝑠 ∧ 𝑦𝑝𝑗 = 𝑡 だった人数を表しています。 こうして計算された 𝜒2 統計量は,自由が (行数-1) × (列数-1) の 𝜒2 分布に従うことが知られ ています。二値データの場合は 1 × 1 = 1 ということです。そしてここまでのプロセスはそのま ま多値項目に対しても拡張可能です。例えば dat の項目のように 6 件法どうしの場合は,自由度 は 5 × 5 = 25 となるわけです。ということで,項目のペア毎に 𝜒2 統計量を計算して統計的仮説 検定を行い,有意だったペアの間は局所独立ではない可能性が疑われる,ということになります。
8.9 適合度 49 𝑄3 統計量 𝑄3 統計量(Yen, 1984)では,回帰分析でいうところの偏相関係数にあたるものを利用します。 偏相関係数は「変数 𝑧 の影響を取り除いた時の 𝑥 と 𝑦 の相関係数」で,疑似相関の影響を検討す る際などに利用されるものです。これを IRT の文脈に当てはめると「変数 𝜃𝑝 の影響を取り除い た時の 𝑦𝑝𝑖 と 𝑦𝑝𝑗 の相関係数」ということで,これが 0 であれば局所独立だろうと言えそうです。 回帰分析では,𝑥 と 𝑦 をそれぞれ 𝑧 で回帰した時の予測値 𝑥,̂ 𝑦 ̂ に対して,残差の相関 𝑟(𝑥−𝑥,𝑦− ̂ 𝑦)̂ のことを偏相関係数と呼んでいました。IRT(ロジスティックモデル)もその中身はロジスティッ ク回帰なので,同じようにして 𝑦𝑝𝑖 と 𝑦𝑝𝑗 をそれぞれ 𝜃𝑝 で予測した時の予測値(期待正答率)を 𝑃 ̂ (𝑦𝑝𝑖 = 1) = 1 (8.40) 1 + exp [−𝛼𝑖̂ (𝜃𝑝̂ − 𝛽𝑖̂ )] というように,𝛼𝑖 , 𝛽𝑖 , 𝜃𝑝 にそれぞれ推定値 𝛼𝑖̂ , 𝛽𝑖̂ , 𝜃𝑝̂ を入れることで計算が出来ます。すると, 残差を 𝑑𝑝𝑖 = 𝑦𝑝𝑖 − 𝑃 ̂ (𝑦𝑝𝑖 = 1) と求めることができます。図 8.29 に 𝑑𝑝𝑖 を表してみました。 1.00 A B ypi 0.75 dpi 0.50 0.25 0.00 -4 -2 0 ^ q p 2 4 図 8.29: 残差 この項目には A さんも B さんも正解しているため,ふたりとも 𝑦 軸の値は 1 です。ですが他 の全項目から推定された 2 人の 𝜃𝑝̂ はそれぞれ −2, 2 だったとします。この場合,A さんは他の 項目の出来から考えるとこの問題にはほぼ間違えるはずなのに正解しています。ということで 𝑑𝑝𝑖 の値が大きくなっています。このようになる理由としては,やはり 𝜃𝑝 以外の別の要因によっ て正解できた,と考えるのが妥当でしょう。これを全体に広げてみた時に,𝑑𝑝𝑖 が高い人ほど別 の項目でも同様に 𝑑𝑝𝑗 の値が大きくなっているとしたら,この 2 つの項目には何か正解に寄与す る 𝜃𝑝 以外の別の要因がある,と考えることができるわけです。 ということで 𝑄3 統計量は 𝑄3 = 𝑟𝐝𝑖 ,𝐝𝑗 (8.41) という形になります。𝑄3 統計量は相関係数なので,サンプルサイズに依存しない効果量の指標 として見ることができ,絶対値が 0.2 を超えてくると怪しいと判断できるようです。
8.9 適合度 50 R でやってみる mirt では residuals という関数でいま紹介した指標を出してくれます。引数 type で,出 してほしい統計量を指定できます。ちょっと厄介なのですが,𝜒2 統計量を出してほしい場合には type = "LD"と指定してください。LD は local dependence の頭文字ですが,そういったら 𝑄3 だって LD だろ,と思われるかもしれません。これは Chen & Thissen(1997)の書き方に則っ ているのだと思います。我慢してください。 X2 統計量を出す 1 2 # 引数df.p = TRUEとすると,自由度とp値を出してくれる residuals(result_2plm, type = "LD", df.p = TRUE) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 Degrees of freedom (lower triangle) and p-values: Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_1 NA 0 0 0.019 0 0.000 0.000 0.093 0.824 0.004 Q1_2 1 NA 0 0.000 0 0.172 0.000 0.208 0.000 0.000 Q1_3 1 1 NA 0.000 0 0.006 0.003 0.000 0.000 0.000 Q1_4 1 1 1 NA 0 0.001 0.193 0.000 0.001 0.485 Q1_5 1 1 1 1.000 NA 0.003 0.000 0.004 0.000 0.019 Q1_6 1 1 1 1.000 1 NA 0.000 0.002 0.099 0.477 Q1_7 1 1 1 1.000 1 1.000 NA 0.032 0.019 0.550 Q1_8 1 1 1 1.000 1 1.000 1.000 NA 0.045 0.121 Q1_9 1 1 1 1.000 1 1.000 1.000 1.000 NA 0.000 Q1_10 1 1 1 1.000 1 1.000 1.000 1.000 1.000 NA LD matrix (lower triangle) and standardized values. Upper triangle summary: Min. 1st Qu. Median -0.109 -0.067 -0.026 Mean 3rd Qu. 0.002 0.063 Max. 0.219 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_1 NA 0.144 0.111 0.048 0.082 -0.098 -0.090 -0.034 0.005 -0.059 Q1_2 50.489 NA 0.171 0.077 0.130 -0.028 -0.071 -0.026 -0.100 -0.103 Q1_3 29.999 71.332 NA 0.102 0.219 -0.055 -0.061 -0.078 -0.109 -0.083 Q1_4 5.504 14.503 25.105 NA 0.094 -0.067 -0.026 -0.079 -0.065 0.014 Q1_5 16.380 41.085 116.151 21.292 NA -0.059 -0.080 -0.058 -0.104 -0.048 Q1_6 23.577 1.863 7.479 10.768 8.539 NA 0.101 0.063 0.033 -0.014 Q1_7 19.698 12.128 8.934 1.695 15.588 24.903 NA 0.044 0.048 0.012 Q1_8 2.820 1.583 14.878 15.199 8.112 9.698 4.621 NA 0.041 0.031 Q1_9 0.049 24.351 29.035 10.195 26.099 2.723 5.519 4.012 NA 0.123 Q1_10 8.358 25.794 16.877 0.487 5.500 0.507 0.357 2.411 37.068 NA 上の行列の上三角が 𝑝 値です。したがってこれが 0.05 を下回っている場合は「局所独立では
8.9 適合度 51 ない」という判断になるのですが,見た感じかなり多くのペアで有意になっています。これは, 統計的仮説検定なのでサンプルサイズの影響を受けているためです。また上の行列の下三角は自 由度です。今回は二値に変換したデータに対して行っているので自由度は 1 になっています。も しも 6 件法のまま (GRM や GPCM で) 推定した結果に対して出した場合には,自由度は先程説 明した通り 5 × 5 = 25 となります。 そして下の行列は実際に計算された 𝜒2 統計量(下三角)および標準化した値(上三角;たぶ んクラメールの連関係数)です。ということで,この行列の上三角について大きいペアが,局所 独立からより強く離れているといえます。 続いて 𝑄3 統計量です。 Q3 統計量を出す 1 2 # 仮説検定は無いので引数df.pは書いても無視される residuals(result_2plm, type = "Q3") 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Q3 summary statistics: Min. 1st Qu. Median -0.235 -0.166 -0.101 Mean 3rd Qu. -0.067 0.033 Max. 0.244 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_1 1.000 0.150 0.109 0.033 0.074 -0.149 -0.148 -0.072 -0.027 -0.115 Q1_2 0.150 1.000 0.188 0.056 0.131 -0.101 -0.173 -0.097 -0.209 -0.193 Q1_3 0.109 0.188 1.000 0.082 0.244 -0.149 -0.170 -0.182 -0.235 -0.187 Q1_4 0.033 0.056 0.082 1.000 0.071 -0.166 -0.122 -0.186 -0.175 -0.068 Q1_5 0.074 0.131 0.244 0.071 1.000 -0.152 -0.197 -0.152 -0.226 -0.145 Q1_6 -0.149 -0.101 -0.149 -0.166 -0.152 1.000 0.054 0.005 -0.048 -0.114 Q1_7 -0.148 -0.173 -0.170 -0.122 -0.197 0.054 1.000 -0.041 -0.047 -0.105 Q1_8 -0.072 -0.097 -0.182 -0.186 -0.152 0.005 -0.041 1.000 -0.045 -0.069 Q1_9 -0.027 -0.209 -0.235 -0.175 -0.226 -0.048 -0.047 -0.045 1.000 0.040 Q1_10 -0.115 -0.193 -0.187 -0.068 -0.145 -0.114 -0.105 -0.069 0.040 1.000 ところどころ高めの相関が出ており,完全な局所独立ではないといえそうです。というのも今 回のデータは,本来 2 因子に分かれるはずの 10 項目を無理やり 1 因子とみなして分析していま す。そもそもデータが一次元性を持っていないはずなので,局所独立性も保たれてはいない,と いうことになります。(本来こういう時には一旦立ち止まり,項目の内容を吟味したりモデルを 変更したりといった手段を考えるべきですが,ここでは気にせず進めます。 ) ちなみに,引数 suppress を与えると,絶対値が suppress 以下の値が NA に置き換わるた め,項目がたくさんある場合でも「閾値を超えているペア」を探しやすくなるかもしれません。 表示する閾値を設定 1 residuals(result_2plm, type = "Q3", suppress = 0.2)
8.9 適合度 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 52 Q3 summary statistics: Min. 1st Qu. Median -0.235 -0.166 -0.101 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_1 1 Q1_2 Mean 3rd Qu. -0.067 0.033 Q1_3 Q1_4 Max. 0.244 Q1_5 Q1_6 Q1_7 Q1_8 1.000 1.000 Q1_9 Q1_10 -0.209 -0.235 0.244 1 0.244 1.000 -0.226 1 1 1 -0.209 -0.235 -0.226 1.000 1 8.9.2 個人適合度 続いて個人の適合度です。IRT における個人適合度の基本的な考え方としては「𝜃𝑝 が高いのに 困難度が低い項目に間違えるのはおかしい」「𝜃𝑝 が低いのに困難度が高い項目に正解するのはお かしい」というものです。するとこれは,図 8.29 的な考え方になります。A さんのほうが予測 値と実測値の間の差 𝑑𝑝𝑖 が大きいことは,前述の通り「𝜃𝑝 以外に回答に影響を与える要因がある (のでモデルが間違っている)」という可能性の他に,「 (もしモデルは正しいと仮定したら)回答 者の方に,何か項目 𝑖 の正解に導いた要因があるのでは」と考える,というわけです。 回答者の方の原因としては,心理尺度でいえば以下のようなものが考えられます(Meijer et al., 2016)。 • 回答者の言語能力が不十分なために項目の内容を十分に理解できていない可能性。低学齢 への実施や第一言語以外での実施で起こる可能性がある。 • 心理尺度が測定しようとしている心理特性に対する認識がない状態(Traitedness)。例え ば「そんなこと考えたこともない」といった場合,項目の意図を他の人とは異なるニュア ンスで捉えてしまい,結果的に「大多数が当てはまると回答する項目に当てはまらないと 回答する」といった反応につながる。 • 回答のモチベーションが低い可能性。やる気が無いので適当に回答している? • 特定の心理特性を持っている人は平均的な人より一貫性の低い回答をする可能性がある。 例えば精神的な問題を抱えていたり,希死念慮を持ち合わせている人はそうでない人と比 べて適合度が低いという先行研究もあるらしい。 ということで,個人適合度が低いイコール適当な回答者(だから除外してもよい)とも言い切 れないのですが,特徴的な回答者をあぶり出すなどの目的でも個人適合度は使い道がありそう です。 実際に個人適合度の指標として提案されているものは相当数ある(Karabatsos, 2003; Meijer
8.9 適合度 53 & Sijtsma, 2001)のですが,わかりやすい&統計的仮説検定のフレームワークが整備されてい ると言った理由から,𝑧ℎ(またの名を 𝑙𝑧 )という統計量(Drasgow et al., 1984)が用いられるこ とが多いようです(Meijer et al., 2016) 。 𝑧ℎ (𝑙𝑧 ) 統計量 𝑧ℎ (𝑙𝑧 ) 統計量では,ある回答者の実際の反応パタン 𝐲𝑝 = [𝑦𝑝1 𝑦𝑝2 ⋯ 𝑦𝑝𝐼 ] と「考えられ る全反応パタン」での尤度を比較することにより,その回答者の反応パタンの「普通さ」を評価 します。 まず,回答者 𝑝 の特性値の最尤推定値 𝜃𝑝̂ に関して,その時の尤度を計算します。局所独立性 が満たされているならば,尤度は単純に全ての項目をかけ合わせたら良いので, 𝐼 𝑙0 = ∏ 𝑃 (𝑦𝑝𝑖 = 1)𝑦𝑝𝑖 𝑃 (𝑦𝑝𝑖 = 0)1−𝑦𝑝𝑖 (8.42) 𝑖=1 となります。ただ掛け算は(コンピュータ的に)大変なので,対数尤度 𝐼 𝑙𝑙0 = ∑ 𝑦𝑝𝑖 ln 𝑃 (𝑦𝑝𝑖 = 1) + (1 − 𝑦𝑝𝑖 ) ln 𝑃 (𝑦𝑝𝑖 = 0) (8.43) 𝑖=1 を使います。この尤度 𝑙𝑙0 は,「特性値が 𝜃𝑝̂ の人が 𝐲𝑝 という回答をする確率」のようなもの*21 です。ここで 𝐲𝑝 以外の回答パタンに目を向けてみると,もしかしたら「𝐲𝑝 よりも 𝜃𝑝̂ と整合的 な回答パタン」があるかもしれません。もちろん反対に「𝜃𝑝̂ と全く合わない回答パタン」という ものもあるかもしれません。 そこで,考えられる全回答パタン(2 件法 10 項目ならば 210 通り)の中で 𝐲𝑝 という回答パタ ンが 𝜃𝑝̂ とどの程度整合的かを相対的にチェックすることを考えます。全回答パタン内で 𝐲𝑝 とい う回答パタンと 𝜃𝑝̂ の整合性(尤度)が平均くらいであれば,その回答パタンはまあ起こりうる ものだろうとみなせる一方で,もし整合性が低い場合には,𝐲𝑝 という回答パタンは起こりにく い(ズレが大きい)だろうと判断できる,ということです。 実際に「全回答パタンでの対数尤度の期待値」は(各項目への回答は離散変数なので) 𝐼 𝐸ℎ = ∑ 𝑃 (𝑦𝑝𝑖 = 0) ln 𝑃 (𝑦𝑝𝑖 = 0) + 𝑃 (𝑦𝑝𝑖 = 1) ln 𝑃 (𝑦𝑝𝑖 = 1) (8.44) 𝑖=1 として求めることができます。つまり 𝑙𝑙0 − 𝐸ℎ が「平均的な回答パタンでの尤度と実際の回答パ タン 𝐲𝑝 での尤度の差」であり,これがマイナス方向に大きいほどモデルとデータの適合度が悪 いということを意味します。 ここで全回答パタンにおける対数尤度の分散を 𝐼 𝜎ℎ2 = ∑ {𝑃 (𝑦𝑝𝑖 = 1)𝑃 (𝑦𝑝𝑖 = 0) [ln 𝑖=1 *21 正確に言えば, 「𝐲 𝑃 (𝑦𝑝𝑖 = 1) 2 ] } 𝑃 (𝑦𝑝𝑖 = 0) (8.45) 𝑝 という回答をした人の特性値として 𝜃𝑝 という値がどの程度もっともらしいか」という感じ ですが,あえて曲解させています。 ̂
8.9 適合度 54 として求めることができるのですが,これを用いると, 𝑙𝑙0 を全回答パタン内で標準化した 値 𝑧ℎ = 𝑙𝑙0 −𝐸ℎ が経験的に標準正規分布に従うということが知られています。ということで, 𝜎ℎ |𝑧ℎ | > 1.96 となった回答者はどうやら「特性値が 𝜃𝑝̂ のもとでは平均的ではない回答パタンであ る」と判断することができるわけです。 図 8.30 に,𝑧ℎ 統計量の考え方を表しました。ヒストグラムは母集団からサンプリングされた 様々な 𝜃𝑝̂ の人と,その人の(ランダムサンプリングされた)回答ベクトル 𝐲𝑝 のもとで計算した 𝑧ℎ を集めたものというイメージです。これが近似的に正規分布に従っているということで,こ の外側 5% を棄却域とみなして検定を行うことができます。B さんの 𝑧ℎ は,この分布の中で割 と平均的なところにあるため,B さんの回答は「 (B さんの推定値 𝜃B̂ における)ふつう」に近い とみなせます。一方 A さんの 𝑧ℎ は,この分布の中でもかなり小さい方にあります。ということ ̂ からすると当てはまりが悪く(つまり困難度の低 は,A さんの回答パタンは A さんの推定値 𝜃A い項目に「当てはまらない」と回答し,困難度の高い項目で「当てはまる」と回答している可能 性が高く),なんなら母集団全体の中でも相当悪いほうだと判断でき,何らかの問題があるので はないかと疑いをかけられるわけです。 図 8.30: 個人適合度の考え方 ちなみに,正規分布的には上側の外れ値も「ふつうではない」ということになりますが,こち らはどう考えたら良いかよくわかりません。理論的には「特性値 𝜃𝑝̂ に対してあまりにも適合し すぎている」という状態で,これは例えば困難度が 𝜃 ̂ 以下の項目では全て「当てはまる」と回答 し,𝜃 ̂ 以上の項目ではすべて「当てはまらない」と回答しているような感じなのですが,だから といって即座に問題とは言えなさそうです。ということで,基本的には 𝑧ℎ < −1.96 の回答者に 着目しておけば良いと思います。 R でやってみる ※個人適合度と回答パタンの関係がわかりやすいので,ここは GRM で推定した結果を使用し ます。 mirt では,personfit() という関数によって 𝑧ℎ 統計量を出すことが出来ます。
8.9 適合度 個人適合度を出す 1 personfit(result_grm) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 outfit z.outfit infit z.infit Zh 1 0.3043410 -2.51418387 0.2901436 -2.630507534 1.11739286 2 0.4319025 -1.53744653 0.4991522 -1.313386687 1.54458648 3 0.5417319 -1.17103185 0.6082584 -0.964496321 1.30478354 4 0.6961604 -0.71567303 0.7021298 -0.707746942 0.38159507 5 0.5979698 -0.95977207 0.5796844 -1.045595986 0.06132109 6 1.6611227 0.99612493 2.2816175 1.585552919 -0.57702558 7 0.2629479 -2.06106332 0.2970130 -1.974928887 1.15924498 8 1.2820317 0.82619048 1.2682515 0.803602305 -1.12899488 9 0.4296002 -0.79327083 0.4483505 -0.806286109 0.88382209 10 0.6018210 -0.83208882 0.6646485 -0.682369844 0.51748259 11 0.7508775 -0.39482320 0.8851335 -0.095554151 0.66704225 12 0.9285353 0.04354098 1.1524578 0.453682964 0.30674141 13 1.9537886 1.83531480 1.9490471 1.865766458 -1.36629763 14 1.0385754 0.23709522 1.0081517 0.168761287 -0.60996496 15 0.8936027 -0.09915012 0.9407077 0.005928737 0.05273411 16 0.6633546 -0.62999115 0.7489752 -0.432539405 0.73371137 17 0.5941524 -1.07615247 0.6027536 -1.058618799 0.98045292 18 1.7516559 1.78554390 1.7275073 1.788745281 -1.75844255 19 1.4892405 1.27332749 1.5781245 1.468740378 -1.18728452 20 1.3794907 0.81105690 1.5238636 1.050442554 -0.07950879 [ reached 'max' / getOption("max.print") -- omitted 2412 rows ] outfit や infit というものは,𝑧ℎ よりもシンプルに予測値と実測値の間の差 𝑑𝑝𝑖 の二乗和 をもとにした適合度指標(Smith et al., 1995)です*22 。ということでこれらについても大きい ほど適合度が悪いと判断できます。実際にこれらの指標の相関を見てみると,いずれもかなり高 い相関になっていることがわかります。 指標間の相関 1 2 # z_hだけは方向が違う(小さいほど悪い)ので他の指標と負の相関 cor(personfit(result_grm)) 1 2 3 4 5 6 outfit z.outfit infit z.infit Zh outfit 1.0000000 0.9029153 0.9767846 0.8895739 -0.8593484 z.outfit 0.9029153 1.0000000 0.8962045 0.9898224 -0.8803762 infit 0.9767846 0.8962045 1.0000000 0.9083434 -0.8420095 z.infit 0.8895739 0.9898224 0.9083434 1.0000000 -0.8617816 Zh -0.8593484 -0.8803762 -0.8420095 -0.8617816 1.0000000 *22 outfit は重み付けない二乗和,infit は情報量(後述)で重み付けた二乗和のようです。 55
8.9 適合度 では実際に,当てはまりが悪い人の回答パタンを見てみましょう。 個人適合度が最も悪い人 1 2 3 4 1 2 3 4 # 引数stats.only = FALSEにすると,統計量に加えて回答パタンも残してくれる pf <- personfit(result_grm, stats.only = FALSE) # Zhの値が最小の人を見てみる pf[which.min(pf[, "Zh"]), ] Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 outfit z.outfit 1 1 1 1 1 6 6 6 6 6 5.804625 5.127515 infit z.infit Zh 581 5.552154 5.122718 -7.275079 581 いま扱っているモデルは 10 項目 1 因子のモデルです。そして 10 項目の識別力は,いずれも正 の値になっています。つまり,基本的にはすべての項目間には正の相関があるため,この回答者 のように「ある項目では 1,別の項目では 6」という回答の付け方には違和感があります。 個人適合度が悪い人 1 subset(pf, Zh < -1.96) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 outfit z.outfit 45 6 5 3 2 3 3 6 3 6 4 1.362101 0.8901956 75 6 6 3 3 1 6 6 5 6 1 2.777980 2.4882677 96 6 6 4 1 1 6 4 5 5 6 2.782816 2.5103934 98 1 1 1 1 2 3 6 3 4 3 2.006426 2.2242279 122 6 6 6 6 6 6 6 3 1 2 3.089210 2.3332034 123 6 6 6 5 6 6 1 3 3 2 2.065644 1.7906200 126 4 6 6 4 6 2 3 2 2 2 1.459869 1.1442225 172 5 2 1 6 3 6 6 6 6 6 4.695997 3.3697866 186 6 5 2 2 2 2 2 6 5 4 1.817030 1.8103247 226 1 1 1 4 6 5 6 5 6 6 4.784881 3.9067830 infit z.infit Zh 45 1.401457 0.9775542 -2.193202 75 2.751967 2.5542871 -3.279082 96 2.855421 2.6749604 -2.107685 98 2.099394 2.4577322 -3.172528 122 3.034305 2.3832158 -2.030530 123 1.955780 1.7068100 -1.971380 126 1.314690 0.8576219 -1.980124 172 3.651095 2.8513846 -3.188176 186 1.791534 1.7810377 -2.201609 226 4.276645 3.7176248 -4.057956 [ reached 'max' / getOption("max.print") -- omitted 121 rows ] 56
8.9 適合度
同様にして,𝑧ℎ < −1.96 の人たちを見てみると,項目ごとにかなりバラバラな回答をしてい
る人たちがあぶり出されています。ということで,𝑧ℎ 統計量によって怪しい回答者をあぶり出す
ことができそうだということがわかりました。
ちなみに反対に 𝑧ℎ > 1.96 の人たちを見たところですべての項目で同じような値をつけてお
り,やはり「適合度が悪い」という判断はなかなか難しそうです(ストレートライン気味な感じ
はありますが)
。
個人適合度が良すぎる? 人
1
subset(pf, Zh > 1.96)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10
outfit z.outfit
119
5
4
4
4
4
4
4
4
3
2 0.1059937 -3.933211
136
6
5
5
6
5
5
5
5
5
5 0.1749942 -2.032175
221
2
2
1
1
2
2
1
2
2
1 0.3892825 -1.326634
674
5
5
5
4
5
5
4
4
5
2 0.2410515 -2.270639
684
4
4
5
5
4
4
4
4
3
2 0.1893306 -3.039619
694
5
5
5
4
5
5
4
4
3
2 0.2993084 -2.123527
739
5
5
5
6
4
4
4
4
3
3 0.2791157 -2.191891
804
5
5
5
5
4
5
5
5
5
3 0.1527831 -2.557928
826
5
5
4
4
4
4
4
3
3
1 0.2959759 -2.472784
834
5
5
4
5
5
4
4
4
3
3 0.1765065 -2.919233
infit
z.infit
Zh
119 0.1122620 -3.899088 2.556666
136 0.2246921 -1.890473 2.065611
221 0.4313157 -1.407809 2.117805
674 0.2776608 -2.150842 1.970796
684 0.2021914 -2.982029 2.001422
694 0.3096588 -2.123376 2.017472
739 0.3059683 -2.111215 1.971686
804 0.1559106 -2.653563 1.963533
826 0.3014073 -2.463963 2.015447
834 0.1741386 -3.000314 2.105461
[ reached 'max' / getOption("max.print") -- omitted 9 rows ]
ここで紹介した 𝑧ℎ 統計量のような個人適合度は,モデルベースで適当回答者を除外するのに
使うこともできます。ただし,個人適合度はあくまでも(ロジスティックなどの IRT)モデルに
基づいて推定されたパラメータをもとに算出される点には少し注意が必要です。そもそもモデル
の式や因子構造の設定が間違っている場合には,正しくない検出をしてしまうかもしれません。
そんなわけで,もちろん適合度の値だけで機械的に除外するのは危ないですが,例えば項目パラ
メータの推定値がなんだかおかしい場合には,適合度指標を参考にしたデータクリーニングに
よって,推定をより安定させられるかもしれません。
57
8.9 適合度 8.9.3 項目適合度 IRT では項目一つ一つの性能(識別力や困難度)に関心があることが多いため,項目レベル でも適合度を考えます*23 。基本的な考え方としては,推定されたパラメータに基づく ICC が, データとどの程度一致するかを確認していきます。適合度が悪い項目が見つかった場合には,項 目の内容を吟味して削除するかを検討したり,場合によってはモデル自体を修正(因子数や因子 構造の修正など)する必要があるかもしれません。項目適合度は,具体的には以下のような手順 によって算出していきます(Orlando & Thissen, 2000)。 1. 項目パラメータと 𝜃𝑝 を推定する 2. 回答者を 𝜃𝑝 の値で並び替える 3. 回答者を 𝜃𝑝 の値によっていくつかのグループに分ける 4. 各グループでの反応確率を計算する 5. 4. で求めた反応確率と「モデル上の期待反応確率」を 𝜒2 的な統計量や視覚的方法を用い て比較する 視覚的なチェック ということで,まずは視覚的に比較してみましょう。mirt では,項目適合度を見るための関数 として itemfit() が用意されています。これに引数 empirical.plot を与えることでプロット を見ることができます。 項目適合度を視覚的にチェック 1 itemfit(result_2plm, empirical.plot = 8) 図 8.31 を見ると,各カテゴリごとに ICC が表示されています。これに重なっている小さいマ ルの Y 座標は,手順 3 で分割した各グループ(ここでは 10 グループ)における 𝑃 (𝑦𝑝𝑖 = 1) を 表しています(多値型の場合はカテゴリ選択確率が表示される)。また,X 座標はそのグループ における 𝜃𝑝 の平均値を表しています。ということで,各カテゴリにおいて,小さいマルと実線 が近いほどデータとモデルの当てはまりが良い,と判断することができるわけです。 ただカテゴリ数や項目数が多くなってくるとすべてを目視で確認するのはかなり大変です。ま た視覚的なチェックだけではどの項目が最も乖離しているのかを判断するのも大変でしょう。と いうことで,怪しい項目にアタリをつけるための統計量を使っていきましょう。 カイ二乗的な統計量 𝜒2 統計量は独立性の検証のところ(Section 8.9.1 )で登場したものです。IRT モデルに基づ いて項目および回答者のパラメータが推定されると,各項目について,指定されたサンプルサイ ズのうち何人がどの選択肢を選ぶと期待されるか,を計算することができます。これを期待度数 *23 もちろんモデル全体での適合度を考えることも可能です。例えば 1 パラメータと 2 パラメータの両モデルを lavaan で実行して,compareFit() することでどちらのモデルが良さそうかを判断したりというプロセスが考 えられます。 58
8.9 適合度 59 Empirical plot for item 8 1.0 P (q ) 0.8 0.6 0.4 0.2 0.0 -4 -2 0 2 4 q 図 8.31: 反応確率の比較プロット 𝐸,また実際にその項目に対する反応の分布を観測度数 𝑂 として 𝜒2 = ∑ (𝑂 − 𝐸)2 𝐸 (8.46) と計算すると,「帰無仮説:モデルがデータに適合している」が正しいならば小さい値におさま るはずだと考えられます。 Yen(1981)は,回答者を 𝜃𝑝 の値の高低によってグループ分けし,グループごとに 𝜒2 の計算 をすることにしました。特にグループ数を 10 に固定して算出した指標は 𝑄1 と呼ばれています。 グループ 𝑔 (𝑔 = 1, 2, ⋯ , 10) における特性値 𝜃𝑝̂ の平均値を 𝜃𝑔̄ とおくと,期待正答率は(2 パラ メータロジスティックモデルならば) 𝐸𝑖𝑔 = 1 (8.47) 1 + exp [−𝛼𝑖̂ (𝜃𝑔̄ − 𝛽𝑖̂ )] と求めることができます。すると,二値項目における 𝑄1 統計量は 2 10 (𝑂𝑖𝑔 − 𝐸𝑖𝑔 )2 [(1 − 𝑂𝑖𝑔 ) − (1 − 𝐸𝑖𝑔 )] 𝑄1 = ∑ 𝑁𝑔 ( + ) 𝐸𝑖𝑔 1 − 𝐸𝑖𝑔 𝑔=1 10 = ∑ 𝑁𝑔 ( 𝑔=1 10 = ∑ 𝑁𝑔 ( 𝑔=1 (1 − 𝐸𝑖𝑔 )(𝑂𝑖𝑔 − 𝐸𝑖𝑔 )2 + 𝐸𝑖𝑔 (𝑂𝑖𝑔 − 𝐸𝑖𝑔 )2 ) 𝐸𝑖𝑔 (1 − 𝐸𝑖𝑔 ) (8.48) (𝑂𝑖𝑔 − 𝐸𝑖𝑔 )2 ) 𝐸𝑖𝑔 (1 − 𝐸𝑖𝑔 ) となります。なお 1 行目の右辺は,第 1 項が正答セルにおける乖離度を表しています。同様に, 第 2 項は (1 − 𝑂𝑖𝑔 ), (1 − 𝐸𝑖𝑔 ) に置き換わっていることから,誤答セルにおける乖離度を表して いることがわかります。したがって,𝑄1 統計量は「グループ」×「項目 𝑖 の正誤」という 10 × 2 クロス表における 𝜒2 統計量を表しているのです。そしてこの 𝑄1 統計量は自由度が 10− パラ
8.9 適合度 60 メータ数の 𝜒2 分布に従うため,これを用いて検定ができるわけです。 同様に,差ではなく比を用いた統計量として 𝐺2 統計量(McKinley & Mills, 1985)というも のもあります。 10 𝐺2 = ∑ 𝑁𝑔 [𝑂𝑖𝑔 ln 𝑔=1 𝑂𝑖𝑔 1 − 𝑂𝑖𝑔 + (1 − 𝑂𝑖𝑔 ) ln ] 𝐸𝑖𝑔 1 − 𝐸𝑖𝑔 (8.49) ほかにもいくつか適合度指標は提案されているのですが,基本的には前述の通り「期待反応確 率と実際の反応確率の差」を何らかの形で統計量に落とし込んでいます。細かな設定としては 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 2 3 4 5 6 7 8 9 10 11 item X2 df.X2 RMSEA.X2 p.X2 G2 df.G2 RMSEA.G2 p.G2 1 Q1_1 375.251 8 0.137 0 458.736 6 0.176 0 2 Q1_2 73.800 8 0.058 0 99.817 5 0.088 0 3 Q1_3 92.664 8 0.066 0 135.273 5 0.104 0 4 Q1_4 102.569 8 0.070 0 153.430 5 0.111 0 5 Q1_5 126.263 8 0.078 0 169.111 5 0.116 0 6 Q1_6 157.177 8 0.088 0 191.087 5 0.124 0 7 Q1_7 151.358 8 0.086 0 219.612 4 0.149 0 8 Q1_8 125.809 8 0.078 0 154.268 6 0.101 0 9 Q1_9 151.495 8 0.086 0 207.812 5 0.129 0 10 Q1_10 572.082 8 0.170 0 754.816 5 0.248 0 G2 の自由度 ( df.G2) がいろいろな値になっていますがあまり気にしなくて OK です。今回 も例によって検定の結果 ( p.*) は全て有意となりました。このように,サンプルサイズが大き い場合これらの統計量はあまり参考にならない可能性が高いです。 …じゃあどうするんだ,という話になるわけですが,一つの方法としては効果量的な考え方を 𝑂 とるために「実際の差 (𝑂𝑖𝑔 − 𝐸𝑖𝑔 ) や比 𝐸𝑖𝑔 を計算する」という方法が考えられます。例えば 𝑖𝑔
8.10 項目情報量・テスト情報量 61 10 1 10 ∑𝑔=1 |𝑂𝑖𝑔 − 𝐸𝑖𝑔 | が大きい項目は,基本的に視覚的にプロットを見ても差があるでしょう。た だ (𝑂𝑖𝑔 − 𝐸𝑖𝑔 ) などを直接計算してくれる関数は用意されていないので,多少頑張る必要があり ます。 8.10 項目情報量・テスト情報量 IRT にかぎらず,統計的推測では推測の精度が問われます。例えば正規分布の母平均の推測の 場合,標準誤差は母分散 𝜎2 およびサンプルサイズ 𝑛 を用いて √ 𝜎𝑛 と表されました。この式は, 2 データの数が多いほど推測の精度が高くなるという自然な考えを表しています。 テストや心理尺度では,標準誤差を真値と誤差の関係から考えます。Chapter 4 で紹介しまし たが,潜在特性の測定(古典的テスト理論)では,観測得点 𝑥 と真値 𝑡 と誤差 𝑒 の関係を 𝑥=𝑡+𝑒 (8.50) 𝜎𝑒2 𝜎𝑥2 (8.51) と表し,信頼性係数を 𝜌 =1− と定義していました。標準誤差は「真値が 𝑡 の人が測定を繰り返した場合に, (平均的に)どの程 度観測得点 𝑥 がばらつくか」を表しています。したがって,𝑥 ∼ 𝑁 (𝑡, 𝜎𝑒2 ) とすると,心理測定に おける標準誤差は 𝜎𝑒 にあたります。信頼性係数 𝜌 に関しては,上の式を変形させると 𝜎𝑒 = 𝜎𝑥2 √1 − 𝜌2 (8.52) となることから,これが間接的に測定の精度を表す指標だということが言えます。なお実際に標 準誤差を計算する場合には,𝜌 の推定値として 𝛼 係数などを代入します。 さて,信頼性係数に基づく上記の標準誤差の式には真値 𝑡 は含まれていません。つまり真値の 値が何であれ標準誤差は同じということを意味しています。ですがこの考え方は割と不自然なも のです。特性値の真値が平均くらいの人では,尺度の項目はほとんどすべてがその人の特性値を 測定するのに意味のある(正解であるほど真値 𝑡 は高いだろうという判断につながる)項目です。 一方で,特性値の真値が非常に高い人では,実はたいていの項目はその人の特性値を測定するの にほぼ無意味です。というのも,ある程度難易度の高い項目に正解する人は,難易度の低い項目 にはほぼ確実に正解するだろうと言えるためです。つまり, 「ある程度難易度の高い項目」の情報 さえあれば「難易度の低い項目」の情報は無くても特に問題ない(真値の推定には影響しない) と言えます。先ほど,標準誤差はサンプルサイズが多くなるほど小さくなる,ということを説明 しました。本来,特性値の真値によって「その人の特性値を測定するのに意味のある」項目の数 は変わるはずだと考えると,標準誤差も個人(の特性値の値)ごとに変わると考えるのが自然な 気がしてきます。 IRT では,この問題を克服するために項目情報量 (item information) およびテスト情報量 (test information) というものを用います。この項目情報量・テスト情報量は,後述するテス トの設計などにおいても大きな力を発揮してくれます。
8.10 項目情報量・テスト情報量 8.10.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 倍(1/52 から 1/4)になっていました。 ここまでの話を IRT に置き換えてみます。𝜃𝑝 の予測に際して,項目がもつ情報量を考えてみ ましょう。図 8.32 には困難度がそれぞれ (0, 3) の 2 つの項目の ICC が描かれています。 ある回答者 𝑝 は,それまでの項目への解答状況から 𝜃𝑝̂ = 0 だと予測されているとします。す ると,この回答者 𝑝 にとって,緑の項目(𝛽𝑖 = 0)は正解するかどうかがかなり不確実(50%) です。そしてこの項目への反応によって 𝜃𝑝̂ = 0 の予測は変化し,多少なりともより正確なもの になるでしょう。そう考えると,緑の項目は 𝜃𝑝̂ = 0 の人にとってそれなりの情報量を持ってい ると考えられます。 62
8.10 項目情報量・テスト情報量 63 P(ypi = 1) bi = 0 bi = 3 qp 0 図 8.32: IRT における情報の考え方 一方で赤い項目(𝛽𝑖 = 3)を見ると,𝜃𝑝̂ = 0 の人はほぼ確実に正解しないだろうという予想が つきます。つまり,出題時点で赤い項目に対する不確実性はほぼゼロ,ということです。実際に この回答者 𝑝 が赤い項目に不正解だったとしても,「そうでしょうね」という話なので 𝜃𝑝 の予 測に対して特に新規の情報をもたらすものとは言えません*24 。同様に,緑の項目でも 𝜃𝑝̂ = 3 の 人にとってはほぼ情報量のない項目,ということになります。このように,ある項目は,その項 目の困難度が回答者の 𝜃 に対してちょうどよいほど,𝜃 の推定に対して多くの情報量を持ってい るのです。ということである項目が持つ項目情報量は,実際には 𝜃 の値によって変動する関数に なっています。この関数のことを項目情報関数 (item information function [IIF]) と呼びま す*25 。この節のはじめにお話したように,IRT では真値によって標準誤差が変わるようにする ことを考えます。この後具体的な関数の形を示しますが,項目情報関数を 𝜃 の関数として表すこ とによってその目的を達成しようというわけです。 8.10.2 項目情報関数 IRT における情報量の考え方のイメージを確認したので,具体的な項目情報関数の式を見てみ ましょう。二値型モデルの場合,項目情報関数 𝐼𝑖 (𝜃) は 𝐼𝑖 (𝜃) = 𝑃 ′ (𝑦𝑖 = 1)2 𝑃 (𝑦𝑖 = 1)(1 − 𝑃 (𝑦𝑖 = 1)) (8.53) という式になります。ここで 𝑃 ′ (𝑦𝑖 = 1) は,𝑃 (𝑦𝑖 = 1) を 𝜃 で微分した導関数です。…といっ てもよくわからないと思うので,実際に見てみましょう。図 8.33 は,図 8.32 に示した ICC を 持つ 2 項目のロジスティックモデルにおける項目情報関数です。どちらの線も先程説明したよう *24 もちろんこの回答者 𝑝 が赤い項目に正解した場合には,𝜃 𝑝 の予測を大きく動かすため,情報量が多くなります。 ですがそもそも「回答者 𝑝 が赤い項目に正解する」という事象の発生確率がほぼゼロなので,「この項目への回 答が 𝜃𝑝 の予測にもたらす情報量の期待値」という観点ではやはりほぼゼロなのです。 *25 項目情報「量」という場合には特定の 𝜃 における値を指す一方で,項目情報「関数」という場合には 𝜃 全域を対 象とすることが多い気がするので,使い分けに気をつけてください。
8.10 項目情報量・テスト情報量 64 に,正解するかが最も不確実な 𝜃𝑝 = 𝛽𝑖 のところで情報量が最大になっています。 Ij(q ) bi = 0 bi = 3 qp 0 3 図 8.33: 項目情報関数 IIF の式 𝐼𝑖 (𝜃) = 𝑃 ′ (𝑦𝑖 = 1)2 𝑃 (𝑦𝑖 = 1)(1 − 𝑃 (𝑦𝑖 = 1)) (8.54) を求めるためには,𝑃 (𝑦𝑖 = 1) の導関数を求める必要があります。正規累積モデルは項目反応関 数の中に積分が含まれているためこれを求めるのはかなり大変なのですが,それと比べるとロジ スティックモデルでは解析的に IIF の式を求めることができます。例えば 2 パラメータロジス ティックモデルの場合,IIF は 𝐼𝑖 (𝜃) = 𝛼2𝑖 𝑃 (𝑦𝑖 = 1)(1 − 𝑃 (𝑦𝑖 = 1)) (8.55) となります。この式からも,IIF は 𝑃 (𝑦𝑖 = 1) = 0.5 となる 𝜃 = 𝛽𝑖 の点において最大値をとる ことがわかります。同時に,識別力 𝛼𝑖 が大きいほど項目情報量も多くなることがわかります。 図 8.34 に,識別力が異なる 2 つの項目 (𝛽𝑖 = 0; ICC は 図 8.8 ) のロジスティックモデルでの IIF を示しました。識別力が高い項目ほど,項目特性関数の傾きが大きくなっていました。その ため 𝜃𝑝 が高い人は正解し 𝜃𝑝 が低い人は不正解する,というメリハリがついていました。𝜃𝑝 の 推定という視点から言えば,𝛼𝑖 が大きい項目ほどその項目への正誤が特性値の推定により確か な情報を与えるという意味で,𝛼 と項目情報量には密接な関係があるわけです。 しかし改めて 図 8.34 をよく見ると,𝜃𝑝 = 0 付近では識別力の高い緑の項目のほうが高い情報 量を持っている一方で,𝜃𝑝 が 0 から離れたところではむしろ識別力の低い青い項目のほうが高 い情報量を持っています。そして対応する 図 8.8 を見ると,識別力 𝛼𝑖 が高くなるほど,𝜃𝑝 が 𝛽𝑖 から離れると急激に 𝑃 (𝑦𝑖 = 1) が 0.5 から離れていることが分かります。つまり「正解するかが 不確実」な(≒情報量を持ちうる)𝜃𝑝 の範囲は識別力の高さと反比例しており,その結果 𝜃𝑝 が 𝛽𝑖 から離れたところでの項目情報量も少なくなってしまうのです。 この事実は「識別力は高ければ高いほど良いわけではない」ことを教えてくれています。極端 な例として「𝜃𝑝 < 𝛽𝑖 の人は確実に間違えるが,𝜃𝑝 ≥ 𝛽𝑖 の人は確実に正解する」項目を考えてみ ましょう。その項目は識別力 𝛼𝑖 → ∞ であり,図 8.35 のような ICC をもつ項目です。
8.10 項目情報量・テスト情報量 65 Ij(q ) ai = 1 ai = 0.5 qp 0 図 8.34: 異なる識別力を持つ項目の IIF P(ypi = 1) ai ®¥ qp bi 図 8.35: 完全な識別力を持つ項目があったら この項目は「𝜃𝑝 が 𝛽𝑖 より上か下か」に関しては完全な情報を教えてくれます (𝐼𝑖 (𝜃𝑝 = 𝛽𝑖 ) = ∞)。その一方で 𝜃𝑝 > 𝛽𝑖 の範囲では, 𝜃𝑝 が 𝛽𝑖 よりどの程度高かろうと 𝑃 (𝑦𝑖 = 1) は 100%(逆もまた然り)なので,「𝜃𝑝 が 𝛽𝑖 よりどの程度上か下か」に関しては全くわかりません (𝐼𝑖 (𝜃𝑝 ≠ 𝛽𝑖 ) ≃ 0)。このように,識別力が高すぎる項目は有効な 𝜃𝑝 の範囲が狭すぎるため,多く の回答者にとって意味のない項目になってしまう可能性があるわけです。 INFO 多値型モデルでの項目情報量 (いつか書きたいという気持ちだけは常に持っている) 8.10.3 テスト情報関数 ここまで,項目情報関数の定義およびその性質を見てきました。この「情報」は,𝜃𝑝 を推定す るための「情報」を表しているわけですが,𝜃𝑝 の推定は尺度・テスト内のすべての項目への回答
8.10 項目情報量・テスト情報量 66 を総合して行われるものです。そのため,次は項目情報関数を拡張して尺度全体での情報関数を 考えてみます。 情報量には「独立な情報の情報量は加法的である」という性質がありました。IRT でもこの性 質は健在です。もしテスト全体が局所独立の仮定を満たしているならば,テスト情報関数 (test information function [TIF]) は単純に 𝐼 𝐼(𝜃) = ∑ 𝐼𝑖 (𝜃) (8.56) 𝑖=1 と表せます。そしてテスト情報量の重要な性質として,真の特性値が 𝜃𝑝 の人における最尤推定 量 𝜃𝑝̂ の誤差分散は漸近的に 𝜎(2𝜃 ̂ |𝜃 ) = 𝑝 𝑝 1 𝐼(𝜃𝑝 ) (8.57) となることが知られています。ということで,標準誤差は 𝜎(𝜃 ̂ |𝜃 ) = 𝑝 𝑝 1 √𝐼(𝜃𝑝 ) (8.58) となります。これは古典的テスト理論の表記に合わせると,最尤推定値 𝜃𝑝̂ = 𝜃𝑝 + 𝜀𝑝 とした場合 の 𝜎𝜀 に相当するものです。このように,IRT では測定の精度を表す標準誤差を,𝜃𝑝 の関数とし て求めることができるわけです。 8.10.4 R で情報量を確認する 項目・テスト情報関数の考え方がわかったところで,実際の項目での関数を確認してみましょ う。これまで ICC を出すために使っていた plot() 関数によって,IIF も出すことができます。 項目情報関数を出す場合には type = "infotrace"とします。 項目情報関数 1 plot(result_2plm, type = "infotrace", facet_items = FALSE) 先ほど説明した通り,低識別力の項目 1 は他の項目と比べて圧倒的に情報量が少ないことがわ かります。また、𝜃𝑝 > 2 のエリアでは、Q1_10 以外のほぼどの項目も情報量が無いため、この 10 項目では特性値が高い人についてはまともな推定はできなさそうだということがわかります。こ の項目情報関数を見れば,例えばどうしても 1 項目削除しないといけないとしたら,基本的には Q1_1 を削除したら良いが,もしも 𝜃𝑝 = −6 付近がメインターゲットだとしたらむしろ Q1_1 は 残したほうが良いとか、𝜃𝑝 > 2 がメインターゲットならば現状の項目セットではどうしようもな いので新規項目を追加する必要がある、というように状況に応じた決定をサポートしてくれるで しょう。 具体的に特定の 𝜃𝑝 の値における項目情報量を出したい場合には iteminfo() 関数を使います。
8.10 項目情報量・テスト情報量 67 Item Information Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 0.4 I (q ) 0.3 0.2 0.1 0.0 -6 -4 -2 0 2 4 6 q 図 8.36: 項目情報関数 項目情報量 1 2 3 4 5 # はじめに項目パラメータの情報だけ抽出する必要がある # 1番目の項目を抽出するなら item1 <- extract.item(result_grm, item = 1) # 引数Thetaには,どの値での特性値を出したいかを指定する iteminfo(item1, Theta = c(-4, -2, 0, 2, 4)) 1 [1] 0.10354648 0.10516020 0.10114364 0.08337954 0.04723855 続いてテスト情報関数です。こちらは引数を type = "info"とすることで表示できます。 テスト情報関数 1 2 # 引数typeは同じ # type = "infotrace"とすると前述の項目情報関数になるので注意 3 plot(result_2plm, type = "info") 尺度(10 項目)全体を見ても,項目困難度がまんべんなく存在している-3 から 0 付近の情報 量が高い一方で,すべての項目の困難度が 0 以上なので 𝜃𝑝 が正のエリアでは情報量がさほど高 くないことがわかります。 プロットする際に引数 type を "infoSE"にすると,標準誤差を合わせて出してくれるように なります。
8.10 項目情報量・テスト情報量 68 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.37: テスト情報関数 テスト情報関数に標準誤差を添えて # 引数typeは同じ plot(result_2plm,type="infoSE") Test Information and Standard Errors 12 2.5 8 1.5 6 1.0 4 0.5 2 0.0 -6 -4 -2 0 2 4 q 図 8.38: テスト情報関数と標準誤差 6 SE(q ) 10 2.0 I (q ) 1 2
8.11 (おまけ)多母集団モデルと特異項目機能 Exclamation-Triangle うまくいかない場合 もしかしたら,以下のようなエラーが出てプロットが表示されないかもしれません。 1 Error in eval(expr, envir, enclos): latticeExtra package is not available. Please install. その場合は,指示に従って install.packages("latticeExtra") をしてから再度実行 してください。 1 最も情報量が高い 𝜃𝑝 = −1.4 付近では,標準誤差はだいたい √2.82 ≃ 0.60 くらいになってい ます。仮に 𝜃𝑝 の真値が −1.4 の人が大量にいた場合,その人達の推定値は平均で 0.6 くらいは上 下するだろう,ということです。裏を返せばその程度の精度でしか推定できていないならば,例 えば 𝜃𝑝̂ = −1.3 の人と 𝜃𝑝̂ = −1.5 の人がいたとしても,前者のほうが真の特性値が高いとはな かなか言い切れないわけですね。 具体的に特定の 𝜃𝑝 の値におけるテスト情報量を出したい場合には testinfo() 関数を使い ます。 テスト情報量 1 2 # 引数Thetaには,どの値での特性値を出したいかを指定する 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 と非常に大きくなって います。つまりこの 10 項目で推定を行っても,𝜃 の真値が 4 の人に対する推定値は平均で ±5 く らいには変動するわけです。そう考えると,この 10 項目ではどうあがいても 𝜃𝑝 が高い層におい てまともな推定値は得られないということが分かります。 8.11 (おまけ)多母集団モデルと特異項目機能 Chapter 7 の 多母集団同時分析のところでは,多母集団同時分析という枠組みを紹介しまし た。当然ながら IRT においても,多母集団同時分析を実行することができます。そもそも SEM や IRT において多母集団モデルを行う理由が何だったのかを思い出してみると,その一つには 69
8.11 (おまけ)多母集団モデルと特異項目機能 「集団ごとに異なる回答傾向の原因を探る」という点がありました。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 を指定しましたが,その時とは指定の仕方が異なるのでご注意ください。 多母集団 IRT モデル(用途 1) 1 2 # modelに1を与えておけば普通の一次元モデルになる # 引数groupにはcharacter型かfactor型のベクトルを与える必要がある 3 result_group <- multipleGroup(dat_bin, model = 1, itemtype = "2PL", group = as.character(dat$gender)) 結果を見るための関数は概ね mirt() のときと同じように使えます。とりあえず推定された項 目パラメータを見てみましょう。 70
8.11 (おまけ)多母集団モデルと特異項目機能 項目パラメータのチェック 1 1 coef(result_group, simplify = TRUE, IRTpars = TRUE) 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 $`1` $items a b g u Q1_1 0.447 -2.061 0 1 Q1_2 1.057 -1.751 0 1 Q1_3 0.971 -1.545 0 1 Q1_4 1.169 -1.271 0 1 Q1_5 1.073 -1.385 0 1 Q1_6 1.249 -1.455 0 1 Q1_7 1.281 -1.086 0 1 Q1_8 1.273 -1.111 0 1 Q1_9 1.045 -0.930 0 1 Q1_10 1.156 0.203 0 1 $means F1 0 $cov F1 F1 1 $`2` $items a b g u Q1_1 0.294 -4.870 0 1 Q1_2 0.853 -3.029 0 1 Q1_3 0.962 -2.178 0 1 Q1_4 0.858 -2.118 0 1 Q1_5 0.871 -2.184 0 1 Q1_6 1.234 -1.613 0 1 Q1_7 1.449 -1.268 0 1 Q1_8 1.228 -1.369 0 1 Q1_9 1.574 -1.041 0 1 Q1_10 1.322 -0.127 0 1 $means F1 0 $cov 71
8.11 (おまけ)多母集団モデルと特異項目機能 43 44 F1 72 F1 1 結果を見ると,$`1`と $`2`という 2 つに分かれて,それぞれ項目パラメータの推定値が出て います。これは,引数 group で指定したグループの値ごとに項目パラメータをそれぞれ推定し た結果です。 (2 つのグループで 𝜃 の分布が同じという前提のもとで)もしも推定された項目パラ メータが大きく異なる場合には,DIF が発生していると考えられるわけです。また,各グループ において $means および $cov はそれぞれ 0,1 となっています。multipleGroup() 関数は,デ フォルトではこのように「各グループの 𝜃 の母集団分布がそれぞれ標準正規分布である」という 設定のもとで項目パラメータを推定します。 特性値 𝜃 の分布のグループごとの差異を見たい場合には, SEM の多母集団モデルと同じ要 領で等値制約を表す引数を与える必要があります。lavaan での多母集団モデルを実行したとき には,因子負荷の等値制約などを引数 group.equal で指定しました。multipleGroup() では, invariance という引数が用意されています。この引数には,表 8.4 に示す値を入れることが出 来ます。free_means および free_vars という値があるということは,これらを設定しない限 り,全てのグループの 𝜃 は標準正規分布に固定されてしまうということなのでご注意ください。 表 8.4: 引数 invariance に指定できるもの 指定 制約されるもの free_means グループ 1 以外の 𝜃 の平均値を自由推定する free_vars グループ 1 以外の 𝜃 の分散を自由推定する slopes 項目の識別力パラメータをグループ間で同じとする intercepts 項目の切片パラメータをグループ間で同じとする 多母集団 IRT モデル(用途 2) 1 2 # 引数invarianceで等値制約を指定する # "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 1 2 coef(result_group2, simplify = TRUE, IRTpars = TRUE) $`1` $items
8.11 (おまけ)多母集団モデルと特異項目機能 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 a b g u Q1_1 0.445 -2.584 0 1 Q1_2 1.072 -1.978 0 1 Q1_3 1.067 -1.555 0 1 Q1_4 1.041 -1.414 0 1 Q1_5 1.038 -1.482 0 1 Q1_6 1.133 -1.396 0 1 Q1_7 1.294 -1.003 0 1 Q1_8 1.162 -1.078 0 1 Q1_9 1.272 -0.792 0 1 Q1_10 1.208 0.242 0 1 $means F1 0 $cov F1 F1 1 $`2` $items a b g u Q1_1 0.445 -2.584 0 1 Q1_2 1.072 -1.978 0 1 Q1_3 1.067 -1.555 0 1 Q1_4 1.041 -1.414 0 1 Q1_5 1.038 -1.482 0 1 Q1_6 1.133 -1.396 0 1 Q1_7 1.294 -1.003 0 1 Q1_8 1.162 -1.078 0 1 Q1_9 1.272 -0.792 0 1 Q1_10 1.208 0.242 0 1 $means F1 0.4 $cov F1 F1 1.013 確かに項目パラメータの推定値がグループ間で同じになり,$`2`の $means および $cov の値 が変わりました。この結果は(2 つのグループで DIF が発生している項目が無いという前提の もとで)グループ 2 のほうが 𝜃 の平均値が高いことを示しているわけです。 73
8.11 (おまけ)多母集団モデルと特異項目機能 74 表 8.5: 第 𝑘 層 (𝑘 = 1, 2, ⋯ , 𝐾) におけるクロス表 項目 𝑗 項目 𝑖 当てはまらない 当てはまる 計 グループ A 𝑥𝐴0 (𝑘) 𝑥𝐴1 (𝑘) 𝑛𝐴 グループ B 𝑥𝐵0 (𝑘) 𝑥𝐵1 (𝑘) 𝑛𝐵 計 𝑛0 (𝑘) 𝑛1 (𝑘) 𝑁 (𝑘) (𝑘) 8.11.1 DIF の評価 ここからは,DIF が発生しているかを評価するいくつかの方法を紹介します。とりあえず二値 IRT モデルについての方法を紹介していきますが,多くの方法はすでに多値モデルに対する拡張 も提案されているはずなので,興味があれば探してみてください。 ノンパラメトリックな方法 DIF の検出法は,大きく分けると IRT モデルに基づく,つまり項目パラメータの推定値によっ て行われるパラメトリックな方法と,IRT を使用しないノンパラメトリックな方法に分けられま す。ノンパラメトリックな方法として最も有名なのが,2 × 2 クロス表に対する 𝜒2 検定を拡張し た Mantel-Haenszel 検定(マンテル・ヘンツェル検定: Mantel & Haenszel, 1959)に基づく方 法です。手法自体は非常に古いものですが,その使いやすさから近年でも DIF 検出の方法とし て非常に多く用いられているようです(Berrío et al., 2020) 。 MH 検定では,まず各集団を何らかの変数の値によって 𝐾 層に層化します。DIF 検出の場合 には,IRT で推定された特性値 𝜃 や合計点などによって層化するのが一般的です。表 8.5 は,そ うして層化されたうちの第 𝑘 層について,ある二値項目に対する回答をグループごとに集計した ものです。MH 検定では,表 8.5 のような 2 × 2 クロス表が 𝐾 個あるときに,それらを全部ひっ くるめて独立性の検定を行う手法と言えます。 まずは一つの層から考えることで,具体的な方法を紹介していきましょう。例えば 𝜃 が最も低 い人達が集められた第 1 層において,グループ A の人たちのほうが「当てはまる」と回答した割 合が高かったとします。このように,特定の層に絞ったときに回答傾向に違いが見られるという のも立派な DIF です。そして DIF が発生している場合には,2 つのグループにおける回答の割 合つまりオッズ比 OR(𝑘) = (𝑘) (𝑘) (𝑘) (𝑘) 𝑥𝐴1 /𝑛𝐴 (8.59) 𝑥𝐵1 /𝑛𝐵 が 1 ではない,という見方ができます。MH 検定では,このオッズ比を全ての層について統合し たもの OR (MH) = 𝐾 (𝑘) (𝑘) 𝐾 (𝑘) (𝑘) ∑𝑘=1 𝑥𝐴1 /𝑛𝐴 ∑𝑘=1 𝑥𝐵1 /𝑛𝐵 (8.60)
8.11 (おまけ)多母集団モデルと特異項目機能 75 が 1 であるかどうかを検定します。 MH 検定の結果が有意であれば,データ全体として一方のグループのほうが「当てはまる」と 回答した割合が高かったと言えます。ICC 的に言えば 図 8.39 のような状況であると言えるわけ です。 P(ypi = 1) Group A Group B qp 図 8.39: MH 検定で検出可能な(均一)DIF 一方で,MH 検定では検出できない DIF も存在します。図 8.40 に示されている DIF では,𝜃𝑝 の低い層ではグループ A のほうが「当てはまる」割合が高い一方で,𝜃𝑝 の高い層では逆転が起 こっています。OR(𝑘) を一つ一つ見れば 1 ではないわけですが,これを統合した OR(MH) は 1 に なってしまうわけです。MH 検定による方法は,図 8.39 に示した均一 DIF は検出可能な一方で, 図 8.40 のような不均一 DIF は検出できないかもしれない方法であることに気をつけましょう。 P(ypi = 1) Group A Group B qp 図 8.40: 不均一 DIF 項目パラメータの推定値を比べる方法 IRT に基づくパラメトリックな方法として,まずは項目パラメータの推定値の差を直接検定す る方法を紹介します。Wald 検定は,パラメータの推定値に対して用いられる検定の一種です。 一般的に用いられる一変量 Wald 検定では,パラメータ 𝜉 の推定値 𝜉 ̂ に対して 𝜉 ̂ − 𝜉0 𝑊 =( ) 𝜎𝜉 ̂ 2 (8.61)
8.11 (おまけ)多母集団モデルと特異項目機能 76 という検定統計量を考えます。ただし 𝜉0 は帰無仮説における値,𝜎𝜉 ̂ は 𝜉 ̂ の標準誤差です。検定 統計量 𝑊 は自由度 1 の 𝜒2 分布に従うことが知られているため,これを用いて「𝜉 ̂ が 𝜉0 である」 という帰無仮説を検定します。 Lord(1980)はこの Wald 検定によって DIF を検出する方法を提案しました。DIF が発生して いる場合には,グループごとに推定された項目パラメータに差があると考えられるため,(8.61) 式の中身を「2 つのグループでの推定値」に置き換えることにより,例えば項目 𝑖 の困難度 𝛽𝑖 で あれば (𝐴) (𝐵) 𝛽 − 𝛽𝑖 𝑊=( 𝑖 ) 𝜎𝛽(𝐴) −𝛽(𝐵) 𝑖 =( 2 (8.62) 𝑖 (𝐴) (𝐵) 𝛽𝑖 − 𝛽 𝑖 𝜎𝛽(𝐴) + 𝜎𝛽(𝐵) 𝑖 2 ) 𝑖 という形で,同じように Wald 検定を実行することができます。 Wald 検定で用いる Wald 統計量 𝑊 は「推定値の差」を「標準誤差」で割った形をしているた め,パラメータの数が複数であっても同じように定義することができます。例えば項目 𝑖 の識別 力と困難度を並べたベクトル 𝝃𝑖 = [𝛼𝑖 𝛽𝑖 ] に対して (8.63) 𝑊 = 𝝃𝑖 𝚺−1 𝝃𝑖⊤ という統計量(𝚺 は標準誤差行列)を考えてやれば,これが自由度 2(=パラメータの数)の 𝜒2 分布に従うため,同じようにして検定ができるのです。 項目特性曲線を比べる方法 IRT では,推定された項目パラメータをもとに項目特性曲線(ICC)を描くことが出来ました。 ICC は,いわば「得点の期待値」を表したグラフです。したがって「グループ間で ICC のずれ が大きい」場合には DIF が発生していると考えることが出来そうです(Raju, 1988) 。ICC のず れとは,図 8.41 でオレンジ色に塗りつぶされた部分の面積のことです。これが大きい場合には, 確かに DIF が発生していると言えそうです。 P(ypi = 1) Group A Group B qp 図 8.41: グループ間での ICC のズレ
8.11 (おまけ)多母集団モデルと特異項目機能 77 グループ A における ICC(正答確率)は,2 パラメータロジスティックモデルならば (𝐴) 𝑃𝑖 (𝜃) = 𝑃 (𝑦𝑝𝑖 = 1|𝜃) = 1 (𝐴) (𝐴) 1 + exp[𝛼𝑖 (𝜃 − 𝛽𝑖 )] (8.64) と表せるため,ICC のズレの面積に基づく DIF の指標として,例えば DIF𝑖 = ∫ ∞ (𝐴) ∣𝑃𝑖 (𝐵) (𝜃) − 𝑃𝑖 (𝜃)∣ 𝑓(𝜃)𝑑𝜃 (8.65) −∞ というものを考えることができるわけです(e.g., Chalmers, 2018; Raju et al., 1995) 。ただし 𝑓(𝜃) は 𝜃 の母集団分布(普通は標準正規分布)を指しています。 R でやってみる DIF の検出法は相当たくさんあります(レビューとして Berrío et al., 2020; Kleinman & Teresi, 2016 など)。ということで検出法の紹介はこのあたりでおしまいにして,mirt パッケー ジに実装されている DIF 検出法を試してみましょう*26 。これまでに紹介した方法と異なるので すが,DIF() 関数を使うと,モデル比較および尤度比検定の観点から DIF の有無を評価してく れます。DIF() 関数のデフォルトの挙動では,第一引数で与えたモデルをベースラインに特定の 1 項目のみに等値制約をかけた(=その項目には DIF が無いと仮定した)ときの変化を計算し ます。ここで等値制約がかかるのは引数 which.par で指定したパラメータです。したがって,2 パラメータモデルの場合は which.par = c("a1", "d") と指定することで「その項目の全ての 項目パラメータの DIF を丸ごと評価する」ことを意味します。(あまり無いケースですが)例え ば識別力には興味がなく,困難度(切片)の DIF だけに関心がある場合には which.par = "d" としたら良いわけです。 DIF の検出(1) 1 DIF(result_group, which.par = c("a1", "d")) 1 2 3 4 5 6 7 8 9 10 11 groups converged AIC SABIC HQ BIC X2 df p Q1_1 1,2 TRUE -24.578 -19.339 -20.363 -12.985 28.578 2 0 Q1_2 1,2 TRUE -32.846 -27.608 -28.632 -21.253 36.846 2 0 Q1_3 1,2 TRUE -18.162 -12.924 -13.948 -6.569 22.162 2 0 Q1_4 1,2 TRUE -13.108 -7.870 -8.894 -1.515 17.108 2 0 Q1_5 1,2 TRUE -13.712 -8.473 -9.497 -2.119 17.712 2 0 Q1_6 1,2 TRUE 2.217 7.456 6.432 13.810 1.783 2 0.41 Q1_7 1,2 TRUE -4.187 1.052 0.028 7.406 8.187 2 0.017 Q1_8 1,2 TRUE -1.149 4.090 3.066 10.444 5.149 2 0.076 Q1_9 1,2 TRUE -15.366 -10.128 -11.151 -3.773 19.366 2 0 Q1_10 1,2 TRUE -9.558 -4.320 -5.344 2.035 13.558 2 0.001 *26 MH 検定は mirt パッケージには用意されていないようです。 R では difR::difMH() などの関数によって MH 検定を行うことができます。
8.11 (おまけ)多母集団モデルと特異項目機能 78 結果を見てみると,SEM の多母集団モデルのときにも登場した AIC や BIC などの情報量規準 および X2(𝜒2 )統計量とこれに関連した自由度・𝑝 値が表示されています。例えば Q1_1 行の AIC の値 -24.578 というのは,「 Q1_1 のみグループごとに項目パラメータが同じ= DIF が無 いという制約をかけることで,AIC が 24.578 悪化した」ということを意味します。したがって, AIC のみで判断するならば Q1_1 には DIF が生じている,と言えるわけです。同様に,𝜒2 検定 の結果も「その項目には DIF が発生していないとする(グループ間で等値制約を置く)ことで モデル適合度がどの程度悪化するか」という観点に基づく尤度比検定が行われています*27 。 実際に DIF がありそうと判断された項目について ICC を見るには,引数 plotdif = TRUE を与えてください。 DIF がありそうな項目の ICC 1 # デフォルトではSABICがマイナスになっている項目のプロットを表示する 2 DIF(result_group, which.par = c("a1", "d"), plotdif = TRUE) 1 2 3 4 5 6 7 8 9 10 11 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 groups converged AIC SABIC HQ BIC X2 df p 1,2 TRUE -24.578 -19.339 -20.363 -12.985 28.578 2 0 1,2 TRUE -32.846 -27.608 -28.632 -21.253 36.846 2 0 1,2 TRUE -18.162 -12.924 -13.948 -6.569 22.162 2 0 1,2 TRUE -13.108 -7.870 -8.894 -1.515 17.108 2 0 1,2 TRUE -13.712 -8.473 -9.497 -2.119 17.712 2 0 1,2 TRUE 2.217 7.456 6.432 13.810 1.783 2 0.41 1,2 TRUE -4.187 1.052 0.028 7.406 8.187 2 0.017 1,2 TRUE -1.149 4.090 3.066 10.444 5.149 2 0.076 1,2 TRUE -15.366 -10.128 -11.151 -3.773 19.366 2 0 1,2 TRUE -9.558 -4.320 -5.344 2.035 13.558 2 0.001 確かに ICC にはそこそこのずれがありそうです。 「項目 1 つ 1 つについて DIF の有無をチェックする」という意味では,全く逆方向からのアプ ローチを取ることも可能です。 パラメータ推定時に invariance = c("free_means","free_vars","slopes","intercepts") を与えると,全ての項目パラメータが同じ,すなわち全ての項目に DIF が無いことを仮定した上 での推定を行うことが出来ました。このときの結果( result_group2)をベースラインとして 考えると,DIF() 関数には特定の 1 項目のみから等値制約を外した(=その項目だけ DIF があ ると仮定した)ときの変化を計算してもらいたいところです。これを実現するためには,DIF() 関数の引数 scheme を指定する必要があります。 *27 項目の数だけ検定を繰り返しているので,実はここで表示されている結果には多重検定の問題が生じています。 これを解消するためには,引数 p.adjust を指定すると良いです。例えば p.adjust = "holm"とすると Holm の方法によって調整された 𝑝 値を表示してくれるようになります。
8.11 (おまけ)多母集団モデルと特異項目機能 79 Item Probability Functions -6 -20 2 4 6 P (q ) Q1_5 Q1_9 Q1_10 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 -20 2 4 6 cat1:1 cat1:2 -6 -20 2 4 6 q 図 8.42: DIF がありそうな項目の ICC DIF の検出(2) 1 # 等値制約を"drop"して検証してほしい 2 DIF(result_group2, which.par = c("a1", "d"), scheme = "drop") 1 2 3 4 5 6 7 8 9 10 11 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 groups converged AIC SABIC HQ BIC X2 df p 1,2 TRUE -10.691 -5.452 -6.476 0.902 14.691 2 0.001 1,2 TRUE -8.581 -3.343 -4.367 3.012 12.581 2 0.002 1,2 TRUE -0.172 5.066 4.042 11.420 4.172 2 0.124 1,2 TRUE 1.756 6.995 5.971 13.349 2.244 2 0.326 1,2 TRUE 2.772 8.010 6.986 14.365 1.228 2 0.541 1,2 TRUE -2.267 2.972 1.948 9.326 6.267 2 0.044 1,2 TRUE 1.468 6.706 5.683 13.061 2.532 2 0.282 1,2 TRUE 0.590 5.829 4.805 12.183 3.41 2 0.182 1,2 TRUE -1.618 3.620 2.596 9.975 5.618 2 0.06 1,2 TRUE 3.055 8.293 7.270 14.648 0.945 2 0.623 結果の見方は先程と同じです。情報量規準の値がマイナスになっている・検定が有意になって いる(𝑝 < .05)項目は等値制約を外したほうが適合度が良くなっている= DIF の疑いあり,と 判断することができます。 続いては Lord の Wald 検定の結果を出してみます。ただ,Wald 検定を行うためには「推定値 の標準誤差行列」が必要になるため,まずは標準誤差を推定するように引数 SE = TRUE を付け
8.11 (おまけ)多母集団モデルと特異項目機能 加えてやり直す必要があります。 標準誤差を推定してほしい 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 を与えることで実行可能です。 DIF の検出(3):Wald 検定 1 DIF(result_group, which.par = c("a1", "d"), Wald = TRUE) 1 2 3 4 5 6 7 8 9 10 11 groups W df p Q1_1 1,2 28.652 2 0 Q1_2 1,2 36.341 2 0 Q1_3 1,2 22.179 2 0 Q1_4 1,2 16.473 2 0 Q1_5 1,2 17.345 2 0 Q1_6 1,2 1.786 2 0.409 Q1_7 1,2 8.373 2 0.015 Q1_8 1,2 5.129 2 0.077 Q1_9 1,2 19.703 2 0 Q1_10 1,2 13.783 2 0.001 結果の見方はこれまでと基本的に同じです。つまり検定結果が有意になった(𝑝 < .05)項目は グループ間でパラメータの推定値が有意に異なる= DIF の疑いあり,ということになります。 最後は ICC のズレに基づく評価方法です。こちらは Chalmers(2018)による方法が DRF() という関数に実装されています。DRF() 関数は,デフォルトでは項目特性曲線(ICC)の代わり にテスト特性曲線(TCC)の差について評価を行います*28 。項目ごとの DIF を見たい場合に は,引数 DIF = TRUE を与える必要があるのでご注意ください。また引数 plot = TRUE を与え ると, 「ICC の差のプロット」を出してくれるようになります。 DIF の検出(4):ICC のズレをチェック 1 DRF(result_group, DIF = TRUE, plot = TRUE) *28 この「テスト全体でのグループ間の挙動の違い」は DIF と同じように特異テスト機能 (Differentian Test Functioning [DTF]) と呼ばれたりします。 80
8.11 (おまけ)多母集団モデルと特異項目機能 81 Signed DIF -6 -4 -2 0 2 4 6 Q1_9 Q1_10 Q1_5 Q1_6 Q1_7 Q1_8 Q1_1 Q1_2 Q1_3 Q1_4 sDIFq 0.3 0.2 0.1 0.0 0.3 0.2 0.1 0.0 -6 -4 -2 0 2 4 6 0.3 0.2 0.1 0.0 -6 -4 -2 0 2 4 6 q 図 8.43: ICC のズレのプロット 結果を見ると,3 種類の DIF が表示されています。 sDIF 符合つき(signed)のズレの面積(単純に差を積分したもの)を表しています。図 8.41 の ようにズレが上側と下側にある場合,これらが打ち消し合って sDIF の値は小さくなりま す。つまり sDIF は「総合的に見てどちらのグループのほうが平均点が高いか」を表して いると言えるでしょう。 uDIF 符号なし(unsigned)のズレの面積(差の絶対値を積分したもの: 8.65 式)を表していま す。したがって,sDIF よりも大きな値が表示されているはずです。 dDIF (8.65) 式とは別のやり方で符号なしのズレを求めるために,絶対値の代わりに差を二乗す ることを考えた DIF です。二乗した分を戻すため,最後に全体にルートをかけます。「二 乗平均のルートをとる」という考え方は,いわゆる逸脱度(deviance)的な考え方なので dDIF と名付けられているのだと思います。多分。uDIF と比べると,ズレの面積が同じ だとしても,大きくズレている箇所があるほど,二乗されている分だけ dDIF のほうが大 きな値になると思われます。 いずれにしても,この (s|u|d)DIF の値が大きい項目は DIF の疑いあり,と言えるわけです。 ICC のズレに基づく検定を行うためには「ICC のズレの標準誤差」のようなものが必要になりま すが,これはパラメータの標準誤差だけでは対応できないようです。そのため,DRF() 関数では 代わりにブートストラップ標準誤差を用いて検定を行います。これを行うためには,引数 draws に適当な数を指定してください。
8.11 (おまけ)多母集団モデルと特異項目機能 82 ブートストラップを使って検定 1 # drawsの値は大きい方が正確になるが,その分時間がかかる 2 DRF(result_group, DIF = TRUE, draws = 500) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 $sDIF groups item sDIF CI_2.5 CI_97.5 X2 df p 1 1,2 Q1_1 0.096 0.060 0.138 25.445 1 0 2 1,2 Q1_2 0.086 0.057 0.118 34.72 1 0 3 1,2 Q1_3 0.077 0.042 0.108 21.433 1 0 4 1,2 Q1_4 0.066 0.032 0.100 14.415 1 0 5 1,2 Q1_5 0.069 0.029 0.105 13.596 1 0 6 1,2 Q1_6 0.022 -0.006 0.058 1.734 1 0.188 7 1,2 Q1_7 0.050 0.012 0.083 7.216 1 0.007 8 1,2 Q1_8 0.041 -0.001 0.078 4.433 1 0.035 9 1,2 Q1_9 0.070 0.029 0.108 13.188 1 0 10 1,2 Q1_10 0.078 0.034 0.123 12.889 1 0 $uDIF groups item uDIF CI_2.5 CI_97.5 X2 df p 1 1,2 Q1_1 0.096 0.064 0.138 26.895 2 0 2 1,2 Q1_2 0.086 0.058 0.118 28.254 2 0 3 1,2 Q1_3 0.077 0.043 0.108 17.54 1 0 4 1,2 Q1_4 0.067 0.038 0.101 17.209 2 0 5 1,2 Q1_5 0.069 0.035 0.105 12.659 2 0.002 6 1,2 Q1_6 0.022 0.006 0.059 2.151 1 0.142 7 1,2 Q1_7 0.050 0.020 0.087 7.638 2 0.022 8 1,2 Q1_8 0.041 0.010 0.079 5.147 2 0.076 9 1,2 Q1_9 0.079 0.053 0.121 21.053 2 0 10 1,2 Q1_10 0.078 0.040 0.125 13.303 2 0.001 $dDIF groups item dDIF CI_2.5 CI_97.5 1 1,2 Q1_1 0.105 0.071 0.150 2 1,2 Q1_2 0.112 0.074 0.155 3 1,2 Q1_3 0.087 0.048 0.132 4 1,2 Q1_4 0.094 0.048 0.137 5 1,2 Q1_5 0.089 0.042 0.139 6 1,2 Q1_6 0.027 0.008 0.076 7 1,2 Q1_7 0.054 0.022 0.097 8 1,2 Q1_8 0.048 0.013 0.095 9 1,2 Q1_9 0.084 0.058 0.129 10 1,2 Q1_10 0.084 0.045 0.136
8.11 (おまけ)多母集団モデルと特異項目機能 83 DIF の検定のプロット 1 2 # plot = TRUEをつけると検定はせずに図だけ表示する 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 図 8.44: DIF の検定のプロット 結果の見方はやはり同じです。つまり検定結果が有意になった(𝑝 < .05)項目は IRF がグ ループ間でズレている= DIF の疑いあり,ということになります。 DIF 検出の実際 ここまで DIF を検出するいくつかの方法を紹介しました。色々な観点からの評価がある,とい う意味では,やはり複数の指標を組み合わせて総合的に判断するのが良いのだろうと思います。 それはそれとして,DIF の有無を評価しようとする場合,そもそも関心のあるグループの間で 特性値 𝜃 の分布に差が無いということはあまり考えられません(無かったら非常にラクなのです が…) 。したがって,そもそも「𝜃 の分布が同じだとして DIF を検出する」というスキームは現実 的に結構厳しいのです。この問題を克服するためには,「グループ間で挙動が変わらない= DIF が無いことが事前にわかっている項目」であるアンカー項目 (anchor items) を加えておくこと が考えられます。アンカー項目の項目パラメータを基準とすることで,𝜃 の分布がグループで異 なる場合であっても特定の項目に DIF が生じているかを評価することが出来るわけです。具体 的にアンカー項目をどのように決定したら良いかはまた難しい話(レビューとして Kopf et al., 2015)なので,ここではとりあえずアンカー項目がある場合のやり方だけ紹介しておきます。 multipleGroup() の引数 invariance には,実は項目の名前(データの列名)を与えること
8.11 (おまけ)多母集団モデルと特異項目機能 84 も出来ます。 アンカー項目を用いた多母集団モデル 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 2 1 # Q1_1からQ1_3のパラメータだけグループ間で同じ coef(result_anchor, simplify = TRUE, IRTpars = TRUE) 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 $`1` $items a b g u Q1_1 0.444 -2.360 0 1 Q1_2 1.013 -1.831 0 1 Q1_3 0.962 -1.456 0 1 Q1_4 1.165 -1.274 0 1 Q1_5 1.061 -1.396 0 1 Q1_6 1.255 -1.450 0 1 Q1_7 1.291 -1.082 0 1 Q1_8 1.279 -1.108 0 1 Q1_9 1.056 -0.924 0 1 Q1_10 1.167 0.202 0 1 $means F1 0 $cov F1 F1 1 $`2` $items Q1_1 Q1_2 Q1_3 Q1_4 a b g u 0.444 -2.360 0 1 1.013 -1.831 0 1 0.962 -1.456 0 1 0.949 -1.177 0 1
8.12 いろいろな IRT モデル 31 32 33 34 35 36 37 38 39 40 41 42 43 44 85 Q1_5 0.957 -1.247 0 1 Q1_6 1.344 -0.734 0 1 Q1_7 1.578 -0.418 0 1 Q1_8 1.353 -0.502 0 1 Q1_9 1.735 -0.204 0 1 Q1_10 1.441 0.625 0 1 $means F1 0.741 $cov F1 F1 0.828 DIF の検出 1 2 # もともと等値制約が置かれているQ1_1からQ1_3については結果が出ない DIF(result_anchor, which.par = c("a1", "d")) 1 2 3 4 5 6 7 8 9 10 11 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 groups converged AIC SABIC HQ BIC X2 df p 1,2 TRUE 0.000 0.000 0.000 0.000 0 0 NaN 1,2 TRUE 0.000 0.000 0.000 0.000 0 0 NaN 1,2 TRUE 0.000 0.000 0.000 0.000 0 0 NaN 1,2 TRUE -1.404 3.835 2.811 10.189 5.404 2 0.067 1,2 TRUE 0.570 5.809 4.785 12.163 3.43 2 0.18 1,2 TRUE -18.972 -13.733 -14.757 -7.379 22.972 2 0 1,2 TRUE -15.059 -9.821 -10.845 -3.467 19.059 2 0 1,2 TRUE -14.729 -9.491 -10.515 -3.136 18.729 2 0 1,2 TRUE -14.361 -9.122 -10.146 -2.768 18.361 2 0 1,2 TRUE -10.458 -5.220 -6.243 1.135 14.458 2 0.001 結果を見て分かるように,アンカー項目を適切に設定しないと,その他の項目の DIF の評価 も適切に行うことができません。そういった意味では,DIF の検出というのは結構職人技だった りするのです。 8.12 いろいろな IRT モデル IRT の本質は, (8.1) 式に集約されていると個人的には思っています。それは項目への反応 (データ)は回答者と項目の交互作用によって確率的に決定するという点です。そのような挙動 をとるモデルの形の一つとして,因子分析と同じようなモデル(正規累積・ロジスティック)や 多値反応に対するモデル(GRM・GPCM など)を紹介してきました。しかし世の中にはもっと 様々な IRT モデルが存在しています。全てを紹介するのは不可能ですが,ここではいくつかの
8.12 いろいろな IRT モデル 86 モデルのさわりだけを(個人的な関心に基づいて)紹介していきたいと思います*29 。 8.12.1 多値反応に対するその他の IRT モデル 多値反応に対する IRT モデルとして,Section 8.6 では GRM や GPCM といったモデルを紹 介してきました。多値反応に対する IRT モデルには他にも様々なものが提案されており,近年注 目を集めている(気がする)のが IRTree(Böckenholt, 2017)と呼ばれるモデル群です。これ らのモデルでは,回答決定のプロセスとして名前通り「ツリー構造」を仮定します。例えば「と ても嫌い (1)」 「嫌い (2)」 「どちらでもない (3)」 「好き (4)」 「とても好き (5)」という 5 件法の項 目があったとしましょう。この項目への回答決定プロセスの考え方の一例としては,図 8.45 の ように 1. まず「どちらでもない (3)」を選択(態度保留)するかを二択で考える 2. 態度を表明する場合には「好き (4-5)」か「嫌い (1-2)」かを二択で考える 3.「とても (1,5)」の選択肢を選ぶほど好き or 嫌いかを二択で考える という感じで3つの二択を順番に行う,というものがあるでしょう。 図 8.45: IRTree の一例 するとこの多値項目の背後には,実は3つの二値項目(query)が存在していると考えること ができ,それぞれが 1.「どちらでもない」を選ぶ傾向(中心傾向):𝜃𝑝 (𝑀) が強い人かどうか (𝐴) 2. 実際にその対象物が好きかどうか:𝜃𝑝 *29 ここで紹介するモデルを始め,多くの(特殊な)IRT モデルは mirt パッケージに限らず R のパッケージとし ては実装されていなかったりします。そういったモデルについては,自分で尤度関数をプログラムして optim() 関数による数値最適化を利用する,あるいは stan などを利用してマルコフ連鎖モンテカルロ法(MCMC)に よって推定する,といった方法を取る必要があります。stan を使うためにはまずベイズ統計学を理解する必要 があるのがちょっと大変ですが,慣れるとほぼなんでもできるようになります。
8.12 いろいろな IRT モデル 87 3. 極端な選択肢を選ぶ傾向(極端傾向):𝜃𝑝 (𝐸) が強い人かどうか という異なる個人特性によって動いていると仮定することができそうです(もちろんプロセス 単位で完全に独立しているとは思いませんが,ある程度は)*30 。各 query に対する反応確率は, 普通の正規累積・ロジスティックモデルなどで表現できます。例えば中心傾向についての query だけで表現できる「どちらでもない」を選ぶ確率は (𝑀) (𝑀) 𝑃 (𝑦𝑝𝑖 = 3|𝜃𝑝 𝛼𝑖 (𝑀) (𝜃𝑝 (𝑀) −𝛽𝑖 ) )=∫ −∞ 1 𝑧2 √ exp (− ) 𝑑𝑧 2 2𝜋 (8.66) (𝑀) (𝑀) (𝑀) = 𝚽 [𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] などとなるわけです。IRTree モデルでは,項目パラメータにも上付き添字がついています。した がって,項目ごとに「どちらでもない」および極端な選択肢の選ばれやすさが異なる,という設 定がされています。 同様に, 「好き (4)」を選ぶ確率は, 「 『どちらでもない』を選ばない」かつ「好き」かつ「 『とて も』ではない」の積なので (𝑀) 𝑃 (𝑦𝑝𝑖 = 4|𝜃𝑝 (𝑀) (1 − 𝚽 [𝛼𝑖 (𝐴) (𝐸) , 𝜃𝑝 𝜃𝑝 ) = (𝑀) (𝜃𝑝 (𝑀) − 𝛽𝑖 (𝐴) )]) 𝚽 [𝛼𝑖 (𝐴) (𝜃𝑝 (𝐴) (𝐸) − 𝛽𝑖 )] (1 − 𝚽 [𝛼𝑖 (𝐸) (𝜃𝑝 (𝐸) − 𝛽𝑖 )]) (8.67) という形で表現できるわけです。 IRTree の目的あるいは利点は,「測定したい回答者の特性値」と「回答傾向」を分離して考え られる点です。一般的に,リッカート尺度への回答には個人のクセのようなものが混ざっている ことが知られています。そのせいで, 「どちらでもない (3)」を選んだ人が「好きと嫌いのちょう ど中間」なのか「本当は『どちらかといえば好き』くらいなんだけど,その程度では『好き』は 選びづらい」のかは区別がつきにくいものです。IRTree モデルならば,上の例で言えば中心傾向 (𝐴) と極端傾向の強さを,個人特性 𝜃𝑝 の一種として考えることができます。そのため,残った 𝜃𝑝 はより純粋な「測定したい回答者の特性値」として解釈できるかもしれないというわけです。当 然ながら,IRTree モデルでは 図 8.45 に示したような回答決定プロセスの仮定がかなり重要に なってきます。逆に言えば,複数のツリー構造間でモデル比較をすることで,どのプロセスが最 も当てはまりそうか,ということを評価したりも出来るかもしれません…。 8.12.2 強制選択形式に対する IRT モデル IRTree モデルはいわば,得られた回答データから回答傾向を統計的に分離しようと頑張って いるモデルの一種です。そもそもリッカート尺度では回答傾向の影響を逃れることができないわ けですが,その対抗策としては他にも「そもそも中間選択肢をなくそう」という方向を考えるこ とができます。また,入社試験などの重要な局面では「ウソをついて本当の自分より良い人に 見せよう」という思惑が働いて回答が歪められる可能性があったりします(フェイキング: e.g., Jackson et al., 2000)。これを防ぐためには,例えば「友だちが多い」「真面目である」といった *30 𝜃 の上付き添字の M, A, E はそれぞれ Middle/Mid-point, Attitude/Agree, Extreme などの頭文字です。 𝑝
8.12 いろいろな IRT モデル 88 複数の短文を同時に提示し,その中から一つだけを選択させたり,ランキングを付けさせるとい う回答方法が考えられます。これならば,前述の「中間選択肢を選びがち」といった問題も防げ ますし,フェイキングの問題もかなり抑制できるでしょう。 複数の選択肢から一つの選択肢を選ばせる回答方法は,強制選択 (Forced Choice) 型項目など と呼ばれます。また,特に選択肢数が2個の場合は一対比較 (Paired Comparison) と呼ばれる こともあります。そしてこれらの出題形式に対しては,やはり普通の IRT とは少し異なるモデ ルが考えられています。最も代表的なモデルが,Thurstonian IRT (TIRT) モデル(Brown & Maydeu-Olivares, 2011)でしょう。TIRT モデルでは,各選択肢がそれぞれ独立に「効用」を 持っていると考えます。この「効用」とは,経済学などで出てくるものと同じようなものだと考 えて OK です。あるいは単純に「各選択肢が自分に当てはまっていると感じる程度」というイ メージでも構いません。とにかく強制選択型項目への回答では,回答者は効用を最大化する選択 肢を選ぶのだ,と考えているわけです。ここでは,回答者 𝑝 が項目 𝑖 の2つの選択肢 𝐴, 𝐵 に対 (𝐴) (𝐵) してそれぞれ効用 𝑢𝑝𝑖 , 𝑢𝑝𝑖 を持っているとしましょう。すると,回答者 𝑝 は効用を最大化する 選択肢を選ぶため,観測されるデータとしては 𝑦𝑝𝑖 = { if 𝑢𝑝𝑖 − 𝑢𝑝𝑖 ≥ 0 (𝐴) (𝐵) if 𝑢𝑝𝑖 − 𝑢𝑝𝑖 < 0 (𝐴) 𝐴 𝐵 (𝐵) (8.68) と表わせます。 TIRT モデルでは,選択肢の効用 𝑢𝑝𝑖 を,よくある因子分析モデルの形 (𝐴) (𝐴) (𝐴) 𝑢𝑝𝑖 = 𝜏𝑖 (𝐴) (𝐴) + 𝑓 𝑝 𝑏𝑖 (𝐴) と定義します。つまり,平均的な効用 𝜏𝑖 (𝐴) (8.69) (𝐴) + 𝜀𝑝𝑖 , 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀(𝐴) ) 𝑝𝑖 (𝐴) と,その選択肢が反映している因子の得点 𝑓𝑖 と因 (𝐴) 子負荷 𝑏𝑝 の和によって表される効用の期待値に,そこに正規分布に従う誤差(独自因子得点) (𝐴) 𝜀𝑝𝑖 が加わって効用の実現値が確率的に決定すると考えているわけです。すると,2つの選択肢 𝐴, 𝐵 の中から 𝐴 を選ぶ確率は最終的に (𝐴) (𝐵) (𝐴𝐵) + (𝑓𝑖 𝑃 (𝑦𝑝𝑖 = 𝐴) = 𝑃 (𝑢𝑝𝑖 − 𝑢𝑝𝑖 ≥ 0) = 𝚽 [𝛼𝑖 (𝐴) (𝐴) (𝐵) (𝐵) 𝛽𝑝 − 𝑓𝑖 𝛽𝑝 )] (8.70)
8.12 いろいろな IRT モデル 89 と表されます*31 。ただし表記を簡略化するため (𝐴𝐵) 𝛼𝑖 (𝐴) = (𝐴) 𝜏𝑖 − 𝜏 𝑖 , √𝜎𝜀(𝐴) + 𝜎𝜀(𝐵) 𝑝𝑖 (𝐴) 𝛽𝑖 (𝐵) 𝛽𝑖 𝑝𝑖 (𝐴) = 𝑏𝑖 + 𝜎𝜀(𝐵) √𝜎𝜀(𝐴) 𝑝𝑖 𝑝𝑖 , (8.71) (𝐵) = 𝑏𝑖 + 𝜎𝜀(𝐵) √𝜎𝜀(𝐵) 𝑝𝑖 𝑝𝑖 という形で別の記号を用いて表現しています。つまり TIRT モデルは,一対比較の場合には,各 項目が2つの因子のみに負荷量を持っていると設定した CFA と実質的には同じことをしていま す。そして一つの項目に3つ以上の選択肢がある場合は,対ごとに異なる項目を構成していると みなして分析を行います。例えば項目 𝑖 が 𝐴, 𝐵, 𝐶, 𝐷 の4択だとすると,この項目は「A と B」 「A と C」 「A と D」「B と C」「B と D」「C と D」という6つの二値項目があると考えます。そ して回答者 𝑝 が選択肢 𝐴 を選んだとすると, 「A と B」 「A と C」 「A と D」については(いずれ も A を選んだという)結果が観測され,残りの「B と C」 「B と D」 「C と D」については欠測扱 いにします。そして観測された3項目についてはそれぞれ (𝐴𝐵) + (𝑓𝑖 (𝐴𝐶) + (𝑓𝑖 (𝐴𝐷) + (𝑓𝑖 𝑃 (𝑦𝑝𝑖 = 𝐴|𝐴, 𝐵) = 𝚽 [𝛼𝑖 𝑃 (𝑦𝑝𝑖 = 𝐴|𝐴, 𝐶) = 𝚽 [𝛼𝑖 𝑃 (𝑦𝑝𝑖 = 𝐴|𝐴, 𝐷) = 𝚽 [𝛼𝑖 (𝐴) (𝐴) (𝐵) (𝐵) 𝛽𝑝 − 𝑓𝑖 𝛽𝑝 )] (𝐴) (𝐴) (𝐶) (𝐶) 𝛽𝑝 − 𝑓𝑖 𝛽𝑝 )] (8.72) (𝐴) (𝐴) (𝐷) (𝐷) 𝛽𝑝 − 𝑓𝑖 𝛽𝑝 )] という形でそれぞれ選択肢 A に関する部分が共通しているという設定のもとで定式化されるわ けです。 8.12.3 回答時間を利用する IRT モデル コンピュータによる測定 (Computer Based Testing [CBT]) では,項目反応以外の情報を得る ことが容易になりました。もっとも代表的な情報として,回答時間 (Response Time [RT]) が挙 げられるでしょう。Chapter 3 では,回答時間を「適当な回答をしている人を検出する」といっ た目的で使用しましたが,その一方で回答者の特性値推定に利用するという方向性も相当に研究 されています。そういうわけで,回答時間を利用する IRT モデルを挙げていくとキリがないの ですが,ここではシンプルなモデルとして,「当てはまる」「当てはまらない」の二択の項目に対 するモデル(Ferrando & Lorenzo-Seva, 2007)を紹介します。 回答時間を利用する IRT モデルでは,そもそも回答者の特性値と回答時間にはどのような関 係があるかを考えるところから始まります。性格特性を尋ねる項目の場合には,「特性の強さが 中途半端な人ほど回答時間が長い」という逆 U 字の関係になることが経験的に知られています。 *31 2つの選択肢の比較のときに,その背後にある効用を比較していて,しかも効用が確率的に変動する,という考 え方をサーストンの比較判断の原則(law of comparative judgment: Thurstone, 1927)と呼ぶため,このモ デルは Thurstonian IRT と名付けられています。
8.12 いろいろな IRT モデル 90 例えば「友だちが多い」という項目に対して,友だちが多い/少ない自信がある人,つまり「外 向性」といった特性が極端な人は,きっとすぐに選択肢を選ぶことでしょう。これに対して, 『友 だちはそこそこ居るけど』という人は,きっと「当てはまる」 「当てはまらない」のどちらを選ぶ かに迷いが生じるはずです。したがって,同じ「当てはまる」を選んだ人同士であっても,素早 く「当てはまる」を選んだ人のほうがより強い特性を持っていると考えるのが自然です。 ということで,この関係を尤度関数に組み込んで IRT モデルを構築していきます。と言って も考え方は簡単で,回答時間 𝑡𝑝𝑖 についても項目反応 𝑦𝑝𝑖 の (8.1) 式と同じように 𝑡𝑝𝑖 = 𝑔(𝜔𝑝 , 𝛾𝑖 ) (8.73) という形で,回答者の要因 𝜔𝑝 と項目の要因 𝛾𝑖 の交互作用として考えます。Ferrando & Lorenzo- Seva(2007)では,これを具体的に ln(𝑡𝑝𝑖 ) = 𝜇 + 𝜔𝑝 + 𝛾𝑖 + 𝑏𝛿𝑝𝑖 + 𝜀𝑝𝑖 , 𝜀𝑝𝑖 ∼ 𝑁 (0, 𝜎𝜀2 ) 𝛿𝑝𝑖 = √𝛼2𝑖 (𝜃𝑝 − 𝛽𝑖 )2 (8.74) と置きました。ここで 𝜇 は全体的な切片項,𝜔𝑝 は回答者 𝑝 の回答スピード,𝛾𝑖 は項目 𝑖 の所要 時間を表現したパラメータです。また,𝛿𝑝𝑖 の中身には (𝜃𝑝 − 𝛽𝑖 )2 というものが入っています。 この項は,𝜃𝑝 = 𝛽𝑖 のとき,つまり反応確率が五分五分のときに最小値を取ります。係数 𝑏 が負 の値であれば,回答者の特性値 𝜃𝑝 が項目困難度 𝛽𝑖 に近いほど回答時間の期待値が長くなること で,逆 U 字の関係を表現しているわけです。 回答時間などの追加情報を用いて特性値を推定する IRT モデルでは,うまく行けばより精度 の高い,あるいは妥当性の高い推定を行うことが出来ると期待されています。ちなみに,IRT に 限らず社会調査などの分野では,このように回答内容以外の情報(位置情報や回答デバイスの情 報などを含む)をパラデータ(Kreuter, 2013)と呼んでいます。パラデータは,やはり適当な回 答の検出に使われるだけでなく,よりノイズの少ないデータを収集するための調査改善など様々 な用途で利用されていたりします。 8.12.4 展開型 IRT モデル これまでに紹介した IRT モデルはすべて, 𝜃 の値が高ければ高いほど高い点数をとる(i.e., 「当てはまる」寄りの回答をする,正解する)という仮定がされていました。項目特性曲線にして も, (𝛼𝑖 > 0 ならば)右上がり(単調増加)でした。ですが質問のタイプによっては「過ぎたるは なお及ばざるが如し」的なものが存在しています。よくあるのは「態度」に関する項目です。例 えば「消費税を 5% にするべきだ」という項目は,「消費税は完全に撤廃してほしい」人も「消 費税はこのままで良い」人も「そう思わない」を選択するでしょう。同様に,心理特性に関して も「カフェで友人と静かに会話をするのは楽しい(因子:外向性)」 (Stark et al., 2006)という 項目に対しては「そもそも他人と話したくない人」も「もっとワイワイ騒ぎたい人」も「当ては まらない」を選ぶと思われます。つまりこうした項目における項目特性曲線は,図 8.46 のよう に,あるところを境に単調減少に転じると考えられるわけです。 ロジスティック・正規累積モデルでは,「項目 (𝛽𝑖 )」と「回答者 (𝜃𝑝 )」は共通の軸の上でマッ
8.12 いろいろな IRT モデル 91 図 8.46: 中間にピークが来る ICC ピングされていると考えることが出来ます。すごく簡単にいえば「困難度が 𝛽𝑖 の項目は 𝜃𝑝 = 𝛽𝑖 の回答者にとって(易しすぎず難しすぎず)ちょうどよい難易度(50%)の項目である」という ことです。これを先程の項目に置き換えて考えると, 「𝛽𝑖 と 𝜃𝑝 の距離が近いほど,その人によっ てちょうどよい程度の(理想の)意見・態度を表している」ということができそうです。このよ うな「理想点 (ideal point)」の考え方に基づく IRT モデルを展開型 (unfolding) IRT モデル と呼びます*32 。 展開型 IRT モデルでは,反応確率 𝑃 (𝑦𝑝𝑖 = 1) が「𝛽𝑖 と 𝜃𝑝 の距離」によって規定されます。し たがって,最も基本的な一次元・二値モデルの場合は 1 𝑃 (𝑦𝑝𝑖 = 1) = exp [− (𝛽𝑖 − 𝛼𝑖 𝜃𝑝 )2 ] 2 (8.75) と表されます(Maydeu-olivares et al., 2006) 。ここでも 2 つの項目パラメータの意味は概ね通 常の場合と同じです。𝛽𝑖 はその項目の「位置」を表しており,𝜃𝑝 がこの値に近い人ほど「当ては まる」と回答する確率が高くなることがわかります。𝛼𝑖 は項目識別力を表していますが,展開型 モデルの場合は,正規分布の分散パラメータのように山の狭さを規定します(図 8.47 ) 。ですが 結局「𝜃𝑝 が理想点 𝛽𝑖 付近か否かを区別する強さ」という意味では,これまでに出てきた項目識 別力と同じような役割を果たしているわけです。 mirt では,itemtype="ideal"と指定することで展開型モデルを実行可能です。 mirt で展開型モデル 1 result_ideal <- mirt(dat_bin, model = 1, itemtype = "ideal") 通常の IRT モデルと同じように coef() でパラメータの推定値を見ることもできるのですが, どうやら IRTpars = TRUE には対応していないようなのでご注意ください。 *32 これに対して通常の IRT モデルは dominance model と呼びます(日本語訳がわからない…) 。
8.12 いろいろな IRT モデル 92 P(ypi = 1) ai = 0.5 ai = 1 ai = 2 qp q p = bi 図 8.47: 識別力の役割(展開型) 展開型モデルの推定値のチェック 1 2 1 2 3 4 5 6 7 8 9 10 11 12 13 # IRTpars = TRUEとするとエラーが出る coef(result_ideal, simplify = TRUE) 14 15 16 17 18 19 20 $items a1 d Q1_1 0.118 -0.716 Q1_2 0.214 -0.468 Q1_3 0.250 -0.571 Q1_4 0.287 -0.609 Q1_5 0.255 -0.600 Q1_6 0.618 -0.312 Q1_7 0.651 -0.455 Q1_8 0.621 -0.489 Q1_9 0.571 -0.653 Q1_10 0.667 -1.187 $means F1 0 $cov F1 F1 1 また,ICC のプロットも同じように実行可能です。 項目特性曲線(展開型) 1 plot(result_ideal, type = "trace", facet_items = FALSE)
8.12 いろいろな IRT モデル 93 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.48: 項目特性曲線(展開型) 例えば 1 番目の(一番山が緩い)項目は, (𝛼𝑖 , 𝛽𝑖 ) = (0.118, −0.716) ということで, 𝜃 = 0.716/0.118 ≃ 6.068 のあたりで最大値 𝑃 (𝑦𝑝𝑖 = 1) = 1 となります。そのため,𝜃 が 6 までしか 表示されていない図 8.48 では右上がりのように見えています。 INFO 多値の展開型 IRT モデル 多値の展開型 IRT モデルの代表的なものは,一般化段階展開型モデル(Generalized graded unfolding model[GGUM]: Roberts, 2019)です。このモデルは,GPCM をベースに展 開型 IRT モデルを構成しています。 例えば 4 カテゴリの項目がある場合には,「全く当てはまらない(下すぎて)」「あまり当 てはまらない(下すぎて) 」~「あまり当てはまらない(上すぎて) 」 「全く当てはまらない (上すぎて)」という要領で 8 個のカテゴリがあると考えます。こうすることで,隣同士の カテゴリの比較は今まで通りできるので GPCM 的な方法が適用できる,というわけです。 そして「全く当てはまらない(下すぎて)」と「全く当てはまらない(上すぎて)」の確率 を合わせたものを「全く当てはまらない」の選択確率として表します。このままだとパラ メータ数が多すぎる(4 カテゴリに対してカテゴリ困難度が 7 個必要になってしまう)の で,理想点を軸にして各カテゴリ閾値間の距離は左右で等しいという仮定をおいた上で定 式化しています。 なお mirt では,itemtype = "ggum"と指定することで実行可能です。
8.13 IRT の使い道 94 8.13 IRT の使い道 最後に,IRT の考え方によって出来るようになることをいくつか紹介します*33 。 8.13.1 等化 Section 8.1 では,IRT を項目依存性と集団依存性をクリアするための理論だと紹介しました。 ですが実際には,IRT モデル(2 パラメータの GRM)はカテゴリカル因子分析と同じなので 𝜃𝑝 はデータをとった集団ごとに標準化されています。つまり最初に例示したような,異なる集団に 対して実施した異なる尺度の比較はこのままではできません。もちろん反対に,異なる尺度を用 いて実施した異なる集団同士の比較もこのままではできません。こうした比較を可能にするため に,IRT では等化 (equating) またはリンキング (linking) という手続き(Kolen & Brennan, 2014)を行います*34 。等化の技術は,TOEIC のような大規模テストでも実際に利用されていま す。そのおかげで,毎回異なるはずの難易度によらず「TOEIC ***点」というように点数を履歴 書に書けるわけです。 (正直言うと,等化=異なる集団や尺度の比較を考えない限りは,特に心理尺度においてわざ わざ因子分析や SEM ではなく IRT を選ぶ理由はほぼ無い気すらします) 等化の基本的な考え方を説明しましょう。ここからは IRT で推定されるパラメータが集団に よらない項目パラメータと項目によらない個人特性だという仮定のもとで話を進めます。 共通回答者による等化 まず,個人特性値が「項目によらない」のであれば,どんな尺度で推定した場合でも 𝜃𝑝 の値 は同じになるはずです。いま,𝑝 さんが同じ構成概念を測定する尺度 1 と尺度 2 による 2 回の調 査の両方に回答したとします。この 2 回のデータをもとにパラメータを推定する場合には,いま までの手順にのっとってそれぞれ 𝜃𝑝 ∼ 𝑁 (0, 1) と仮定して推定を行います。その結果 • 尺度 1 での推定値 𝜃𝑝̂ = 0.5 (2) • 尺度 2 での推定値 𝜃𝑝̂ = −0.5 (1) だったとしましょう。仮に尺度 1 による推定が「項目によらない個人特性」の値とすると,尺 度 2 での推定値は本来の値より 1 小さい値になっているようです。𝑝 さんの特性値が尺度 1 回答 時と尺度 2 回答時で変わっていないとしたら,この 1 のズレは「尺度 2 は尺度 1 よりも特性値を (2) 小さめに推定する尺度である」ことが原因と考えられます。そこで,𝜃𝑝̂ = −0.5 がどうにかし (1) て 𝜃𝑝̂ と同じ 0.5 になるように調整する(ゲタを履かせる)必要がありそうです。 もう少し規模を広げて,同じようにこの尺度 1 と尺度 2 による 2 回の調査の両方に回答した人 が 100 人いたとします。その結果,2 回の推定値の散布図が 図 8.49 のようになったとします。 *33 R での実践例までやると膨大になってしまうので,概念の紹介だけにとどめます。気になる人は自分で調べたり 直接質問してください。 *34 R では equateIRT や plink といったパッケージによって等化のメソッドが提供されています。
8.13 IRT の使い道 95 そしてこの 2 回の推定値の回帰直線を求めたら (2) (8.76) (1) 𝜃𝑝̂ = 𝐴𝜃𝑝̂ + 𝐵 (2) (1) (2−1) であったとします。この場合,全員の 𝜃𝑝̂ をなるべく 𝜃𝑝̂ に近づくように調整した値 𝜃𝑝̂ (2−1) 𝜃𝑝̂ (2) = 𝜃𝑝̂ − 𝐵 𝐴 は, (8.77) という形で求められそうだとわかります。 3 qp (2 ) 2 1 0 -1 -2 -1 0 1 (1 ) 2 qp 図 8.49: 2 回の測定における推定値 2 パラメータロジスティックモデルでは,項目反応関数は 𝑃 (𝑦𝑝𝑖 = 1) = でした。ここで, 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.78) 𝜃𝑝 − 𝐵 𝐴 𝛽 −𝐵 𝑖 𝛽𝑖∗ = 𝐴 𝛼∗𝑖 = 𝐴𝛼𝑖 (8.79) 𝜃𝑝∗ = とおくと, 𝑃 (𝑦𝑝𝑖 = 1) = = = 1 1 + exp [−𝛼∗𝑖 (𝜃𝑝∗ − 𝛽𝑖∗ )] 1 𝜃 −𝐵 1 + exp [−𝐴𝛼𝑖 ([ 𝑝𝐴 ] − [ 𝛽𝑖𝐴−𝐵 ])] 1 1 + exp [−𝛼𝑖 (𝜃𝑝 − 𝛽𝑖 )] (8.80)
8.13 IRT の使い道 96 となり,もとの項目反応関数と同じになります。探索的因子分析では,潜在変数のスケールの不 定性の問題から,全ての観測変数を標準化していました。等化では,このスケールの不定性をあ る種逆手に取ることで 𝛼, 𝛽 や 𝜃 を適当なスケールに変換しているわけです。具体的には,2 つの (2) (1) 尺度に両方とも回答した全員の 𝜃𝑝̂ がなるべく 𝜃𝑝̂ に近づくようにスケールを変換し,これに応 じて項目パラメータにも変換をかけます。こうして得られた変換後の項目パラメータならば,尺 度 A と尺度 B でも問題なく比較ができるようになります。 共通項目による等化 先程説明した等化の流れは,共通の項目がある場合でも同じように考えられます。つまり,2 つの集団に対して実施する際に,それぞれの尺度に共通の項目をいくつか入れておけば, (2) 𝛽𝑖 − 𝐵 𝐴 (2−1) (2) 𝛼𝑖 = 𝐴𝛼𝑖 (2−1) 𝛽𝑖 = (8.81) となるような等化係数 𝐴, 𝐵 を求めることができ,それに応じた変換が可能となります。等化係 数 𝐴, 𝐵 を求める方法としては,推定値の要約統計量を用いる方法 (mean-mean 法, mean-sigma 法) や,「変換後の ICC がなるべく重なるようにする」方法(Stocking-Lord 法: Stocking & Lord, 1983; Haebara 法: Haebara, 1980)が一般的です。 等化は IRT の可能性を大きく広げてくれる重要な技術ですが,その適用には様々なハードル があります。何より 2 つの尺度が測定している 𝜃 が等質であるという大前提が重要です。もしも 尺度 1 と尺度 2 がそもそも異なる構成概念を測定しているならば,パラメータを調整したとこ ろで同じモノサシの上で比較することはナンセンスです。本来は全部まとめてデータを取った上 で一次元性の検証ができれば良いのですが,等化をしたい場合にはそんなことができないので, 理論的側面などから測定内容の等質性を担保してあげる必要があります。また,等化係数はパラ メータの推定値をもとに計算されます。つまり,等化の係数には少なからず誤差がある,という ことです。厄介なのは,尺度 2 での項目パラメータの推定値それ自体にも誤差があるという点 です。つまり等化後のパラメータには「等化前の誤差+等化係数の誤差」という 2 つの誤差が重 なっているので,等化をしない場合と比べるとどうしても推定値は不安定になりやすいです。理 論的には,等化は繰り返し行うことも可能です(尺度 4 →尺度 3 →尺度 2 →尺度 1)が,実際に 等化を繰り返すとなかなか厳しいことになるので注意が必要です。 8.13.2 テストの設計 テストの一番の目的は,回答者の特性値 𝜃 を精度良く推定することです。ですがテストの目的 によって「どの範囲の 𝜃 をどの程度の精度で推定できたら良いか」は変わってきます。例えばふ つうの心理尺度では,比較的広い範囲の特性値をそれなりの精度で推定したいと思うでしょう。 この場合には,困難度が高い項目から低い項目までをまんべんなく用意するのが良さそうです。 一方で,心療内科が扱うような「診断」を目的とした尺度や合格/不合格を決めるテストなどの 場合,ある特性の値が具体的にどの程度かよりも「基準値以上か以下か」に強い関心があったり
8.13 IRT の使い道 97 します。このような場合には,困難度が基準値周辺の,識別力ができるだけ高い項目ばかりを用 意するほうが良さそうです。このように最適な項目の特性はテストの目的に応じて異なります。 IRT では項目・テスト情報関数が 𝜃 の関数になっているおかげで,様々な 𝜃 の値に合わせて オーダーメイドで尺度を構成することができます。幸いなことにテスト情報関数は局所独立の仮 定が満たされていれば項目情報関数の和なので,他の項目との関係を気にすることなく,一つ一 つの項目に対して「その項目を追加したらテスト情報量がどの程度増加するか」のみを考えたら よいわけです。特に自動テスト構成 (automated test assembly) では,コンピュータの力を使っ て,例えば「𝜃 = (−4, −2, 0, 2, 4) の各点における標準誤差が 0.2 以下=テスト情報量が 25 以上 になる最小項目数のセット」などの条件にあった最適な項目セットを探し出してくれます*35 。 ただこれだけのテストの設計を行うためにはかなりの事前準備が必要となります。というのも 項目情報関数は項目パラメータが既知でないとわかりません。また,目的に応じて最適な項目 セットを作成するために,なるべく候補となる項目の数は多いほうが良いです。もちろんすべて の項目のパラメータを一気に推定するわけにはいかないので,多くの場合は複数回に分けて集め たデータをもとにパラメータを等化したりして項目プールを作成する必要があるでしょう。とい うことで心理尺度の場合にはなかなかそこまで大量の項目を用意すること自体が難しいので,こ こまで理論的な尺度構成を行っている先行事例はたぶんほとんど存在しません。一応用途として は,100-200 項目くらい用意した候補の中から「最高の 20 項目」を選んで心理尺度を作成する, といった方法はありそうな気がします*36 。 8.13.3 適応型テスト テストや尺度の中には通常様々な困難度の項目が混ざっているため,ほぼ確実に「その人には ほぼ意味のない項目」が含まれることになります。この問題は,先ほど紹介したテストの設計の 考え方を持ってしても避けられません。そこで自動テスト構成をもっと局所的に行うと「個人に 応じて最適な項目セットを選び出す」という考えになります。つまり視力検査のように 𝜃𝑝 が高 そうな人にはより困難度の高い項目を,𝜃𝑝 が低そうな人にはより困難度の低い項目を出していけ ば,効率的にテスト情報量を稼ぐことができそうです。コンピュータ上で行われるコンピュータ 適応型テスト (computerized adaptive testing [CAT]) *37 では,この考えに基づいて個人 ごとに次の項目を出し分けていくことができます。 ここでは IRT に基づく最もシンプルな項目選択法を紹介しましょう。推定の標準誤差は 𝜎(𝜃 ̂ |𝜃 ) = 𝑝 𝑝 1 √𝐼(𝜃𝑝 ) で表されていました。したがって推定精度をなるべく上げるためには,𝐼(𝜃𝑝 ) を大きくする=項 目情報量 𝐼𝑖 (𝜃𝑝 ) が大きい項目を選んでいけばよいわけです。ただし 𝐼𝑖 (𝜃𝑝 ) は「真値 𝜃𝑝 における *35 R では lpSolve や glpkAPI,Rsymphony などのパッケージを使うと自動テスト構成ができるようになります。 *36 100-200 項目すべてに真面目に回答してくれる人たちを数百人用意する必要があるのですが,それさえクリアで *37 きれば… R では catR や mirtCAT といったパッケージが役に立つかもしれません。
8.13 IRT の使い道 項目情報量」のため,真値 𝜃𝑝 がわからない状況では 𝐼𝑖 (𝜃𝑝 ) の値もわかりません。そこで,1 問 解答するたびに推定を行い,その時点での暫定的な推定値 𝜃𝑝∗̂ における項目情報量 𝐼𝑖 (𝜃𝑝∗̂ ) が最大 となる項目をとりあえず選び続けたらなんだかうまく行きそうです*38 。 CAT はうまく行けば従来式と同程度の推定精度を半分以下の項目数で達成できたりもする, かなり夢のある技術です。ただすべての回答者にとって最適な項目を用意できないと効果が薄れ てしまうため,単純に自動テスト構成するときよりももっと大量の項目を用意しておく必要があ ります。また,テストを実施しリアルタイムで次の項目を選択するためのシステムの構築など, 実装に向けたハードルはかなり高いものです。とりあえずは「そういう技術も IRT によって実 現できるんだなぁ」くらいに思っておいていただければ良いと思います*39 。 *38 もちろん他にも項目の選び方は色々あります。ざっくりいうと「推定精度をなるべく高める」方向性と「なるべ く人によって異なる項目が出るようにする(項目情報が外部に漏れるのを防ぐため) 」という方向性があります。 *39 IRT によらない CAT もあることにはあります。ただ IRT に基づく CAT のほうが圧倒的マジョリティです。 98
参考文献 参考文献 Berrío, Á. I., Gómez-Benito, J., & 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 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., & 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.p114-127 Brown, A., & 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., & 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 Cho, E. (2023). Interchangeability between factor analysis, logistic IRT, and normal ogive IRT. Frontiers in Psychology, 14, 1267219. https://doi.org/10.3389/fpsyg.2023.1267219 Choi, Y.-J., & 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., & Williams, E. A. (1984). Appropriateness measurement with polychotomous item response models and standardized indices (Issue ADA141365). Model Based Measurement Laboratory, University of Illinois. Ferrando, P. J., & 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 99
100 参考文献 Gessaroli, M. E., & De Champlain, A. F. (1996). Using an approximate chi-square statistic to test the number of dimensions underlying the responses to a set of items. Journal of Educational Measurement, 33(2), 157–179. 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/psycholres19 54.22.144 南風原 朝和(2000) .局所独立性 Retrieved July 4, 2022 from https://www.p.u- tokyo.ac.jp/~haebara/local_ind/ Jackson, D. N., Wroblewski, V. R., & 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., & 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., & Brennan, R. L. (2014). Test equating, scaling, and linking (3rd ed.). Springer. Kopf, J., Zeileis, A., & Strobl, C. (2015). Anchor selection strategies for DIF analysis. Educational and Psychological Measurement, 75(1), 22–56. https://doi.org/10.1177/0013 164414529792 Kreuter, F. (Ed.). (2013). Improving surveys with paradata: analytic use of process information. Wiley series in survey methodology. John Wiley & Sons. Lord, F. M. (1980). Applications of item response theory to practical testing problems. L. Erlbaum Associates. Lord, F. M., & Novick, M. R. (1968). Statistical theories of mental test scores. AddisonWesley series in behavioral science quantitative methods. Addison-Wesley. Mantel, N., & Haenszel, W. (1959). Statistical aspects of the analysis of data from retrospective studies of disease. Journal of the National Cancer Institute. https://doi.or g/10.1093/jnci/22.4.719
参考文献 Maydeu-olivares, A., Hernández, A., & 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 McKinley, R. L., & Mills, C. N. (1985). A comparison of several goodness-of-fit statistics. Applied Psychological Measurement, 9(1), 49–57. https://doi.org/10.1177/014662168500 900105 Meijer, R. R., Niessen, A. S. M., & 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 Meijer, R. R., & Sijtsma, K. (2001). Methodology review: evaluating person fit. Applied Psychological Measurement, 25(2), 107–135. https://doi.org/10.1177/01466210122031957 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/0146621692 01600206 Orlando, M., & 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., & 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 W. J. van der Linden (Ed.), Handbook of item response theory: Volume one, models (pp. 189–209). CRC Press. Roberts, J. S. (2019). Generalized graded unfolding model. In W. J. van der Linden (Ed.), Handbook of item response theory: Volume one, models (pp. 369–392). CRC Press. Roussos, L. A., & Ozbek, O. Y. (2006). Formulation of the DETECT Population Parameter and Evaluation of DETECT Estimator Bias. Journal of Educational Measurement, 43(3), 215–243. https://doi.org/10.1111/j.1745-3984.2006.00014.x 101
102 参考文献 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., & 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, San Francisco, CA. Stark, S., Chernyshenko, O. S., Drasgow, F., & 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., & Lord, F. M. (1983). Developing a common metric in item response theory. Applied Psychological Measurement, 7 (2), 201–210. https://doi.org/10.1177/0146 62168300700208 Stout, W. (1987). A nonparametric approach for assessing latent trait unidimensionality. Psychometrika, 52(4), 589–617. https://doi.org/10.1007/BF02294821 Stout, W., Habing, B., Douglas, J., Hae Rim Kim, Roussos, L., & Jinming Zhang. (1996). Conditional covariance-based nonparametric multidimensionality assessment. Applied Psychological Measurement, 20(4), 331–354. https://doi.org/10.1177/0146621696020004 03 Takane, Y., & de Leeuw, J. (1987). On the relationship between item response theory and factor analysis of discretized variables. Psychometrika, 52(3), 393–408. https: //doi.org/10.1007/BF02294363 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. (1981). Using simulation results to choose a latent trait model. Applied Psychological Measurement, 5(2), 245–262. https://doi.org/10.1177/014662168100500212 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. (2007). Conditional covariance theory and detect for polytomous items. Psychometrika, 72(1), 69–91. https://doi.org/10.1007/s11336-004-1257-7
参考文献 Zhang, J., & 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 103