Tsai, R.C, Bockenholt, U. (2002) Two-level linear paired comparison models: estimation and identifiability issues. Mathematical Social Sciences. 43, 429-449.

　もういいかげんうんざりしているのだけれど、ひきつづき一対比較の論文である。
　Thurstonian一対比較モデルの識別性について論じている Tsai(2003) で、ペア特有誤差の分散が一定なときの識別条件についてこの論文がreferされていたので、読んでみた。なお、Tsai(2003)ではペア特有誤差の分散を制約しない識別条件がこの論文に頼らずに示されているし、ペア特有誤差がないときの識別条件ならTsai(2000) にわかりやすい導出があるので、べつにこれを読む必要はないんだけど、毒も喰らわば皿まで、ということで…

1. イントロダクション
[略]

2. 線形一対比較モデル
　\(r\)個の項目に対する完全一対比較実験を考える。評定者の反応過程についてのモデルは山ほどある。Bossuyt(1990 書籍), David(1988 書籍)をみよ。一番よく使われているのは線形反応モデルで、その代表がBradley-Terry-LuceモデルとThurstonianモデルである。
　評定者が複数いるときは、評定者内分散と評定者間分散の両方を考えないといけないはずだが、伝統的には後者は無視されてきた。そこでご提案します。2レベルのモデルで表現しましょう。

2.1 評定者内モデル
　ヒト\(i\)にとっての項目\(j\)の効用の平均を\(\mu_{ij}\)とする。項目\(j\)と\(k\)の比較を$$ y_{ijk} = \mu_{ij} – \mu_{ik} + \varepsilon_{ijk}$$ とする。\(\varepsilon_{ijk}\)が正規分布するならThurstonianモデル、ロジスティック分布するならBTLモデルになるし、Stern(1990 JASA)のようにガンマ分布を考えれば他のモデルも導出できる。
　\(r\)個の項目の一対比較は\(\mathbf{A}^*_i \mu_i + \varepsilon_i \)と書ける。\(\mu_i\)は\(r\)次元ベクトル、\(\varepsilon_i\)は\(C(r, 2)\)次元ベクトル。\(\mathbf{A}_i^*\)は\(C(r,2) \times r\)の計画行列で、ランクは\(r-1\)。右端列を削ってフルランクにした行列を\(\mathbf{A}\)とする。
　なお、\(\mathbf{A}_i\)はspans only the main effect spaceなわけだけど[なんて訳せばいいの？ 「主効果の空間を張っている」？ まあとにかく主効果しか考えてないってことね]、対応する交互作用空間を得ることもできる。その行列\(\mathbf{B}_i\)は\((r-1)(r-2)/2\)列となる。
　反応はたいてい二値か順序なので、カテゴリ数\(H\)として閾値\(\tau_0 < \tau_1 < \cdots, < \tau_H\)を考え、\(\tau_0 = -\infty, \tau_H=\infty\)として、\(\tau_{h-1} < y_{ijk} < \tau_h\)のときに反応\(h\)が得られると考える。項目の左右の効果がないと仮定すれば\(\tau_h = -\tau_{H-h}\)となる。二値なら\(\tau_1 = 0\)ということになる。

2.2 評定者間モデル
　個人効用を固定要素とランダム要素に分けて、$$ \mu_i = \mu + \nu_i$$ とする。項目の共変量\(\mathbf{X}\)があったら、\(\mu = \mathbf{X} \beta\)とするか、\(\nu_i = \mathbf{X} \zeta_i + \delta_i\)とする。\(\nu_i\)ないし\(\zeta_i, \delta_i\)はなにかパラメトリックな族から選ぶか、データから決める。

2.3 normal-normalの場合
　\(\varepsilon_i \sim N(0, \sigma^2 \mathbf{I}), \nu_i \sim N(0, \Sigma_\nu)\), 両者は独立としよう。このとき$$ \left[ \begin{array}{c} \mathbf{y}_i \\ \nu_i \end{array} \right] \sim N \left( \left[ \begin{array}{c} \mathbf{A}^* \mu \\ \mathbf{0} \end{array} \right], \left[ \begin{array}{cc} \mathbf{A}^* \Sigma_\nu \mathbf{A}^{*\top} + \sigma^2 \mathbf{I} & \mathbf{A}^* \Sigma_\nu \\ \Sigma_\nu \mathbf{A}^{*\top} & \Sigma_\nu \end{array} \right] \right)$$ となる。共変量があったら、\(\zeta_i \sim N(0, \Sigma_\zeta), \delta_i \sim N(0, \sigma^2_\delta \mathbf{I})\), 両者は独立として、\(\Sigma_\nu = \mathbf{X} \Sigma_\zeta \mathbf{X}^\top + \sigma^2_\delta \mathbf{I}\)となる。

