【統計】ロジスティック回帰分析

本記事では、ロジスティック回帰分析について記載します。

【目次】

 

今回はモデル作成までです。モデルの評価は別記事で扱います。

 

計算式等


ロジスティック回帰は二値分類(y\in\{0,1\})で用いられる分析手法です。
モデル式は回帰分析と似ていますが、シグモイド関数を使って0~1の範囲の結果を出すことが特徴です(通常の回帰分析では0~1の範囲を超えてしまい、解釈不能な結果になってしまいます)。
(青線:ロジスティック回帰、緑線:単回帰分析)

f:id:cochineal19:20210112220847p:plain

モデル式、シグモイド関数および損失関数はそれぞれ次のとおりです。

モデル式
\quad \widehat{y}=\beta _{0} + \beta _{1} x_{1} +...+ \beta _{k} x_{k} = \beta ^{T} X

シグモイド関数(活性化関数)
\quad z=\sigma \left( \widehat{y} \right) = \dfrac{1}{1+exp \left( -\widehat{y} \right)} = \left( 1 + exp \left( -\widehat{y} \right)\right)^{-1} = \dfrac{exp \left( \widehat{y} \right)}{1 + exp \left(\widehat{y}\right)}

交差エントロピー誤差(分類問題の損失関数)
\quad E\left(\beta\right)=-\sum \{ y_{i} \cdot ln\left(z_{i}\right) + \left(1-y_{i} \right) \cdot ln \left( 1-z_{i} \right) \} =0
 ※ ln( ) は自然対数。 y_{i} は実測値で y_{i}\in\{0,1\} です(正解ラベルとも言われます)。

交差エントロピー誤差は、対数尤度にマイナスをつけて0を最小値としたものです(負の対数尤度)。式の第一項(  y_{i} \cdot ln\left(z_{i}\right) )で実際にイベントが起きたデータ(  y_{i}=1 )に対する予測誤差を計算します(予測90%なら10%が誤差)。式の第二項(  \left(1-y_{i} \right) \cdot ln \left( 1-z_{i} \right) )でイベントが起きなかったデータ(  y_{i}=0 )に対する予測誤差を計算します(予測10%なら10%が誤差)。
第一項も第二項も  y_{i}=z_{i} の時に  ln\left(1\right)=0 となることから、実測値と予測値が完全一致すれば合計0で誤差なしとなります。
※イメージ的には紫の棒が  y_{i}=1 の時の  z_{i} 、オレンジの棒が  y_{i}=0 の時の  1-z_{i} です。

f:id:cochineal19:20210112230833p:plain

対数尤度関数の計算方法についても掲載します。
対数尤度関数
\quad ln\ L\left( \beta\right)=f \left(y_{i}, z_{i}\right)=\sum \{ y_{i} \cdot ln\left(z_{i}\right) + \left(1-y_{i} \right) \cdot ln \left( 1-z_{i} \right) \}
尤度関数
\quad L \left( \beta \right) =\Pi \{ z_{i}^{y_{i}} \times \left( 1-z_{i} \right)^{ \left( 1 - y_{i} \right) } \}

※尤度関数は0~1の確率を "かけ算" するため、値がかなり小さくなってしまいます(コンピュータ計算に不向き)。そのため、対数変換して確率の "足し算" とした対数尤度関数が広く使われます。


そして、肝心の回帰係数 β の計算方法は次のとおりです。
回帰係数 β の計算方法

\quad \dfrac{\partial E\left( \beta \right) }{\partial \beta }=\dfrac{\partial E\left( \beta\right) }{\partial z}\dfrac{\partial z}{\partial \widehat{y}}\dfrac{\partial \widehat{y}}{\partial \beta }=\sum \left( z_{i}-y_{i}\right) x_{i}=0

式中に  z_{i} があるため、"回帰係数 β を算出するために回帰係数 β が必要" ということになります。よって、本記事ではニュートン・ラフソン法(以降、ニュートン法)による回帰係数の更新方法を取り上げます。

