統計学 学習ノート

はじめに

直線を当てはめると、確率が1を超えてしまう

「勉強時間で合格率を予測したら、確率が120%になったんだけど？」

例えば、勉強時間 $x$ と合否のデータに、これまでの回帰直線を当てはめたとしましょう。合格を1、不合格を0として $\hat{y} = 0.1x$ という式が得られたとします。12時間勉強した人の予測値は1.2です。確率が1を超えてしまいました。逆に、切片が負なら確率がマイナスになることもあります。

事故の件数でも同じことが起きます。件数は0以上の整数なのに、直線はマイナスの件数を平気で予測します。さらに、件数のばらつきは平均が大きいほど大きくなるので、誤差の分散が一定という仮定も成り立ちません。

実は、応答の種類に合わせて線形モデルを少しだけ改造すれば、この問題は解決します。

そこで登場するのが、今回の主役 一般化線形モデル（generalized linear model、GLM）です。一般化線形モデルとは、説明変数の線形和はそのまま使い、応答の分布と、線形和との結びつけ方を取り替えられるようにしたモデルの枠組みです。

前回の記事では、回帰診断法でモデルを点検しました。今回の記事では、二値の応答を扱うロジスティック回帰モデルとプロビットモデル、カウントの応答を扱うポアソン回帰モデルを学び、最後に3つを一般化線形モデルとしてまとめて解説します。

ロジスティック回帰モデル：確率をオッズの対数で表す

確率をロジットへ両替すると、線形和の世界で直線になる

はじめに、最もよく使われるロジスティック回帰モデルについて解説します。

確率 $p$ は0から1の間に閉じ込められています。一方で、線形和 $\beta_0 + \beta_1 x$ は、どんな値にもなります。住む世界が違うので、そのままでは等号で結べません。

そこで、海外旅行の両替を思い浮かべてください。日本円のままでは現地で買い物ができないので、両替所で現地のお金に替えますよね。同じように、確率を線形和の世界へ両替します。

両替は2段階です。まず、確率を オッズ $p/(1-p)$ に替えます。オッズは「起こる確率が起こらない確率の何倍か」を表し、値の範囲は0から無限大です。次に対数をとると、範囲はマイナス無限大から無限大になります。これが ロジット です。

ロジスティック回帰モデル は、このロジットを線形和で表すモデルです。

$$\log\frac{p}{1-p} = \beta_0 + \beta_1 x \quad\Longleftrightarrow\quad p = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x)}}$$

$p$ は $y = 1$ となる確率、$x$ は説明変数、$\beta_0, \beta_1$ は係数です。右の式は $x$ に対してS字の曲線を描き、0と1の間からはみ出しません。

具体的に見てみましょう。合格のモデルが $\log\{p/(1-p)\} = -4 + 0.5x$ だったとします。

勉強時間 $x$	4	8	10	12
ロジット	$-2$	0	1	2
合格確率 $p$	0.119	0.500	0.731	0.881

12時間でも確率は0.881で、1を超えません。

ここがポイントです。係数の解釈は「確率」ではなく「オッズ」で行います。$x$ が1増えるとロジットは $\beta_1$ 増えるので、オッズは $e^{\beta_1}$ 倍になります。この倍率が オッズ比 です。上の例では $e^{0.5} \approx 1.65$ で、「勉強時間が1時間増えると、合格のオッズが約1.65倍になる」と読みます。確率が1.65倍になるのではありません。実際、$x$ が8から9に増えるときと、12から13に増えるときでは、確率の増え方は違います。

推定には、前に学んだ最尤法を使います。応答 $y_i$ は0か1なので、ベルヌーイ分布の尤度から、対数尤度は次の式になります。

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \Big\{ y_i \log p_i + (1 - y_i)\log(1 - p_i) \Big\} = \sum_{i=1}^{n} \Big\{ y_i \eta_i - \log(1 + e^{\eta_i}) \Big\}$$

ここで $\eta_i = \boldsymbol{x}_i^\top \boldsymbol{\beta}$ です。これを $\boldsymbol{\beta}$ で微分して0とおくと、次の 尤度方程式 が得られます。

$$\sum_{i=1}^{n} (y_i - p_i)\,\boldsymbol{x}_i = \boldsymbol{0}$$

言葉で言うと、「実測と予測の差が、説明変数に対して打ち消し合う」という条件です。形は最小二乗法の正規方程式に似ていますが、$p_i$ が $\boldsymbol{\beta}$ の非線形な関数なので、解を式では書けません。そこで、ニュートン法のような反復計算で解きます。係数の標準誤差はフィッシャー情報行列 $X^\top W X$（$W$ は $p_i(1-p_i)$ を並べた対角行列）の逆行列から求まり、$\hat{\beta}/\text{SE}$ を標準正規分布と比べるのがワルド型検定です。

