読了:Maydeu-Olivares & Bockenholt (2005) 順位付け・一対比較データをSEMで分析する(というかMplusで分析する)

Maydeu-Olivares, A., Bockenholt, U. (2005) Structural Equation Modeling of Paired-Comparison and Ranking Data. Psychological Methods, 10(3), 285-304.

 仕事の都合で読んだ奴。
 順位付けや一対比較のデータをSEMのソフトで分析しましょう、という解説論文。supplementary materialとしてMplusのコードも用意されている。

 いわく。
 選好・態度・価値の測定に一対比較や順位付けが広く使われている。理由:

  • 回答行動への制約が小さいから[←?]。特に、選択肢間の差が小さいときにもちゃんと情報が得られる。
  • 内的一貫性をチェックできるから。
  • 個人差とか選択肢間類似性について豊かな情報が得られるから。

 一対比較・順位付けについての統計モデルのほとんどはThurstone(1927)に遡る。彼は決定の結果の不整合性を捉えるため、決定とは確率的なものだと考えた。選択課題において、各選択肢はなにかの潜在的な連続的効用についての判断を引き出す。対象者はその瞬間に効用が最大である選択肢を選ぶ。効用は対象者を通じて正規分布する。彼のこの発想は、現在でいう潜在変数モデルである。
 Thurstoneはこのアプローチの下でいくつかのモデルを考えた。Case IIIモデルは、潜在効用が無相関だと仮定する奴。Case Vモデルはこれに加えて、潜在効用の分散が共通だと仮定する奴。
 現代では、Thurstone流の選択モデルはSEMの枠組みで捉えられるし、そうすることでモデルをいろいろ拡張できる。またSEMのソフトで簡単に推定できる。

 というわけで、本論文の目的は2つ。

  • Thurstonian選択モデルについてやさしく解説してあげます。
  • Thurstonian選択モデルをMplusでどうやって推定するか教えてあげます。

 先行研究。SEMでの一対比較・順位づけの扱いについてはMaydeu-Olivares(1999 Psychometrika; 2001 Psychometrika)をみよ。[…ほかにもいくつか挙げてるけど省略]
 本論文の新規性:

  • 一対比較だろうが順位付けだろうがほぼ同じ方法で推定しちゃいます。
  • パラメータ解釈とモデル選択について解説します。
  • SEMで推定します。選択肢が多くても大丈夫。

1. 比較判断の二値コーディング
 まずはコーディングについて。
 一対比較の場合、項目数を\(n\), 可能なペアの数を\(\tilde{n} = n(n-1)/2\)とする。ペア\(l\)が項目\(\{i,k\}\)だとして、\(i\)が選ばれたら\(y_l = 1\), \(k\)が選ばれたら\(y_l = 0\)とする。
 順位付けの場合も同様とする。つまり、\(n\)項目の順位付けの結果も、\(\tilde{n}\)個の項目ペアに対する一対比較として表現する。

2. 順位付けのThurstonian選択モデル
2.1 順位付けの反応モデル
 項目\(i\)の効用を潜在確率変数とし、これを\(t_i\)とする。\(t_i \geq t_k\)のとき\(y_l = 1\)、そうでないとき\(y_l = 0\)とする(\(t_i\)は連続なので等号はどっちにつけてもよい)。\(y_l\)には誤差項がないことに注意。
 これは、潜在変数\(y^*_l = t_i – t_k\)が0より上か下かで\(y_l\)が決まる、と書き変えてもよい。以下、潜在効用のベクトルを\(\mathbf{t}\)とし、長さ\(\tilde{n}\)の潜在ベクトル \( \mathbf{y}^* = \mathbf{A t} \)を考える。\(\mathbf{A}\)は\(\tilde{n} \times n\)の計画行列で、各行に1がひとつ、-1がひとつはいっているわけね。

 さて、Thurstoneのモデルでは、$$ \mathbf{t} \sim N(\mathbf{\mu}_t, \mathbf{\Sigma}_t)$$ と考える。すると\( \mathbf{y}^{*} \)もまたMVNに従い$$\mathbf{\mu}_{y^*} = \mathbf{A} \mathbf{\mu}_t $$ $$ \mathbf{\Sigma}_{y^*} = \mathbf{A} \mathbf{\Sigma}_t \mathbf{A}^\top$$ である。
 [ここまでの共分散構造が図にしてある。はいはい、ここまではわかります。指標数\(\tilde{n}\), 因子数\(n\), 因子の分散共分散を自由推定、指標にかならずパスが2本刺さる、パス係数は固定、というCFAになるわけね。指標に誤差項がないというところがポイントである]

