ロジスティック回帰の回帰係数を解釈したい

はじめに

今年の6月ごろから、GoogleがCourseraにて提供しているGoogle Advanced Data Analytics Professional Certificateを受講しており（*1）、最近全7コースのうちの5つ目に当たる以下のコースを修了することができました。

Regression Analysis: Simplify Complex Data Relationships

一方、コースを進めていく中で、現在進行形でデータサイエンティストとしてお金をもらって仕事をする上では、こちらのコースで紹介されている以上の数理的バックグラウンドを理解している必要性を感じました。
そういう意味で「ロジスティック回帰わからん」となったので、本記事では、ロジスティック回帰で算出した回帰係数を解釈してインサイトを導くことをゴールとして、ロジスティック回帰について解説していきます。なお、回帰分析一般の話はカバーしませんのでご了承ください。

ロジスティック回帰のモチベーション

基本形とオッズについて

ロジスティック回帰というと、目的変数が2値である場合に適応されるモデ
ルであり、以下の式で表されます。

$${log (\frac{p}{1-p}) = \beta_0 + \beta_1 x_1 + \cdots + \beta_n x_n}$$

この式の左辺にあらわれているある事象の生起確率のオッズ（*2）がこのロジスティック回帰で必要な前提知識と言えるでしょう。このオッズは、ある事象が発生しない確率（1-p）に対する発生確率（p）の割合を表しているといえ、両者がランダムに発生する（=50%）時オッズは1を取り、値の差が広がるにつれてオッズ自身の大小も大きく変わります。

試しに1%単位で取りうるすべての確率pに対してこのオッズを求め、プロットしてみます。

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

# 任意の確率pを与えてそのオッズを計算する関数
def calculate_odds(probability):
    if probability <= 0 or probability >= 1:
        return "Invalid probability. Probability should be between 0 and 1 (exclusive)."
    odds = probability / (1 - probability)
    return odds

# 0から100まで1%単位でprobabilityを変化させてオッズを計算
odds_list = [calculate_odds(p / 100) for p in range(1,100)]

# オッズの傾向をプロット
sns.lineplot(odds_list)

# オッズの度数分布をプロット
sns.histplot(odds_list)

発生確率が100%に近づくにつれて指数関数的に値が大きくなっていますね。
もう一つ試しに、オッズのヒストグラムを描いてみます。

かなりテールヘビーな分布になっていますね。これはすなわち、発生確率が低い時（P<0.5）のオッズは0から1の範囲に収まる一方、発生確率が高い時（P>0.5）のオッズはより大きな幅の値を取ることを意味します。このようにオッズを導入したことによるスケールの不均衡を解消する、というのが続いて解説するロジット変換のモチベーションになっているようです。

ロジット変換

上記の通り、オッズに対数変換を施すことでスケールの修正をすることがロジット変換のモチベーションと捉えられそうです。
では先ほど同様実際に1%単位で取りうるすべての確率pに対してオッズ（今回は対数を取ります）を求め、その値をプロットしてみます。

# 0から100まで1%単位でprobabilityを変化させて対数オッズを計算
log_odds_list = [np.log(calculate_odds(p / 100)) for p in range(1,100)]

# 対数オッズの傾向をプロット
sns.lineplot(log_odds_list)

# 対数オッズの度数分布をプロット
sns.histplot(log_odds_list)

いわゆるロジット関数と呼ばれるやつですね。続いてヒストグラムを見てみます。

p=0.5の時を0として、p<0.5の時とp>0.5の時でスケールが揃っているのがわかります。

ロジスティック回帰モデルの解釈方法

では続いて、実データをロジスティック回帰モデルに当てはめた時に結果をどのように解釈するのかについてみていきます。
今回のデータはseabornのirisデータを使用します。

# データセットの読み込み
iris_df = sns.load_dataset('iris')
# 2品種に絞る
iris_df = iris_df[(iris_df['species']=='versicolor') | (iris_df['species']=='virginica')].reset_index(drop=True)

データセットの中身はこんな感じ。目的変数はカテゴリを表す speciesで、説明変数となるそれ以外の項目は連続値になっています。 

データの中身が確認できたところで、今回構築するロジスティック回帰モデルを式で表現してみます。

