ロジスティック回帰モデルとは?

2022.11.16

要約

　ロジスティック回帰モデルとは、目的変数Yが0か1を取るモデルの一つです。説明変数Xの値を変えると、Y＝1になる確率が定まります。あくまで確率が定まるだけなので、Yが0か1は不確定です。

　例えば、曲線がY=1を取る確率、点が観測された現象とすると、次のようなロジスティック回帰モデルがありえます。

全体像

（1）問題の構造

　「ロジスティック回帰モデルとは何か？」という問いを設定します。この問いを

①ロジスティック回帰モデルの定義は何か

②ロジスティック回帰モデルの定義はどう解釈できるか

③ロジスティック回帰モデルをどう結果の予測に使えるか

④ロジスティック回帰モデルをどう原因の説明に使えるか

に分解します。

（2）前提の選択

　ロジスティック回帰モデルの定義（①）については、前提に含めてしまいます。

（3）論点の選択

　定義の解釈（②）については「モデル式をわかりやすく書き換えるには？」を、結果の予測（③）については「条件付き期待値はいくらか？」を、原因の説明（④）については「限界効果はいくらか？」を論点とします。

$$条件付き期待値E(Y|X_1,X_2 \cdots X_k)$$

$$限界効果\frac{\partial E(Y|X_1,X_2 \cdots X_k)}{\partial X_1}$$

$$Y：目的変数、X：説明変数、k：説明変数の種類数$$

（4）付録一覧

　冗長さを避けて、可読性を上げるために、以下の内容は付録に回します。

「潜在変数による定義は何か？ 」「ロジット関数による定義は何か？」「オッズ比はいくらか？」「平均的限界効果の定義は何か？」「限界効果をどう計算したか？ 」「Rコード」

前提

　ロジスティック回帰モデルとは、目的変数Yが0か1かをとる回帰モデルです。

　Y=1になる確率をpとすると、

$$p=P(Y=1|X_1,X_2 \cdots X_k)$$

$$=\frac{e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}{1+e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}$$

$$=\frac{1}{1+e^{-(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}$$

となるのをロジスティック回帰モデルと定義します。

　例えば、β0=0.1、β1=0.2、k=1で、さきほどの発生させたデータを用いると、pと観察データは次のように図示できます。

結果

（1）モデル式の解釈

　ロジスティック回帰モデルの式は、次のように作ることが出来ます。

　第一に、Y=0の確率を1としたときの確率比率をP*と定義します。Y=0は基準カテゴリ(reference category)です。

$$P^*(Y=0)=1$$

　第二に、Y=1の確率比率を次のように定義します。

$$P^*(Y=1)=e ^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}$$

　第三に、確率比率 P*を確率Pに変換します。Y=0とY=1は互いに排反であることに注目すると

$$P(Y=0または1)=P(Y=0)+P(Y=1)$$

$$P(Y=1)=\frac{P^*(Y=1)}{P^*(Y=0)+P^*(Y=1)}$$

$$=\frac{e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}{1+e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}$$

となります。ロジスティック回帰モデルができました。

（2）Yの条件付き期待値

　説明変数Xが決まったときの目的変数Yの期待値は

$$E(Y|X_1,X_2 \cdots X_k)$$

$$=\frac{e^{\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k}}{1+e^{\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k}}$$

です。なぜなら、Yは1か0しかとらないので

$$E(Y|X_1,X_2 \cdots X_k)$$

$$=1 \times P(Y=1) + 0 \times P(Y=0)$$

$$= 1\times p +0 \times (1-p)$$

$$=\frac{e^{\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k}}{1+e^{\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k}}$$

だからです。

（3）Xの限界効果

数式

　X1が微小1単位大きいと、Yが平均的に大きい量を、限界効果といいます。ロジスティック回帰モデルの限界効果は

$$限界効果＝\frac{\partial P(Y=1|X1,X_2 \cdots X_k)}{\partial X_1}$$

$$=\frac{e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k }}{ \Big(1+e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k } \Big)^2}\beta_1$$

