仕事の都合でGelmanらのMRP(マルチレベル回帰・事後層別)についてちょっとだけお話しする機会があって、資料を準備していてふと気がついたんだけど、Gelmanさんが前に書いた”Struggles…”という面白い論文があった。あれはMRPの萌芽だったんじゃないだろうか。あの論文ではそうは呼ばれていなかったけど。

Gelman, A. (2007) Struggles with survey weighting and regression modeling. (with commentaries.) Statistical Science. 22(2), 153-164.

　私は2014年にメモをとりながら読み終えている。メモを読み返してみると、手法提案論文というよりはもっと大きな枠組みを提示した論文なのだが、確かに、MRPの基本アイデア(非確率標本に対する階層回帰と事後層別)が提示されている。

　というわけで、今回ざっと再読したのだが、当時と現在とでは関心の持ち方も少し変わっているし、昔のブログ記事を書き換えるのもなんなので、もう一度メモをとってみた。

1. 背景
　調査ウェイティング、それはゴミ屋敷だ(“Survey weighting is a mess”)。単純な平均や割合の推定を別にすれば、ウェイトをどう使えばいいのかはっきりしないことが多い。
　魅力的な代替案として、ウェイティングのかわりに回帰モデリングをするという手があるが、膨大な交互作用をどう扱うかによって結果はどうにでも変わってしまう。
　本論文は、標本と母集団の間の差を調整するひとつの戦略として、事後層別と階層回帰の併用を提案する。

1.1 標本から母集団の量を推定する
　[数値例。New York City Social Indicators Surveyという縦断調査には、電話回線数、家族構成、エスニシティ、年齢、教育の分布が母集団に合うようにしたウェイトをがついている。同一の質問に対する2時点の反応率について、(a)各時点のweightedの平均の差(そのロジットスケールでの表現), (b)2時点を縦に積み、ウェイトに使った共変量と時点を説明変数として線形回帰(ロジスティック回帰)したときの時点の係数、を比べている。たとえば「ニューヨークの成人は健康な状態にあると思う」という回答の割合でいうと、(a)では、1999年の割合は75%, 2001年の割合は78%、差は3.4%。(b)の線形回帰では時点の係数は6.6%。ちょっと違うわけです。
　8年前に読んだときは「うーんなんで違うんだろう」と思ったけれど、今読んでみると当然だという気がする。(a)は事後層別セル別の時点差を平均しているが、(b)は事後層別セル別の時点差が同じだと仮定しているわけだから]

1.2 事後層別とウェイティング
　通常、モデルに基づく推測というものは、データ収集時のデザインが「無視可能」だと暗黙のうちに想定している(BDA2の7章をみよ)。つまり、回帰の文脈でいえば、標本抽出なり無回答なりに影響するすべての変数が含まれているという想定であり、標本抽出の文脈でいえば、事後層別セルのなかでは抽出確率が等しいという想定である。

　ここで、ウェイティングと事後層別を統一的に扱う枠組みを導入しておこう。
　\(X\)は標本において観察され、母集団について同時分布が既知だとする。母集団は十分に大きいとする。

• 事後層別。\(X\)が離散的だとして、その可能なカテゴリを事後層別セルとよび、\(j\)番目のセルの母集団サイズを\(N_j\), 標本サイズを\(n_j\)とする。どの事後層別セルでもデータは単純無作為抽出 (SRS) だと考える。標本サイズの割り当て方はこの話とは無関係である (古典的な層別抽出も事後層別の一種と考えるわけだ)。\(N_j\)は既知だとしよう (未知な場合も多いが、その推定の話は脇に置いておく)。
  　任意の変数の母平均 $$ \theta = \frac{\sum_j N_j \theta_j}{ \sum_j N_j} $$ の推定量は$$ \hat{\theta}^{PS} = \frac{\sum_j N_j \hat{\theta}_j}{ \sum_j N_j}$$ である。
• ウェイティング。個体ウェイトを\(w_i\)として、weightedの平均は$$ \bar{y} = \frac{\sum_i^n w_i y_i}{\sum_i^n w_i} $$ である。事後層別セルの中で個体ウェイトは個体を通して等しい。セルウェイトを\(W_j = n_j w_i\)として $$ \bar{y} = \frac{\sum_j W_j \bar{y}_j}{\sum_j W_j} $$ とも書ける。一般に、ウェイトはデザインと収集したデータの両方に依存する。
• 抽出確率に基づくウェイティング。事後層別ではない状況で抽出確率の逆数をウェイトとすることがある。たとえば、電話調査で世帯当たり電話回線数の逆数をウェイトにするような場合がそれだ(ウェイト値はデータと無関係に決まる)。でも、こういう固定ウェイトを使っていると、世帯当たり電話回線数と無回答の間に関係があったときにバイアスが生じる。本論文ではそういうのもみんな事後層別に組み込んだ場合を考える(つまり、固定ウェイトを使わず、世帯当たり電話回線数で事後層別する場合について考える)。

