今回はAbadie, Athey Imbens and Wooldridge (2017)（以下AAIW）についてお話します。とはいっても、ほとんど自分の勉強の為の備忘録といった内容で、すごくつまらないと思うので先に謝っておきます。

この論文はNBER Working Paperで、論文というより計量経済手法のユーザー向けガイダンスといった感じです。少し前に回帰分析をしていて、ふと「固定効果モデルを使う時ってクラスターはどう考えればいいのかな」と思い見つけた論文です。

因みに、この論文に関しては世界銀行で開発経済学者をしているDavid McKenzieさんがブログを書いていますので、もし興味がある方は是非そちらも読んでみてください。

問題定義：計量経済理論と実証研究のギャップ

AAIWが本論文を書いた一番の動機は実証経済文献内でクラスターの使い方に関して混乱があり、 その大きな理由に「理論ベース」と「デザインベース」なクランスタリングへの見方にギャップがあるから。

ここで著者がいう「デザインベース」とは、サンプリングデザイン及び実験デザインの事をいう。

彼らはクラスタリングは、根本的にサンプル抽出及び実験デザインに関する問題として扱われるべきであるとする。

実証経済文献内でクラスターの使い方に関して混乱は、大きく２つの問題に分けられる。

問題１：クラスター調整の動機

計量経済理論に基づいたクラスター調整の動機としては、従属変数において見えない部分がクラスター内（例えば、村）で相関性を含んでいる可能性があげられる。

この動機づけだと、以下にあげる３つの点において曖昧なままである。

1. どのクラスターを使うべきなのか明確でない。
2. ランダム化実験でクラスターが使われないのか説明できない。
3. 固定効果モデルにおいてクラスター調整を動機づけできない。

１点目は、例えば、教育年数が賃金に及ぼす影響を見る回帰式を推定する際、都道府県クラスターまたは年齢グループ毎にクラスター調整するべきなのか、 上述の動機づけだと明確でないことを言及する。

問題２：クラスターレベルの設定方法

一般的には、クラスターレベルは一番大きなグループを使うのが良いとされ、クラスターの数が少なくなりすぎないようにするとされている。

本論文では、大きなクラスターは必ずしも良い方法ではない事も明確にする。

問題の実演

AAIWは実際クラスターがどう標準誤差に影響するのか見るために、ロバスト分散とクラスター分散を架空データを使って比較する。 （標準誤差 \(= \sqrt{\frac{分散}{観測数}}\)なので、分散だけに注目。）

セットアップ

推定したい回帰式は以下の通り。

