« 読了：Allenby & Rossi (2006) 階層ベイズモデルとはなにか (feat.「SPL 狼よ静かに死ね」) | メイン | 読了：Bass(1969) 諸君、俺の名をとってこれをBassモデルと呼び未来永劫語り継ぐがよい (とは書いてないけど) »
2016年4月18日 (月)
Vanhonacker, W.R., Lehmann, D.R., Sultan, F. (1990) Combining Related and Sparse Data in Linear Regression Model. Journal of Business and Economic Statistics. 8(3), 327-335.
Sultan, et al.(1990)の後続研究としてSultanさんが紹介していた奴。実はこっちのほうがわかりやすいんじゃないかと疑って読んでみた。その予想はあたっていて、細かいところまできちんと説明してくれている論文であったが(というか、Sultan, et al.(1990)が説明を端折りすぎなのである)、しっかし、ねえ...
イントロ部分に早速わからない箇所があるので訳出しておく。
事前情報を使うことによって得られる効率性は、ベイズ統計手法が幅広く用いられ発展してきた理由でもある。Lindley & Smith (1972)は階層事前分布の使用を奨めており、多くの適用例の基礎となっている。独立に生成された知見を事前分布に統合しようという初期の試みはこのアプローチを採用している。しかし、階層事前分布の根底には交換可能性という非常に強い想定があり、これは過度な単純化だとして批判されてきた(Kass, 1983 JASA)。この想定は、通常は正規事前分布として捉えられているのだが(Leamer, 1978書籍)、知見はお互いに関連しており、データはその関連性についての情報を含んでいる、ということを述べている。この関連性には、研究における差異のせいで、ほぼ確実に異質性が含まれている。異質性についての事前情報はスパースであることが多く、グルーピングや条件付き事前分布の使用は現実的でないので、事前の研究のすべてに等しい情報が含まれているという厳密な想定を伴う伝統的アプローチが、妥当な批判にも関わらず生き延びている。
体系的変動をあきらかにするという試みにおいて、メタ分析の結果は独立に実行された研究群を通じたパラメータ異質性を明示的に捉える。従って、その結果を事前情報として用いれば、階層ベイズ事前分布の根底にある、厳密な交換可能性という想定を回避し、こうした事前分布に頼ることに対して実証研究の中でしばしば表明されてきた懸念に応えることができる。新しい研究を、メタ分析モデルで指定されたモデレータ変数について先行研究とマッチングさせることによって、事前分布の情報内容を差異化できる。いいかえれば、異質性が事前情報において明示的に捉えられていることになる。
ううむ。ひらたくいえばLenk&Raoのような階層ベイズモデルを使うアプローチをディスっているのだと思うのだが、階層事前分布における交換可能性の想定はきつすぎる、という批判の意味がよくわからない。もしかすると、階層事前分布には共変量を入れられないけどメタ分析モデルなら入れられるよ、という意味なのだろうか。だとしたら、Lenk&Raoが提案しているように、共変量をいれた階層ベイズモデルならいいんでしょ、という話になるはずで、交換可能性という概念を持ち出す必要はないだろう...
まあいいや。以下、本文についてのメモ。
いま、オリジナルのm個の研究で重回帰が使われているとする。
$y_i = X_i \beta_i + u_i$
$i$は研究番号。$y_i, X_i$のデータサイズを$n$、$X_i$の変数数を$p$とする[$n$と$p$に添字$i$がついてないけど、ほんとはつけるべきだろう]。すべての$i$において$E(u_i) = 0$, 異なる2つの研究$i,j$の間で $E(u_i u_j')=\sigma^2_i I$とする。
ここから「メタ分析事前分布」を導出する。
いま、研究の属性(製品カテゴリとか推定手法とか)によって、研究を$k$個のグループに分類できるとする。研究$j$がグループ$l$に属しているとき、ランダム係数を想定して
$\beta_j = \bar{\beta}_l + e_j$
ただし$E(e_j) = 0$、グループ$l$に属する異なる2つの研究$i,j$の間で$E(e_i e_j') = \Sigma_l$とする。
階層ベイズ事前分布は、すべての$i$について、パラメータ・ベクトル$\beta_i$がある平均ベクトルと共分散行列を持つある事前分布からのドローだと想定する。従って$\beta_i$は交換可能であり、この事前分布はすべての研究のあいだの「つながりの程度」を捉えていることになる。[中略]この強い想定について検証する方法もある。この想定によって、ある研究グループを分離して分析できるようになる。[中略。いっぽう]メタ分析によって異質性を明示的に捉えて記述すれば、すべての研究グループを同時に分析できるようになる。こちらのほうが効率が良いし、理解も促進される。[やっぱりそうだ。階層ベイズモデルには共変量がいれられないというのが話の前提になっているみたいだ]
まとめて表現しよう。$m$個の研究について、パラメータ、平均、誤差をそれぞれ縦積みする。
$\beta' = (\beta'_1, \ldots, \beta'_m)$
$\bar{\beta}' = (\bar{\beta}'_1, \ldots, \bar{\beta}'_k)$
$e' = (e'_1, \ldots, e'_m)$
とする。2本目の式に注意。たとえば3つ研究があって、所属グループが$(1,1,2)$なら、$\bar{\beta}' = (\bar{\beta}'_1, \bar{\beta}'_1, \bar{\beta}'_2)$である。
以上を使って
$\beta = \bar{\beta} + e$
さて、$\bar\beta$を各研究のデザインを列とする$Z$の線形和とみて
$\beta = Z c + e$
とする。これがメタ分析モデルだ。
実際には真のパラメータベクトル$\beta$は観察されないわけで、そのかわりに
$\hat\beta = Z c + \epsilon$
を推定する。 誤差項$\epsilon$について考えよう。これは真のパラメータに含まれる誤差項$e$と、パラメータの推定誤差$\hat{u} = \hat\beta - \beta$の和である。$E(\epsilon)=0$となる。では$E(\epsilon \epsilon’)$はどうなるか。
[さあ、この辺から話がややこしくなってくるんだけど]
$E(\epsilon \epsilon’)$を$\Omega$としよう。それは$m \times m$のブロック対角行列で、グループ$l$の対角ブロックが$\Sigma_l + \sigma^2_i (X'_i X_i)^{-1}$となる。[なぜか理解していないが、信じることにしよう]
で、ここから、$\hat{c}$のBLUEな推定量は下式となるのだそうだ[これも信じることにしよう]。
$\hat{c} = D\hat\beta$
$D = (Z' \Omega^{-1} Z)^{-1} Z' \Omega^{-1}$
さて、新しい研究がやってまいりました。そこでは
$y_o = X_o \beta_o + u_o, \ \ E(u_o)=0, \ \ E(u_o u'_o)=\sigma^2_o I$
が得られているとしましょう。この研究のデザインは$Z_o$で、これをメタ分析モデルに放り込み、 $\beta_o$の事前推定値
$b_o = Z_o \hat{c}$
が得られたとしましょう。
この$b_o$は、こういう性質を持っている、ということが示せるのだそうだ。
$b_o - Z_o D d = Z_o D M \beta_o + Z_o D \gamma$
$D$はさっき定義した行列。$d$は、すべての研究について、それが属するグループの平均パラメータのベクトルと、新研究が属するグループ平均パラメータのベクトルの差を縦積みしたもの。$\gamma$は$\hat{u}$と、「それが属するグループの$e$と新研究が属するグループの$e$の差」を足しあげたもので... ま、これは実際には使わないので、あんまし気にしなくてよさそうだ。
ってことはですね。
$y_o = X_o \beta_o + u_o$
$b_o - Z_o D d = Z_o D M \beta_o + Z_o D \gamma$
この2つを縦に積み上げて、一気に推定すればいいや、ということになる。縦に積んだのをこうやって表現することにします。
$y = X \beta_o + v$
ただし、$y$は$y_o$と$b_o - Z_o D d$を縦に積んだもの。$X$は$X_o$と$Z_o D M$を縦に積んだもの[本文中には $X' = \left[X'_o (Z_o D M)'\right]$とあるが、$X' = \left[X'_o, (Z_o D M)'\right]$の誤植であろう]。$v$は$u_o$と$Z_o D \gamma$なる謎のベクトルを縦に積んだもの。
さあ、この$\beta_o$の推定量はどうなるか。こうなることが示せるのだそうだ。
まず、誤差項$v$の共分散行列$\Omega_v$について考えよう。それはブロック対角行列となり、左上には$E(u_o u'_o)=\sigma^2_o I$が、右下には$Z_o D E(\gamma \gamma') D' Z_o$がくる。
では、$E(\gamma \gamma')$はどうなるのか。それは次の2つの行列の和になるのだそうだ。ひとつめ、$\Omega$。ふたつめ、新しい研究がグループ1に属するとして、すべての要素が$\Sigma_1$である行列。
さて、$\beta_o$のBLUEな推定量は
$\hat\beta = (X'\Omega_v^{-1}X)^{-1}X'\Omega_v^{-1}y$
これをrecursive estimatorという。
実際の推定は次のように行う。[一応メモするけど、要するにtake it easyという感じ。原文から少し添え字を変えた]。
グループ$l$の研究数を$m_l$とする。新研究はグループ1に属するということにする。
- $\hat{\sigma}^2_i = \frac{(y_i - X_i \hat{\beta}_i)'(y_i - X_i \hat{\beta}_i)}{n-1}$
- $\hat{\sigma}^2_o = \frac{1}{m_1} \sum_{i=1}{m_1} \hat{\sigma}^2_i$
- $\hat{\bar{\beta}}_l = \frac{1}{m_l} \sum_{j=1}{m_l} \hat{\beta}_j$
- $\hat{\Sigma}_l = \frac{1}{m_l} \sum_{j=1}{m_l} (\hat{\beta}_j - \hat{\bar{\beta}}_l)(\hat{\beta}_j - \hat{\bar{\beta}}_l)'$
- $d$は単純に、$\hat{\bar{\beta}}_l$同士を引き算して求める。
事例が2つ。1つめのマーケティング・ミクス・モデルはとばして、二つ目の新製品普及のほうのみ読んだ。
ご存知Bassモデルは、時点$t$において未受容者が受容する条件付き確率を$P(t)$、既受容者を$Y(t)$として
$P(t) = p + (q/m) Y(t)$
離散データの場合はラグをとって
$S_t = a + b Y_{t-1} + c Y^2_{t-1}$
とすればOLS推定できる。ほんとはNLSのほうがいいんだけど、話の都合上OLSで推定する。
過去の19個の普及曲線について$a, b, c$を推定した結果を集めてきた。これを消費財と産業財の2群にわける。
ほんというと、過去の普及曲線について$a, b, c$の共分散がほしいんだけど、手に入らないのであきらめて、共分散はみんな0だということにした[$\omega$を単なる対角行列にし、対角要素は当該パラメータの分散にした、ということであろうか]。
誤差分散の推定値$\hat\sigma^2_i$も手に入らないので、$\hat\sigma^2_o$のかわりに一年目の売上をいれた[いいのか、そんなことで。というか1年間の売上ってなんのことだ、予測対象であるカラーテレビのことか]
新研究としてカラーテレビの普及曲線を使い、1期目のみ、2期目まで、の2つについて事後パラメータを出す。2期でかなり予測できるようになった。3期だともうぴったり予測できちゃう。云々。
いやあ... メモしてみたはいいけれど... これ、死ぬほどややこしいわ... そのわりに、線形回帰に限定されている、という... Bassモデルの推定だったら、やっぱし階層ベイズモデルをつかったほうがいいんじゃなかろうか。
話を逆向きにたどり、推定手順をメモしておこう。
- それぞれの研究について、誤差分散の推定値を求める。すなわち$\hat\sigma^2_i = \frac{(y_i - X_i \hat\beta_i)'(y_i - X_i \hat\beta_i)}{n-1}$。
- 新研究について、誤差分散の推定値を求める。それがグループ1に属するとして、$\hat\sigma^2_o = \frac{1}{m_1} \sum_{i=1}{m_1} \hat\sigma^2_i$ 。
- それぞれのグループについて、パラメータ平均の推定値を求める。すなわち$\hat{\bar{\beta}}_l = \frac{1}{m_l} \sum_{j=1}{m_l} \hat{\beta}_j$。
- それぞれのグループについて、パラメータの分散を求める。すなわち$\hat{\Sigma}_l = \frac{1}{m_l} \sum_{j=1}{m_l} (\hat{\beta}_j - \hat{\bar{\beta}}_l)(\hat{\beta}_j - \hat{\bar{\beta}}_l)'$。
- すべての研究について、「そのグループの平均パラメータ」と「グループ1の平均パラメータ」の差を求める。これが$d$。
- 過去の研究のパラメータ推定値$\hat\beta$について、その誤差項$\epsilon$の分散行列$\Omega$を求める。それは$m \times m$のブロック対角行列で、グループ$l$の対角ブロックが$\Sigma_l + \sigma^2_i (X'_i X_i)^{-1}$。
- 話の都合上登場した謎の変数$\gamma$の共分散行列$E(\gamma \gamma')$を求める。それは、新研究がグループ1に属するなら、$\Omega$の全要素に$\Sigma_1$を足したもの。
- $\hat\beta$から事前パラメータを求める際に使う係数行列$D=(Z' \Omega^{-1} Z)^{-1} Z' \Omega^{-1}$を求める。
- 事前パラメータ$b_o=Z_o D\hat\beta$。
- 縦積みモデルの従属変数ベクトル$y$をつくる。それは新研究の従属変数$y_o$と、$b_o - Z_o D d$を縦積みしたベクトル。
- 縦積みモデルの誤差項の共分散行列$\Omega_v$を求める。それはブロック対角行列で、左上に$\sigma^2_o I$が、右下に$Z_o D E(\gamma \gamma') D' Z_o$がくる。
- 事後パラメータの推定量は、$(X'\Omega_v^{-1}X)^{-1}X'\Omega_v^{-1}y$。
つ、疲れる...
論文：データ解析(2015-) - 読了：Vanhonacker, Lehmann, Sultan (1990) 重回帰のメタ分析で事前分布をつくってベイズ推定