です。理由は付録5に載せました。

限界効果は変化する

　限界効果の数式からわかるのは、説明変数Xの値によって

・プラスとマイナスは、変わらない

・絶対値は、変わる

ということです。限界効果は、β1だけでは決まりません。ここが重回帰モデルと異なるところです。

　例えば、β0=0.1、β1=0.2、k=1の場合の限界効果は、下図です。

※レポートでは、平均的限界効果の数値を報告します。平均的限界効果の定義には、2つあります。詳しくは、記事最後の付録4へ。

考察

（1）結論

　まとめると、ロジスティック回帰モデルとは、次のような回帰モデルです。

定義

　Y=1になる確率をpとすると、

$$p=P(Y=1|X_1,X_2 \cdots X_k)$$

$$=\frac{e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}{1+e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}$$

モデル式の解釈

$$確率比率P^*(Y=0)=1とすると$$

$$P^*(Y=1)=e ^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}$$

である回帰モデルと解釈できます。

結果の予測

　Xが決まったときにYの期待値は、以下のとおりです。

$$E(Y|X_1,X_2 \cdots X_k)$$

$$=\frac{e^{\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k}}{1+e^{\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k}}$$

原因の説明

　X1がYに与える影響は、次で表せます。

$$限界効果＝\frac{\partial P(Y=1|X1,X_2 \cdots X_k)}{\partial X_1}$$

$$=\frac{e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k }}{ \Big(1+e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k } \Big)^2}\beta_1$$

（2）妥当性評価

　前提が合っていれば、以上の内容は正しいでしょう。

（3）意義

　2つの事象がランダムに現れる現象をモデルに出来ます。まず

Y=1：テストに合格する

Y=0：テストに不合格になる

のように変数を定義してやると、ロジスティック回帰モデルに変換できます。すると、次のようなデータ形式についてロジスティック回帰分析を用いて、パラメーターβを推定できます。

合格=1
不合格=0		1週間勉強時間テスト前勉強時間
Aさん115
Bさん12.510
Cさん003
・・・・・・・・・・・・
Hさん1010
Iさん001

　他にも

・シュートが決まる/シュートが決まらない

・残業する/残業しない

・プロジェクトが成功する/プロジェクトが失敗する

といった場合のモデル化ができます。

付録1：潜在変数による定義は何か？

　ロジスティック回帰モデルを潜在変数モデルで捉えることもできます。

　Yが0か1のどちらを取るのかは、次の潜在変数で決定されます。

　潜在変数は↓で決まります。

$$潜在変数Y^*=\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k + U$$

$$Y^*：潜在変数、X：説明変数$$

$$\beta：母回帰係数、U：誤差項$$

$$U|X_1,X_2 \cdots X_k はロジスティック分布に従う$$

　潜在変数が決まると、次のように目的変数が決定されます。

$$潜在変数Y^* >0 のとき　目的変数Y=1$$

$$潜在変数Y^* ≦0 のとき　目的変数Y=0$$

　説明変数X1,X2・・・Xkが決まると、潜在変数Y*がおおまかに決まります。最終的には、ロジスティック分布にしたがってランダムに決定される誤差項Uによって、目的変数が1になるか0になるか決まります。

　例えば、潜在変数が

$$潜在変数Y^*=0.1+0.2 X + U$$

$$U|Xはロジスティック分布に従う$$

で決まる場合を考えてみましょう。Rにて乱数を生成すると、下図のように潜在変数Y*が決まりました。（付録4にRコード有）

　この場合、潜在変数Y*と0との関係で、↓のように目的変数が1か0に振り分けられます。

付録2：ロジット関数による定義は何か？

　ロジスティック回帰モデルには、ロジット関数を用いて線形回帰モデルにした定義式もあります。

　ロジット関数を次のように定義すると