\[Y_{i} = \alpha + \tau W_{i} + \epsilon = \beta'X_{i} + \epsilon\]

\(W_{i}\)は特定のクラスター\(C_{i} \in \{1, ...,C\}\)に所属するかどうかを表すダミー変数。

３つ目のイコールは\(\beta' = (\alpha, \tau)\)と\(X_{i} = \begin{pmatrix} 1 \\ W_{i} \end{pmatrix}\)と設定して書き換えられたもの。

\(\beta\)の推定は最小二乗法を使うと、お馴染みの\(\hat{\beta} = (\boldsymbol{X}'\boldsymbol{X})^{-1}(\boldsymbol{X}'\boldsymbol{Y})\)が得られる。 （\(\boldsymbol{X}\)はNx2マトリックス、\(\boldsymbol{Y}\)はN-ベクター。）

ここで重要なのは、理論ベースで考えると、\(\epsilon_{i}\)の塊であるN-ベクター\(\boldsymbol{\epsilon}\)はランダムと考えられるが、 クラスター内での相関性がある場合、エラー\(\boldsymbol{\epsilon}\)と\(X\)独立ではなく、クラスターも条件に入れなくては行けないところ。

なので\(\boldsymbol{C}\)をクラスターのN-ベクターとすると、 \(E[\boldsymbol{\epsilon} \vert \boldsymbol{X}, \boldsymbol{C}] = 0\)と\(E[\boldsymbol{\epsilon \epsilon'} \vert \boldsymbol{X}, \boldsymbol{C}] = \Omega\) という条件をつける必要がある。

更に大切なのは、クラスターがない場合、\(\Omega\)は対角成分が分散で、残りは0の対角マトリックスであり、 分散均等性を仮定できる場合、\(E[\boldsymbol{\epsilon \epsilon'} \vert \boldsymbol{X}] = \sigma^2\)。

よって、\(\hat{\beta}\)の分散は

\[V(\hat{\beta}) = \sigma^2(\boldsymbol{X}'\boldsymbol{X})^{-1}\]

となる。

分散非均一性を仮定下でのロバスト分散値は

\[V_{robust}(\hat{\beta}) = (\boldsymbol{X}'\boldsymbol{X})^{-1} \bigg ( \sum_{i=1}^{N} \Omega_{ii}X_{i}X'_{i} \bigg ) (\boldsymbol{X}'\boldsymbol{X})^{-1}\]

であり、本文ではこの分散の推定値がベンチマークとして使われる。

クラスターがある場合、観測単位のエラーが独立せず（\(E[\epsilon_{i}\epsilon_{j} \vert \boldsymbol{X}, \boldsymbol{C}] \neq 0 , i \neq j\)）、対角成分以外が０ではなくなってしまうので、\(\Omega\)が対角マトリックスであると仮定できない。

この状況下で分散共分散行列は

\[\Omega_{ij} = \left\{ \begin{array}{ll} 0 & \text{if } C_{i} \neq C_{j} \\ \rho \sigma^2 & \text{if } C_{i} = C_{j}, i \neq j \\ \sigma^2 & \text{if } i = j \\ \end{array} \right.\]

１つ目は\(i\)が\(j\)違うグループに所属する場合で、共分散はない。２つ目は\(i\)が\(j\)が同じグループに所属している場合、共分散に何らかの値がつく。 この値の大きさはグループ内の相関性の強さ\(\rho\)による。３つ目は対角成分をしめる分散値。

よって、クラスターのサイズは均等と仮定した場合の分散値は

\[V_{cluster}(\hat{\tau}) = V(\hat{\beta})\bigg( 1 + \rho_{\epsilon} \rho_{W} \frac{N}{C} \bigg )\]

\(\rho_{\epsilon}\)はクラスター内におけるエラーの相関性。\(\rho_{W}\)はクラスター内における独立変数の相関性。 クラスター相関性が標準誤差に影響を与えるかどうかが、この２つの数値にかかっている。

この式からすると、\(\rho_{\epsilon}\)か\(\rho_{W}\)のどちらかが０であれば、\(V_{cluster}(\hat{\tau}) = V(\hat{\beta})\)でクラスター調整する必要はないと考えられる。 ということは、\(W\)がランダムに割り当てられていれば、\(\rho_{W}\)は必然的に０となり、クラスター調整不要と結論づけられる。 さらに、固定効果モデルを使えば、エラーにおける相関性は排除されるので、これまた調整不要と考えられる。 以下の比較で、AAIWはこれがいかに間違った考え方か実演してくれる。

さらに、パネルデータを使用する場合等、時間依存性に対応するクラスター分散値は

\[V_{robustCluster}(\hat{\beta}) = (\boldsymbol{X}'\boldsymbol{X})^{-1} \bigg ( \sum_{c=1}^{C} X'_{c} \Omega_{c} X_{c} \bigg ) (\boldsymbol{X}'\boldsymbol{X})^{-1}\]

本文ではこの分散の推定値（\(\hat{V}_{robustCluster}(\hat{\beta})\)）をベンチマークである\(\hat{V}_{robust}(\hat{\beta})\)と比較する。

（\(\hat{V}_{robustCluster}(\hat{\beta})\)と\(\hat{V}_{robust}(\hat{\beta})\)の式は本文をご参照ください。）

比較

AAIWは２つの誤信を実演によって明かしている。

誤信１：\(\rho_{\epsilon}\)と\(\rho_{W}\)の両方が＞０でないかぎり、クラスター調整は不要。

シミュレーション環境は、観測数N=100,323、クラスター数C=100、クラスター一つの大きさは平均1,000ほど。 OLSを使って\(Y_{i} = \alpha + \tau W_{i} + \epsilon\)を推定する。

まずは\(\rho_{\epsilon}\)と\(\rho_{W}\)の推定。\(\rho_{.}\)はクラスター内相関性の強さなので、

\[\frac{\vert クラスター有り分散値ークラスター無し分散値 \vert}{クラスター無し分散値}\]

として計算される。その結果、\(\hat{\rho_{\epsilon}} = 0.001\)と\(\hat{\rho_{W}}=0.001\)といずれの相関性も低いことがわかる。

驚きなのが、相関性が低いにも関わらず\(\hat{V}_{robustCluster}(\hat{\beta})\)と\(\hat{V}_{robust}(\hat{\beta})\)を使った標準誤差にかなり大きな差がある。係数が\(\hat{\tau}^{OLS} = -0.120\)に対して、\(SE_{robust} = 0.004\)、\(SE_{robustCluster} = 0.100\)という結果だった。

この誤信の理由としてAAIWが挙げているのは、\(\rho_{\epsilon}\)と\(\rho_{W}\)の個々が大事という誤解。著者によると、本当に大切なのは\(\rho_{\epsilon W}\)、つまりクラスター内における残差と独立変数の積の相関性の強さとのこと。実際に、本シミュレーションでは\(\rho_{\epsilon W}=0.5\)。確かに後者の数値を見た場合の結論は大きく変わる。

更に、上の回帰式を固定効果を入れてOLSで推定すると、\(\hat{\tau}^{FE} = -0.120\)に対して、\(SE_{robust} = 0.003\)、\(SE_{robustCluster} = 0.243\)という結果で、いずれにしても大きな差がでた。

誤信２：クラスター内相関性に懸念がある場合はクラスター調整するべき。

クラスター内相関性がありそうであれば、サンプルサイズが大きければ害はないので、とりあえずクラスター調整しておいたほうがいいという誤信。

この誤信を解くために使われるシュミレーション環境は、母集団の大きさが10,000,000個体、\(C=100\), 個々のクラスターのサイズはきっちり100,000個体。個体\(i\)は50\%の確率で\(W_{i}\in \{0,1\}\)をランダムに割り当てられる。\(i\)の結果変数は\(Y_{i} = \tau_{C_{i}}W_{i} + \nu_{i}\)と定められている。（\(\nu \sim N(0,\sigma^2)\)）一つ注目なのは、係数の\(\tau_{C_{i}}\)はクラスターごとの係数である。加えて、50\%の確率で\(W_{i}\)の数値を与えられているため、平均治療効果は０。10,000回シミュレーションして結果は以下の通り。

固定効果無					固定効果有
\(\sqrt{\hat{V}_{robust}}\)	Coverage rate	\(\sqrt{\hat{V}_{robustCluster}}\)	CR		\(\sqrt{\hat{V}_{robust}}\)	CR	\(\sqrt{\hat{V}_{robustCluster}}\)	CR
0.007	0.950	0.051	1.00		0.007	0.950	0.131	0.986

Coverage rateとは、10,000回シュミレーションした内、真の平均治療効果（０）が95\%信頼区間に収まるシュミレーションの割合。つまり、この数値は0.95に近くあるべき。

結果を見ると、ロバストクラスター推定値は固定効果があってもなくても、ロバスト推定値と比べるとかなり大きい。ロバストクラスター推定値が大きい理由は、この計算方法ではサンプルされたクラスターの数（100）よりも多くのクラスターが母集団に存在すると仮定されており、今回のシミュレーションではこの想定が間違っていることである（母集団のクラスター数も100とされている）。

解決案：「デザインベース」な考え方

よって理論だけではクラスター調整するべきなのか、及びどんなふうに調整するべきなのか明確でない。 このジレンマにAAIWは以下の実務的なアドバイスをだしている。

1. 研究デザインを考慮
   • 研究者はまずサンプリングの時点でクラスターされているか考えるべき。
   • その次は治療の割当がクラスターされていたか考えること。

この２つの質問への答えがノーであれば、クラスター調整するべきではない。

加えて、

1. ロバストクラスター標準誤差は以下の３つの条件の内の一つが当てはまらない限り基本的に保守的な数値である。
   • 治療効果が均一
   • 実際にあるクラスターの一部しか観察できない
   • 最低でも各クラスターから１個体はサンプルされていること
2. クラスタリングによる上向きバイアスはサンプルにおけるすべてのクラスターを含め、且つ、各クラスター内で治療割当に変化があれば修正される。
3. 固定効果を回帰式入れる場合、クラスター調整するには治療効果における不均一性がなくてはならない。

意見・感想：

何が問題で、何が解決策なのか直観的ロジックと実演でとても分かりやすく説明してあり、すごく勉強になりました。 論文の後半にはフォーマルな証明がされているので、興味がある人は是非！ （前半はMostly Harmlessチックで、後半は大変。笑）

この論文の勉強のお供に使っていたのが、私の学校の計量経済の先生でBruce Hansenという人が書いた教科書。教科書というよりもはや百科事典という感じですが、分散共分散行列の話等昨年度習ったのにもう忘れちゃったの復習に使いやすいので興味のある方は是非。（この内容で如何せんタダですし、目次にハイパーリンクがついてないから使いづらいって言ったら、きっと怒られるだろうなー。）

最後に

今週は私の勉強に付き合っていただきありがとうございました！なにかおすすめの計量経済、特に実務に密着した論文等ありましたら、是非教えてください。

本ブログ記事に対するご感想や、本ブログ全体に関わるご意見などもありましたら、下にコメントを残すか、econ.blog.japan@gmail.comまでご連絡ください。

また、Twitterアカウントの@EconJapanのフォローもよければお願いします。Twitterでの本ブログのコメント・拡散も歓迎です。その際は、＃econjapanblogをお使いください！

渋谷

コメント

元橋:

本ペーパーで、（１）サンプリングデザインのためのクラスタリングと、（２）実験デザインのためのクラスタリング（Treatmentをどのレベルで設定するかに関わる）の大きく２点について議論されているが、（１）については馴染みがなかったので勉強になった。

（２）については、TreamentのAssignmentがなされたレベルで、SEをクラスタリングすれば良いと、学んできたので、あまり新規性はなかった（このペーパーに基づいて授業で教えられていただけかもしれないが…）。他方、固定効果がある場合に、SEのクラスタリングの影響があるのが、Treatment effectに異質性がある場合のみというのは勉強になった。

鈴木:

「ちゃんと読まなければ…」と2年半くらい思っていた論文を解説していただいて助かりました笑 記事にある通り、論文内で「こういうときにこういう clustered SE を使いましょう」「こういうときは使ってはいけません」という指針があって、これ自体はとても有用だと思うのですが、いざ自分のしている研究にあてはめてみようとすると「ん？あれ？これどうすりゃいいんだ？」みたいになりました。 論文で行われているものに加え、特に自然実験を想定したケースのシミュレーションをしてみると、どういうときにどういうバイアスが生じるかわかりやすいかなと思いました（次の私の pet project になりそうである…）。

標準偏差の上方バイアスは、「正しい標準偏差を用いれば『推定値が統計的に有意に0と異なる』という結果が得られたであろう場合にも『推定値が統計的に有意に0と異ならない』という結果を得てしまう」ことにつながり、統計的に有意な結果が重視される昨今のクソみたいな事情から、そのような研究は学術誌に投稿されず、研究プロジェクトがお蔵入りしてしまうことになるかもしれません（file drawer problem）。 ちょっと関連する話で、この記事ではStataの areg コマンドで計算される clustered SE は上方バイアスがかかっており、プロジェクトがお蔵入りする原因の一つになっているという話がされています。 もちろん、今回の論文にある方法で正しい標準偏差を計算しても、「統計的に有意な結果が得られなかったプロジェクトがお蔵入りする」という file drawer problem が改善されるわけではありませんが、より正しいと思われる標準偏差を使うことで、file drawer にしまわれそうだったプロジェクトが世に出て、それが社会の役に立つのであればそれはいいことなのかなとぼんやり思いました。

文献：

Abadie, A., Athey, S., Imbens, G.W. and Wooldridge, J., 2017. When should you adjust standard errors for clustering? (No. w24003). National Bureau of Economic Research.