はじめに
ロジスティック回帰は、説明変数の情報にもとづいて
- データがどのクラスに属するかを予測・分類する（例：ある顧客が商品を買うか買わないかを識別する）
- 注目している出来事が発生する確率を予測する（例：ある顧客が何％の確率で商品を買うか予測する）
ために利用されるモデルです。
この記事では、Scikit-learnライブラリを使い、ロジスティック回帰により確率を予測する方法を備忘録として書いておきます。
Scikit-learn について
Scikit-learnは、Pythonの機械学習ライブラリの一つです。
ロジスティック回帰について
ロジスティック回帰は、確率（例：機器の故障率）を予測するためのアルゴリズムです。機械学習系の文献では、分類のためのアルゴリズムとして紹介されていることが多いロジスティック回帰ですが、予測した確率（スコア）をもとにクラス分類を行う（例：機器を故障するクラス or 故障しないクラス に分類する）ものであるため、当然確率そのものを予測するためにも利用できます。
ロジスティック回帰では、対数オッズを説明変数 $x_i$ の線形和で表現します。予測したいこと（正事象）の確率を $p$ としたとき、オッズは $p/(1−p)$ と書くことができ、正事象の起こりやすさを表します。オッズの対数をとったものが、対数オッズです。
\log(\frac{p}{1−p}) = w_0x_0 + w_1x_1 + \cdots + w_mx_m = \sum_{i=0}^m w_ix_i
ここで、重み $w_0 $は $x_0=1$ として切片を表します。ロジスティック回帰は、対数オッズと複数の説明変数の関係を表すモデルの重み $w_i$ を学習することが目的です。
ただ、ロジスティック回帰を利用するときに関心があるのは、説明変数の値を与えたときの正事象の確率 $p$ です。そこで、上式を左辺が $p$ になるように変形すると、
p=\frac{1}{1+\exp(−\sum_{i=0}^m w_ix_i)}
となります。モデルの重み $w_i$ を学習後、この式を利用して説明変数の値が与えられたときの正事象の確率を求めることができます。
ロジスティック回帰モデル
scikit-learnでロジスティック回帰をするには、linear_modelのLogisticRegressionモデル（公式ドキュメント：https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html ）を使います。主に利用するメソッドは以下の通りです。
- fitメソッド：ロジスティック回帰モデルの重みを学習
- predict_probaメソッド：説明変数の値から各クラスに属する確率を予測
ここでは、sklearn.datasets ライブラリの make_classification メソッドで生成したデータに対して確率を予測します。以下のコードでは、200,000サンプルのデータを生成しており、説明変数（特徴量）の数は20、目的変数は0か1のクラスで、クラス0とクラス1のデータ数の比率はおよそ1:1です。（以降のコードの動作環境は、Python 3.7.3, scikit-learn 0.20.3 です。）
from sklearn.datasets import make_classification
X, Y = make_classification(n_samples=200000, n_features=20, n_informative=2, n_redundant=2)
以降では、説明変数Xを利用してロジスティック回帰モデルを構築し、目的変数が1となる確率を予測します。
ロジスティック回帰モデルの構築
ロジスティック回帰モデルのインスタンスを作成し、fitメソッドで説明変数の重みを学習することで、ロジスティック回帰モデルを構築します。ここでは、scikit-learnライブラリのmodel_selection.train_test_splitメソッドで、データをモデル構築用データ（学習データ）と予測精度検証用データ（検証データ）に分割し、学習データをfitメソッドの引数として与え、ロジスティック回帰モデルを構築しています。
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.5, random_state=0) # 50%のデータを学習データに、50%を検証データにする
lr = LogisticRegression() # ロジスティック回帰モデルのインスタンスを作成
lr.fit(X_train, Y_train) # ロジスティック回帰モデルの重みを学習
ロジスティック回帰モデルのインスタンスを生成する際、モデルの重みを調整する正則化などの設定をするハイパーパラメータの設定ができます。詳細は公式ドキュメント、もしくは Scikit-learn でロジスティック回帰（クラス分類編） をご覧ください。
学習により得られた、ロジスティック回帰モデルの切片 $w_0$ はintercept_属性に、説明変数の係数 $w_1, \cdots, w_{20}$ はcoef_属性に格納されます。
print("coefficient = ", lr.coef_)
print("intercept = ", lr.intercept_)
coefficient = [[ 1.07367390e-03 5.17036366e-01 5.21447708e-03 -9.98281654e-03
3.44933503e-01 2.66695180e+00 -8.91173651e-03 2.27093901e-02
5.33371028e-03 -1.55647553e-02 -1.26011741e-02 -1.06324472e-03
3.09072395e+00 -4.44766398e-03 7.88774312e-03 2.99896065e-02
-1.35581106e-02 -1.96393314e-02 2.37722077e-03 2.48242512e-02]]
intercept = [-0.89110432]
このモデルに対し、検証データの説明変数の値を引数としてpredict_probaメソッドを実行すると、それぞれの検証データの目的変数（クラス）が0である確率と、1である確率が予測できます。出力結果の1列目が目的変数が0である確率、2列目が目的変数が1である確率です。
probs = lr.predict_proba(X_test)
print(probs)
[[7.17584367e-02 9.28241563e-01]
[9.99301416e-01 6.98583960e-04]
[1.80200179e-03 9.98197998e-01]
...
[3.34485484e-02 9.66551452e-01]
[6.32885474e-02 9.36711453e-01]
[4.45884313e-04 9.99554116e-01]]
ロジスティック回帰モデルの性能評価
確率予測の性能評価には、主に以下のような指標が用いられます。
- キャリブレーションプロット
- ロジスティック損失（logistic loss）
- brier score（ブライヤースコア）
キャリブレーションプロットは、予測した確率が実際の確率（割合）に一致するかを確認するものであり、ロジスティック損失やbrier scoreは、モデル間の予測精度の比較に用いることが多いです。
キャリブレーションプロットは、横軸をモデルで予測した目的変数が1である確率（スコア）を複数の区間（ビン）に区切ったもの、縦軸を各ビンに該当するデータのうち目的変数が実際に1であるデータの割合をとったグラフです。スコアと、実際に目的変数が1のデータの割合（分布）が一致するかを確認するために利用します。スコアと実際の割合が一致するとき、キャリブレーションプロットは45度線に乗るため、キャリブレーションプロットが45度線に近いほど、モデルの予測性能が高いといえます。
たとえば、以下のようなコードでキャリブレーションプロットを作成できます。ここでは、検証データに対して目的変数が1である確率を予測し、sklearn.calibrationライブラリのcalibration_curveメソッドを利用してキャリプレーションプロットを作成しています。今回の検証データで作成したキャリブレーションプロットは、45度線より若干下ぶれていますが、おおむね45度線に近いかな、という印象です。
import matplotlib.pyplot as plt
from sklearn.calibration import calibration_curve
prob = lr.predict_proba(X_test)[:, 1] # 目的変数が1である確率を予測
prob_true, prob_pred = calibration_curve(y_true=Y_test, y_prob=prob, n_bins=20)
fig, ax1 = plt.subplots()
ax1.plot(prob_pred, prob_true, marker='s', label='calibration plot', color='skyblue') # キャリプレーションプロットを作成
ax1.plot([0, 1], [0, 1], linestyle='--', label='ideal', color='limegreen') # 45度線をプロット
ax1.legend(bbox_to_anchor=(1.12, 1), loc='upper left')
ax2 = ax1.twinx() # 2軸を追加
ax2.hist(prob, bins=20, histtype='step', color='orangered') # スコアのヒストグラムも併せてプロット
plt.show()
ロジスティック損失は、確率予測の正確さを測るための指標のひとつで、誤差が大きいほど損失が指数関数的に大きくなる特徴があります。式で書くと、
logistic \ \ loss = -\frac{1}{n}\sum_{i=1}^n \{y_i \log p_i + (1-y_i)\log(1-p_i)\}
です。ここで、$p_i$ は $i$ 番目のデータの目的変数が1であると予測した確率（スコア）、$y_i$ は $i$ 番目のデータの実際の目的変数の値です。ちなみに、ロジスティック回帰は、ロジスティック損失が最小となるようにモデルの重みを学習しています。
以下のコードでは、sklearn.metrics ライブラリの log_loss を利用して、検証データに対するロジスティック損失を算出しています。
from sklearn.metrics import log_loss
print('logistic loss = ', log_loss(y_true=Y_test, y_pred=prob))
logistic loss = 0.13207794568463985
brier socre（ブライヤースコア）は、ロジスティック損失と同様にスコアと実際の目的変数の値の誤差を表した指標で、スコアと実際の目的変数の値の平均二乗誤差で定義されます。式で書くと、
brier \ \ score = \frac{1}{n} \sum_{i=1}^n (p_i - y_i)
です。$p_i$ と $y_i$ の定義は、ロジスティック損失の説明に記載した定義と同様です。
以下のコードでは、sklearn.metrics ライブラリの brier_score_loss を利用して、検証データに対するbrier score を算出しています。
from sklearn.metrics import brier_score_loss
print('brier score = ', brier_score_loss(y_true=Y_test, y_prob=prob))
brier score = 0.031796375907874444
参考までに、スコアに対してロジスティック損失とbrier scoreがどのような値をとるかを見てみます。下図では、目的変数の実績値が1のときのロジスティック損失（青実線）とbrier score（緑実線）、目的変数の実績値が0のときのロジスティック損失（黄点線）とbrier score（赤点線）をプロットしています。横軸が、モデルで予測した目的変数が1となる確率、縦軸がロジスティック損失およびbrier scoreの値です。この図を見ると、スコアと実績値の差が大きくなったときの値は、ロジスティック損失のほうが大きくなる傾向があることがわかります。
おわりに
この記事では、scikit-learnライブラリでロジスティック回帰モデルを構築し、確率予測の性能評価をする方法について簡単に触れました。
参考
- [第2版]Python機械学習プログラミング 達人データサイエンティストによる理論と実装（https://www.amazon.co.jp/dp/B07BF5QZ41/ref=dp-kindle-redirect?_encoding=UTF8&btkr=1 ）