高校生のためのデータ分析入門 (23)：離散型アウトカムの回帰モデル（前編）

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回：高校生のためのデータ分析入門 (22)：稀なイベントの回数ならポアソン分布！ - ねこすたっと

おさらい：線形回帰モデルでは上手くいかないケース
2値変数をモデル化する
係数βはどう解釈したらいいの？
おわりに

おさらい：線形回帰モデルでは上手くいかないケース

以前、線形回帰モデルを説明しました。

necostat.hatenablog.jp

線形回帰モデルは、

観測値 = 期待値＋誤差
期待値は説明変数の関数
誤差は正規分布に従う

というモデルでした（線形性や誤差に関する他の仮定など、他にも大事な前提条件があるので、忘れていたら前の記事を振り返ってください）。

誤差が、正規分布という連続型確率分布に従うので、応答変数も連続変数でなければなりませんでしたね。次の例に挙げるような、誤差の分布が正規分布とみなせない離散型の応答変数の場合は、線形回帰モデルはうまく当てはまりません。

試験に合格するかどうか（2値変数*1）
学校の売店でたまにしか売れない商品の月販売個数（平均値が大きくない個数・回数）

2値変数をモデル化する

例えば、模擬試験の結果をもとに、3か月後の国家試験に合格するかどうかを予測したいとしましょう。この国家試験は得点開示がなく、合格・不合格しか分かりません。これまでに模擬試験を受けて、本番の結果を教えてくれた1000人のデータを使います。

説明変数X：模擬試験の点数
応答変数Y：1=合格, 0=不合格

期待値を説明変数の関数で表す

線形回帰モデルでは、期待値を説明変数の関数で表しました。今回も、

$\begin{aligned} E[Y_i|X_i] &= \beta_0 + \beta_1 X_i \end{aligned}$

としたいところですが、ここで一つ問題が...。

応答変数が2値変数の場合、期待値はY=1のカテゴリーを取る割合です。つまり、左辺の E[Y|X] は「模試でX点だった人が合格する割合」なので、0から1の範囲しか取りません。これに対して、右辺 $\beta_0 + \beta_1 X_i$ は、βやXの値によって-∞から+∞まであらゆる値を取る可能性があります。 0から1の範囲から外れた値は、合格割合として解釈できません（下図）。

もちろん、応答変数が連続変数であっても、データによっては期待値がYの取りうる範囲外になることはあります。下の図では、XもYも全て0以上の値ですが、回帰直線の一部は（Xの取りうる範囲内にも関わらず）0未満になっています。

ただ、このような状況は2値変数のときより遥かに少ないです（2値変数ではほぼ毎回起こる！）。

ロジット変換でつなぐ！

そこで、期待値と $\beta_0 + \beta_1 X$ を結びつける間に、一工夫はさみます。見やすさのため、E[Y|X]=pと書き直しました。

$\begin{aligned} \log \left( \frac{p_i}{1-p_i} \right) &= \beta_0 + \beta_1 X_i \end{aligned}$

左辺の変換、つまり確率を「あることが起こる確率を起こらない確率で割ったものの自然対数*2」にする変換をロジット変換（logit transformation）と呼びます。下は、 $f(p) = \log \left( \frac{p_i}{1-p_i} \right)$ （ロジット関数）のグラフです。確率pが0に近づくと、ロジットlog(p/1-p)はp=0の線に近付きながら-∞に向かい、pが1に近づくと、ロジットはp=1に漸近して+∞に向かいます。単調増加関数なので、ロジットが1つに決まれば、対応する確率も1つです。

ロジット変換を挟めば、 $\beta_0 + \beta_1 X$ がどんな値になっても、確率として適切な範囲にある値が1つに決まりますね（すばらしい！）

さっきのデータに当てはめると、下のグラフのようになります（確率の軸が縦方向なので、90°向きが変わっていることに注意）。

上の図のS字カーブは、ロジット関数の逆関数（pについて解いたもの）で、ロジスティック関数（logistic function）と呼びます。期待値の式を書き換えると、

$\begin{aligned} E[Y_i|X_i] &= \frac{e^{\beta_0 + \beta_1 X_i}}{e^{\beta_0 + \beta_1 X_i}+1} \\ &= \frac{1}{1+e^{-(\beta_0 + \beta_1 X_i)}} \\ \end{aligned}$

となります。

この例のように、「2値変数の期待値にロジスティック関数を当てはめるモデル」をロジスティック回帰モデル（logistic regression model）と呼びます。

誤差は？