2.2 順位付けモデルが含意する閾値とテトラコリック相関
 ここで指標はカテゴリカルである。SEMの場合、閾値とテトラコリック相関を推定することになる。
 \(\mathbf{y}^*\)を標準化した得点を \(\mathbf{z}^* = \mathbf{D}(\mathbf{y}^* – \mathbf{\mu}_{y^*}) \)とする。ただし\(\mathbf{D} = [Diag(\Sigma_{y^*})]^{-1/2}\)である。\(\mathbf{z}^*\)もMVNに従い、平均0、相関行列は$$ \mathbf{P}_{z^*} = \mathbf{D} ( \mathbf{\Sigma}_{y^*} ) \mathbf{D} = \mathbf{D}(\mathbf{A} \mathbf{\Sigma}_t \mathbf{A}^\top) \mathbf{D}$$ となる。これがテトラコリック相関行列という。
 \(z^*_l \geq \tau_l\)のとき\(y_l = 1\), そうでないとき\(y_l = 0\)とする。\(\tilde{n}\)個の閾値\(\tau_l\)を縦に並べたベクトル\(\mathbf{\tau}\)について$$ \mathbf{\tau} = – \mathbf{DA\mu}_t $$である。

 推定のプロセスは以下の通りである。(1)順序付けデータを一対比較の二値データに変換する。(2)閾値とテトラコリック相関行列を求める。(3)関心あるパラメータである\(\mathbf{\mu}_t, \mathbf{\Sigma}_t\)を推定する。

2.3 共分散構造
 潜在効用の共分散行列\(\mathbf{\Sigma}_t\)をいろいろ制約したり検定したりできる。Thurstoneの古典的な定義では次の3つがある:

  • 無制約。
  • Case III。非対角は0。
  • Case V。\( \mathbf{\Sigma}_t = \sigma^2 \mathbf{I} \)。

2.4 識別制約
 無制約なThurstonianモデルの場合、3つの制約が必要になる。

  1. \(\mu_t\)のうちひとつを0にする。
  2. 最後の要素に関する共分散を0にする(つまり、\(\mathbf{\Sigma}_t\)の非対角成分のうち一番下の行にあるのを0にする)。
  3. 最初と最後の潜在効用の分散を1にする(つまり、\(\mathbf{\Sigma}_t\)の非対角成分のうち右上の端と左下の端を1にする)。

 たとえば4項目の場合、固定したパラメータには上添字に\(*\)をつけるとして、$$ \mu_t = \left( \begin{array}{c} \mu_1 \\ \mu_2 \\ \mu_3 \\ 0^* \end{array} \right) $$ $$ \Sigma_t = \left( \begin{array}{cccc} 1^* & \sigma_{21} & \sigma_{31} & 0^* \\ \sigma_{21} & \sigma^2_2 & \sigma_{32} & 0^* \\ \sigma_{31} & \sigma_{32} & \sigma^2_3 & 0^* \\ 0^* & 0^* & 0^* & 1^* \end{array} \right) $$ とする。

 Case IIIの場合は、\(\mu_t\)のうちひとつを0, \( \mathbf{\Sigma}_t\)の対角要素のうちひとつを1にする。非対角はすべて0。
 Case Vの場合は、\(\mu_t\)のうちひとつを0, \( \mathbf{\Sigma}_t\)の対角要素をすべて1にする(つまり\(\sigma^2 = 1\))。非対角はすべて0。

3. 順位付けデータのSEM
 では、これをSEMにどう組み込むか。

3.1 SEMモデルとしてのThurstonian順位付けモデル
 Mplusやlisrelに従い、SEMのモデルを以下とする。
