Gambino, J. G. (2009) Design Effect Caveats. The Americal Statistician, 63(2), 141-146.

　調査データ解析に関してときどき出てくる話題である「デザイン効果」についての資料を探していたら、2009年のAmerican Statisticianに解説が載っているのを見つけた。ありがたや。考えてみれば2009年は最近とはいえないが、話題自体がわりかし古いので、この話題の解説としては新しめであるといえよう。

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

2. 不均一確率抽出のもとでの単純平均と分散の性質
　有限母集団\(\{y_1, y_2, \ldots, y_N\}\)があるとしますね。合計を\(Y = \sum_{i=1}^N\), 平均を\(\bar{Y} = Y/N\), 分散を\(S^2 = \frac{1}{N-1} \sum_{i=1}^N(y_i – \bar{Y})^2\) とする。
　サイズ\(n\)の確率標本を得たとしよう。標本包含確率と同時標本包含確率を\(\pi_i, \pi_{ij}\)とする。標本包含インジケータを\(\delta_i\)とする。$$ \bar{y} = \frac{1}{n} \sum_{i=1}^n y_i = \frac{1}{n} \sum_{i=1}^N \delta_i y_i$$ とする。ここから$$ E(\bar{y}) = \frac{1}{n} \sum_{i=1}^N \pi_i y_i$$ $$ E[\sum^n y^2_i] = \sum_{i=1}^N \pi_i y_i^2$$ $$ E[(\sum^n y_i)^2] = \sum_{i=1}^N \sum_{j=1}^N \pi_{ij} y_i y_j$$ なので、標本分散$$ s^2 = \frac{1}{n-1} \sum^n (y_i – \bar{y})^2 = \frac{1}{n-1} \left[ \sum^n y_i^2 – \frac{1}{n} (\sum^n y_i)^2 \right]$$の期待値は $$ E(s^2) = \frac{1}{n-1} \sum_{i=1}^N \pi_i y^2_i – \frac{1}{n(n-1)} \sum_{i=1}^N \sum_{j=1}^N \pi_{ij} y_i y_j$$となる。SRSの場合は\(\pi_i = \pi_{ii} = \frac{n}{N}, \pi_{ij} = \frac{n(n-1)}{N(N-1)} \) であるから、\(E(s^2) = S^2\)である。
　ここまではいいですね？ 文句ないですね？

3. デザイン効果の適用
　SRSの下では、母合計の推定量は\(N \bar{y}\)、その分散は\(V_{SRS} = N^2(1-\frac{n}{N}) \frac{S^2}{n}\)である。
　複雑なデザインの下では、ふつうHT推定量 \(\hat{Y} = \sum_{i=1}^n \frac{y_i}{\pi_i} \)が用いられるけれど、そこはここでは問わない。なんだかしらんが、その分散を\(V_c\)としよう。
　\(\hat{Y}\)のデザイン効果とは二つの分散の比、すなわち$$ d = \frac{ V_c(\sum_{i=1}^n \frac{y_i}{\pi_i}) }{ V_{SRS}(\frac{N}{n} \sum_{i=1}^n y_i) } = \frac{ V_c }{ N^2(1-\frac{n}{N}) \frac{S^2}{n} }$$である。HT推定量でなくてもこのかたちになる。
　逆に、\(V_c\)がデザイン効果を使って推定されることも多い。なかには、$$ \hat{V}_E = d N^2 (1-\frac{n}{N}) \frac{s^2}{n} $$ としている場合もある。しかしその期待値は$$ E(\hat{V}_E) = V_c E(s^2) / S^2$$ であって\(V_c\)ではない。一般には\(E(s^2) \neq S^2\)だからである。\(\hat{V}_E\)の式のなかの\(s^2\)を不偏推定量\(\hat{S}^2\)に差し替えないといけないのである。