ニュートン法
\quad \beta ^{new}=\beta^{old} - H^{-1}\nabla E\left(\beta \right)

ニュートン法では既存の回帰係数  \beta ^{old} と演算後の回帰係数  \beta ^{new} との差がなくなるまで更新を繰り返します。最適解の算出には交差エントロピー誤差のヘッセ行列(H)と重み付き行列(R)を用います。このRはシグモイド関数微分の対角行列です。
\quad \nabla E \left( \beta \right) = \sum \left( z_{i} - y_{i} \right) x_{i} = \left( Z - Y \right) X^{T} = \sum x_{i} y_{i} - \sum x_{i} z_{i}
\quad H=\nabla\nabla E \left( \beta \right)=\sum z_{i}\left(1-z_{i}\right)x_{i}x_{i}^{T}=X^{T}RX
\quad R = z_{i} \left( 1 - z_{i} \right) = \dfrac{1}{1+exp \left( -y \right)} \times \dfrac{exp \left( -y \right) }{1+exp \left( -y \right)} = \dfrac{ exp \left( -y \right)}{\left( 1+exp \left( -y \right)\right)^{2}}

これらをまとめると次の図解になります。
f:id:cochineal19:20210117010158p:plain

近似的な方法

簡易的にロジスティック回帰分析を行う方法として、ロジット変換した値を目的変数にした回帰分析があります。二値を Yes / No の数で表すとロジット変換は次の式です。
 \quad logit = ln\left( \dfrac{P_{yes}}{1-P_{yes}}\right) = ln\left( \dfrac{N_{yes}/N_{all}}{N_{no}/N_{all}}\right) = ln\left( \dfrac{N_{yes}}{N_{no}}\right)

なお、0の対数変換 ln(0) は計算できないため、0.5 を加算してロジット変換する方法があります。これを経験ロジットと言います。 
 \quad 経験logit = ln\left( \dfrac{N_{yes}+0.5}{N_{no}+0.5}\right)

また、この回帰分析に重みを加えたい場合、ロジットの分散の逆数を重みとします(重み付き回帰分析)。経験ロジットを使う場合は3つ目の式です。
 \quad w_{i}=p_{i}\left( 1-p_{i}\right) n_{i}
 \qquad p_{i}=\dfrac{N_{yes}}{N_{all}}
 \qquad p_{i}=\dfrac{\left(N_{yes}+0.5 \right)}{\left(N_{all}+1\right)}

 

最後に、覚えておくと良さそうなシグモイド関数の性質をメモします。

シグモイド関数の性質
 \quad z=\dfrac{1}{1+exp\left(-\widehat{y}\right)}
 \quad 1-z=\dfrac{exp\left(-\widehat{y}\right)}{1+exp\left(-\widehat{y}\right)}
 \quad \dfrac{z}{1-z}=exp\left(-\widehat{y}\right) ←オッズ比

 

ノート


今回も簡単なデータで計算手順を見てみます。

■ サンプルデータ

X1 X2 Nyes Nall Pyes Logit 経験Logit
5 0 1 20 0.050 -2.944 -2.565
5 1 6 22 0.273 -0.981 -0.932
10 0 9 25 0.360 -0.575 -0.552
10 1 16 30 0.533 0.134 0.129
15 0 19 30 0.633 0.547 0.528
15 1 28 30 0.933 2.639 2.434


今回は近似的な方法で算出した回帰係数 β を初期値にしたいと思います。
その後、ニュートン法を用いて回帰係数 β を更新し、最適解を探していきます。

 

初期値設定 

f:id:cochineal19:20210130204334p:plain

LINEST関数で計算してみます。目的変数は経験ロジットとしました。今回は重み付き回帰分析はしていません(単に初期値を得るためなので)。
回帰係数 β はそれぞれ  \beta_{0}=-4.09208,\ \beta_{1}=0.32291,\ \beta_{2}=1.40674 でした。
※LINEST関数の使い方は単回帰分析で触れていますので参照ください。