$$ \mathbf{y}^* = \mathbf{\nu} + \mathbf{\Lambda \eta} + \mathbf{\epsilon} $$ $$ \mathbf{\eta} = \mathbf{\alpha} + \mathbf{B \eta} + \mathbf{\xi} $$ Thurstonian選択モデルは、\(\mathbf{\nu} = \mathbf{0}\)、\(\mathbf{\Lambda} = \mathbf{A}\)、\(\mathbf{\eta} = \mathbf{t}\)、\(\mathbf{\epsilon} = \mathbf{0}\)、\(\mathbf{\alpha} = \mathbf{\mu_t}\)、\(\mathbf{B} = \mathbf{I}\)である。

 SEMにおける \(\mathbf{y}^*\)の平均と分散は $$ \mathbf{\mu}_{y^*} + \mathbf{v} + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\alpha}$$ $$ \mathbf{\Sigma}_{y^*} = \mathbf{\Lambda}(\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B})^{-1\top} \mathbf{\Lambda}^\top + \mathbf{\Theta} $$ Thurstonian選択モデルでは、\(\mathbf{\Psi} = \mathbf{\Sigma}_t\)、\(\mathbf{\Theta} = \mathbf{0}\)である。

3.2 二値観察変数の場合のSEMモデルの推定
 LisrelなりMplusなりで推定する場合は、まず閾値とテトラコリック相関を推定し、次にモデルパラメータを推定する。
 Mplusの場合、モデルパラメータの推定の際には以下を最小化する: $$ F = [\hat{\kappa} – \kappa(\theta)]^\top \hat{\mathbf{W}} [\hat{\kappa} – \kappa(\theta)] $$ ここで\(\hat{\kappa}\)はすべての閾値とテトラコリック相関の推定値を縦に積んだもの、\(\kappa(\theta)\)はモデルパラメータ\(\theta\)による閾値・テトラコリック相関への制約を表す。\(\theta\)とはここでは\(\mathbf{\mu}_t\)と\(\mathbf{\Sigma}_t\)である。
 では\(\hat{\mathbf{W}}\)はなにか。2つの方法がある。
 標本の閾値とテトラコリック相関の漸近共分散行列を\(\mathbf{\Xi}\)としよう。ひとつはWLS推定で、\(\hat{\mathbf{W}} = \hat{\mathbf{\Xi}}^{-1} \)とする。もうひとつはDWLSで、\(\hat{\mathbf{W}} = [Diag(\hat{\mathbf{\Xi}})]^{-1} \)とする。どちらの方法でも、一致性のある漸近正規推定量が手に入る。
 標本サイズを\(N\), カイ二乗推定量を\(T = N \hat{F}\)とする。閾値・テトラコリック相関に課した制約について検定するためには…[以下略]

 [ここでSupplementary materialsに目を通したのだけど、Mplusに即していえばこういうことらしい。
 Mplusで二値データを扱うとき、データからわかるのは閾値とテトラコリック相関であって、ここから共分散構造を推定しないといけないわけだ。測定誤差\(\epsilon\)の共分散\( \Theta \)について2つの考え方がある。

  • \( \Theta \)の対角要素をパラメータではなく、元の共分散行列の関数とみて、$$ \Theta = I – diag( \Lambda(I-B)^{-1} \Psi (I-B)^{-1\top} \Lambda^\top) $$ とするやりかた。(よくわからんがdelta approachのことであろう)
  • モデルの下での\(y^*\)のSDの逆数を持つ対角行列を\(\mathbf{D} = (Diag(\Sigma))^{-1/2} \) として、$$ \mathbf{D} (\nu + \Lambda (\mathbf{I} – \mathbf{B})^{-1} \alpha), \ \ \mathbf{D} ( \Lambda (\mathbf{I}-\mathbf{B})^{-1} \Psi (\mathbf{I} – \mathbf{B})^{-1\top} \Lambda^\top + \Theta) \mathbf{D}$$を使う方法。(theta アプローチについて説明しているのだと思うが、使うという意味がよくわからない…)

Thurstonianモデルでは後者を使わないといけない。Mplusではこれをthetaアプローチという。
 推定量としてはDWLS(MplusでいうWLSM)がお勧め。切片\(\nu\)は0に固定し、因子平均\(\alpha\)を最後のひとつを除き自由推定すること。因子負荷は正しく固定すること。\(\Theta\)の対角は自由推定すること(←これは一対比較の場合の話だろう)。]