2.4 ランダム効果のその他の分布
　normal-normalじゃないモデルもつくれる。logistic-normalとか(\(\mathbf{y}, \nu\)の同時分布が閉形式にならないけど)。

3. 識別
　一般に、潜在反応\(\mathbf{y}\)はスケールも位置も識別できない。ここでは、位置を固定し(例, \(\mu_r = 0\))、個人内反応のスケールを固定することにする(例, \(\sigma^2 = 1\))。
　ランダム効果の分布の指定によるが、それでも不定性が残る場合もある。ここでは効用が正規分布する場合のいくつかの必要条件について論じる。
　データは離散値的な構造を持っているので、\(\mathbf{V} = \mathbf{A}^* \Sigma_\nu \mathbf{A}^{*\top} + \Sigma^2 \mathbf{I}\)は識別できず、\(\mathrm{Corr}(\mathbf{V})\)のみが識別できる。以下では\(\mathrm{Corr}(\mathbf{V})\)の識別可能な構造について論じる。
　説明のため、\((r-1) \times r\)の対比行列\(\mathbf{C} = [\mathbf{I}_{r-1} | \mathbf{-1}]\)を考える。

系1. \(\Sigma_1, \Sigma_2\)を2つの\(r \times r\)共分散行列とする。\(\mathbf{C} \Sigma_1 \mathbf{C}^\top = \mathbf{C} \Sigma_2 \mathbf{C}^\top\) が成り立つ必要十分条件は、なんらかの実数ベクトルを\(\mathbf{d}\)として\(\Sigma_1 = \Sigma_2 + \mathbf{d} \mathbf{1}^\top + \mathbf{1}\mathbf{d}^\top\)が成り立つことである。

証明 Tsai(2000)をみよ。

系2. \(\mathbf{A}\)を\(C(r, 2) \times (r-1)\)の列フルランクな一対比較計画行列とする。\(\mathbf{A} \Sigma_1 \mathbf{A}^\top = \mathbf{A} \Sigma_2 \mathbf{A}^\top\)が成り立つ必要十分条件は、\(\Sigma_1 = \Sigma_2\)が成り立つことである。

証明 \(\mathbf{A}\)は列フルランク行列、すなわち、\(\mathbf{M}\)の列空間を\(C(\mathbf{M})\)と書くとして\(C(\mathbf{I}_{r-1}) \subseteq C(\mathbf{A})\)であるから、\(\mathbf{A} \Sigma_1 \mathbf{A}^\top = \mathbf{A} \Sigma_2 \mathbf{A}^\top\)なのは\(\Sigma_1 = \Sigma_2\)であるとき、そのときに限る。[うへえ… いわれてみりゃそうなんだけど、よく思いつくね…]

補題1 \(\Sigma_1, \Sigma_2\)を2つの\(r \times r\)共分散行列とする。 \(\mathbf{A}\)を\(C(r, 2) \times (r-1)\)の列フルランクな一対比較計画行列とする。\(\mathrm{Corr}(\mathbf{A} \mathbf{C} \Sigma_1 \mathbf{C}^\top \mathbf{A}^\top + \mathbf{I}) = \mathrm{Corr}(\mathbf{A} \mathbf{C} \Sigma_2 \mathbf{C}^\top \mathbf{A}^\top + \mathbf{I}) \)が成り立つのは、なんらかのベクトル\(\mathbf{d}\)について\(\Sigma_1 = \Sigma_2 + \mathbf{d} \mathbf{1}^\top \mathbf{1} \mathbf{d}^\top \)が成り立つとき、そのときに限る。

証明　Appendixをみよ。[みましたけど、後件⇒前件は系1のせいで一瞬でおわるが、前件⇒後件はめっちゃ面倒な話になっている… ついていけない…]