$${log( \frac{P(species='virginica')}{1-P(species='virginica')} )= intercept + \beta_1 x_{sepallength} + \beta_2 x_{sepalwidth} + \beta_3 x_{petal_length}+ \beta_4x_{petalwidth}}$$

さらに、左辺の対数を外すと

$${\begin{array}{} \frac{P(species='virginica')}{1-P(species='virginica')} &=& e^{intercept + \beta_1 x_{sepallength} + \beta_2 x_{sepalwidth} + \beta_3 x_{petal_length}+\beta_4x_{petalwidth}}\\ &=& e^{intercept} e^{\beta_1 x_{sepallength} }e^{\beta_2 x_{sepalwidth}}e^{\beta_3 x_{petal_length}}e^{\beta_4 x_{petalwidth}}\end{array}}$$

となります。こちらを念頭に、データをロジスティック回帰モデルに当てはめて結果を出してみます。今回のモデリングにおいては、speciesのうちのverginicaが発生する確率を対象とします。ここで、verginicaが発生しない確率はversicolorが発生する確率と言い換えることもできます。

from sklearn.linear_model import LogisticRegression

# 説明変数・目的変数の設定
X = iris_df.drop("species",axis=1)# 説明変数
y = iris_df['species'].map({'versicolor': 0, 'virginica': 1}) # versicolorをクラス0, virginicaをクラス1とする

# ロジスティック回帰モデルへのフィッティング
clf = LogisticRegression(random_state=0).fit(X, y)

モデリングの結果

$$
\begin{array}{|c|c|c|} 
説明変数 & exp(回帰係数) & 回帰係数 \\ \hline
sepal length & -0.39 & 0.67  \\ \hline
sepalwidth & -0.51 & 0.60 \\ \hline
petallength & 2.93 & 18.74 \\ \hline
petalwidth & 2.42 & 11.21 \\ \hline
\end{array}
$$

sepal_length / sepal_widthとpetal_length / petal_widthで傾向が異なりますね。また、上記のlogを外した式から見て取れるように、オッズは各説明変数の指数変換の乗法モデルといえそうです。したがって、LogisticRegression().fit()が出力するcoef_のべき乗により、目的変数のオッズの変化を解釈することができ、その値を上の表の一番左に示しています。
これを見ると、sepal_lengthはその長さが1単位変化するごとにオッズが0.67倍に変化するということになり、それはすなわちオッズが縮小していくことを意味しており、verginicaが発生しない＝versicolorが発生する確率を高めると言えそうです。sepal_widthも同様で、その程度が1単位の変化に対して0.6倍のオッズの縮小をもたらすといえるでしょう。
一方、petal_length / petal_widthについては1単位変化するごとにオッズをそれぞれ18.74倍、11.21倍拡大すると解釈することができ、verginicaが発生する方向に影響していると言えるでしょう。

さいごに

正直調べれば調べるほど、どれだけ調べてもロジスティック回帰の数理的バックグラウンドがわからんとなる一方で、上記のようなモデル化を実現する神秘というか美しさを感じました。

今回の記事は、筆者の勉強のための備忘録という意味合いが強いため、王道の説明と外れたストーリーや少々乱暴な論理展開も見られるかと思います。（それも含めてロジスティック回帰理解への挑戦ということにします。笑）
とはいえいくらなんでもな表現は明らかな間違いがありましたら本当に教えてくださいmm

今回はこの辺で。ありがとうございました！

（*1）コースを受けてみての所感は、全7コース修了後にでもしたためようかと思っています。

（*2）元々は、発生確率pを以下のロジスティック関数で表すことができる。

$${p = \frac{exp( \beta_0 + \beta_1 x_1 + \cdots + \beta_n x_n)}{1+exp( \beta_0 + \beta_1 x_1 + \cdots + \beta_n x_n)}}$$

これを展開していくとオッズが出てくるため、オッズはロジスティック回帰にとって重要な前提となっているが、そもそもなぜ上記のロジスティック関数で発生確率pを表現することができるのか今回は理解しきれませんでした。。。誰か分かる方いましたらぜひ教えてください🙇

参考

統計モデルと推測

ロジスティック回帰を理解する(1)_オッズとロジット変換について - Qiita

7-1. ロジスティック回帰分析1 | 統計学の時間 | 統計WEB

7-2. ロジスティック回帰分析2 | 統計学の時間 | 統計WEB

Scikit-learn でロジスティック回帰（クラス分類編） - Qiita