1.3 さまざまな推定方法
　母平均の推定のためにウェイトつき平均を用いるのは標準的だが、回帰のような複雑な分析の場合にどうすべきかは明確でない。たとえば\(y\)の\(z\)への回帰で、重みつき最小二乗法を使うという方法もあるし、ウェイティングに使った\(X\)をコントロールして重みなしの回帰をするという方法もある。方法によってSEの推定値も変わってくる。[数値例が示されている]

1.4 交互作用が果たす重大な役割
　[ここ、8年前は文意がとりにくく戸惑ったところだ。細かくメモする]
　\(y\)の\(z\)への回帰を考えよう。標本包含確率が\(X\)に依存しているとする。\(y\)は\(z\)と\(X\)の両方に依存しているかもしれない。そのときは\(y\)を\(X\)と\(z\)に回帰して、\(X\)の母集団分布を通じて平均するのがよい。このとき、\(z\)と\(X\)の交互作用もモデルに含めるのが大事だ。
　モデルが交互作用を含むとき、事後層別が必要になる。たとえば母集団における、収入(対数) \(y\) の身長 \(z\) への回帰に関心があるとしよう。データは男女比が母集団に一致するように調整されている。交互作用を含む回帰式を推定した。$$ y = 8.4 + 0.017 z – 0.079 \mathrm{male} + 0.007 z \cdot \mathrm{male} + \mathrm{error} $$ 所与の\(z\)について、$$ E(y|z) = 8.4 + 0.017z – 0.079 E(\mathrm{male} | z) + 0.007 z \cdot E(\mathrm{male} | z) $$ となる。\( E(\mathrm{male} | z) \)はこの調査自体から推定できるだろう。
　\(E(y|z)\)は\(z\)の線形関数でない。つまり、私たちは収入の身長に対する回帰を定義できるけれど、それがなぜ関心の対象となるのかがはっきりしない。[原文: it is not clear why it should be of any interest.]
　調査において調整があるとき、回帰係数の解釈はこのように困難になる。これが、1.1節での母平均の単純比較において注意が必要な理由である。例として、白人と非白人の間の年収(対数)の差の平均を推定するという問題を考えよう。交互作用を含めた回帰式は$$ y = 9.5 – 0.02 \mathrm{white} + 0.20 \mathrm{male} + 0.41 \mathrm{white} \cdot \mathrm{male} + \mathrm{error} $$ と推定された。従って $$ E(y | \mathrm{white}) – E(y | \mathrm
{not white}) $$ $$ = -0.02 + 0.20 \{E(\mathrm{male} | \mathrm{white}) – E(\mathrm{male} | \mathrm{not white})\} + 0.41 E(\mathrm{male} | \mathrm{white}) $$ となる。\( E(\mathrm{male} | \mathrm{white}), E(\mathrm{male} | \mathrm{not white})\)はこの調査自体から推定できるだろう。

　このように、調査に調整があるとき、私たちは交互作用モデルをあてはめるので、単純な回帰係数(ここでは\(\mathrm{white}\)の回帰係数)だけみているわけにはいかなくなり、交互作用項を事後層別セルを通して平均しなければならなくなる。
　本論文の焦点は、調査における反応のモデルと、それに対応するweightedの平均推定値との関係である。本論文の究極の目標は、ウェイト構築のためのモデル・ベースの手続きを得ることである。別の言い方をすれば、調査に調整があるときに、回帰モデルから有効かつほぼ不偏な推定値を得るための枠組みをつくることである。
　[ああ、そういう話だったのか… 8年前に読んだときは、回帰モデルに交互作用項がはいると係数の解釈が難しくなるよというごく一般的な指摘なのか？ と戸惑ったのだが、話の重点は、データに\(X\)による調整が入っているからには\(X\)は\(y\)と\(z\)の共変量なんだから、当然\(X\)と\(z\)の交互作用項もモデルにいれないといけないよね、だから母集団における\(z\)の係数を推定するために交互作用項を事後層別セルを通して平均しなければならなくなるよね、という点なのだ]