4. SRS分散とデザイン効果
　デザイン効果の式には\(S^2\)が入っている、という問題について。
　関心あるのが母合計\(Y\)だとしよう。仮に標本がSRSだったとして、推定量$$ \hat{Y}_{SRS} = \frac{N}{n} \sum_{i=1}^n y_i = N\bar{y}$$ の分散は$$ V_{SRS} = N^2 \left( 1-\frac{n}{N} \right) \frac{S^2}{n}$$ である。このように、現実のデザインのもとでのデザイン効果を推定したいのに、やっぱり\(S^2\)の推定値が必要になってしまう。2節で示したように、標本からふつうの式で求めた\(s^2\)は(SRSでない限り)\(S^2\)の推定値ではない。
　では、\(S^2\)をどうやって推定するか？
　[ここから、SRSでない標本からどうにかして\(S^2\)の推定値を得る方法について説明がはじまるんだけど、そもそも\(S^2\)の推定値もソフトが出してくれたりしませんかね？ そんなことはないか？ うーん、調べてみないとわかんないな…]

　思うに、$$ (N-1) S^2 = \sum_{i=1}^N y_i^2 – \frac{Y^2}{N} $$ ですよね。\(u_i = y^2_i\)とし、第1項を\( \sum_{i=1}^N y^2_i = \sum_{i=1}^N u_i = U \)と書こう。\(Y\)なり\(U\)なりの推定値はおそらくソフトが出力してくれるだろう。残るは\(Y^2\)だ。いま\(Y\)の不偏推定量\(\hat{Y}\)があれば、それがなんであれ、\( V_c(\hat{Y}) = E(\hat{Y}^2) – Y^2 \)と書ける。つまり\(Y^2 = E(\hat{Y}^2) – V_c(\hat{Y})\)だ。おそらく、ソフトは\(V_c(\hat{Y})\)の(近似的な)不偏推定量\(v_c(\hat{Y})\)を用意してくれるだろう。$$ Y^2 = E(\hat{Y}^2) – E[v_c(\hat{Y})] = E[\hat{Y}^2 – v_c(\hat{Y})]$$ より、\( \hat{Y}^2 – v_c(\hat{Y}) \)は\(Y^2\)の不偏推定量だ。
　まとめると、$$ \hat{V}_{SRS} = N^2 \left( 1-\frac{n}{N} \right) \frac{1}{n(N-1)} \left[ \hat{U} – \frac{\hat{Y}^2 – v_c(\hat{Y})}{N} \right] $$ やったぜ、推定できるぜ。

　この特殊ケースとして、\(\hat{Y} = \sum_{i=1}^n w_i y_i\)の場合について考えよう。\( \hat{N} = \sum_{i=1}^n w_i \)は\(N\)と等しいとする(でないとしても近似してますよね)。\( \hat{\bar{Y}} = \hat{Y} / \hat{N} \)と書く。
　このとき次の式が成り立つ。Kish(1965)にもCochran(1977)にもSarndal et al.(1992)にも書いてあるから導出はそっちをみてね。$$ \hat{V}_{SRS} = N^2 \left( 1 – \frac{n}{N} \right) \frac{1}{n} \times \left[ \frac{1}{N-1} \sum_{i=1}^n w_i(y_i – \hat{\bar{Y}})^2 + \frac{v_c(\hat{Y})}{N(N-1)} \right] $$
　さて、この式の\( \frac{v_c(\hat{Y})}{N(N-1)} \approx v_c(\hat{\bar{Y}}) \) というのは\(d S^2/n\)の推定値で、\(n\)が大きければ無視できる。とっぱらうと$$ \hat{\hat{V}}_{SRS} = N^2 \left( 1 – \frac{n}{N} \right) \frac{1}{n} \times \left[ \frac{1}{N-1} \sum_{i=1}^n w_i(y_i – \hat{\bar{Y}})^2 \right] $$ $$ = N^2 \left( 1 – \frac{n}{N} \right) \frac{1}{n(N-1)} \left[ \hat{U} – \frac{\hat{Y}^2}{N} \right] $$ かくして、\( \hat{d} = v_c / \hat{\hat{V}}_{SRS}\)が手に入りますね。
　では、\(d S^2/n\)が十分に小さくない場合はどうなるか。\(\hat{\hat{V}}_{SRS}\)は過小評価になっちゃうけど、でも\(n/N\)が十分に小さければ$$ \hat{V}_{SRS} \approx \hat{\hat{V}}_{SRS} + \frac{v_c(\hat{Y})}{n} $$ だから、$$ \hat{d}’ = \frac{v_c}{ \hat{\hat{V}}_{SRS} + \frac{v_c}{n}} = \frac{\hat{d}}{1+\frac{\hat{d}}{n}} $$ と近似できる。[な・る・ほ・どーーーー]