「期待値」に関しては説明しましたが、「誤差」についてはどうでしょうか。ここで、二項分布の特殊形、ベルヌーイ分布の登場です。

例えば、模試の点数が60点だった人の期待値（=合格確率）が90%と推定されたとしましよまう。模試の点数が60点だった人のなかには、無事合格する人もいれば、残念な結果になる人もいます。模試の点数が同じなのに、結果が違うのは全くの運です。あたかも神様が、当りが90%のクジを引いて、受験者の合否を決めているみたいですね。

線形回帰モデルでは、期待値と実際に観察された値の違いは正規分布で説明されました。ロジスティック回帰モデルでは、その違いを二項分布（ベルヌーイ分布）で説明しています。

ロジスティック回帰モデルのまとめ

まず、ロジスティック回帰モデルを式で表してみます。

$\begin{aligned} \log \left( \frac{p_i}{1-p_i} \right) &= \beta_0 + \beta_1 X_i \\ Y_i &\sim Binom(1, p_i) \end{aligned}$

これを日本語で噛み砕くと、

説明変数 $X_i$ の線形関数として確率 $p_i$ が決まる。つまり、それぞれの人は $X_i$ に応じて $p_i$ が与えられる。
説明変数と確率を直接はつなげないので、ロジット変換を介してつなぐ。
応答変数 $Y_i$ は、確率 $p_i$ をパラメータとするベルヌーイ分布に従う。

図で表すとこんな感じでしょうか。

係数βはどう解釈したらいいの？

線形回帰モデルでは、ある説明変数の係数βは「その説明変数が1単位変化したときの期待値の変化量」として解釈できました。

necostat.hatenablog.jp

ロジスティック回帰モデルではどうでしょうか。線形回帰モデルのときと同じように、E[Y|X=x] と E[Y|X=x+1] の差を計算してみましょう。読みやすくなるように、E[Y|X=x+1] = P1, E[Y|X=x] = P0と書くことにします。

$\begin{aligned} P_1 &= \frac{1}{1+e^{-(\beta_0 + \beta_1 (x+1))}} \\ P_0 &= \frac{1}{1+e^{-(\beta_0 + \beta_1 x)}} \\ \end{aligned}$

数式は整理できましたか？全然キレイな形にできませんね...。

ロジスティック回帰モデルでは、係数βと期待値の変化量を結びつけることが簡単ではありません。代わりに、下のように期待値のロジットを計算してみましょう。

$\begin{aligned} \log \left( \frac{P_1}{1-P_1} \right) &= \beta_0 + \beta_1 (x+1) \\ \log \left( \frac{P_0}{1-P_0} \right) &= \beta_0 + \beta_1 x \\ \end{aligned}$

第1式から第2式を引くと、下のようになります。

$\begin{aligned} \log \left( \frac{P_1}{1-P_1} \right) - \log \left( \frac{P_0}{1-P_0} \right)&= \beta_1 \\ \end{aligned}$

右辺はすっきりしましたね。左辺は、log A-log B= log(A/B) を使って変形します。

$\begin{aligned} \log \left( \frac{P_1/(1-P_1)}{P_0/(1-P_0)} \right) &= \beta_1 \\ \frac{P_1/(1-P_1)}{P_0/(1-P_0)} &= e^{\beta_1} \end{aligned}$

さて、ここで左辺の分母と分子に登場した P/(1-P)は「あることが起こる確率を、起こらない確率で割ったもの」で、オッズ（odds）と呼ばれます。オッズ起こる確率Pが小さいとき（目安：<10%）は、起こる確率Pとオッズは近い値になります*3。

左辺は「X=x+1のときにY=1が起こる確率」と「X=xのときにY=1が起こる確率」の比で、オッズ比（odds ratio）と呼ばれます。オッズ比が1であれば、2つの発生確率は同じであり、オッズ比が1から離れるほど発生確率に違いがある、ということになります。

係数βを指数変換したもの $e^{\beta}$ を使えば、「説明変数が1単位変化したときに、オッズが何倍になるか」、すごく平たく言えば「何倍起こりやすくなるか*4」を推定できるわけです。

おわりに

線形回帰モデルで係数を推定するために使った「最小二乗法」は使えません。仮定したモデルでデータを説明するときに、最も尤もらしい*5 パラメータを求める方法（最尤法）を使います。
ロジスティック回帰モデルでも、モデルの仮定に無理がないか確認する方法はあります。
次回：高校生のためのデータ分析入門 (24)：離散型アウトカムの回帰モデル（後編） - ねこすたっと