例題：2群の比較をロジスティック回帰で解く

新しい薬を40人に使うと30人が回復し、対照の40人では15人が回復しました。薬の有無を $x$（薬あり1、なし0）として、ロジスティック回帰の係数とオッズ比の95%信頼区間を求めましょう。

【解き方】

1. 説明変数が0と1だけなので、最尤推定値は群ごとの比率に一致する
2. 対照群のロジット： $\hat{\beta}_0 = \log(15/25) \approx -0.511$
3. 薬の群のオッズは $30/10 = 3$、対照群は $15/25 = 0.6$ なので、オッズ比は $3 \div 0.6 = 5$、$\hat{\beta}_1 = \log 5 \approx 1.609$
4. 標準誤差： $\sqrt{1/30 + 1/10 + 1/15 + 1/25} = \sqrt{0.24} \approx 0.490$
5. ワルド統計量： $1.609 \div 0.490 \approx 3.29 > 1.96$ なので有意
6. 信頼区間： $\exp(1.609 \pm 1.96 \times 0.490)$ より、約1.91から約13.1

答えはオッズ比5、95%信頼区間は約1.91から約13.1です。区間は係数の世界で作ってから指数をとるので、オッズ比の5に対して左右対称になりません。

試験のポイント：ロジスティック回帰の係数 $\beta$ は対数オッズ比です。$e^{\beta}$ がオッズ比で、$\beta$ そのものを確率の増分と読むのは誤りです。

以上が、ロジスティック回帰モデルについての解説です。

プロビットモデル：正規分布で確率を表す

見えない水位が堤防を超えた確率が、プロビットの確率

次に、もう1つの二値モデルについて説明します。

両替所は1か所ではありません。確率を線形和の世界へ替える方法は、他にもあります。プロビットモデル は、標準正規分布の分布関数 $\Phi$ を使った両替です。

$$p = \Phi(\beta_0 + \beta_1 x) \quad\Longleftrightarrow\quad \Phi^{-1}(p) = \beta_0 + \beta_1 x$$

なぜ正規分布なのでしょうか？ ここで、堤防を思い浮かべてください。川の水位は見えませんが、堤防を超えた瞬間に「あふれた」という結果だけが分かります。

二値の応答も同じです。目に見えない連続の量 $y^*$ を考えます。

$$y^* = \beta_0 + \beta_1 x + \varepsilon, \qquad y = \begin{cases} 1 & (y^* > 0) \\ 0 & (y^* \le 0) \end{cases}$$

これを 潜在変数 と呼びます。$\varepsilon$ が標準正規分布なら、$P(y = 1) = P(\varepsilon > -\eta) = \Phi(\eta)$ となり、プロビットモデルが出てきます。$\varepsilon$ がロジスティック分布なら、ロジスティック回帰モデルになります。つまり、2つのモデルの違いは誤差の分布だけです。

例えば、殺虫剤の濃度と虫の死亡を考えます。虫1匹ごとに耐えられる限界の濃度があり、それが正規分布に従うとすれば、プロビットモデルが自然です。

具体的に見てみましょう。$\eta = -1 + 0.5x$ のとき、$x = 3$ なら $\eta = 0.5$ で、$p = \Phi(0.5) \approx 0.691$ です。$x = 2$ なら $\eta = 0$ で、$p = 0.5$ です。

2つのモデルは、見た目がよく似ています。標準ロジスティック分布の標準偏差は $\pi/\sqrt{3} \approx 1.81$ で、標準正規分布の1より大きいです。そのため、同じデータでは、ロジスティック回帰の係数はプロビットの係数の1.6から1.8倍ほどになります。

$z$	0.5	1.0	1.5
プロビット $\Phi(z)$	0.691	0.841	0.933
ロジスティック（係数を1.7倍）	0.701	0.846	0.928

ほとんど同じですね。違いが出るのは、確率が0や1に近い裾の部分だけです。実務では、オッズ比で解釈できるロジスティック回帰がよく使われます。

なお、説明変数の効き目は、どちらのモデルでも $x$ の位置で変わります。$x$ を1増やしたときの確率の変化率は、プロビットでは $\phi(\eta)\beta_1$、ロジスティックでは $p(1-p)\beta_1$ です。$\phi$ は標準正規分布の密度です。

以上が、プロビットモデルについての説明です。

ポアソン回帰モデル：件数を対数で表す