$$logit(p)=\log_e \left( \frac{p}{1-p} \right)$$

　ロジスティック回帰モデルは

$$logit(p)=\beta_0+\beta_1 X_1 +\beta_2 X_2 +\cdots + \beta_k X_k$$

と定義できます。なお、ロジット関数は対数オッズともみなせます。オッズについては付録をご覧ください。

$$logit(p)=\log_e \left( \frac{p}{1-p} \right)=\log_e \left( オッズ \right)$$

付録3：オッズ比はいくらか？

　ロジスティック回帰モデルでは、限界効果は変化します。しかし、限界効果に関係する概念であるオッズ比は一定です。そこで、ロジスティック回帰ではオッズ比を報告する場合が多いです。

結論

　ロジスティック回帰モデルにおけるX1のオッズ比は

$$オッズ比=e^{\beta_1}$$

です。

前提1：オッズとは？

　オッズとは、「Y=1が起きる確率をp」としたときに

$$オッズ＝\frac{起きる確率}{起きない確率}=\frac{p}{1-p}$$

となる数値のことです。

前提2：オッズ比とは

　オッズ比とは、確率pのオッズと確率qのオッズの比

$$オッズ比=\Bigg(\frac{p}{1-p} \Bigg) \Bigg/ \Bigg(\frac{q}{1-q} \Bigg) $$

です。ここで

$$p=P(Y=1|X_1=a,X_2 \cdots X_k)$$

$$q=P(Y=1|X_1=a-1,X_2 \cdots X_k)$$

と定義すると、オッズ比はX1が1単位増えたときに、Y=1になる確率がどれくらい増えるかに関わる概念になります。

前提3：ロジスティック回帰モデルにおけるp

　定義よりロジスティック回帰モデルにおけるpは

$$p=\frac{e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}{1+e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}$$

です。

結果

　ロジスティック回帰モデルの場合、オッズは

$$オッズ\frac{p}{1-p}$$

$$=\frac{ \Bigg(\frac{e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}{1+e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}} \Bigg)} { \Bigg(1-\frac{e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}{1+e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}} \Bigg)}$$

$$=e^{(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}$$

$$=e^{\beta_1 a +他}　←X_1=a$$

となります。つまり、オッズ比は、X1以外の変数が不変とすると

$$オッズ比＝\Bigg(e^{\beta_1 a+他} \Bigg) \Bigg/ \Bigg(e^{\beta_1 (a-1) +他}\Bigg) $$

$$=e^{\beta_1}$$

　これはX1の値によって変化しないです。

　したがって、ロジスティック回帰モデルでは、オッズ比が一定かつ回帰係数β1から簡単に計算できます。

付録4：平均的限界効果の定義は何か？

　平均的限界効果の定義には、2つあります。サンプル・サイズがnだとすると

1．個々の説明変数を用いて限界効果を計算→全ての限界効果を平均を計算

$$\frac{1}{n} \sum_{i=1}^n \frac{e^{ \beta_0+\beta_1 X_{1i}+\beta_2 X_{2i}+\cdots + \beta_k X_{ki} }}{ \Big(1+e^{ \beta_0+\beta_1 X_{1i}+\beta_2 X_{2i}+\cdots + \beta_k X_{ki} } \Big)^2}\beta_1$$

$$ただし、X_{1i}はデータ番号iの説明変数X_1の数値$$

2．説明変数の平均を計算→平均的個人の限界効果を計算

$$\frac{e^{ \beta_0+\beta_1 \overline{X_1}+\beta_2 \overline{X_2}+\cdots + \beta_k \overline{X_k} }}{ \Big(1+e^{ \beta_0+\beta_1 \overline{X_1}+\beta_2 \overline{X_2}+\cdots + \beta_k \overline{X_k} } \Big)^2}\beta_1$$

$$ただし、\overline{X_1}は平均　\frac{1}{n} \sum_{i=1}^n X_{1i}$$

の二つが平均的限界効果の定義です。

付録5：限界効果をどう計算したか？

結論

$$限界効果=\frac{e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k }}{ \Big(1+e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k } \Big)^2}\beta_1$$