　さらなる特殊ケースとして、\(y_i\)が二値の場合について考えよう。標本における割合を\(\hat{p}\)とする。もし標本がSRSならば$$ v_{SRS} = N^2 (1-\frac{n}{N}) \frac{\hat{p}(1-\hat{p})}{n} $$ と書ける(\(\frac{n}{n-1}\)は端折った)。
　いっぽう実際に求めたのは$$ \hat{p}_w = \frac{1}{N} \sum_i w_i y_i = \frac{\hat{Y}}{N} $$ だとする(実際には\(\frac{\hat{Y}}{\hat{N}}\)かもだけど)。これを使って$$ \hat{\hat{V}}_{SRS} = N^2 \left( 1-\frac{n}{N} \right) \frac{\hat{p}_w (1-\hat{p}_w)}{n}$$ が求まる(例によって\(\frac{N}{N-1}\)は端折った)。
　デザイン効果を求める際に、分母を\( \hat{\hat{V}} \)じゃなくて、手元の標本からうっかり求めた\(v_{SRS}\)にしてしまうと、\(\hat{p}\)と\(\hat{p}_w\)が違っていてしかも1/2から離れているときに、変な結果になっちゃうわけである。
　
5. デクラスタリングの影響が過大評価されるのはなぜか
　本章では、クラスタ抽出デザインをSRSに置き換えたとき(デクラスタリング)のベネフィットが過大評価されていることと、その理由を示す。問題は、良く用いられているデザイン効果の式の誤りにある。

　次の式が有名である。$$ d = 1 + (m-1) \rho_I $$ \(m\)はクラスタから抽出された最終要素数のクラスタ当たり平均。\(\rho_I\)は級内相関係数である。この式が、標本平均のような単純な推定量における、はたまた回帰係数のような複雑な推定量における、デザイン効果を表している… と思っていませんか。