線形和が1増えるごとに、平均が同じ倍率で増える

続いて、カウントの応答を扱うモデルについて解説します。

事故の件数や来店客の数は、0, 1, 2, … の整数で、平均 $\mu$ は正の値です。そこで、平均の対数を線形和で表します。これが ポアソン回帰モデル です。

$$Y_i \sim \text{Po}(\mu_i), \qquad \log \mu_i = \beta_0 + \beta_1 x_i$$

対数の世界で足し算をすることは、元の世界で掛け算をすることです。貯金の複利のように、$x$ が1増えるごとに平均が同じ倍率で増えます。その倍率が $e^{\beta_1}$ で、率比 と呼ばれます。

例えば、交差点の1年間の事故件数が $\log \mu = 0.2 + 0.3x$ だったとします。$x$ は1日の交通量（万台）です。$x = 2$ なら $\mu = e^{0.8} \approx 2.23$ 件、$x = 3$ なら $\mu = e^{1.1} \approx 3.00$ 件です。交通量が1万台増えるごとに、事故は $e^{0.3} \approx 1.35$ 倍になります。

現実には、観測の期間や人数がデータごとに違うことがあります。5年分の件数と1年分の件数をそのまま比べるのは不公平ですよね。そこで、期間 $t_i$ を線形和の外から足します。

$$\log \mu_i = \log t_i + \beta_0 + \beta_1 x_i$$

$\log t_i$ は係数を推定しない項で、オフセット と呼びます。両辺を整理すると $\log(\mu_i / t_i) = \beta_0 + \beta_1 x_i$ で、1年あたりの率をモデル化していることになります。例えば、上の交差点を5年間観測すると、$x = 2$ の期待件数は $2.23 \times 5 \approx 11.1$ 件です。

ただし、弱点もあります。ポアソン分布は平均と分散が等しい分布でした。実際のデータでは、分散のほうが平均より大きいことがよくあります。これを 過分散 と呼びます。放っておくと標準誤差が小さく出すぎて、効果を過大評価します。対処には、分散を定数倍する準ポアソンモデルや、負の二項分布を使います。

例題：2つの地区の事故率を比べる

地区Aでは10年間に30件、地区Bでは12年間に18件の事故が起きました。オフセットつきのポアソン回帰で、地区Aの率比と検定統計量を求めましょう。

【解き方】

1. 年あたりの率： Aは $30/10 = 3$、Bは $18/12 = 1.5$
2. 率比： $3 \div 1.5 = 2$ なので、$\hat{\beta}_1 = \log 2 \approx 0.693$
3. 標準誤差： $\sqrt{1/30 + 1/18} \approx 0.298$
4. ワルド統計量： $0.693 \div 0.298 \approx 2.33 > 1.96$

答えは率比2で、有意水準5%で差があると判断します。件数の逆数の和で標準誤差が書ける点は、さきほどのオッズ比と同じ形です。

以上が、ポアソン回帰モデルについての解説です。

一般化線形モデル：3つの部品で全部をまとめる

分布とリンク関数を取り替えれば、同じ枠でモデルが作れる

最後に、ここまでのモデルを1つの枠組みに整理していきましょう。

パソコンの組み立てを想像してください。部品を取り替えれば、同じ枠で事務用にもゲーム用にもなります。一般化線形モデルも、3つの部品の組み合わせでできています。

【一般化線形モデルの3つの部品】

1. 確率分布： 応答 $Y$ が従う分布。指数型分布族から選ぶ
2. 線形予測子： $\eta = \boldsymbol{x}^\top \boldsymbol{\beta}$
3. リンク関数： 平均 $\mu$ と線形予測子を結ぶ単調な関数 $g(\mu) = \eta$

リンク関数が、ここまで「両替所」と呼んできたものです。部品の組み合わせは、次のとおりです。

モデル	分布	リンク関数 $g(\mu)$	分散関数 $V(\mu)$
線形回帰	正規分布	恒等 $\mu$	1
ロジスティック回帰	二項分布	ロジット $\log\{\mu/(1-\mu)\}$	$\mu(1-\mu)$
プロビットモデル	二項分布	プロビット $\Phi^{-1}(\mu)$	$\mu(1-\mu)$
ポアソン回帰	ポアソン分布	対数 $\log \mu$	$\mu$

ふつうの線形回帰も、この枠組みの特別な場合です。

では、なぜこの分布たちは同じ枠に入るのでしょうか？ どれも 指数型分布族 という共通の形に書けるからです。

$$f(y; \theta, \phi) = \exp\left\{ \frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi) \right\}$$