cochineal19.hatenablog.com

 

回帰係数 β の更新1回目 

f:id:cochineal19:20210130205729p:plain

初期値の回帰係数 β をもとに予測値 \widehat{y},\ z を算出します。
サンプルデータの1行目は  x_{1}=5,\ x_{2}=0 なので次のとおりです。※上図Excel A列の  x_{0} は切片項用に 1 を設定しています。
\quad \widehat{y}=-4.09208+0.32291\times5+1.40674\times0\fallingdotseq-2.478
\quad z=\dfrac{1}{1+exp\left(-\widehat{y}\right)}=\dfrac{1}{exp\left(-1\times-2.478\right)}\fallingdotseq 0.077

次にヘッセ行列を求めます。
<式の再掲>
\quad H=\nabla\nabla E \left( \beta \right)=\sum z_{i}\left(1-z_{i}\right)x_{i}x_{i}^{T}=X^{T}RX
\quad R = z_{i} \left( 1 - z_{i} \right) = \dfrac{1}{1+exp \left( -\widehat{y} \right)} \times \dfrac{exp \left( -\widehat{y} \right) }{1+exp \left( -\widehat{y} \right)} = \dfrac{ exp \left( -\widehat{y} \right)}{\left( 1+exp \left( -\widehat{y} \right)\right)^{2}}

サンプルデータの1行目を取り上げると R は  R=z_{1}\left(1-z_{1}\right)=0.077\left(1-0.077\right)\fallingdotseq 0.071 です。
ヘッセ行列は次のとおりで、今回は切片+説明変数2つの3次元なので3×3の行列になります。

\quad H=\begin{bmatrix} \sum H_{00} \quad \sum H_{01} \quad \sum H_{02} \\ \sum H_{10} \quad \sum H_{11} \quad \sum H_{12} \\ \sum H_{20} \quad \sum H_{21} \quad \sum H_{22} \end{bmatrix}=\begin{bmatrix} -27.1181 \quad -289.7161 \quad -13.9390 \\ -289.7161 \quad -3455.7462 \quad -132.4044 \\ -13.9390 \quad -132.4044 \quad -13.9390 \end{bmatrix}

\qquad \sum H_{00}=\sum x_{0i}\times x_{0i}\times R_{i}\times N_{all\ i},\quad \sum H_{01}=\sum x_{0i}\times x_{1i}\times R_{i}\times N_{all\ i}, \quad \sum H_{02}=\sum x_{0i}\times x_{2i}\times R_{i}\times N_{all\ i}
\qquad \sum H_{10}=\sum x_{1i}\times x_{0i}\times R_{i}\times N_{all\ i},\quad \sum H_{11}=\sum x_{1i}\times x_{1i}\times R_{i}\times N_{all\ i},\quad \sum H_{12}=\sum x_{1i}\times x_{2i}\times R_{i}\times N_{all\ i}
\qquad \sum H_{20}=\sum x_{2i}\times x_{0i}\times R_{i}\times N_{all\ i},\quad \sum H_{21}=\sum x_{2i}\times x_{1i}\times R_{i}\times N_{all\ i},\quad \sum H_{22}=\sum x_{2i}\times x_{2i}\times R_{i}\times N_{all\ i}
 

※ 今回の集計表はレコードをまとめているため、計算式に「  \times N_{all\ i} 」を入れています。

さらに  \nabla E \left( \beta \right) を求めます。