4. 順位付けへの適用例: スペインの心理学専攻学部生のキャリア選好のモデリング
[メモ省略]

5. Thurstonian一対比較モデル
 一対比較の場合は非整合的な回答がありうるから、誤差項をつける。$$ y_l^* = t_i – t_{i’} + e_l $$
\(e_l \sim N(0, \omega^2_l) \)とし、他のペアとも潜在効用とも無相関とする。行列で書くと $$ \mathbf{y}^+ + \mathbf{At} + \mathbf{e} $$ \(\mathbf{t}, \mathbf{e}\)がMVNなんだから\(\mathbf{y}^*\)もMVNであり、$$ \mathbf{y}_{y^*} = \mathbf{A} \mathbf{\mu}_t $$ $$ \mathbf{\Sigma}_{y^*} = \mathbf{A} \mathbf{\Sigma}_t \mathbf{A}^\top + \mathbf{\Omega}^2 $$ \(\mathbf{\Omega}\)の対角は自由推定してもいいし、\(\mathbf{\Omega}^2 = \omega^2 \mathbf{I}\)と制約してもよい。
 潜在的な差の相関行列は$$ \mathbf{P}_{z^*} = \mathbf{D} ( \mathbf{\Sigma}_{y^*} ) \mathbf{D} = \mathbf{D}(\mathbf{A} \mathbf{\Sigma}_t \mathbf{A}^\top + \mathbf{\Omega}^2) \mathbf{D}$$ となる。
 … というのが順位付けとのちがい。あとは同じである。

 [ここから大事なところなので、逐語訳に近い形でメモを取り直しました]
 いかなるSEMモデルも、全体的な適合度という観点からは区別できない複数の等価なモデルを持っている。しかしそれらのモデルは実質的解釈が異なるかもしれない(MacCallum et al., 1993 Psych.Bull.)。長年の研究を通じて、等価なSEMモデルをみつけるための一連のルールが開発されている(Lee &amp Hershberger, 1990 MBR; Stelzl, 1986 MBR)。Thurstonian選択モデルの場合、所与の推定されたモデルと等価なモデルの集合を見つけるためのルールをTsai(2003 Psychometrika)が提案している。
 一対比較のThurstonianモデルを推定したとしよう。このモデルにおける潜在効用の共分散行列を\(\Sigma_1\), 誤差共分散行列を\(\Omega^2_1\)とする。どちらも正則である。このときそのモデルは、別のモデル$$ \Sigma_2 = c \Sigma_1 + \mathbf{d}\mathbf{1}^\top + \mathbf{1} \mathbf{d}^\top$$ $$ \Omega^2_2 = c \Omega^2_1 $$ と等価である(\(c\)は任意の正の定数、\(\mathbf{d}\)は任意の定数ベクトル。ただし\(\Sigma_2, \Omega^2_2\)はどちらも正則でないといけない)。同じことは順位付けでも生じる(ただし\(\Omega^2_1=\mathbf{0}\))。
 [数値例。中略]
 等価なモデルが存在するということは、Thurstonian選択モデルの実質的解釈に対して重要な意味を持っている。
 まず、Case IIIモデルとCase Vモデルの解釈に際して問題が生じる。いまCase Vモデルの適合が良かったとしても、それと等価でありかつ潜在効用がすべて相関しているモデルが存在するわけだから、潜在効用が無相関だと主張することはできない。主張できるのは、このデータのもっとも倹約的な表現がCase Vモデルでしたということ、そしてそのもでるによればデータは「潜在効用が無相関だ」という仮説と整合しているということ、だけである。この解釈を検証するためには追加データを集める必要がある。Bockenholt(2004 Psych.Methods)をみよ。
 等価なモデルの存在は、無制約モデルの解釈にも影響する。モデルの共分散の推定値の大きさを解釈することはできない。なぜなら、適合度が同じだが異なる推定値を持つ等価なモデルが常にありうるからである。推定された共分散は、最初と最後の選択肢の間の共分散 \( \sigma_{1n} = 0 \)に対して相対的に解釈することしかできない。この共分散は、推定された共分散のスケールをセットしている。正の共分散は「\( \sigma_{1n}\)より関連が強い」、負の共分散は「\( \sigma_{1n}\)より関連が弱い」と解釈すべきである。
 この点を理解するために、共分散行列を潜在効用の平方距離の行列 \(\Delta\)に変換してみよう。その要素は$$ \delta^2_{ik} = \sigma^2_i + \sigma^2_k – 2 \sigma_{ik} $$ である。さっきの変換に出てきた任意の定数\(c\)について\(\Delta / c\)を求めると、これは変換に対して不変である。従って、共分散は相対的に解釈すべきだが、相対的距離は一意に解釈できる。正の共分散は効用間の距離が小さいこと、負の共分散は効用間の距離が大きいことを示す。
 このように、異なるThurstonianモデルのパラメータ解釈と比較に際しては十分な注意が必要である。推定されたパラメータ(平均, 分散, 共分散)は相対的にしか解釈できない。Case IIIモデルやCase Vモデルがデータの倹約的表現であることがわかっても、潜在効用が独立だと推論することはできない。しかし、無制約モデル・Case IIIモデル・Case Vモデルのあいだの区別は依然として有用である。それは効用の共分散行列がどのような倹約的な形式を持っているかを調べる際の出発点になる。

