Pritikin, J.N. (2020) An exploratory factor model for ordinal paired comparison indicators. Heliyon, 6.
一対比較データを扱うRパッケージのひとつにpcFactorStanというのがあって、その元論文。このパッケージを使うつもりはあまりないんだけど、順序尺度の一対比較データから刺激の個人効用を推定する、それもStanで… という点には大いに関心があるので、試しに目を通してみた。
掲載誌は初めて見たけど、オープン・アクセス・ジャーナル。版元はCell Pressだから、そんなに変な雑誌ではない…と思いたい。
1. イントロダクション
相対的指標に基づく新しいタイプの因子分析モデルを提案する。
たとえば、「キャロブ・ミルクの味はココアミルクに比べて{良い, 悪い, 良くも悪くもない}」, 「ココアソイの飲み心地はココアミルクに比べて{良い, 悪い, 良くも悪くもない}」「キャロブ・ソイの色はキャロブ・ミルクに比べて{良い, 悪い, 良くも悪くもない}」…というような一連の質問があるとしよう。モノのペアを「味」「飲み心地」「色」といった項目(ファセット)について比較する質問群である。
ファセットがひとつならBradley-Terry-Luceモデルを使って分析するところだ。ファセットが複数あるとき、それらの相関をモデリングするという研究もある(Davidson & Bradley, 1969 Biometrika)。ここでやりたいのはファセットについての因子分析モデルである。それも相関行列を一旦推定してから因子分析モデルを推定するんじゃなくて(それだと不確実性が失われてしまう)、直接モデリングしたい。
一対比較データに対する因子分析という研究はすでにある。モノの価値を潜在変数とし、それが因子の指標であるとみなすとか[Thurstoninan IRTのことね]。モノ間の比較に関しては単一項目分析に留まっている。
社会的望ましさバイアスに対抗するため、あるブロックのなかに異なる潜在次元(ビッグファイブとか)の項目をいれて、ブロック内で順序づけする、という研究もある(Cheung, 2006 Edu.Psych.Meas.)。しかし、5つの特性の分散を潜在変数で説明しようとしているわけではない。[待って待って、頭が混乱してきた… この研究は、ヒトに項目ブロックの順序づけをさせて、項目を指標にしたCFAモデルを推定して、ヒトのビッグファイブ得点を出す、ってことですよね。これって多次元のThurstonian IRTですよね？]
以降の説明に当たっての注意。数式の表記はStanにあわせます。行列のコレスキー因子は上三角じゃなくて下三角です。正規分布のふたつめの引数は分散じゃなくてSDです。
2. モデル
2.1 項目モデル
モノの数を\(N\)とする。モノのインデクスを\(m, n\)とする(\(m < n\))。モノの潜在的価値を\(\mu_1, \ldots, \mu_N\)とする。
反応段階数\(H\)は3以上の奇数とする。閾値を\(\tau_1 < \cdots < \tau_{H-1}\)とし、\(\tau_0 = -\infty, \tau_H = \infty\)とする。
観察データは\(y_{mn} \in \{1, \ldots, H\}\)とする。もっとも\(m\)寄りの回答が\(1\), もっとも\(n\)寄りの回答が\(H\)である。その裏の潜在変数を\(Z_{mn}\)とする。累積分布関数を\(P\)として$$ \pi(Y_{mn} \leq y_{mn}) = P(\tau_{y_{mn}} + \mu_n – \mu_m)$$ である[\(\pi\)とは母比率のことかな]。順序効果はないものとする。
閾値は \(\Delta_d (d \in \{1, \ldots, D\}, D \equiv (H-1)/2) \)によって完全にパラメータ化できる。\(\delta_d \equiv \sum_{q=1}^d \Delta_q\)を\( (\delta_D, \ldots, \delta_1, -\delta_1, \ldots, -\delta_D)\)という順に並べたのを\(T\)として、$$ \pi(y_{mn} \leq h | \alpha, s) = \frac{1}{1+\exp(-\alpha(T_h + s(\mu_n – \mu_m)))}$$ とする。\(\alpha > 0\)は識別力パラメータ、\(s > 0\)はスケールである。すべての\(\Delta\)について\(\Delta > 0\)とする。
[原文では\(T\)ではなくて縦長の長方形みたいな記号が使われている。
それにしても、わっかりにくい説明だな… 読み解いておこう。
話を簡単にするために、回答は5件法だとしよう。段階反応モデルについておさらいする。5件法回答\(y\)の背後に潜在的な連続変数\(y^*\)があると考え、回答カテゴリが1となる確率のロジットは\(\tau_1 + y_{mn}^*\), 2以下となる確率のロジットは\(\tau_2 + y_{mn}^*\), 3以下となる確率のロジットは\(\tau_3 + y_{mn}^*\), 4以下となる確率のロジットは\(\tau_4 + y_{mn}^*\)だとする( \( \tau_1 < \tau_2 < \tau_3 < \tau_4)\)。言い換えると、回答カテゴリ1, 2, 3, 4, 5の間に閾値\(-\tau_1, -\tau_2, -\tau_3, -\tau_4\)があり、\(y^*\)が\(-\tau_h\)を超えると、回答がカテゴリ1からhまでに落ちる確率が0.5を超えるわけである。
この論文の場合でいうと、$$ logit(\pi(y_{mn} \leq h)) = \tau_h + y_{mn}^* $$ $$ \pi(y_{mn} \leq h) = logit^{-1} (\tau_h + y_{mn}^*) = \frac{1}{1+\exp(- (\tau_h + y_{mn}^*))} $$ となる。
一対比較回答の裏にある潜在変数\(y_{mn}^*\)は、ふたつのモノのあいだの価値の差が大きいときに大きくなるはずである。著者は、潜在変数は価値の差を正の定数倍したものだと考えている。その定数を\(\beta(> 0)\)とすれば、$$ \tau_h + y_{mn}^* = \tau_h + \beta(\mu_n – \mu_m) $$である。
この論文では、\(\tau_h + y_{mn}^*\)をわざわざ\(\alpha(T_h + s(\mu_n – \mu_m))\)と書き換えている(\(\alpha > 0, s > 0\))。本来なら定数\(\beta\)ひとつで済むものを、わざわざ\(\alpha\)と\(s\)に分けているわけだ。次の節で述べられるように、これは\(\mu_1, \ldots, \mu_N\)をスケーリングしたいからである。まあとにかく、\(T_h = \tau_h/\alpha\)であり、\(T_h\)もまた閾値に相当する定数である。
さて、著者によれば、5件法の場合\((T_1, T_2, T_3, T_4)\)は\((\delta_2, \delta_1, -\delta_1, -\delta_2)\)と書ける。つまり著者は、4つの閾値は0を中心として左右対称に並んでいると考えている。なるほど、一対比較だからね、自然な考え方だ。
で、著者は\(\delta_1\) (つまり0から\(T_2\)と\(T_3\)までの距離, 言い換えると回答カテゴリ3の幅の1/2) を\(\Delta_1\)と呼び、\(\delta_2 – \delta_1\) (回答カテゴリ2と4の幅)を\(\Delta_2\)と呼んでいる。5件法の場合、\(\Delta_1, \Delta_2\)さえ決まれば閾値が決まるわけである。なるほど。
著者いわく、\(\Delta_1 > 0, \Delta_2 > 0\)とする。あれれー？ \((\delta_2, \delta_1, -\delta_1, -\delta_2)\)は小さい順に並んでるはずだから、\(\Delta_1, \Delta_2\)は負値じゃないの？ 私がどこかで勘違いしているか、どこかに誤植があるのだろう(前者の確率が圧倒的に高い)。]
2.2 識別
\(\mu_n – \mu_m\)についての情報しかないわけだから、\(\mu_1, \ldots, \mu_N\)のうち推定できるのは\(N-1\)個だけである。\(\mu_1 = – \sum_{m=1}^N \mu_m\)という風にハードな制約をかけるか、事前分布\(\mu \sim N(0, \sigma)\)というようにソフトな制約をかけるか、である。StanでIRTモデルを推定するときはふつう後者である。ここでも後者を採用する。
ソフトは変数の分散が1に近くなるようにスケーリングするとうまく推定できることが多い。というわけで、識別力パラメータ\(\alpha\)とは別にスケーリングパラメータ\(s\)を用意する。同時に識別はできない。決め方は後述する。
2.3 相関と因子モデル
項目数を\(I\)とする。\(\mu_1, \ldots, \mu_N\)も\(\alpha\)も\(\Delta_1, \ldots\)も、カテゴリ数\(H\)までも、項目によって異なりうるとしよう。[まじか、項目間で閾値も変えるのか。なるほどね、その代わり\(\mu\)の分散を項目間で揃えるつもりなのだろう]
価値の行列を\(\mu\) (サイズ\((I, N)\)), 因子負荷行列を \(\lambda\)(サイズ\((I, F)\)), 因子得点行列を\(\eta\) (サイズ\((F, N)\)), 正規分布に従う残差を\(\epsilon\) (サイズ\((I, N)\))として $$ \mu = \lambda \eta + \epsilon$$ とする。
反応\(\pi(y_{i,mn})\)のモデル、\(I\)個の項目の相関のモデル、\(I\)個の項目の因子構造\(\mu\)のモデルのそれぞれについて尤度を求めることができる。以下では反応の対数尤度、つまり\(\log\pi(y_{i,mn})\)の合計を用いる。
2.4 Stanによるベイジアンな統計的推論
RパッケージpcFactorStanを作りました。なお、適合度診断については付録を見て下さい。
2.5 項目閾値
最初は閾値の事前分布について指数分布を考えてたんですが、どうもうまくいきませんでした。結局、すべての\(\Delta\)について\(\Delta = \iota (max(\mu) – min(\mu)), \iota \sim \beta(1.1, 2)\)としました。
2.6 尺度
まず、項目間の関係を無視し、項目ごとに\(s_i\)を決める。
[ここからよくわからんのでほぼ逐語訳]
分散\(\sigma \sim\) inv_gamma\((1,1)\)と標準化価値\(W_n = N(0, 1), n \in \{1, \ldots, N\}\)を推定する。これにより、スケールと位置タイプパラメータ\(\mu \leftarrow \sigma^{1/2} W\)を分離する、中心化されていないパラメタ化 が可能になる。なんらかの定数 \(j \geq 1\)について\(s \leftarrow sd(\mu)^j\)と設定することで、適応的事後分布は1.0において\(sd(\mu)\)を目指すようになる。ロジスティック関数を正規累積分布にあわせるため、\(\alpha=1.479\)と固定した。
[なに言ってんだかさっぱりわからないので、勝手に推測する。5件法だとして、ある項目だけに注目し、$$ \pi(y_{mn} \leq h) = logit^{-1}(1.479(T_h + s(\mu_n – \mu_m)) \ \ \mathrm{for} \ h=(1,2,3,4)$$ $$ T_4 = -T_1 = \Delta_1 + \Delta_2, \ T_3 = -T_2 = \Delta_1$$ $$ s = \mathrm{sd}(\mu_1, \ldots, \mu_N)^j $$ $$\Delta_d = \iota_d (\mathrm{max}(\mu_1, \ldots, \mu_N) – \mathrm{min}(\mu_1, \ldots, \mu_N)) \ \ \mathrm{for} \ d=(1,2)$$ $$ \iota_d \sim \beta(1.1, 2) \ \mathrm{for} \ d = (1,2)$$ $$ \mu_n = \sigma^{1/2} W_n \ \ \mathrm{for} \ n = (1,\ldots,N) $$ $$ \sigma \sim IG(1,1)$$ $$ W_n \sim N(0, 1) \ \mathrm{for} \ n = (1, \ldots, N)$$ というモデルを推定しました(\(j\)は定数)、以降では\(s\)の点推定値を使います、という話ではないだろうか？]
定数\(j\)は適応の強さをコントロールする。30人, 総当たり比較数450のシミュレーションをやった。\(\mu\)を正規分布からドローし(分散を10通り)、\(\alpha\)を適当に決め(10通り)、計100通りについて、\(j\)を\((1,\ldots,9)\)にして\(\mu\)の推定のRMSEを調べた。\(j\)は5くらいまで大きくすると良いようである。
[結局, \(j\)をどうセッティングしたときの\(s\)の点推定を使うのだろうか。\(j=5\)だろうか。わからん…なにもわからん…]
2.7 項目識別
今度は\(s\)を点推定値に固定し、\(\alpha\)を推定する。事前分布は\(N^+(1.749, 0.2)\)とする。
2.8 相関モデル
[ここ、わからんので原文をメモ]
Worths \(\mu_{N \times I} \leftarrow W_{N \times I} \Sigma^{T/2}\) where \(I\) dimension correlation matrix Cholesky factor \(\Sigma^{1/2} \sim lkj(2.5)\), and \(W_{N \times I}\) are element-wise \(N(0, 1)\).
[さすがにここはなにか誤植か脱字があるのではなかろうか。勝手に推測すると、各要素を\(N(0, 1)\)からドローて行列\(W\) (サイズ\(N, I\))をつくり、相関行列をLKJ分布 \(lkj(2.5)\)からドローして\(W\)の後ろから掛け、これを価値の行列\(\mu\)とする、ということではないかと思う。
わ・か・ら・ん。なんでここで項目間の相関行列をドローしないといけないのさ。その後ろに因子モデルを考えるんでしょ？ それとも、因子モデルを考えないときの話をしているの？]
2.9 因子モデル
行列\(\eta’\)(サイズ\((I, F)\))は各要素が\(N(0, 1)\)に従うとする。因子数が1ならばこれを因子負荷行列\(\eta\)とする。因数が2以上ならば相関行列を\(ljk(2.5)\)からドローして\(\eta’\)の後ろから掛け、これを\(\eta\)とする。
[因子得点は項目じゃなくてモノに与えられるんだろうから、サイズは\((F, N\)じゃないの？ これはさすがに誤植だろうなあ]
因子負荷行列\(\lambda\)(サイズ\(I, F\))は、各要素を\(\lambda\)として\( (\lambda+1)/2 \sim \beta(3,3)\)とする。
\(\epsilon\)は要素ごとに\(N(0,1)\)を事前分布とする、というわけにはいかない。\(\lambda \eta\)のパラメータ化でスケールと位置を掛けたからだ。そこで[…なんか書いてあるんだけどよくわからん。メモ省略]
因子負荷行列はそれだけでは解釈できないので […解釈できるように標準化するというような話。パス]
2.10 モデル適合度の評価
LOO交差検証をお勧めする。[…略…]
3. シミュレーション研究
[2頁ちょっと。読まずにとばした]
4. 適用例
[読まずにとばした]
5. 考察
今後の課題: 群間の測定不偏性の検討。測定誤差のもっと精緻な取り扱い(現状、比較反応が推移律を破っていてもモデルはそのことに敏感でない)。
有望な応用場面として、痛みの評価があげられる。痛みは多次元的なのに、従来の一対比較による推定では一次元性が仮定されていた。
云々。
——-
正直、私にとってはものすごおおおおくわかりにくううううい論文で、途中で心が折れてしまった。というか、これあまりに不親切でない？ 誤植も多くない？ 素人にもわかってもらおうという気持ちが足りなくない？(やつあたり)
ま、どういうデータを想定したどういうモデルなのかがわかったので良しとしよう。たくさんのヒトに、複数の評定項目についてモノとモノの一対比較(3件法とか5件法とか)を繰り返してもらったデータについて、評定項目を指標にしたEFAをやる、ってことなのね。当面このモデルを使うことはなさそうだし、必要になったらまた勉強する、ということで…
それにしても、ベイジアン・モデリングの際の事前分布の設定って、ときどき「シェフの秘伝のソース」みたいな奴が出てくることがあって、げんなりしてしまう。この論文でいうと項目閾値の事前分布とか。なんで\(\Delta = \iota (max(\mu) – min(\mu)), \iota \sim \beta(1.1, 2)\)にしようと思いついたのさ？ (やつあたり)