　以上を踏まえるとこういうことになる。\(\sigma^2\)を定数に固定すれば、\(\Sigma_\nu\)と観察上等価な\(\Sigma^*_\nu\)は\(\Sigma^*_\nu = \Sigma_\nu + \mathbf{d} \mathbf{1}^\top + \mathbf{1} \mathbf{d}^\top\)である。つまり、識別できるのは\(\alpha = \mathbf{C}\mu, \upsilon = \mathbf{C}\nu, \Sigma_\nu = \mathbf{C} \Sigma_\upsilon \mathbf{C}^\top\)だけである。たとえば\(\Sigma_\nu = \mathbf{I}\)という仮説をほかの\(\Sigma^*_\nu = \Sigma_\nu + e \mathbf{11}^\top\)という仮説と区別することはできないわけだ。
　項目間類似性に関心があるんだったら、\(\Sigma_\nu\)の対角をなんらかの定数\(c\)に固定するのがお勧めである。

　いっぽう、共変量がある場合は話がかわってきて[…中略…]。

4. 推定
　反応\(W_{ijk}\)が\(h\)である確率は、レベル1のリンク関数を\(\Psi\)として$$ p_{ijk}(h; \theta_i) = Pr(W_{ijk} = h | \theta_i) = \Psi[\tau_h – (\mu_{ij} – \mu_{ik})] – \Psi[\tau_{h-1} – (\mu_{ij} – \mu_{ik})]$$ である。もし\(\varepsilon_i\)が正規分布なら\(\Psi\)は標準正規累積分布関数ね。
　反応が\(\theta_i\)のもとで条件付き独立だとすれば、\(z_{ijkh}\)を反応のダミー変数として$$ Pr(W_{ijk} = w_{ijk}; j < k) = \int^\infty_{-\infty} \cdots \prod_{j=1}^{r-1} \prod_{k > j}^{r} \prod_{h=1}^{H} [p_{ijk}(h; \theta_i)]^{z_{ijkh}} g(\nu_i\; \Sigma_\nu) d \nu_i$$ となる。こう書くとややこしいけど、normal-normalならこうである。$$ Pr(W_{ijk} = w_{ijk}; j < k | \theta_i) = \Psi([\mu_{w_{i12}}, \mu_{w_{i13}}, \ldots, \mu_{w_{i(r-1)r}}], \mathbf{V}) $$ ただし$$ \mu_{w_{ijk}} = \frac{\tau_{w_{ijk}} – (\mu_j – \mu_k)}{\sqrt{\nu_{jk}}} $$ であり、\(\mathbf{V}\)は一対比較の共分散行列で、その対角は\( \upsilon_{jk} = \sigma^2 + \sigma^2_{\nu_j} – 2 \sigma_{\nu_{jk}} + \sigma^2_{\nu_{k}} \) である。[記号の使い方のせいでなんだかめんどくさくなっておる…]

　推定について。項目数が6までなら数値積分で解ける。それより大きいときの推定方法についていろいろ提案があるんだけど、ここではモンテカルロEMアプローチ(MCEM)を使う。

5. 事例
　[きちんと読んでないけど…]
　Schwqartz & Bilsky(1987)の基本的価値の10次元ってのがある。powerとかachiavementとかなんとか。この10次元が2次元空間上で放射線状になっている。あくまで理論上はという話だけど。
　各次元あたり1項目、総当たりで45ペアの一対比較課題をやった。7件法評定。次の3つのモデルを仮定した。

• (a) \(\Sigma_\nu = \Sigma_\nu^2 \mathbf{I}\)
• (b) \(\Sigma_\nu = \mathbf{X} \Sigma_{\zeta} \mathbf{X}^\top + \Sigma^2_\delta \mathbf{I} \)と仮定し、\(\mathbf{X}\)を「単位行列の対角よりも4つ右側まで1を埋めた」行列とし、\(\Sigma_{\zeta}\)は「対角要素のうち5つだけが自由推定であとは0」とする。[んんん？ これ、SEMでいうとどういうこと？ 10個の1次因子{F1, …, F10}の後ろに5個の2次因子{G1, …, G5}を考えて、G1はF1,F2,F3,F4,F5に負荷1, …, G5はF5, F6, F7, F8, F9に負荷1, … あれ、辻褄が合わない…???]
• (c) 無制約の共分散行列