$\theta$ は 自然パラメータ、$\phi$ は分散の大きさを決める母数です。この形から、次の性質が出ます。

$$E[Y] = b'(\theta), \qquad V[Y] = b''(\theta)\,a(\phi)$$

ポアソン分布で確かめてみましょう。確率関数は $e^{-\mu}\mu^y/y! = \exp(y\log\mu - \mu - \log y!)$ なので、$\theta = \log\mu$、$b(\theta) = e^{\theta}$、$a(\phi) = 1$ です。すると、$b'(\theta) = b''(\theta) = e^\theta = \mu$ で、平均と分散が等しいという性質が出てきます。同じように、ベルヌーイ分布では $\theta = \log\{p/(1-p)\}$、$b(\theta) = \log(1 + e^\theta)$ となります。

$\theta$ を線形予測子に等しくするリンク関数を、標準リンク関数 と呼びます。ロジットと対数は、自然パラメータそのものだから選ばれていたのです。プロビットは標準リンクではありません。

モデルの当てはまりは、逸脱度 で測ります。データを完全に再現する飽和モデルの対数尤度を $\ell_{\text{sat}}$ として、次のように定義します。

$$D = 2\big(\ell_{\text{sat}} - \ell(\hat{\boldsymbol{\beta}})\big)$$

正規分布では、逸脱度は残差平方和を $\sigma^2$ で割った量になります。つまり、逸脱度は残差平方和の親戚です。入れ子の2つのモデルの逸脱度の差は尤度比検定の統計量で、帰無仮説のもとで近似的に、増やした係数の個数を自由度とするカイ二乗分布に従います。

例題：逸脱度でモデルを比べる

ロジスティック回帰で、説明変数を2個増やしたところ、逸脱度が65.3から52.9に下がりました。増やした2個の係数がともに0という帰無仮説を、有意水準5%で検定しましょう。

【解き方】

1. 検定統計量： $65.3 - 52.9 = 12.4$
2. 自由度： 増やした係数の個数なので2
3. 棄却限界： $\chi^2_{0.05}(2) = 5.99$
4. 判定： $12.4 > 5.99$ なので棄却する（p値は約0.002）

答えは、増やした変数は有意に当てはまりを改善する、です。

試験のポイント：GLM は「分布・線形予測子・リンク関数」の3点セットです。逸脱度の差の自由度は、増やした係数の個数です。

以上が、一般化線形モデルについての解説です。

細かい指数型分布族の式を忘れてしまっても問題ありません。ここでは、「応答に合った分布を選び、リンク関数で線形和と結びつける」という点だけ押さえておきましょう。

今回のまとめ

二値も件数も、部品を替えればGLMの枠に入る

それでは、今回の内容を振り返ります。

今回は、二値やカウントの応答を扱う回帰について紹介しました。直線がはみ出すなら、応答の世界を線形和の世界へ両替すればよい。一般化線形モデルの発想は、この1つです。

今回のポイントを振り返りましょう。

振り返り

• ロジスティック回帰は、確率をオッズの対数へ両替する。$e^{\beta}$ はオッズ比で、確率の倍率ではない。
• プロビットは正規分布の分布関数で両替する。堤防を超える潜在変数の見方で、ロジスティックと結びつく。
• ポアソン回帰は平均の対数を線形和で表す。$e^{\beta}$ は率比で、期間の違いはオフセットで調整する。
• 分散が平均より大きい過分散では、標準誤差が小さく出すぎる。
• 一般化線形モデルは、分布・線形予測子・リンク関数の3つの部品で、これらをまとめて扱う。

キーワードのおさらい

• ロジスティック回帰モデル: 確率のロジットを線形和で表す二値応答のモデル
• プロビットモデル: 確率を標準正規分布の分布関数で線形和と結ぶ二値応答のモデル
• ポアソン回帰モデル: カウントの平均の対数を線形和で表すモデル
• 一般化線形モデル: 指数型分布族の分布と、リンク関数と、線形予測子からなるモデルの枠組み
• オッズ比: オッズの倍率。ロジスティック回帰では $e^{\beta}$
• 逸脱度: 飽和モデルとの対数尤度の差の2倍。当てはまりの悪さを測る

さて、応答の種類はこれでだいぶ広がりました。でも、「いつ壊れるか」のような時間の応答には、別の難しさがあります。調査が終わっても、まだ壊れていないものがあるからです。

次回は、途中までしか観測できないデータを扱う「生存時間解析」と、非線形なモデルである「ニューラルネットワークモデル」について紹介します。

それでは、次の記事へ進んでいきましょう。