　次の超母集団モデルを考えよう。クラスタ数を\(C\)、各クラスタに属する最終要素数を\(M\)とする。クラスタ平均は $$ (y_i, x_i) \sim N \left( 0^\top, \left[ \begin{array}{cc} \sigma^2_y & \rho_m \sigma_x \sigma_y \\ \rho_m \sigma_x \sigma_y & \sigma^2_x \end{array} \right] \right) $$ 最終要素は $$ (y_{ij}, x_{ij}) = (y_i, x_i) + e^\top_{ij} $$ $$ e_{ij} \sim N \left(0, \left[ \begin{array}{cc} \tau^2_y & \rho_c \tau_x \tau_y \\ \rho_c \tau_x \tau_y & \tau^2_x \end{array} \right] \right) $$ とする。$$ V(y_{ij} = \sigma^2_y + \tau^2_y$$ $$ V(x_{ij}) = \sigma^2_x + \tau^2_x$$ $$ cov(x_{ij}, y_{ij}) = \rho_m \sigma_x \sigma_y + \rho_c \tau_x \tau_y$$ $$ \mathrm{corr}(x_{ij}, y_{ij}) = \rho_u = \frac{\rho_m \sigma_x \sigma_y + \rho_c \tau_x \tau_y}{\sqrt{(\sigma^2_x + \tau^2_x)(\sigma^2_y + \tau^2_y)}} $$ となる[なんでいちいち書いているかというと、あとで補足変数を使った回帰推定量の話をするから]。級内相関は$$ \mathrm{corr}(y_{ij}, y_{ik}) = \rho_I = \frac{\sigma^2_y}{\sigma^2_y + \tau^2_y} $$ となる。

　まずは、クラスタを無視して、サイズ\(n\)のSRS標本を得たとしよう。$$ V(\bar{y}) = (1-\frac{n}{N}) \frac{S^2}{n} $$ である。\(S^2\)とは超母集団の量\(V(y_{ij}) = \sigma^2_y + \tau^2_y\)の推定量であるとみることができる。これを$$ V(\bar{y}) \sim (1-\frac{n}{N}) \frac{\sigma^2_y + \tau^2_y}{n} $$ と書こう。

　では、\(c=n/M\)個のクラスタの抽出だけをやったら？ 選んだクラスタの集合を\(s\)として、$$ \bar{y}_m = \frac{1}{c} \sum_{i \in s} \frac{1}{M} \sum_{j=1}^M y_ij = \frac{1}{c} \sum_{i \in s} \bar{Y}_i$$ の分散は、$$ V(\bar{y}_m) = (1-\frac{c}{C}) \frac{S^2_m}{c} $$ ここで\(S^2_m\)は\( V(\bar{Y}_i = \sum^2_y + \tau^2_y / M\)の推定量とみなせるから、$$ V(\bar{y}_m) \sim (1-\frac{n}{N}) \frac{\sigma^2_y + \tau^2_y / M}{c}$$ と書ける。デザイン効果は$$ d = \frac{M \sigma^2_y + \tau^2_y}{\sigma^2_y + \tau^2_y} = 1+(M-1)\rho_I$$ となる。[あながち嘘でもないってわけね]

　こんどは、サイズ\(n\)のSRS標本で、推定量が回帰推定量$$\bar{y}_r = \bar{y} + \beta_u (\bar{X} – \bar{x}) $$である場合を考えよう。\(\bar{X}\)と\(\beta_u = \rho_u \frac{S_y}{S_x}\)は既知とする。$$ V(\bar{y}_r = (1 – \frac{n}{N}) \frac{S^2_y}{n}(1-\rho^2_n) $$ $$ \sim (1 – \frac{n}{N}) \frac{\sigma^2_y + \tau^2_y}{n}(1-\rho^2_n) $$ と書ける。つまり\(1-\rho^2_n\)が標本平均と比べたゲインなわけね。「推定量効果」といってもよい。
　最後に、\(c\)このクラスタを抽出して回帰推定量$$ \bar{y}’_m = \bar{y}_m + \beta_m(\bar{X} – \bar{x}_m), \ \ \beta_m = \rho_m \frac{S_{my}}{S_{mx}} $$ を使う場合。$$ V(\bar{y}’_r) = (1-\frac{c}{C}) \frac{S^2_{my}}{c}(1-\rho^2_m)$$ $$ \sim (1-\frac{n}{N} \frac{M \sigma^2_y + \tau^2_y}{c} (1-\rho^2_m) $$ 「推定量効果」は\(1-\rho^2_m\)である。

　まとめよう。上記の超母集団モデルのもとで、分散を\( (1-\frac{n}{N})\frac{1}{N}\)で割った量は、

• SRSの平均: \(\sigma^2_y + \tau^2_y\)
• SRSの回帰推定量: \( (\sigma^2_y + \tau^2_y) (1-\rho^2_u)\)
• クラスタ抽出の平均: \(M\sigma^2_y + \tau^2_y\)
• クラスタ抽出の回帰推定量: \( (M\sigma^2_y + \tau^2_y) (1-\rho^2_m)\)

となるわけである。\(\rho_m\)はクラスタレベルの相関、\(\rho_u\)はクラスタを無視したときの相関である。

　推定量が平均だったら、デザイン効果は$$ d = \frac{M \sigma^2_y + \tau^2_y}{\sigma^2_y + \tau^2_y} = 1 + (M-1) \rho_I $$ となる。しかし回帰推定量だったら(つまり、どちらのデザインでも回帰推定量を使っていたら)、デザイン効果は$$ d_r = d + \frac{1-\rho_m^2}{1-\rho^2_u} $$ となるのである。
　ね？ \(d\)は過大評価でしょ？

　クラスタ抽出の回帰推定量と、SRSの平均の比は、デザイン効果というよりも、「デザインx推定量効果」なのである。より一般的にいえば、SRS分散\( (1-\frac{n}{N} \frac{S^2}{n}\)は、クラスタ当たり最終要素数を\(m\)として\(d = 1 + (m-1) \rho_I\)だけ増えるけれど、しかし\(1 – \rho^2_{xy}\)だけ減るわけである。

　以上、デザインと推定量の組み合わせが、デクラスタリングのベネフィットの過大評価につながるという話をお届けした。
　実は、過大評価の原因はこれだけではない。原因は全部で3つある。

• 推定量。
• インシデンス。ある疾患に10%の人が罹患していて(クラスタとは無関係に)、クラスタあたり10人抽出していたら、クラスタ当たりの罹患者の期待値は1である。母罹患者数のデザイン効果は1に近くなる。
• 抽出単位。世帯調査では人についてSRSを仮定することが多いが、現実にはそうでない(そしてもちろん世帯とはクラスタである)。[んんん？ そのせいでデザイン効果が過大になるの？　過小じゃなくて？ 頭がこんがらがってきた]

6. デザイン効果の算出におけるそのほかの問題
3つの問題を指摘する。

• SRS仮定の水準。たとえば、\(V_{SRS}\)を国レベルで求め、国レベルのデザイン効果を求めることができる。つまり、現実のデザインと、国レベルでのSRSとの比較である。いっぽう、実際の層\(h\)ごとに\(V_{h,SRS}\)を求めてこれを合計し、デザイン効果の分母とすることもできる。後者のほうがフェアだと思うかもしれないけれど、それぞれの層について\(V_{SRS}\)を求めるのは難しいかもしれないし、国レベルなら\(v_c(\hat{Y})N\)を無視できるけど層レベルではそうでないかもしれない。
• 推定量の選択。さっきのよりも単純な例を挙げると、\(y\)の推定に比推定量\(\frac{\hat{Y}_c}{\hat{X}_c} X\)を使っているとき、比べるべきは\(\hat{Y}_{SRS}\)ではなく\(\frac{\hat{Y}_{SRS}}{\hat{X}_SRS} X\)であろう。
  　分母側を単純すぎる推定量の分散にしてしまう理由の一つは、複雑な推定量を扱うのが難しいからである。[…比推定量の分散の推定の話とか、まあ分母と分子で推定量が違っていてもしょうがないかもね、とか、そういう話。パス]
• 最終単位の選択。多段デザインでは、デザイン効果を求める際につかう単位をどこにするかがはっきりしないことがある。特に最後の段階が全数だった場合がそうだ。もしSRSが実現可能ならどうしていたか、という観点から考えるという路線とか、そもそも知りたい変数はなんなのかと考える路線(世帯特性が知りたいのなら世帯のSRSだと考えるのが自然だ)とかがある。

——————-
　いやー、これは勉強になった…！
　正直なところ、知りたいこととはちょっと主旨が違っていて、最大のポイントである5章については「ふーん」という感じだったんだけど(クラスタ抽出って仕事の中であまり出てこないし)、その手前の4章、デザイン効果の分母とすべきSRS分散を推定する方法についての説明に感心しきりで、やたらに細かくメモをとってしまった。
　非専門家向けの啓蒙論文とはこのようなものであってほしい。Kish先生もちょっと見習ってほしいね。(すいません冗談です)