尤度比検定すると(c)の勝ち。
　解釈の際には、生の\(\hat{\Sigma}_\nu\)ではなくて、対角成分が1だと仮定したときの\(\hat{\Sigma}_\nu\)をみるのがよい[うーん、これはそういう制約をつけて推定し直すってこと？ それとも単に\(\hat{\Sigma}_\nu\)を相関行列に変換するってこと？]。理論的に想定された円状の相関に近かったけど、ちょっと例外があった。
　[実は7件法は多すぎるようだという話。中略]

6. 結論
[略]
—————–
　うーむ。
　この話題もいいかげん飽きてきたので、一対比較モデルの識別性について、昨年秋ごろから延々と悩んでいる点をメモしておく。恥をさらしているような気もするけれど。
　Thurstonian一対比較モデル(本論文でいうnormal-normalな線形一対比較モデル)を想定する。対象者\(i\)の項目\(j\)への平均効用を\(u_{ij}\)として、項目ペア\((j, k)\)への潜在反応\(y_{ijk}\)が$$ y_{ijk} = u_{ij} – u_{ik} + \varepsilon_{ijk}$$ \(\varepsilon_i = (\varepsilon_{i11}, \ldots, \varepsilon_{iJ(J-1)})^\top, \ \mathbf{u}_i = (u_1, \ldots, u_J)^\top\)として$$\varepsilon_i \sim_{iid} N(\mathbf{0}, \sigma^2 \mathbf{I}), \ \ \mathbf{u}_i \sim_{iid} N(\mu, \Sigma)$$ だとしますね。
　疑問点:

• 私の場合、関心があるのは\(\sigma^2, \mu, \Sigma\)ではなく、むしろ個人効用の行列\(\mathbf{U}\)、つまり、行\(i\)が\(\mathbf{u}^\top_i\)であるような行列なのです。いま、なんらかの識別条件を課して、パラメータ推定値\(\hat{\sigma^2}, \hat{\mu}, \hat{\Sigma}\)を得て、そこから個人効用の推定値の行列 \(\hat{\mathbf{U}}\)を得たとしますね。でも、いや個人効用の共分散行列は \( \hat{\Sigma}’ = \hat{\Sigma} + \mathbf{d}\mathbf{1}^\top + \mathbf{1}\mathbf{d}^\top\)だってことにしようよ、そっちのほうが解釈しやすいからさ、という話になったら、個人効用の推定値の行列も\(\hat{\mathbf{U}}’\)に代わる。\(\mathbf{U}, \hat{\mathbf{U}}, \hat{\mathbf{U}}’\)の間にはどういう関係があるのだろうか？ 直感としては、\(\mathbf{u}_i\)と\(\hat{\mathbf{u}}_i\)や\(\hat{\mathbf{u}}’_i\)との関係は線形ではなく、なにか単調な非線形関係になるような気がする。その関数形がわからないと、怖くて\(\hat{\mathbf{u}}_i\)を使えない…
• 反応が3値以上のとき、閾値を固定したら識別条件はガラッと変わらないだろうか？ たとえば3件法なら閾値を(-1, +1), 4件法なら(-1, 0, +1), 5件法なら(-c, -1, +1, +c)と固定するだけで、もはや\(y_{ijk}\)の共分散行列が推定できるので、\(\mu\)には識別制約が要るが、\(\sigma^2, \Sigma\)はそのまま識別できるのでは？ 甘い？ この発想は甘いかなあ？
• 識別制約の代わりに弱情報事前分布を与えてベイズ推定することはできる？ 仮にできるとして、どんな弱情報事前分布ならいいの？ 特に、\(\Sigma\)への事前分布の与え方がわからない。Mplusだとデフォルトで\(\Sigma \sim IW(\mathbf{I}, r+1)\)となり、これは相関行列の非対角成分の周辺分布が\((-1, +1)\)上で一様になるような分布だから、もうちょっと制約しないといけないと思うんだけど、その一方で、この分布の対角成分(つまり個人効用の個人間分散)の周辺分布は\(IG(1, 1/4)\)であり、無情報事前分布のような顔をして実はそれなりの情報を持っているんだから、これでも収束するはずじゃないの？ という気もする。