読了: Brown & Maydeu-Olivares (2012) 強制選択データのThurstonian IRTモデルをMplusで推定する

Brown, A., Maydeu-Olivares, A. (2012) Fitting a Thurstonian IRT model to forced-choice data using Mplus. Behavior Research Methods. 42(2).

 一対比較データからの個人効用推定の問題で延々と悩んでいる関係で、毒もくらわば皿まで、という気分で読んでみた奴。

 いわく。
 強制選択設問のうち、各項目がそれぞれひとつの心理的属性を測っているのを多次元強制選択(MFC)という。
 ふつう強制選択設問の得点化の際には、ブロック内の順位の逆数を足し上げる。ひとりの持ち点は固定になる(イプサティブデータ)。イプサティブデータは、(1)相対的で、(2)構成概念妥当性が歪み(相関が負に歪む)、(3)基準関連妥当性が歪み(外的基準との相関が歪む)、(4)信頼性が歪む。これは得点化の仕方が悪いからである。
 われらがThurstonian IRTモデルはMFCをうまく扱えます。Mplusで簡単に推定できる。Mplusのコードを書くのが厄介だが、コードを生成する簡単なExcelマクロも作ったよ。

Thurstonian IRTモデル

強制選択回答のコーディング
 \(n\)項目からなるブロックについて順位付けを求めるとしよう。反応を\(\tilde{n} = n(n-1)/2\)ペアの一対比較に置き換えてコード化する。ペア\(l\)の結果を二値変数\(y_l\)とする。
 部分的な順序づけ(最良と最悪を選ぶとか)の場合は欠損を持つペアがあることになる。この欠損はMCARではないがMARではある。

選好反応と潜在特性との関係のモデリング
 項目\(i\)の効用を\(t_i\)とする。Thurstoneに従い効用は正規分布に従うとする。ペア\(l = \{i, k\}\)の効用差を\( y^*_l = t_i – t_k\)とし、0以上のときそのときに\(y_l = 1\)とする。
 心理属性が\(d\)個あって $$ t_i = \mu_i + \sum_{a=1}^{d} \lambda_{ia} \eta_a + \epsilon_i$$ 行列で書いて $$ \mathbf{t} = \mathbf{\mu} + \mathbf{\Lambda} \mathbf{\eta} + \mathbf{\epsilon} $$ とする。\(\mathbf{\eta}\)は共通属性ベクトル、\(\mathbf{\Lambda}\)は因子負荷行列, \(\mathbf{\mu}\)は項目切片ベクトル。\(\mathbf{\eta}\)は独自因子ベクトルで互いに無相関とする。
 \(\mathbf{\Phi} = var(\mathbf{\eta})\)とし、識別のため対角を1とする。\(\Psi^2 = var(\mathbf{\epsilon}) \)とする。こっちは対角行列である。[二乗をつけたりつけなかったりすんの気持ち悪いな…]

 効用差は$$ y^*_l = t_i – t_k = -\gamma_l + \sum_a^d (\lambda_{ia} – \lambda_{ka}) \eta_a + \epsilon_i – \epsilon_k$$ ただし\(\gamma_l = -(\mu_i -\mu_k)\)。これを行列で書いて $$ \mathbf{y}^* = -\mathbf{\gamma} + \mathbf{\breve{\Lambda}} \mathbf{\eta} + \mathbf{\breve{\epsilon}} $$ ブロック数を\(p\)として、\(\mathbf{\gamma}\)は長さ\( p \times \tilde{n} \)のベクトル。\( \mathbf{\breve{\Lambda}} \)は\( ( (p \times \tilde{n}), d) \)の因子負荷行列で、計画行列を\( \mathbf{A} \)として\( \mathbf{\breve{\Lambda}} = \mathbf{A \Lambda}\)である。\( \mathbf{\breve{\epsilon}} \) の共分散行列は\( \mathbf{\breve{\Psi}}^2 = \mathbf{A} \mathbf{\Psi}^2 \mathbf{A}^\top \)である。
 [あああ… 因子負荷っていう言葉が何を指しているのか気をつけないといけないな。項目の因子負荷行列\(\mathbf{\Lambda}\)と、項目ペアの因子負荷\( \mathbf{\breve{\Lambda}} \)があるのだ。後者には計画行列がはいっている]