2. 問題
2.1 単純な平均とトレンドの推定
　1.1節の事例に戻ろう。(a)(b)のどちらを信じるべきか。以下では正しい答えを与えてくれて、かつより複雑な推定対象にもスムーズに一般化できるアプローチについて考えよう。

2.2 deep事後層別
　ウェイティングでは事後層別のセル数がすぐに膨大になってしまう。交互作用のモデリングについて、\(X\)と\(z\)の交互作用のうちこれは含めるがこれは含めない、という選択が必要になる。

3. ウェイティングと事後層別を回帰モデリングで結合する
　セル平均がある種の線形回帰で推定されたとき、事後層別推定値[\(\hat{\theta}^{PS}\)のこと。全体の母平均の推定値]は重み付き平均として解釈できる。

3.1 古典的モデル

• 完全な事後層別。セル推定値\(\hat{\theta}_j\)としてセル平均\(\bar{y}_j\)を使って $$ \hat{\theta}^{PS} = \frac{\sum_j N_j \hat{y}_j}{\sum_j N_j}$$ 事後層別セルのインジケータを含む古典的回帰とみることもできる。
• ウェイティングしない。すべての\(i\)について\(w_i = 1\)とする。セル推定値\(\hat{\theta}_j\)として\(\bar{y}\)を使うといってもいいし、定数項のみの古典的回帰といってもいい。
• セル特性に対する古典的回帰。上記の2つの中間に位置する。事後層別に用いた変数を回帰に投入する(交互作用は必ずしも投入しない)。
  　層別変数のデータ行列を\(X\)として、回帰モデルは $$y \sim N(X \beta, \sigma^2_y I) $$ $$\hat\beta = (X^\top X)^{-1} X^\top y $$ 事後層別セルの母集団サイズのベクトルを\(N^{POP}\), 層別変数の行列を \(X^{POP}\)とする。セル平均の推定値は \(X^{POP} \hat\beta\) だ。では母平均の推定値はどうなるか。それはセル平均の推定値の加重平均であるから $$ \hat\theta^{PS} = \frac{1}{N} \sum N_j (X^{POP} \hat\beta) $$ \(\hat\beta\)を代入して $$\hat\theta^{PS} = \frac{1}{N} (N^{POP})^\top X^{POP} (X^\top X)^{-1} X^\top y$$ これを $$\hat\theta^{POP} = \frac{1}{n} \sum w_i y_i$$ と書きなおそう。\(w\) は $$w = \left( \frac{n}{N} (N^{POP})^\top X^{POP} (X^\top X)^{-1} X^\top \right)^\top $$ \(w\)の合計は1になる。つまり、これもまたウェイティングだと捉えることができる。なお、\(w\)はデータとモデルに依存しているが\(y\)には依存していない点に注意。

3.2 階層モデル
　今度は、セル平均\(\hat{\theta}_j\)を階層モデルで推定する場合を考えよう。ここからは\(w_i\)が\(y\)に依存する。

