« 読了:DuMouchel & Duncan (1983) 回帰分析も「ウェイト・バック」すべきでしょうか | メイン | 映画「カンフー・ジャングル」に登場する香港映画人たち »
2016年5月22日 (日)
先週都合により目を通した論文のメモをのせておく。ほんとはそれどころじゃないんだけど...
Rossi, P., Gilura, Z., Allenby, G.M. (2001) Overcoming Scale Usage Heterogeneity: A Bayesian Hierarchical Approach. Journal of the American Statistical Association, 96(453), 20-31.
調査票に5件法の項目をずらっと並べたところ、回答をみると妙に右側にばかりマルしている人もいれば妙に左側にばかりマルしている人もおり、これはいかなる前世の因縁かとくよくよ思い悩む世界1000万人のリサーチャーを憐れみ、ベイジアン・モデリングのビッグ・スターAllenby兄貴とその仲間たちが授ける必殺の階層ベイズモデル。(必殺でもないかも。憐れんでもいないかも。1000万人もいないかも)
仕事の都合で急遽目を通した。兄貴の主著"Bayesian Statistics and Marketing"(以下BSM本)のCase Study 3とほぼ同じ内容だと思うのだが、同じのを二度読むのもなんだか癪にさわったし、本を持ち歩くのは重かったので。
前置きは省略して、いきなりモデルにはいると。。。
対象者 $i$ ($N$人) の項目 $j$ ($M$個) への回答の背後に連続的潜在変数$y_{ij}$を仮定する。全部$K$件法だとして、閾値$c_0, \ldots, c_K$を仮定し、$c_0 = -\infty, c_K=\infty$とする。[ここでは閾値が全員共通だと考えている点に注意]
ベクトル$y_i = [y_{i1}, y_{i2}, \ldots, y_{iM}]'$について$y_i \sim N(\mu^*_i, \Sigma^*_i)$と仮定する。
[えーっと... ここまではまだ、潜在変数が回答に落ちる仕組みの話であって、潜在変数の生成についてはなんら語ってないわけね]
で、$y_i$についてこうモデル化する。
$y_i = \mu + \tau_i \iota + \sigma_i z_i, \ \ z_i \sim N(0, \Sigma)$
よって$\mu^*_i = \mu + \tau_i \iota, \Sigma^*_i = \sigma^2_i \Sigma$である。で、尺度の左端につけやすい人は$\tau_i$が大きく$\sigma$が小さい人、両極につけやすい人は$\tau_i$がゼロで$\sigma$が大きい人、という風に考える。
[注目!ここでいう$y_i$は態度のベクトルではなくて、態度と回答スタイルをコミにしたベクトルなのだ。$M$個の回答の背後にある態度ベクトルの分布としてまず多変量正規分布 $N(\mu, \Sigma)$を考え、$i$さんの潜在ベクトルの分布はそれを$\mu$から$\sigma_i$倍して$\tau_i \iota$ずらしたものだ、と考えるわけね]
あるいはこう考えてもよい。上では閾値を全員共通に$c_1, \ldots, c_K$と捉えたが、$y_i$の分布は全員で共通と考え、閾値の側に個人差があると捉えて
$c^*_i = \tau_i + \sigma_i c$
これでも結局は同じことになる。[そうやね、こっちのほうがわかりやすいわね]
モデルに階層をいれます。
標準的な階層モデルであれば、位置$\tau_i$と尺度は$\sigma_i$はアプリオリに独立で、条件付き共役事前分布に従うと仮定するところだが、ここでは次のように考える。
$[\tau_i, \log \sigma_i]' \sim N(\varphi, \Lambda)$
$\tau_i$と$\sigma_i$が相関を持ちうることに注意。
識別の都合上、$\varphi_1 =0 ,\varphi_2 = -\lambda_{22}$と制約します。$E[\tau_i]=0, E[\sigma_i]=1$としているわけです。
[2018/01/01追記: $\lambda_{22}$とは$\Lambda$の2番目の対角要素であろう。このくだり、この論文を読んだときには疑問に思わなかったのだが、このたび振り返ってみるとよく理解できない。BSM本の該当箇所(p.242)のほうが詳しいし、制約のかけ方が若干変更されているので、訳出しておく]
$\tau_i$を識別するため、制約$E[\tau_i]=0$をかける。$\tau_i$の分布は対称かつユニモーダルなので、この分布の中央値と最頻値もゼロにしたことになる。さて、$\sigma_i$の分布は歪んでいるから、識別制約の選択には十分な注意が必要である。ひとつの論理的な選択は$E[\sigma^2_i]=1$とすることであろう。これは識別制約として$\varphi_2 = -\lambda_{22}$をかけることに対応する。しかし、散らばりのパラメータ($\lambda_{22}$)が大きくなるにつれて、$\sigma_i$の分布は太い右裾を持ちつつ、1より左にある最頻値の周りに集中する。これは事前分布としてあまり望ましいものではない。[...]
もっと理にかなったアプローチは、$\sigma_i$の事前分布の最頻値を1に制約するというものである。これは識別制約として$\varphi_2 = \lambda_{22}$をかけることに対応する。[...] $\lambda_{22}$が大きくなるにつれて、$\sigma_i$の分布は1のまわりに固まるが、右裾が厚くなって散らばりが大きくなる。というわけで、我々は
$\varphi_1 = 0, \varphi_2 = \lambda_{22}$
という識別制約を採用する。なお、Rossi et al.(2001)では平均を1にするという方略を採用していた。
[うううむ。云っている意味はわかったけど、$E[\sigma_i]=1$とするためには$\varphi_2 = -\lambda_{22}$とすればよいという理由がやっぱりわからない。記号のせいでどうも混乱してしまうので、頭の悪い俺のために、$\tau_i$と$\sigma_i$を独立とみなし、$\sigma_i, \varphi_2, \lambda_{22}$をそれぞれ$X, A, B$と書き換えよう。すると
$\log X \sim N(A, B)$
$X$は対数正規分布に従っている。対数正規分布の平均は
$E[X] = \exp(A+B/2)$
だ。$E[X]=1$とするためには$A + B/2 = 0$、つまり$A = -B/2$とすればよい。もとの記号に戻すと$\varphi_2 = -\lambda_{22}/2$となる。あっれーーー? どこで間違えたんだろう?]
閾値のほうは、$\sum_k c_k = m_1, \sum_k c^2_k = m_2$と制約します。さらに、$c_k = a + b k + e k^2$という制約も入れ、自由パラメータを$e$だけにします。[細かいところよくわかんないけど、後者の制約はちょっと勘弁してほしい... 前者の制約だけで自由パラメータは$K-2$個になる模様だ、5件法とかならもうそれで十分じゃん]
最後に、事前分布の導入。
$\pi(\mu, \Sigma, \varphi, \Lambda, e) = \pi(\mu) \pi(\Sigma) \pi(\varphi) \pi(\Lambda) \pi(e)$
として...
2018/01/01追記: ここからは、この論文のメモから離れ、この論文、BSM本、Rのbayesm::rscaleUsage()のマニュアル、の3つを比較してみる。
- 潜在変数の平均ベクトル$\mu$の事前分布$\pi(\mu)$について。
- この論文: $\pi(\mu) \propto constant$
- BSM本: $\pi(\mu) \propto constant$
- bayesm::rscaleUsage(): $\mu \sim N(mubar, Am^{-1})$; mubarのデフォルトはrep(k/2, p) (kは段階数, pは項目数), Amのデフォルトは0.1*I。
- 潜在変数の分散行列$\Sigma$の事前分布$\pi(\Sigma)$について。
- この論文: $\Sigma^{-1} \sim W(\nu_\Sigma, V_\Sigma)$
- BSM本: $\Sigma \sim IW(\nu_\Sigma, V_\Sigma)$
- bayesm::rscaleUsage(): $\Sigma \sim IW(nu, V)$
- $\pi(\Sigma)$のハイパーパラメータについて。
- この論文: $\nu_\Sigma = dim(\Sigma) + 3 = K+3,\ V_\Sigma = \frac{1}{\nu_\Sigma} I$
- BSM本: $\nu_\Sigma = dim(\Sigma) + 3 = K+3,\ V_\Sigma = \nu_\Sigma I$
- bayesm::rscaleUsage(): nuのデフォルトはp+3(pは項目数)、Vのデフォルトはnu*I
仮にそうだとして... $M$次元多変量正規分布の共分散行列の事前分布を$IW(M+3, (M+3)I)$とするというのは、ここに限ったことではなくて、無情報事前分布としては一般的なやりかたらしい。BSM本の2.8.3節。 - 尺度使用の異質性を表すパラメータ$\tau_i, \sigma_i$の平均ベクトル$\varphi$は、この論文では$\varphi = [0 ,-\lambda_{22}]'$, BSM本とbayesmでは$\varphi = [0, \lambda_{22}]'$なので、その事前分布$\pi(\varphi)$については考えなくてよろしい。
- $\tau_i, \sigma_i$の共分散行列$\Lambda$の事前分布$\pi(\Lambda)$について。
- この論文: $\Lambda^{-1} \sim W(\nu_\Lambda, V_\Lambda)$
- BSM本: $\Lambda \sim IW(\nu_\Lambda, V_\Lambda)$
- bayesm::rscaleUsage(): $\Lambda \sim IW(Lambdanu, LambdaV)$
- $\pi(\Lambda)$のハイパーパラメータについて。
- この論文: $\nu_\Lambda = dim(\Lambda) + 3 = 2+3,\ V_\Lambda = \frac{1}{\nu_\Lambda} I$
- BSM本: (長くなるので後述)
- bayesm::rscaleUsage(): Lambdanuのデフォルトは20, LambdaVのデフォルトは(Lambdanu-3)*Lambda, Lambdaのデフォルトはmatrix(c(4,0,0,0.5),ncol=2)
- 閾値のパラメータ$e$の事前分布$\pi(e)$について。
- この論文: $\pi(e) \propto unif[-0.2, +0.2]$
- BSM本: $\pi(e) \propto unif[-0.2, +0.2]$
- bayesm::rscaleUsage(): e $\sim$ unif on a grid. なお、この点についてはBSM本のp.252に解説がある。いわく、Rの実装にあたっては、$e$の事前分布として、間隔(-0.1, 0.1)におけるグリッド上の一様分布を使い、このグリッドのうえでランダムウォークした。つまり、グリッドの点の右か左に等しい確率で移動した(ただし両端においては内側に移動した)。云々。
BSM本における$\Lambda$の事前分布についての説明。まわりくどくて面倒なので、虚心に丸ごと全訳する。
$\Lambda$の事前分布は、$\tau_i, \sigma_i$の推定値がシュリンクする程度に影響する。我々の階層モデルは、$\tau_i, \sigma_i$の分布についてデータが持っている情報に適応できるが、ただしそれはハイパーパラメータ$\Lambda$の事前分布の支配下での話である。$\tau_i, \sigma_i$の推定値の基盤となる設問は、あまり数が多くないのがふつうだろう。このことは、$\Lambda$の事前分布の影響がきわめて大きくなるであろうということを意味している。階層モデルの多くの適用事例では、十分な量の情報を提供してくれるユニットのサブセットが存在し、このサブセットのおかげで、適応的シュリンケージによって$\Lambda$を決定できる。しかし我々の状況では、こうしたサブセットを利用することができず、$\Lambda$の事前分布が通常よりも大きな影響力を持ちうるのである。以上の理由により、$\Lambda$の事前分布の選択には注意が必要である。また、ここで検討・提案する事前分布のセッティングは、階層モデルの文脈で通常用いられるものよりも幾分タイトなものになる。
$\Lambda$の事前分布の役割は、$\tau_i, \sigma_i$の事前分布を導出することである。事前分布のハイパーパラメータの選択が持つ含意について検討宇するために、$\tau_i, \sigma_i$の周辺事前分布をシミュレーションで計算してみよう。周辺事前分は下の式で定義される。
$\pi(\tau, \sigma) = \int p(\tau, \sigma | \Lambda) \pi(\Lambda | \nu_\Lambda, V_\Lambda) d\Lambda$
$\nu_\Lambda, V_\Lambda$について評価するために、$\tau_i, \sigma_i$がとりうる値については幅広く検討するが、事前分布の分散については許容可能な範囲より大きくならないように制約する。$\tau$については範囲を$\pm 5$とする。これは10件尺度の大部分を包含するという意味で非常に幅の広い範囲である。$\sigma$については、このパラメータによって潜在変数が取り得る値の範囲がどのように制約されるかといを考えないといけない。$\sigma$が小さいということは、回答者が尺度のうち狭い部分しか使わないということを表し、$\sigma$が大きいということは、回答者が尺度の全体を使うということを表す。いま、「小」尺度使用者(たとえば、下の3段階とか、上の3段階とかしか使わないような人)、「中」尺度使用者(10件法で5段階しか使わない人)、「大」尺度使用者(尺度の両端を使う人)がいるとしよう。「中」使用者のSDに対する「小」使用者のSDの比は、$\sigma$の「小さい」値に対応するだろう。また、「中」使用者のSDに対する「大」使用者のSDの比は、$\sigma$の「大きい」値に対応するだろう。このように考えると、$\sigma$の範囲は幅を広く取って(0.5, 2)程度だといえるだろう。そこで次のセッティングを採用する。これは$\tau_i, \sigma_i$についての比較的に情報的な事前分布に対応している。
$\nu_\Lambda = 20, \ V_\Lambda = (\nu_\Lambda - 2 - 1) \bar{\Lambda}$
$\bar{\Lambda}$は2x2の対角行列で、対角要素は[4, 0.5]とする。
このセッティングは、$E[\Lambda]$が$\bar{\Lambda}$にほぼ等しいことを保証する。$\tau_i, \sigma_i$の周辺事前分布は、変に大きな値を取ることなく、適切な範囲をカバーする。
やれやれ... 追記はここまでにして、論文のメモに戻る。
....でもって、これをGibbsサンプラーで推定する。[付録に細かいことがいろいろ書いてあるけど、パス]
実例。顧客満足度調査、n=1811, 10件法6項目。このモデルで推定した真値と支出額との相関が高かった、というような話。
... 読みなおしてみて、やっぱり面白い論文であった。Allenby兄貴のチームは、論文は面白いのにね、なんであの本はあんなにわかりにくいんだろうか。
感想が4点。
- こうやって$M$項目の回答データの背後に$M$次元の多変量正規分布を仮定し、尺度使用の異質性を取り除いた真の態度を推定したとして、じゃあどうすんの、と思う。きっと推定したスコアを回帰分析や因子分析や潜在クラス分析にかけることを想定しているのだと思うけれど、そのときにはスコアが多変量正規分布に従うという想定を外すことになるわけだ。ほんとは尺度使用の異質性だけを先にモデリングするんじゃなくて、そのあとで想定するモデルまで全部一発で同時推定したほうがいいんじゃないかろうか。この論文のような二段階アプローチってどこまで許されるんだろうか、その辺のところ、ご専門の方に訊いてみたいものだ。
- scale usage heterogeneityをモデル化する動機はなんだろうか。この論文の理屈では、ほっとくと項目間の共分散行列が底上げされちゃうことが問題なのだろうけれど(そうか、結局この論文はcommon method varianceの話なんだな)、実際の調査データ分析では、CMVそのものを問題にすることより、群間での回答スタイルの差(典型的には国によるちがい)を問題にすることのほうが多いんじゃないかしらん。えーと、そうすると、この論文の枠組みでいえば、上の階層で$\tau_i$と$\log \sigma_i$を群を表す共変量に回帰させればいいんだろうな。正直に言うと、群内の閾値の異質性なんてあんまし関心なくて、Mplusでいえば
MODEL: F1-F3 by u1-u20 (*1);
MODEL US: [u1$1-u20$1 u1$2-u20$2 u1$3-u20$3 u1$4-u20$4];
MODEL JP: [u1$1-u20$1 u1$2-u20$2 u1$3-u20$3 u1$4-u20$4];
...
というような感じのモデルを柔軟に組みたいんだけどなあ(上のコードが正しいかどうかは別にして)。ま、Stanとかで自分でやれってことなんでしょうね。 - いっぽう、尺度使用の記述(?)としてこのモデルを使うのはすごく面白いなあ、と思った次第である。そんなこんなで、k件法尺度の項目をある個数以上聴取したデータを扱うときには、この論文の手法を実装したRのbayesm::rscaleusage()を回してみたりしているんだけど、うーん、どのくらいの項目数があればいいのかしらね。直感としては、項目数が少ないときはうまくいかないだろうと思う。そのへんを教えてくれるとありがたいんだけど。
- ところで、尺度使用の異質性のせいでセグメンテーションにバイアスがかかるという実証研究としてGreenleaf (1992, JMR)というのがあるのだそうだ。回答スタイルの研究については以前原稿の都合でかなり頑張って論文を集めたのに、見落としがあった模様。やになっちゃうなー。[→2018/01/07追記: ちゃんと読んでいた。なんだかなあ]
せっかくなので(なにがだ)、Allenby兄貴の「兄貴」称号の由来となった、香港映画最後のカンフー・スター、ドニー・イェン兄貴の名場面から。日本での最新公開作「カンフー・ジャングル」 (2014, 陳徳森監督)の最後の対決シーンより。
ドニー兄貴の大ヒット作といえば「イップ・マン」シリーズだが、現代劇に限って言えば、本作が兄貴の最高傑作なのではないだろうか。個人的な意見でございますが、回り道をし歳を食い、カンフーこそすごいけれどもちょっとダメなところのある男を演じ始めたとき、ドニーさんは俄然光り輝き始めたように思う。「孫文の義士団」のギャンブル中毒とか。
この映画でいうと、敵役のワン・バオチャンがいかに非人間的なまでに強いかという描写に力を注いだところも勝因だけれど、主演のドニーさんが妹弟子の写真をひそかにスマホの待ち受け画像(?)にしているというあたりが、なんというかその、ぐっとくるわけであります。
(2018/01/01: 大幅に追記した)
論文:データ解析(2015-) - 読了:Rossi, Gilura, Allenby (2001) k 件法項目で高いほうにばかり答える人や低いほうにばかり答える人がいるのをなんとかする