6. 一対比較への適用例: コンパクト・カーの選好のモデリング
[メモ省略]

7. 一対比較・順位付けのためのThurstonian因子モデル
 こんどは因子に構造をいれます。そうすることの利点:

  • 共分散行列を無制約にするよりも倹約的なので、パラメータの推定精度が上がる。
  • 因子間の類似性の構造が視覚されるので解釈しやすい。
  • 項目間の相対的距離は不変なので解釈しやすい。

 探索的因子分析モデルを考えます。$$ \mathbf{t} = \mathbf{\mu}_t + \mathbf{\Lambda}_t \xi + \varsigma$$
\(\xi\)は\(m\)次元の共通因子で、平均0、分散1, 互いに無相関とする。 要は二次因子分析モデルのようなモデルである。
 独自因子\(\varsigma\)は互いに無相関で、平均0, 分散\(\Psi^2_t\)の正規分布に従うとする。潜在効用の共分散は$$ \Sigma_t = \Lambda_t \Lambda_t^\top + \Psi^2_t $$ となる。
 \(\mathbf{P}_{z^*}\)はどうなるかというと… [メモ省略]

 識別制約の入れ方。因子負荷行列\(\Lambda_t\)の右上の三角を0に制約するのが簡単である(つまり, 共有因子数が3だったら、一番右上とその左・その下の3つが0になる)。さらに、\(\Lambda_t\)の一番下の行を0、\(\mu_t\)の一番下を0, \(\Psi^2_t\)の一番右下を1に制約するのがお勧め。
 共通因子数1,2,3,4のモデルの識別のためには項目数が5,6,8,9個必要である。どう考えればよいかというと、$$ \Sigma_t = \Lambda_t \Lambda_t^\top + \Psi^2_t$$で共通因子数が\(m\)なら、パラメータは\(\Psi^2\)に\(n-1\)個、\(\Lambda_t\)に\(2(2n-m-1)/2\)個ある。しかしその合計は無制約モデルにおける\(\Sigma_t\)の識別可能パラメータの数\( ((n(n-1))/2)-1 \)を超えられない。[この辺、考えるのが面倒なので写経している]

 ふつうの因子分析とちがって、一対比較の場合には項目の平均にすごく関心がある。これを共通因子に依存させる定式化もできる。\(\mu_t\)を共通切片\(\mathbf{1}\gamma\)に置き換えて$$ \mathbf{t} = \mathbf{1} \gamma + \mathbf{\Lambda}_t \mathbf{\xi} + \mathbf{\varsigma} $$ とし、\(\xi\)の平均\(\mu_\xi\)を推定するわけである。識別のためには\(1\gamma\)の最後の要素を0に固定すること。

8. 一対比較への適用例・再訪: コンパクトカーの選好を因子モデルでモデリングする
[メモ省略]

