## はじめに 武力行使の政治学 第3章の中でたびたび登場する`マルチノミアル・ロジットモデル`と`ヘックマン・プロビットモデル`に関連して、ロジットモデルとプロビットモデルに関する補足を行うことがこの資料の目的である。なお著者は数式自体を追ってはいるが、その裏にある数理統計的な裏付け(なぜ、〜と与えるとうまくいくのか?)については理解していないので、その辺りには触れないで頂けると幸いです。 # 統計学入門 (東大出版) ## 第13章 回帰分析 ### intro 回帰分析は、2変数$X$,$Y$のデータがあるとき**回帰方程式**と呼ばれる説明の関係を定量的に表す式を求めることを目的としている。 (注意するべきは、$Y$を$X$で説明することを試みるのであって、$X$と$Y$との相関分析ではない。) ### 13.1 回帰分析 説明変数$X$、被説明変数$Y$について、$i$番目を$Y_i, X_i$とし、ばらつきの部分を$\varepsilon_i$とすると、母集団において $$Y_i = \beta_1 + \beta_2 X_i + \varepsilon_i$$ となる。このモデルを母回帰方程式とよび、$\beta_1,\beta_2$を母(偏)回帰係数という。 誤差項 $\varepsilon_i$は以下の3つの性質を満たす確率変数 1. 期待値 $E(\varepsilon_i) = 0$ 2. 分散は一定: $V(\varepsilon_i) = \sigma^2$ 3. 異なった誤差項は無相関: $i\neq j$ならば$Cov(\varepsilon_i, \varepsilon_j) = 0$ これらより、$Y_i = \beta_1 + \beta_2 X_i + \varepsilon_i$のモデルが $$E[Y_i] = \beta_1 + \beta_2 X_i $$ であることを意味する。($E(\varepsilon_i) = 0$と、$E[X_i]=X_i$であることより($X_i$は確率変数ではなく、すでに確定した値をとっていることに注意)) ### 13.2 回帰係数の決定 #### 13.2.1 最小二乗法による回帰係数の推定 (詳しくは、東大出版のp.260-264)(略) #### 13.2.2 決定係数 (略) ### 13.3 偏回帰係数の統計的推測 #### 13.3.1 偏回帰係数の標本分布 最小二乗法の回帰係数の推定において、正規方程式を解いた結果が、以下のようになるので $$\hat \beta_2 = \frac{\sum (X_i - \bar X) (Y_i - \bar Y)}{\sum (X_i- \bar X)^2}, \hat \beta_1 = \hat Y - \hat \beta_2 \bar X$$ 最小二乗推定量$\hat \beta_1, \hat \beta_2$(標本偏回帰係数)は以上のように表せる。 ここで、13.1であげた誤差項 $\varepsilon_i$に対する仮定に加えて、誤差項$\varepsilon_1, \varepsilon_2, \varepsilon_3, ..., \varepsilon_n$が独立で共通の正規分布$N(0.\sigma^2)$に従うとする。すると、標本偏回帰係数$\hat \beta_2$は、正規分布にしたがっている$\varepsilon_i$の線形関数であることより、(ここいまいち分かっていない。)正規分布の再生性より、その標本分布はふたたび正規分布となり、 $$N(\beta_2, \frac{\sigma^2}{\sum (X_i - \bar X)^2})$$ となる。 #### 補足 回帰残差 実測値$Y_i$の回帰方程式に定められた回帰値$Y_i$からのズレ$\hat e_i = Y_i - \hat Y_i = Y_i - \hat \beta_1 - \hat \beta_2 X_i$は、$X$で説明されず残った部分であり、回帰残差と呼ぶ。なお、この回帰残差は、正負両方の値を含むが、全体の和は0となる。 また、誤差項$\varepsilon_i$の分散$\sigma^2$は、回帰方程式のあてはまりの良さを表し、 $$s^2 = \sum \hat e_i^2 / (n-2)$$ で推定される。 (補足終了) このとき、$\sigma$は未知であるので、推定、検定を行う際は、回帰残差$\hat e_i = Y_i - \hat \beta_1 -\hat \beta _2 X_i$を用いて、上記の補足で求めたように、$s^2 = \sum \hat e_i^2 / (n-2)$のように標本から推定を行い、推定値の標準誤差$s.e. = \sqrt{s^2} = \sqrt{\sum \hat e_i^2 / (n-2)}$とする。 以上から、$\hat \beta_2$の標準誤差は $$s.e.(\hat \beta_2) = \frac{s.e.}{\sqrt{\sum (X_i - \bar X)^2}}$$ となる。 ### 13.4 重回帰分析 重回帰方程式は不空数の説明変数を含むもので $$Y_i = \beta_1 + \beta_2 X_{2,i} + \beta_3X_{3,i} + ...+ \beta_k X_{k,i} + \varepsilon_i$$ となるモデル。(ただし、$X_2,X_3,...,X_k$は説明変数であり、$\varepsilon$は誤差項で以前に述べた性質を満たすものとする。) 重回帰方程式は、$k$個の未知の母偏回帰係数$\beta_1,\beta_2,...,\beta_k$を含むがその推定には、最小二乗法が用いられる。(これは単回帰の場合と同様) すなわち、$S = \sum \varepsilon_i^2$を考え、これを最小にする$\beta_1,\beta_2,...,\beta_k$を求める。 これは、 $$\frac{\partial S}{\partial \beta_1} = 0,\frac{\partial S}{\partial \beta_2} = 0,...,\frac{\partial S}{\partial \beta_k} = 0$$ を解くことにほかならない。 (解く過程は略) これによって得られた解を、$\hat \beta_1,\hat \beta_2,...,\hat \beta_k$とすると、 $$Y = \hat \beta_1 + \hat \beta_2 X_2 + \hat \beta_3 X_3 + ...+ \hat \beta_k X_k$$ となる標本重回帰方程式が得られる。 # 自然科学の統計学(東大出版) ## 第8章 質的データの統計的分析 ### 8.1 二値データ #### 8.1.1 質的データ 分析対象となる事象には数値データを得ることができず、標本の質的特性のみがデータとして得られる場合(今回の場合は、アメリカが単独軍事行動したかどうか?など)のように、対象がある状態やカテゴリーの項目に属しているかどうかの身を知ることができるデータを**質的データ**と呼ぶ。 質的データの中でも最も基本的であり、重要なのは、観測結果が成功、失敗のように2つの状態をとる二項的(binary)な場合であり、この場合、成功の場合を1、失敗の場合を0のダミー変数$Y_i$を使って観測結果を表す。また、このようなデータは2値データと呼ばれる。(本文中でも、`1が軍事行動の発動された場合、0がなかった場合`のように観測結果を2値データ化している。) #### 8.1.2 二値データの確率モデル 以下では例として薬物の投与水準に基づく生物の死亡数について考える。 $Y_i$を$i$番目の個体に対する薬物の効果を表す変数として、薬物によって生物が死亡する場合 1、生存する場合 0とする。$Y_i = 1$となる確率(=死亡率)は、一定量を単位とした薬物の対数濃度$X_i$に影響され、濃度が高くなるに従ってその値は大きくなるとする。 このようなとき、ある$X_i$の値に対する$Y_i = 1$となる確率$P_i$は、$Y_i$の$X_i$に関する条件付き確率$P(Y_i=1|X_i) = \frac{P(X_i \cap Y_i = 1)}{P(X_i)}$であるが、これは$X_i$の関数であって、以下のように表せる。 $$P(Y_i=1|X_i) = F^*(X_i)$$ ここで、$X_i$の関数としたが、$F^*$としてどのような関数を使うかが問題になる。 最も簡単なのは$F^*$を線形関数$F^*(X_i) = \alpha_0 + \alpha_1 X_i$とすることであり、これは線形確率モデルと呼ばれるが、$X_i$の小さな値や大きな値の際に$F^*(X_i)$が負になったり1を超えてしまったりするような問題を孕んでいる。(※ $0 \leq F^* \leq 1$が成立していないといけない。(確率なので)) 二値データの分析で一般に使われるモデルは回帰分析を応用した以下の考えであり、それは$Y_i$が0をとるか1をとるか決定する仮想的な因子$Y_i^*$があり、それが $$Y_i^* = \beta_0 + \beta_1 X_i + \varepsilon_i$$ で表すことができるとするものである。($\varepsilon_i$は誤差項と呼ばれる。) 我々は、この仮想的な$Y_i^*$を直接観測できないが、$Y_i$は$Y_i^*$の符号によって $$Y_i = \left\{ \begin{array}{ll} 1 & (Y_i^* > 0) \\ 0 & (Y_i^* \leq 0) \end{array} \right.$$ と決定されるものとする。 いま$F$を$-\varepsilon_i$の累積分布関数とする。すると上記より$Y_i^* = \beta_0 + \beta_1 X_i+\varepsilon_i$なので$Y_i=1$となるのは、$-\varepsilon_i \leq \beta_0 + \beta_1 X_i$のときであるから、$Y_1$となる確率は、この$F$を用いて $$P(Y_i=1|X_i) = F^*(X_i) = F(\beta_0 + \beta_1 X_i)$$ で与えられる。 プロビット・モデルでは$F^*$として標準正規分布が、ロジット・モデルでは、ロジスティック分布が上記のような形で使用される。 ![model](https://i.imgur.com/qA64efN.jpg) ### 8.2 ロジット・モデルとプロビット・モデル #### 8.2.1 プロビットモデル プロビットモデルでは、$\varepsilon_i$が標準正規分布に従うものとし、その累積分布関数 $$\phi(z) = \int_{-\infty}^{z} \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$$ を使い、$Y_i = 1$となる確率を、$F^*(X_i) = \phi(\beta_0 + \beta_1 X_i)$とするものである。 <!-- このとき、$Y^*_i = \beta_0 + \beta_1 X_i + \varepsilon_i$は、 --> #### 8.2.2 ロジットモデル ロジット・モデルでは、誤差項の確率分布がロジスティック分布であるとし、ロジスティック分布の累積分布関数 $$\Lambda(z) = \frac{e^x}{1+e^z}$$ を用いて、$Y_i = 1$となる確率を $$F^*(X_i) = \Lambda (\beta_0 + \beta_1 X_i)$$ とするものである。 ところで、$P_{0,i} = P(Y_i=0|X_i)$と$P_{1,i} = P(Y_i=1|X_i)$は条件付き確率の性質と、二値であることより $$1 = P(Y_i=0|X_i) + P(Y_i=1|X_i)$$が成立し、以上より $$ P(Y_i=0|X_i) = 1 - P(Y_i=1|X_i) = \frac{1}{1 + exp(\beta_0 + \beta_1 X_i)}$$ となる。さらに、 $$P_{1, i} = \frac{exp(\beta_0 + \beta_1 X_i)}{1+ exp(\beta_0 + \beta_1 X_i)}$$であることより、 $$\frac{P_{1,i}}{P_{0,i}} = exp(\beta_0 + \beta_1 X_i)$$ となり、$$\log{\frac{P_{1,i}}{P_{0,i}}} = \beta_0 + \beta_1 X_i$$ であることがわかる。 この結果は、ロジット・モデルは$Y_i$が0となるか1となるかの2つの確率の比(オッズ)の対数が$X_i$の線形関数であることを示す。 #### 8.2.3 モデルの推定方法 ある標本$Y_1,Y_2,...,Y_n$を得る同時確率は$Y_i=1$の場合$F(\beta_0 + \beta_1 X_i)$、$Y_i = 0$であれば、$1- F(\beta_0 + \beta_1 X_i)$を用いて考えるので、以下のようになる。 $$L( \beta_0,\beta_1) = \prod_{Y_i = 1} F(\beta_0 + \beta_1 X_i) \cdot \prod_{Y_i = 0} [1- F(\beta_0 + \beta_1 X_i)]$$ (しれっと最尤推定を使うことにしているが、本当はこれも吟味するべき(ちなみに、数理統計の本を読みましたが何にも分かりませんでした。...)) 最尤推定量$\hat \beta_0$、$\hat \beta_1$は$L(\beta_0, \beta_1)$を最大にするものであるが、まずは、最尤推定量を求める定石通り、両辺対数をとって $$\log {L( \beta_0,\beta_1) }= \sum_{i = 1}^{n} (Y_i)\log{F(\beta_0 + \beta_1 X_i)} + (1-Y_i) \log{(1- F(\beta_0 + \beta_1 X_i))}$$ とする。そして、 $$\frac{\partial \log{L(\beta_0,\beta_1)}}{\partial \beta_0} = 0, \frac{\partial \log{L(\beta_0,\beta_1)}}{\partial \beta_1} = 0$$ となる。$\beta_0, \beta_1$を求めるわけであるが、プロビット、ロジット・モデルのいずれにしても解析的に解を求めることは困難なので、ニュートン法を用いて、数値計算で求めることになる。このとき、$\log{L(\beta_0,\beta_1)}$が$\beta_0, \beta_1$の凸関数であるので、上記の関数が複数の解が存在することはなく、常に一意(計算誤差$\varepsilon$を除き)に定まる。 ### 8.3 省略 ### 8.4 説明変数が2個以上の場合 (武力行使の政治学がまさにこれ) $Y_i$が$p+1$個の変数$X_{1,i}, X_{2,i}, ... X_{p,i}$による場合に拡張ができる。(これは通常の回帰分析と同様) 式で表すと $$Y^*_i = \beta_0 + \beta_1 X_{1,i} + ... + \beta_p X_{p} + \varepsilon_i$$ となる。そしてこれまでと同様に$Y^*_i$の符号によって$Y_i$(=観測結果)が決定されるとする。 8.1のように$F$を$-\varepsilon_i$の累積分布関数とすると、$Y_i = 1$となる確率は、 $$F(\beta_0 + \beta_1 X_{1,i} + ... + \beta_p X_{p,i})$$ で与えられる。 尤度関数についても同様に与えられるので略