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

 ロジスティック回帰モデルとは、なんでしょうか? 定義から使用法まで、総合的に解説します。

要約

(1)問い

 ロジスティック回帰モデルとは、なんでしょうか?

(2)答え

 ロジスティック回帰モデルとは、目的変数Yが0か1かを取る場合についての回帰モデルの一種です。

(3)意義

 統計学、機械学習で使われている手法を知れます。応用範囲は、医学、心理学、マーケティングなど多岐にわたります。

前提

(1)定義

 ロジスティック回帰モデルには3つの定義があります。潜在変数、確率p、ロジット関数を用いたものです。

 それぞれ解説しますが、お好きな定義を見つけてください。個人的には潜在変数が一番好きで、確率pが一番使いやすいです。

ロジスティック回帰モデル(潜在変数Y*)

 ロジスティック回帰モデルは、目的変数Yが0か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に振り分けられます。

ロジスティック回帰モデル(確率p)

 ロジスティック回帰モデルには、Y=1になる確率pによる定義もあります。結論としては

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

$$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{1}{1+e^{-(\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k)}}$$

となります。この式は、次の3段階を経て作られています。

 第一に、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に変換します。

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

$$なおY=0と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)}}$$

となります。ちなみに、Y=0については以下が成り立ちます。

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

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

 ここで、Zを次のようにおくと、確率pは

$$Z=\beta_0+\beta_1 X_1 +\cdots + \beta_k X_k$$

$$p=\frac{1}{1+e^{-Z}}$$

と表せます。これは(標準)シグモイド関数と呼ばれるもので、Zが大きいとpが1に近づき、Zが小さいとpが0に近づきます。この定式化の利点は、

・確率は0以上1未満という確率の定義との整合性がある

点です。

ロジスティック回帰モデル(ロジット関数)

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

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

$$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)$$

(2)仮定する条件

 なお、観測されたデータからロジスティック回帰モデルのパラメーターを求めようとする場合(=ロジスティック回帰分析)、次の仮定が必要になります。

$$1.データ(X,Y)は独立で同一の分布に従い$$

$$ロジスティック回帰モデルに従う$$

$$2.標本が無作為抽出$$

$$3.説明変数Xの分散が0でない$$

 なお、ロジスティック回帰モデルには「潜在変数の誤差項はロジスティック分布に従う」も含まれているので

$$外生性:E(U|X)=0$$

$$誤差項の系列相関はなし$$

$$均一分散:Var(U_i)=\frac{\pi^2}{3}$$

$$誤差項はロジスティック分布に従う$$

も同時に仮定されています。

(3)事実

回帰モデルの目的

 回帰モデルを用いる動機には「Yを予測したい!」「XがYに与える効果を知りたい!」の二つがあります。

方法

(1)問題の構造

 「ロジスティック回帰モデルとは何か?」という問いは

・定義

・使用法(予測、効果)

に分けられます。

(2)論点と判断基準

 定義についてはすでに述べました。

 そこで、使用法について

・説明変数が決定したときに、Yはどう予測できるか(期待値)

・説明変数Xが1変化したときに、Yがどれくらい変化するか(限界効果)

という論点を設定します。

 説明したいものが説明変数と回帰係数だけで表せれればよいでしょう。

結果

(1)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}}$$

だからです。

(2)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$$

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

限界効果は変化する

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

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

・絶対値は、変わる

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

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

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

考察

(1)結論

 ロジスティック回帰分析とは以下のものです。

定義

$$「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)}}$$

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の限界効果

$$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)妥当性評価

定義について

 この記事では、確率pに関する定義をベースに議論しました。ロジット関数の定義は、確率pに関する定義から導けます。

 ただ、潜在変数の定義から、確率pの定義を導くことはしていません。この点に関して検証が必要でしょう。

使用法について

 ロジスティック回帰モデルの回帰係数の解釈としては、以上の通りとなるでしょう。

 ただ、得られたデータをロジスティック回帰モデルで解釈する場合(=ロジスティック回帰分析)、満たさなくてはいけない仮定が強いので、実際の使用の際は注意する必要があります。

(3)意義

 ロジスティック回帰モデルは、機械学習や人工知能にも使われる考え方です。また、統計手法としては、医学、心理学、社会学、経済学、マーケティングと広く使われる手法です。

 ロジスティック回帰モデルを理解することは重要でしょう。

付録1:オッズ比

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

結論

 ロジスティック回帰モデルにおける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から簡単に計算できます。

付録2:平均的限界効果の定義

 平均的限界効果の定義には、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}$$

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

付録3:限界効果の計算

結論

$$限界効果=\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$$

付録4: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))

 出力③↓

コメント欄 お気軽にコメントをお寄せください!

タイトルとURLをコピーしました