independent-clusters Thurstonian IRTモデルのパラメータ
 たいていの場合、各項目はどれかひとつの特性を測っている[EFAでいう単純構造のことかな]。こういうモデルをindependent-cluster basisという。以下ではこのモデルを仮定する。そうでなくても推定方法は同じだが識別性条件が変わってくる。
 項目\(i, k\)がそれぞれ\(\eta_a, \eta_b\)を測っているとしよう。するとこう単純化できる: $$ y^*_l = -\gamma_l + \lambda_i \eta_a – \lambda_k \eta_b + \epsilon_i – \epsilon_k $$ いっぽう、もしどちらも\(\eta_a\)を測っていたら$$ y^*_l = -\gamma_l + (\lambda_i – \lambda_k) \eta_a + \epsilon_i – \epsilon_k $$

 Thurstonian IRTモデルを推定した結果として手に入るのは次の3つである。

  • 閾値パラメータ \(\gamma_l\)。\(p \times \tilde{n}\)個。
  • 因子負荷パラメータ。\(p \times n\)個 [ここでいっているのは\(\mathbf{\Lambda}\)に入ってる要素のことで、\(\mathbf{\Lambda}\)は単純構造だから\(p \times n\)個なのである]。
    • \(n=2\)のとき、それぞれの項目の因子負荷は、\( \mathbf{\breve{\Lambda}} \)に一度だけ登場する[あ… 同一項目の複数回提示は考えてないってことね]。
    • \(n > 2\)のとき、それぞれの項目は\(n-1\)ペアに含まれるから、なので、それぞれの項目の因子負荷は\( \mathbf{\breve{\Lambda}} \)に\(n-1\)回登場する。[たとえば\(n=3\)として、あるブロックは\( \mathbf{\breve{\Lambda}} \)上に部分行列 $$ \left( \begin{array}{ccc} \lambda_1 & -\lambda_2 & 0 \\ \lambda_1 & 0 & -\lambda_3 \\ 0 & \lambda_2 & -\lambda_3 \end{array} \right) $$ を持つ。確かに各項目の因子負荷は2回登場するが、符号は入れ替わることがあるわけだ]。
  • 独自性パラメータ \(\psi^2_i\)。\(p \times n\)個。
    • \(n=2\)のとき、\( \mathbf{\breve{\Psi}}^2 \)は、対角が\(\psi^2_1+\psi^2_2, \psi^2_3+\psi^2_4, \ldots\)となり、非対角は0となる。
    • \(n > 2\)のとき、\( \mathbf{\breve{\Psi}}^2 \)はブロック対角行列になる。ブロックのサイズは\((\tilde{n}, \tilde{n})\)である。[うわー面倒くさい。いちおうメモしておくと、たとえば\(n=3\)として、ブロックの要素(1,1)はペア(1,2)の誤差分散なので\( \psi^2_1 + \psi^2_2 \)となる。要素(1,2)はペア(1,2)とペア(1,3)の誤差分散なので\(\psi^2_1\)となる。要素(1,3)はペア(1,2)とペア(2,3)の誤差分散なので\(-\psi^2_2\)となる。頭痛くなってきた…]

モデル識別
 ある特性だけを測定するように設計されたMFC項目 (IRTではこれをmulti-unidimensional structureという) のためのThurstoninan IRTモデルを識別する場合は…
 潜在特性の分散を1とする。項目誤差分散は、ブロックサイズ\(n > 2\)ならば各ブロックあたり1項目で固定すれば良い。以下では最初の項目の誤差分散を1にする。\(n = 2\)ならば、項目の誤差分散は識別できないので、二値アウトカムの誤差分散を1に固定する。
 [えええ、まじで? \(n (> 2)\)項目の順位付け課題の場合、\(\tilde{n}\)個の二値アウトカムは遷移律を必ず満たすわけで、二値アウトカムは独自性を持たないとみるべきでは? \(n = 2\)の場合、Maydeu-Olivares & Bockenholt (2005) では誤差分散は自由推定でも等値制約でもいいっていってたじゃんー!
 …と錯乱したのだが、考えてみたらここでは項目効用のさらにその裏に共通因子を考えているので、たとえ順位付け課題でも項目誤差を持たせたいわけだ。それにしても、一対比較の場合、アウトカムの測定誤差と項目効用の独自性を分離できないというのは、ちょっと辛くないですか… 別にいいのかなあ…]

 \(n=2\)で\(d=2\)の場合は本質的にEFAとなるので、さらなる識別制約が必要になる。後述する。また、項目が複数の属性を測っていたら[交差負荷があったら]、さらなる識別制約が必要になる。

 あるブロックの因子負荷が等しかったり、データから区別できなかったりすると、モデルが識別できなくなるだろう。たとえば\(\lambda_i=\lambda_k = \lambda\)のとき、効用差は$$ y^*_l = -\gamma_l + \lambda(\eta_a – \eta_b) + (\epsilon_i – \epsilon_k) $$となり、(たとえば)最後の項目の負荷\(\eta_d\)からの差しかわからなくなる。こういう不定性を突き止めるのはなかなか難しい(Mplusがなにか警告してくれたりするわけではない)。出力をよくよくみること。どこかの因子負荷がゼロに近いとか、因子間相関のSEが大きいとか、因子間相関がやけに大きいとか、Mplusに”the latent variable covariance matrix (Psi) is not positive definite”とかいわれたら要注意である[嗚呼、耳が痛い…]。ブロック内の因子負荷に等値制約を置いたり、因子間相関のどれか一つを実質理論的に期待される値に固定したりしてみるべし。