9. 結論
 モデル選択に当たっては次のステップをお勧めします。

  1. まず無制約のモデルを推定する。項目ペアの誤差分散を等値にするモデルと項目ペア別にするモデルを比較して、どっちを採用するか決める。どっちも当てはまりが悪かったら別のモデルを探す。
  2. Case IIIモデルを推定して、無制約モデルと比べる。
  3. 無制約よりCase IIIのほうが適合度が高かったら、次はCase Vに進む。
  4. Case IIIより無制約のほうが適合度が高かったら、因子分析モデルを試す[ああ、なるほど… そういう位置づけなのか]。因子モデルが良かったら、今度は効用平均を因子に依存させるモデルを試す。

[整理しておくと、無制約モデルに対する制約の方向として、選択肢の効用が無相関だと仮定する方向と、選択肢の効用の背後に共通因子があると仮定する方向がある。著者らは、このふたつの方向のどっちに進むべきかを適合度の観点で決めましょう、といっているわけだ。
 うーん、それはちょっと違和感があるなあ。その決定には既有知識や分析の目的が強く効くのではなかろうか。自分の仕事に置き換えると、たとえば n 個の製品への消費者の選好を一対比較で測定したとして、製品の知覚マップを描きたいというような話だったら因子分析モデルに進むだろうけれど、個人効用を層別集計したいというような話だったら、仮に因子分析モデルの適合度が高かったとしても無制約モデルを採用すると思う。個人効用推定に際して余計な仮定が入っているのが嫌だからだ。でもあれかな、「データ生成メカニズムを因子分析モデルでよりよく表現できることがわかっているのに、なぜその知見を活用せず無制約モデルを採用するのか、そのせいでパラメータの推定精度が落ちるではないか」と云われたら、ちょっと困っちゃうかな…
 著者の先生方に訊いてみたいものだ。いま心理学かなんかの研究をしているとして、対象者の選択肢に対する個人効用を一対比較で測定したとしますね。効用の選択肢間類似性とかにはまるきり関心がなくて、単にその次のステップのために効用を推定したいだけ、推定した効用で対象者を層別したり、他の変数との関連を観察したり、特定の選択肢への効用が高い/低い人を選んで別の調査をかけたいだけ、だとしますね。一対比較回答についてモデルを組んでみたら、無制約モデルより因子数3の因子分析モデルのほうが適合度がよかったとしますね。なんでか知らんけど因子数3の適合度が高いわけ。そのとき、因子数3の因子分析モデルで推定した効用を採用します? そこまで強気になれます?]
——————
 なかなか親切な論文であった。誉めて遣わすぞ。
 一対比較課題で選択肢の価値を個人別に推定したいんだけど、個人あたりの試行回数が少なくて古典的なBradley-Terryモデルではうまくいかないとき、要は個々の回答に+1と-1のパスが刺さるカテゴリカルCFAじゃん… と考えて分析していたのだが、その発想が正しいことがわかって嬉しかった。今ふと気が付いたけど、そういう分析方法は豊田先生の本にも載っていたような気がするな… 別に私が自分で思いついたわけではなかったかもしれない。
 いっぽう、識別制約の掛け方については大変勉強になりました。勤務先の若い人に一対比較データをStanでモデリングしてもらったとき、なかなか収束しないんですよとぼやいているので不思議に思っていたのだが、その理由がようやく理解できた。

 恥かしながら、よくわからなかった点:

  • 効用の共分散行列\(\Sigma_t\)にかける識別制約には実質的なインパクトがあるのではないか。無制約の場合、最初と最後の潜在効用の分散を1, 最後の潜在効用の共分散をすべて0にせよというお勧めだが、それって因子の並べ方で結果が変わるということを意味していないだろうか? 共分散行列が変わるのはいいとして、そこから算出される潜在効用の類似性の行列 \( \Delta / c\)も変わっちゃうのは困る。参るなあ。ベイズ推定だと制約をかけなくてもなんとなく収束しちゃうんだけど…
  • MplusではPARAMETERIZATION = THETAと指定しないといけないそうだが、なぜか。一対比較でカテゴリカル指標の誤差分散を推定したいからだというのならわかるんだけど、順位付けのときにもそう指定しないといけない理由がわからない。

 ところで本文中の注で、効用utilityという言葉について短く解説してあって(心理学の論文だからね。ちょっと違和感があるかなと思ったのだろう)、Kahneman(2003, Am.Econ.Rev.)というのが引用されていた。これ面白そうだな。

2022/12/19追記: 原文に当たり直してちょっと追記した。