« 読了：Aguilera, et al. (2006) ロジスティック主成分回帰 | メイン | 読了：Hellevik (2009) 二値の従属変数に対してロジスティック回帰とかじゃなくて線形回帰をやって、いったい何が悪いんだよ »
2016年6月23日 (木)
週末に学会の大会を聴講しに行った空き時間で読んだ奴。自分が発表しない学会というのは気楽なものである。
King, G., Tomz, M., Wittenberg, J. (2000) Making the most of statistical analysis: Improving interpretation and presentation. Amerian Journal of Political Science, 44(2), 341-355.
社会科学者よ、統計学的な分析結果を人に伝えるときはこういう風にしなさい、という啓蒙論文。Google様によれば被引用回数3241。すげえ。
あ、第一著者のKingって、名著と名高い「社会科学のリサーチ・デザイン」の著者だ... ちゃんと読んでませんけど。
冒頭に挙げられている良い伝え方の例：「ほかの点がすべて同一であれば、教育年数が１年増えると、年収は平均して1500ドル(±約500ドル)増えるでしょう」。悪い伝え方の例：「教育の係数は有意水準0.05で統計的に有意でした」。
ポイントは次の３つ。(1)標準的な統計モデルから、関心の対象となる新しい量を抽出すること。(2)その量の不確実性を評価すること。(3)統計的訓練を受けていなくてもわかる結果に変換すること。
以上を実現するための有力な武器がシミュレーションである。
統計モデルの非常に一般的なクラス、すなわち
$Y_i \sim f(\theta_i, \alpha)$
$\theta_i = g(X_i, \beta)$
を考える。一本目は統計モデルの確率的コンポーネントで、従属変数$Y_i$が確率密度$f(\theta_i, \alpha)$からのランダムドローとして生成されている。確率密度関数の特性はオブザベーションによって変動するかもしれないし($\theta_i$)、一定かもしれない($\alpha$)。二本目はモデルのシステマティックなコンポーネントで、$\theta_i$がどう変動するかを示している。$g(\cdot, \cdot)$はリンク関数と呼ばれることが多い。
このクラスのメンバーであるなんらかのモデルを考える。たとえば線形正規回帰モデルなら
$Y_i \sim N(\mu_i, \sigma^2)$
$\mu_i = X_i \beta$
ロジットモデルなら
$Y_i \sim Bernoulli(\pi_i)$
$\pi_i = \frac{1}{1+\exp(-X_i \beta)}$
ですわね。まあとにかく、なんらかのモデルをつくり、結果が得られた、としましょう。
ここからが本題である。モデルから得られた結果をどうやって解釈するか。
多くの研究者は$\hat\alpha$, $\hat\beta$の符号と「統計的有意性」しかみない。でもそれらはふつう、研究の動機となっている実質的問題と間接的にしか関連してない。実質的な関心が直接に持たれるような量を提示すべきだ。
また、つぎの2つの不確実性を無視してはならない。(1)推定の不確実性。$\beta$と$\alpha$は完全にはわからない。(2)根本的な不確実性。仮に$\beta$と$\alpha$が完全にわかったとしても$Y$には不確実性が残る。
そこでだ。諸君、シミュレーションしたまえ。
シミュレーションとは、サーベイ・サンプリングの理屈を使って複雑な数学的計算を近似することだ。たとえば確率分布$P(y)$の平均を計算するために、$E(Y) = \int^{\infty}_{-\infty} y P(y) dy$を求めるんじゃなくて、$P(y)$から$M$個の値をドローしてきて平均するわけである。$M$を増やせば正確になる。
まずはパラメータのシミュレーション。手順は次の通り。
- (1)$\alpha$, $\beta$の点推定値とその分散共分散行列を求める。通常のソフトで出力できる(たぶん最尤法で)。以下、$\hat\beta$と$\hat\alpha$を縦積みしたベクトルを$\hat\gamma$、その共分散行列を$\hat{V}(\hat\gamma)$とする。
- (2)多変量正規分布$N(\hat\gamma, \hat{V}(\hat\gamma))$からベクトルを一本ドローする。これを$\tilde\gamma$とする。
(2)を$M$回繰り返す。たとえば1000回とか。
次に、予測値のシミュレーション。手順は次の通り。
- (1)どんな予測値が得たいのかを決め、説明変数の値を決める。このベクトルを$X_c$としよう。
- (2)ある$\tilde\gamma$をつかって、$\tilde\theta_c = g(X_c, \tilde\beta)$を算出する。
- (3)前項の$\tilde\gamma$を使い、$f(\tilde\theta_c, \tilde\alpha)$から値をひとつドローする。これを$\tilde{Y}_c$とする。
必要ならば、従属変数の期待値についてもシミュレーションするがよい。なお、正確に言うと「従属変数の期待値」と「従属変数の予測値の平均」とはちがうのだが、非線形性がシビアでないかぎり両者はだいたい近くなる。
従属変数の期待値のシミュレーションと、予測値のシミュレーションとはちがうぞ。後者には二種類の不確実性がはいっているが、前者には推定の不確実性しか入っていない。たとえば選挙結果の予測とか為替レートの予測という場面では後者が大事だが、特定の説明変数の平均的な効果に関心がある場合には前者が大事かも。
手順は次の通り。
- (1)説明変数の値$X_c$を決める。
- (2)ある$\tilde\gamma$をつかって、$\tilde\theta_c = g(X_c, \tilde\beta)$を算出する。
- (3)前項の$\tilde\gamma$を使い、$f(\tilde\theta_c, \tilde\alpha)$から値をm回ドローして平均する。これが$\tilde{E}(Y_c)$。$m$を大きくすることで根本的な不確実性を取り除くことができる。
(2)と(3)を$M$回繰り返す。このとき、$M$と$m$は十分に大きくすること。なお、線形正規モデルやロジットモデルでは$E(Y_c)=\theta_c$なので$\tilde\theta_c$をそのまま使えばよろしい。
第一階差のシミュレーション。第一階差とは、2つの期待値の差のこと。上の手順の(1)で、$X_c$を2つ用意する($X_s, X_e$としよう)。で、(5)で$\tilde{E}(Y_s)$ と$\tilde{E}(Y_e)$の差を求める。これを繰り返して平均する。
ところで、たとえば順序プロビットモデルで$P(Y=3)$を求めるというような場合には、期待値を推定するアルゴリズムをちょっと修正しなければならない[←??? なんでだろう...]。そんなときは僕らが作ったCLARIFYというソフトを使うといいよ。
本論文で紹介したのとはちょっと別なアプローチとして以下がある。
- 完全にベイジアンな手法。本論文のように中心極限定理に頼って正規分布に基づく漸近解をもとめるんじゃなくてMCMCでやる。でも収束の判定が難しい。CRARIFYに載っているよ。
- ブートストラッピング。使いやすくて強力。でもある種の量($Y$の最大値とか)については推定値が歪む。CRARIFYに載っているよ。
- シミュレーションじゃなくて、デルタ法。これは確率変数の非線形関数を近似するあるツールを使う方法で... 非線形関数$g(X_c, \beta)$をテイラー展開して...[読んでもよくわからんのでパス]。計算が難しいし、結局は近似なので、あんまし使われていない。
まとめると、シミュレーションは便利だ。分析的な解がないときにも正確な結果が得られる。教育上も良い。ある研究者は「それでも分析的手法を教えるべきだ」という理由を挙げることができた人に5000ドルをあげると宣言しているが、この賞金を受け取った人はまだいない[ははは。Simonという人だそうだ]。ま、一番いいのは両方教えることだけどね。
シミュレーションの際のコツ。
- シミュレーションにはパラメータ推定量のSEだけじゃなくて、完全な分散行列$\hat{V}(\hat{\gamma})$が必要である。ちゃんと出力するよう、ソフトを設定しなければならない。
- MVNからドローするときいくつかのパラメータを除外しちゃうひとがいるけど($\beta$と$\alpha$のどちらかだけとか)、これはよくある間違い。時間はかかるかもしれないけど、必ず全部ドローすること。
- モデルによっては、$\gamma$の要素が直交していて、ソフトが別々の分散行列を出してくることがある。ちゃんと並べてブロック対角行列をつくること。めんどくさいようでも、とにかく$\hat{V}(\hat{\gamma})$を構成しちゃったほうが間違いがない。
- 有限標本で多変量正規近似が維持されるようにするためには、$\gamma$の要素を再パラメータ化しなければならない。正規分布の場合のように、すべてのパラメータに制約がかかってなくて論理的にシンメトリックであれば話は別だが、そうでない限りパラメータは再パラメータ化する必要がある。たとえば、$\sigma^2$のような分散パラメータはゼロより大きくないといけないので、$\sigma^2=\exp(\eta)$というような表現をつかって再パラメータ化するとよい。で、$\gamma$の要素の一つが$\eta$だと考える。こうすると$\gamma$はMVNに従うから$\eta$は$-\infty$と$+\infty$の間の値をとる。で、あとで$\hat\sigma^2 = \exp(\hat\eta)$を求め直すわけである。相関$\rho$だったら、FisherのZ変換$\rho=(\exp(2\eta)-1)/(\exp(2\eta)+1)$を使う。確率$\pi$だったら、ロジスティック変換$\pi = [1+\exp(-\eta)]^{-1}$を使う。
- 実をいうと、いつもいつも$Y$をシミュレートしないといけないわけではない。たとえばロジットモデルだったら、$\tilde{E}(Y)$を求めるためには$\tilde\pi$を求めれば十分なわけで、わざわざ二値の$Y$をドローする必要はない。でも、よくわかんないなと思ったら、とにかく$Y$をシミュレートすること。
後半は事例。線形回帰、ロジットモデル、時系列クロスセクショナルモデル、多項ロジットモデル、打ち切りのあるワイブル回帰モデルについて例を挙げている。最初のふたつだけ読んだ。メモは省略するけど、正直、前半の説明よか事例のほうがはるかにわかりやすいよ... 最初っから事例を使って説明してくれればいいのに...
結局、仕事の役には立たなかったんだけど、勉強になりました。特に、「2種類の不確実性」という言い方が勉強になった。以前、若い友人たちと新宿のカフェで応答曲面モデルの勉強会をやったとき、回帰モデルから得られる、ある条件下での期待値の信頼区間と、その条件下での予測の信頼区間とは全然ちがう問題なのよという話になり、私の説明が下手なせいであんまし納得してもらえなかったんだけど、そうだよな、こういうことなんだよな。
いっぽう、素朴な疑問でこっぱずかしいが、こうやってパラメータについて無理やりMVNを仮定するんじゃなくて、全部ブートストラップ法でいいんじゃない？ というモヤモヤ感がある。
また、この論文はあるモデルから得られる結果をどうやってわかりやすく伝えるかという話けど、モデルの不確実性(変数選択やリンク関数選択の不確実性)をも考慮して、一連のモデルからなる集合から得られる結果について伝えるときに、こういうシミュレーションによる方法をうまく使えないもんかなあ、という疑問もある。
論文：データ解析(2015-) - 読了：King, Tomz, & Wittenberg (2000) 統計モデルから得られる知見について人々にわかりやすく伝えるためのシミュレーションの手引き