Mplusによるパラメータ推定と適合度検定
 \(y^*_l\)は観察できないので、二値アウトカムを扱えるソフトが必要である。Mplusとかね。
 誤差に相関があるのでFIMLは使えない。ULSかDWLS(MplusでいうULSMV, WLSMV)で推定する(どっちでも大差ない)。theta approachで推定すること。\(\mathbf{A}\)も\(\mathbf{\breve{\Psi}}^2 = \mathbf{A} \mathbf{\Psi}^2 \mathbf{A}^\top\)もフルランクでないので、”the residual covariance matrix (theta) is not positive definite”と文句を言われるが、これは仕方ない。
 テトラリック相関行列に対するモデルの適合度を検定できる。\(n > 2\)のときは自由度の修正が必要になる[…中略…]

個人得点の推定
 二値アウトカム\( y_l \)の項目特性関数(ICF)は$$ Pr(y_l = l | \eta_a, \eta_b) = \Psi \left( \frac{-\gamma_l + \lambda_i \eta_a – \lambda_k \eta_b}{\sqrt{\psi^2_i + \psi^2_k}} \right) $$ となる。モデルのパラメータが全部推定できたら、対象者の\(\mathbf{\eta}\)をMAP推定できる。Mplusがやってくれる。実は\(n \> 2\)のときICFは独立でないにもかかわらずMplusは局所独立性を仮定しているのだが、推定値の正確性への影響は小さい[←あ、なるほど。ここは気がつかなかった]。

チュートリアル: MplusのコードをExcelマクロでつくる
[パス。いずれ必要になったら読もう]

まとめ
 このチュートリアルとExcelマクロを使えば強制選択データを簡単に分析できるが、強制選択という形式とそのデータの「得意とするところ」については立ち止まってよく考えるべし。相対的判断から絶対的得点を復元するためには注意深い設計が必要である。たとえば1因子の場合、因子負荷にばらつきが必要である。いっぽう15因子以上の場合、すべての項目がpositively keyedでも(因子の相関が低ければ)復元できるといわれている。
 […後略…]
———-
 前半の定式化のみ読み、後半のチュートリアルをまるごととばしてしまったので、頭の体操になったという程度の感想しか持てないのだが、ま、必要になったら全部読もう。

 しかしなあ… 必要になることがあるのかなあ…
 この論文は、項目の比較課題で各項目の効用の背後にさらに少数因子を想定し、確認的因子分析モデルを置いて分析する、という話なんだけど、自分の仕事ではそういう場面をちょっと想像しにくい。
 なぜかを考えてみたんだけど、探索的場面と確認的場面にわけるとして、

  • 探索的場面ならば(つまり、少数因子に対する項目の負荷がわからない場面では)、たぶん評定課題でやってしまう。共通手法分散が怖いけれど、因子構造を調べるだけならそれほど実害がないのではなかろうか(共通手法因子のようなものが出るだけで)。あ、Maxdiff(Best-Worst法)で各項目の個人効用を無制約に推定し、それを探索的因子分析にかける、ということはあるな。個人効用推定値はipsativeなので(ふつう推定値の個人内平均は固定されている)、ほんとはおかしいんだけど、項目数が多いのでそこには目をつぶるわけだ。そこんところをThurstonian IRTモデルで正しく分析し項目の個人効用と因子得点を同時推定しますといわれたら、それは確かに魅力的だ。でもこの論文の主旨とはちょっと違う。
  • 確認的場面で比較課題を用いるときは、各因子あたり項目1個に絞り込んで、項目の個人効用を無制約に推定すると思う。思い浮かぶ理由は3つだ。
    • 回答負荷。消費者に比較課題を求めるのは大変なので、設問数を減らしたい。
    • 想定する因子数が多い。たとえばブランドイメージなら、8因子とか10因子なんてあたりまえだ。一対比較だと大変な設問数になってしまう。
    • 「因子あたり項目1個だ、測定誤差を分離できないのは以て瞑すべし」と割り切ることへの抵抗が小さい。実際、あるアウトカムに影響しそうな項目を片っ端からあつめて測定したのに(この段階では探索的因子分析を想定している)、アウトカムと個々の項目の関係を定量化したいから全項目を投入した回帰をやってくれ、といわれて困ってしまうことも少なくない。気持ちはわかる。因子分析を通じて各項目の独自性が捨てられてしまうのが怖いのだ。これは結局、関心の対象である個人特性についての実質的理論が乏しいことの反映だと思う。教育測定でも同じことが起きると思いますね。能力についての実質的理論が乏しい場面ほど、IRTへの抵抗感が強くなり、各テスト項目と外的基準の関係を細かーく調べたくなる。

 うーん、どうなんだろうか。ひょっとすると私の経験が乏しく視野が狭いだけかもしれない。この論文には適用例として、Forced-Choice Five Factor Markers(Brown & Maydeu-Orivares, 2011 Edu.Psych.Measurement), IRT-scoreed Occupational Personality Questionnaire (Brown & Bartram, 2009 Conf.), Customer Contact Styles Questionnaire (Brown, 2010 学位論文)というのがあげられている。