• 階層回帰。モデルを$$ y \sim N(X \beta, \Sigma_y) $$ $$ \beta \sim N(0, \Sigma_\beta)$$ とする。事前の精度行列\(\Sigma_\beta^{-1}\)は対角行列で、階層回帰の係数でないところは0になる。たとえば、予測子が定数項, 性別, エスニシティ(黒人/非黒人), 性別xエスニシティ, 年代(4), 学歴(4)、年代x学歴(16)だとしよう。古典的回帰なら、予測子は1(定数項)+1(性別)+1(エスニシティ)+1(性別xエスニシティ)+3(年代)+3(学歴)+9(年代x学歴)=19個。階層回帰の場合は、1+1+1+1+4+4+16=28個と考え、事前精度行列の対角成分は、最初の4個を0, 次の4個を同一, 次の4個を同一, 次の16個を同一とする。[←いまいち腑に落ちないんだけど…まあいいや、先に進もう] $$ \hat\beta = ( X^\top \Sigma^{-1}_y X + \sigma^{-1}_\beta)^{-1} X^\top \Sigma^{-1}_y y$$ $$ \hat{\theta}^{PS} = (1/N) (N^{POP})^\top X^{POP} \times (X^\top \Sigma^{-1}_y X + \Sigma^{-1}_\beta)^{-1} X^\top \Sigma^{-1}_y y$$となる[wの式は省略]。\(\Sigma_y, \Sigma_\beta\)のもとで重み付き平均となる。
• 交換可能な正規モデル。階層回帰モデルについての理解を深めるため、その特殊ケースである、セル平均についての交換可能な正規モデルを考えよう。$$ \bar{y}_j \sim N(\theta_j, \sigma^2_y/n_j) $$ $$ \theta_j \sim N(\mu, \sigma^2_\theta)$$ [原文では一本目の式の分散は\(\sigma^2/n_j\)だが、添字\(y\)を勝手につけた。] 事後層別推定値は\(\sigma_y, \sigma_\theta\)のもとで重み付き平均となる。
  　セル平均の事後平均からはじめよう。ここから一時的にセル番号の添字を\(k\)にする。$$ \hat\mu = \frac{\sum_k \frac{\bar{y}_k}{\sigma^2_y / n_k + \sigma^2_\theta}}{\sum_k \frac{1}{\sigma^2_y / n_k + \sigma^2_\theta}} $$として $$ \hat\theta_k = \frac{(n_k/\sigma^2_y) \bar{y}_k + (1/\sigma^2_\theta) \hat{\mu}}{n_k / \sigma^2_y + 1/\sigma^2_\theta}$$ である。ここから$$ \hat{\theta}_k = \sum_j c_{kj} \bar{y}_j$$ と書ける。\(c_{kj}\)は[…省略…]。$$ \hat{\theta}^{PS} = \sum_k \sum_j \frac{N_k}{N} c_{kj} \bar{y}_j$$ となる。つまり\(\sum_k \frac{N_k}{N} c_{kj}\)を重みにした重み付き平均になっている。個体ウェイトに直すと[…中略…]。このように、階層事後層別は、上の\(\hat{\theta}_k\)の式と同じ因子でウェイトを縮小していることになる。
  　このように、一般に階層回帰モデルでは、ウェイトの縮小の程度が、アウトカム\(y\)の層内分散と層間分散に依存して決まる。[なんだかわからんが、へー、そうなんですか]
• その他の階層モデル。上の2つのモデル(つまり、セルレベルの分散成分を持つ階層回帰)に、予測変数が順序尺度だったときの隣接セル間の相関を組み込んだモデルが提案されている。Lazzeroni & Little(1998 J.OfficialStat.), Elliott & Little(2000 J.OfficialStat.)をみよ。またロジスティック回帰への拡張もあり得る。このときは\(\hat{\theta}^{PS}\)はもはや重み付き平均ではなくなってしまう。

3.3 モデル・ベースの事後層別推定値の諸特性

• 標準誤差。[略]
• 暗黙的ウェイトのアウトカムへの依存性。古典的なウェイトは、\(n_j, N_j, X\)には依存しても\(y\)には依存しない。ただし、\(y\)によってはわざわざウェイティングする必要がない、という意味での依存性はあるけれど。しかし、階層回帰で暗黙的に用いられている個体ウェイトは、ハイパーパラメータ\(\Sigma_y, \Sigma_\beta\)に依存していて、ハイパーパラメータはデータから推定する。従って、\(y\)が変わればウェイトも変わる。

4. 今後の展開
まとめると、標本と母集団の既知の差を調整するふたつの標準的な方法として、ウェイティングと回帰モデリングがあるわけだ。

• ウェイティングの実務的な限界。(1)回帰係数のような複雑なestimandsに対してどうウェイトを使えばいいのかわからない。(2)weightedの推定値は標準誤差の推定が難しい。(3)ウェイトの構築はあなたが思うより難しい。主観的な判断がすごく入ってくる。
• モデリングの実務的な限界。理論的には、すべての事後層別セルで条件付けないといけないわけだが、そうするとモデルが複雑になりすぎる。とはいえ、すごく複雑な回帰モデルによって母平均をうまく推定できることもあるわけで、モデルのパラメータ数だけでなく、関心ある量とモデルとの関係も問題となるわけだ。
• 階層モデル・事後層別による統合。出発点は、表1のような例でどちらの推定値が良いかを理解することだ。[とかなんとか…中略…] 本論文で提案したような、信頼がおけてかつ簡単な統合的アプローチの開発が求められている。云々。

————————-
　やれやれ、疲れた…
　勉強にはなるんですけど、Gelman兄貴の論文って、ちょっとわかりにくいところがあるように思う。天才肌というか、ちょっと説明をすっ飛ばすところがあるというか… いやまあ、私のレベルが低いからでしょうけど。

　この論文には5人の識者によるコメントと著者の返答がついている。8年前は疲れちゃって目を通さなかったのだが、今回はちゃんと読んでみた(むしろそっちが目的であったともいえる)。長くなるので、そちらのメモは別のエントリで。