\quad \nabla E \left( \beta \right) = \sum \left( z_{i} - y_{i} \right) x_{i} = \left( Z - Y \right) X^{T} = \sum x_{i} y_{i} - \sum x_{i} z_{i}
\quad \nabla E \left( \beta \right) = \begin{bmatrix} \sum \left( x_{0i} \times N_{yes \ i} \right) - \sum \left( x_{0i} \times z_{i} \times N_{all \ i} \right) \\ \sum \left( x_{1i} \times N_{yes \ i} \right) - \sum \left( x_{1i} \times z_{i} \times N_{all \ i} \right) \\ \sum \left( x_{2i} \times N_{yes \ i} \right) - \sum \left( x_{2i} \times z_{i} \times N_{all \ i} \right) \end{bmatrix}=\begin{bmatrix} 79 - 80.84 \\ 990 - 1008.99 \\ 50 - 51.49 \end{bmatrix} = \begin{bmatrix} -1.842 \\ -18.990 \\ -1.489 \end{bmatrix}

※ 今回の集計表はレコードをまとめているため、計算式に「  \times N_{all\ i} 」と「  \times N_{yes\ i} 」を入れています。「  \sum x_{i} y_{i} 」は  y=1 の時だけ0以外の値を持つため「  \times N_{yes\ i} 」としています。 

これら計算結果から回帰係数 β を更新すると次のようになります。
\quad \beta ^{new}=\beta^{old} - H^{-1}\nabla E\left(\beta \right)
\quad \beta ^{new}=\begin{bmatrix} -4.09208 \\ 0.32291 \\ 1.40674 \end{bmatrix} - \begin{bmatrix} -27.1181 \quad -289.7161 \quad -13.9390 \\ -289.7161 \quad -3455.7462 \quad -132.4044 \\ -13.9390 \quad -132.4044 \quad -13.9390 \end{bmatrix}^{-1}\begin{bmatrix} -1.842 \\ -18.990 \\ -1.489 \end{bmatrix}=\begin{bmatrix} -4.09514 \\ 0.32093 \\ 1.32180 \end{bmatrix}

 \beta_{0}=-4.09514,\ \beta_{1}=0.32093,\ \beta_{2}=1.32180 それぞれの更新前後の差は 0.00306, 0.00198, 0.08494 で差がありますので更新を続けます。

 

回帰係数 β の更新2回目 

f:id:cochineal19:20210130211424p:plain

更新2回目は  \beta_{0}=-4.10022,\ \beta_{1}=0.32135,\ \beta_{2}=1.32386 で更新前後の差は 0.00508, -0.00043, -0.00206 です。だいぶ良くなりましたが、まだ少し差がありますのでもう一度更新します。

 

回帰係数 β の更新3回目 

f:id:cochineal19:20210130211714p:plain更新3回目は  \beta_{0}=-4.10022,\ \beta_{1}=0.32135,\ \beta_{2}=1.32386 で更新前後の差は 0.000000, 0.000000, 0.000000 です。更新前後の回帰係数 β が同値で差がなくなりました。

よって更新はここまでにして次を最終モデルとします。

\quad \widehat{y}=-4.10022 + 0.32135 x_{1} + 1.32386 x_{2}

 

Rでの実行:

> ads <- data.frame(
+   x1=c(rep(5,20),rep(5,22),rep(10,25),rep(10,30),rep(15,30),rep(15,30))
+   ,x2=c(rep(0,20),rep(1,22),rep(0,25),rep(1,30),rep(0,30),rep(1,30))
+   ,y=c(rep(0,19),rep(1,1),rep(0,16),rep(1,6),rep(0,16),rep(1,9),rep(0,14),rep(1,16),rep(0,11),rep(1,19),rep(0,2),rep(1,28))
+ )
> fit1 <- glm(y ~ x1 + x2, family=binomial, data=ads)
> summary(fit1)
Call:
glm(formula = y ~ x1 + x2, family = binomial, data = ads)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0812  -0.8307   0.4935   0.8906   2.2684  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.10022    0.72055  -5.690 1.27e-08 ***
x1           0.32135    0.05562   5.778 7.58e-09 ***
x2           1.32386    0.40344   3.281  0.00103 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 217.64  on 156  degrees of freedom
Residual deviance: 164.56  on 154  degrees of freedom
AIC: 170.56
Number of Fisher Scoring iterations: 4


SASでの実行:

data ads(drop=i);
  do i=1 to  1; x1= 5; x2=0; y=1; output; end;
  do i=1 to 19; x1= 5; x2=0; y=0; output; end;
  do i=1 to  6; x1= 5; x2=1; y=1; output; end;
  do i=1 to 16; x1= 5; x2=1; y=0; output; end;
  do i=1 to  9; x1=10; x2=0; y=1; output; end;
  do i=1 to 16; x1=10; x2=0; y=0; output; end;
  do i=1 to 16; x1=10; x2=1; y=1; output; end;
  do i=1 to 14; x1=10; x2=1; y=0; output; end;
  do i=1 to 19; x1=15; x2=0; y=1; output; end;
  do i=1 to 11; x1=15; x2=0; y=0; output; end;
  do i=1 to 28; x1=15; x2=1; y=1; output; end;
  do i=1 to  2; x1=15; x2=1; y=0; output; end;
run;

ods output LackFitChiSq=HOSLEM1;
proc logistic data=ads descending;
  model y(event='1') = x1 x2 / outroc=ROC1 lackfit;
  output out=OutDS;
run;

f:id:cochineal19:20210130212212p:plain



 

統計ソフトってすごく便利。。と体感した。

 

プログラムコード


■ Rのコード

fit <- glm(y ~ x1 + x2, family=binomial, data=ads)
summary(fit)

 

SASのコード

proc logistic data=[InputDS] descending;
  model y(event='1') = x1 x2;
  output out = [OutDS];
run; quit;

 

Pythonのコード

整備中

 

雑談


回帰係数 β の計算方法で示した交差エントロピー誤差の偏微分の詳細です。
「連鎖律」という方法で微分します。

交差エントロピー誤差
\quad E\left(\beta\right)=-\sum \{ y_{i} \cdot ln\left(z_{i}\right) + \left(1-y_{i} \right) \cdot ln \left( 1-z_{i} \right) \} =0
偏微分
\quad \dfrac{\partial E\left( \beta \right) }{\partial \beta }=\dfrac{\partial E\left( \beta\right) }{\partial z}\dfrac{\partial z}{\partial \widehat{y}}\dfrac{\partial \widehat{y}}{\partial \beta }
 \qquad \qquad =-\sum \{ y_{i} \times \dfrac{1}{z_{i}} \times z_{i}\left( 1-z_{i}\right) \times x_{i} + \left(1-y_{i}\right) \times \dfrac{1}{1-z_{i}} \times \left( -z_{i}\left( 1-z_{i} \right) \right) \times x_{i}  \}
 \qquad \qquad =-\sum \{ y_{i} \left( 1-z_{i}\right) + \left(1-y_{i}\right) \left( -z_{i} \right)\} x_{i}
 \qquad \qquad =-\sum \{ y_{i} - y_{i} z_{i} - z_{i} + y_{i} z_{i}\} x_{i}
 \qquad \qquad =-\sum \{ y_{i} - z_{i} \} x_{i}
 \qquad \qquad =\sum \left( z_{i}-y_{i}\right) x_{i}=0

図解するとこのような感じ。シグモイド関数微分は公式で z(1-z) になります。

f:id:cochineal19:20210117223042p:plain 

-----

なお、SASやRのロジスティック回帰では「ニュートン法」ではなく、「フィッシャースコアリング法」が用いられているようでした(記事を書き終わってから気づいた)。
この方法は、ニュートン法で使った「ヘッセ行列」を「フィッシャー情報行列」に置き換えたものです。計算コストが良い様子。
今回は更新の仕組みが分かれば良いと思っていたため、この記事は終了にします。

フィッシャースコアリング法について Wikipedia を参照。

尤度方程式 - Wikipedia

 

参考


統計学(出版:東京図書), 日本統計学会編

多変量解析実務講座テキスト, 実務教育研究所

統計学入門−第10章

PRML4.3.3

 


サイトマップ

cochineal19.hatenablog.com

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。