前提：合成関数の微分

$$y=F(x)をxで微分したい$$

が関数Fが複雑で直接的には求められない場合、合成関数の微分を使うとうまくいくときがあります。

　合成関数の微分とは

$$y=f(u)　u=g(x)　つまりF=f(g(x))のとき$$

$$\frac{ \partial y}{\partial x}=\frac{ \partial y}{ \partial u} \frac{ \partial u}{ \partial x} $$

を利用して、微分を繰り返し、yをxで微分した結果を求めることです。

結果

　E、A、Bを次のように

$$B=-(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)$$

$$A=1+e^{-(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}=1+e^B$$

$$E=\frac{1}{1+e^{-(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}=\frac{1}{A}$$

おいた上で、合成関数の微分

$$限界効果=\frac{ \partial E}{\partial X_1}$$

$$=\frac{ \partial E}{ \partial A} \frac{ \partial A}{ \partial B} \frac{ \partial B}{ \partial X_1}$$

$$= 　(- A^{-2} )\times e^B \times (- \beta_1)$$

$$=\frac{e^{ – ( \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k ) }}{ \Big(1+e^{ -( \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k )} \Big)^2}\beta_1$$

　分子と分母に、分子の2乗をかけると

$$=\frac{e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k }}{ \Big(1+e^{ \beta_0+\beta_1 X_1+\beta_2 X_2+\cdots + \beta_k X_k } \Big)^2}\beta_1$$

付録6：Rコード

作図には、R言語を用いて行いました。R言語については「しまうまのRでデータ分析入門」をご覧ください。

　美しい図示のため、パッケージとしてggplot2を用いました。

#パッケージ
library(ggplot2)

　ちなみに、Rでロジスティック回帰分析をするコード

#ロジスティック回帰分析のコード
glm(y~x1 + x2 + x3,family = binomial("logit"))
　
#lmは最尤法による一般化線形モデルの推定。yが目的変数、xが説明変数、data1が使用するデータの名前。ロジスティック回帰にはbinomial("logit")が必要です

　入力①↓

#確率密度関数を出力。横軸は-10から10
ggplot(data.frame(x = c(-10, 10)), aes(x))+
  stat_function(fun=dlogis)+ #ロジスティック分布
  stat_function(fun=dlogis)　#標準正規分布

　出力①↓

　②と③のための関数作成とデータ作成↓

#確率pの関数作成
p <- function(x) (exp(1)^(0.1+0.2*x))/(1+(exp(1)^(0.1+0.2*x)))
　
#サンプル数sample_numberを決めます。今回は100です。
sample_number <- 100
　
#sample_number個の正規分布に従う乱数を発生。平均mean=0.5、標準偏差sd=5
x <- rnorm(n=sample_number, mean =-0.5, sd=5)
　
#sample_number個のロジスティック分布に従う乱数を発生
u <- rlogis(sample_number)
　
#  yy = -0.1 + 0.2*x+ u　に従う潜在変数yyを作る
yy <- 0.1 + 0.2*x+ u
　
# 0を以上なら1、未満なら0の二値変数yを作る
y <- ifelse(yy>=0,1,0)
　
#データフレーム化
data_L <- data.frame(x,y,yy)

　出力②↓

ggplot(data_L,aes(x =x, y=yy))+
  geom_point()+
  scale_x_continuous(limits = c(-20, 20))+
  scale_y_continuous(limits = c(-5, 7))

　出力②↓

　入力③↓

ggplot(data_L,aes(x =x, y=y))+
  geom_point()+
  stat_function(fun=p)+
  scale_x_continuous(limits = c(-20, 20))+
  scale_y_continuous(limits = c(0, 1))

　出力③↓