はじめに

ロジスティック回帰は、説明変数の情報にもとづいて

データがどのクラスに属するかを予測・分類する（例：ある顧客が商品を買うか買わないかを識別する）
注目している出来事が発生する確率を予測する（例：ある顧客が何％の確率で商品を買うか予測する）

ために利用されるモデルです。
この記事では、Scikit-learnライブラリを使い、ロジスティック回帰によりクラス分類を行う方法を備忘録として書いておきます。

Scikit-learn について

Scikit-learnは、Pythonの機械学習ライブラリの一つです。

公式ドキュメント：http://scikit-learn.org/stable/index.html

ロジスティック回帰について

ロジスティック回帰は、分類のためのアルゴリズムです（例：電子機器を故障するクラス or 故障しないクラスに分ける）。単純なアルゴリズムなので実装しやすい一方で、特徴量空間（説明変数となるデータがある空間）が線形分離可能な場合のみ高い性能を発揮します。線形分離可能とは、説明変数が二次元平面上にあるとすると、あるクラスに対応する説明変数（の値の集合）と、他のクラスに対応する説明変数（の値の集合）を、一本の直線を挟んで隔てることができることを言います。

ロジスティック回帰では、対数オッズを説明変数 $x_i$ の線形和で表現します。予測したいこと（正事象）の確率を $p$ としたとき、オッズは $p/(1-p)$ と書くことができ、正事象の起こりやすさを表します。オッズの対数をとったものが、対数オッズです。

\log(\frac{p}{1-p}) = w_0x_0 + w_1x_1 + \cdots +w_mx_m = \sum^m_{i=0}w_ix_i

ここで、重み $w_0$ は $x_0=1$ として切片を表します。ロジスティック回帰は、対数オッズと複数の説明変数の関係を表すモデルの重み $w_i$ を学習することが目的です。

ただ、ロジスティック回帰を利用するときに関心があるのは、説明変数の値を与えたときの正事象の確率$p$です。そこで、上式を左辺が$p$になるように変形すると、

p = \frac{1}{1+\exp(-\sum^m_{i=0}w_ix_i)}

となります。モデルの重み $w_i$ を学習後、この式を利用して説明変数の値が与えられたときの正事象の確率を求めることができ、その確率が一定値以上であれば正事象のクラスに分類することができます。

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

scikit-learnでロジスティック回帰をするには、linear_modelのLogisticRegressionモデル（公式ドキュメント：https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html ）を使います。主に利用するメソッドは以下の通りです。

fitメソッド：ロジスティック回帰モデルの重みを学習
predictメソッド：説明変数の値からクラスを予測

ここでは、UCI Machine Learning Repository (http://archive.ics.uci.edu/ml/datasets/Iris ) で公開されている、アヤメの品種データを使います。以下のコードでは、seaborn ライブラリに付属のデータセットを読み込んでいます。（以降のコードの動作環境は、Python 3.7.3, scikit-learn 0.20.3 です。）

import seaborn as sns
iris_df = sns.load_dataset('iris') # データセットの読み込み
iris_df = iris_df[(iris_df['species']=='versicolor') | (iris_df['species']=='virginica')] # 簡単のため、2品種に絞る

データを覗いてみると、こんな感じです。

iris_df.head()

各変数（データ項目）の説明は以下の通りです。

変数	説明
sepal_length	がく片の長さ
sepal_width	がく片の幅
petal_length	花弁の長さ
petal_width	花弁の幅
species	品種

各変数の値の関係を見てみましょう。以下のコードでは、各変数の関係を可視化する散布図行列をプロットしています。品種ごとの、がく片と花弁の長さおよび幅の値と品種の関係が見て取れます。特に、3行3列目のグラフに着目すると、花弁の長さ（petal_length）が小さい場合の品種はversicolorであり、大きい場合の品種はvirginicaに分類できそうです。

import matplotlib.pyplot as plt
sns.pairplot(iris_df, hue='species')
plt.show()

以降では、花弁の長さ（petal_length）を使って品種を分類するロジスティック回帰モデルを構築します。

ロジスティック回帰モデルの構築

ロジスティック回帰モデルのインスタンスを作成し、fitメソッドで説明変数の重みを学習することで、ロジスティック回帰モデルを構築します。ここでは、品種versicolorをクラス0に、品種virginicaをクラス1としており、scikit-learnライブラリのmodel_selection.train_test_splitメソッドで、データをモデル構築用データ（学習データ）と予測精度検証用データ（検証データ）に分割し、学習データをfitメソッドの引数として与え、ロジスティック回帰モデルを構築しています。

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression

X = iris_df[['petal_length']] # 説明変数
Y = iris_df['species'].map({'versicolor': 0, 'virginica': 1}) # versicolorをクラス0, virginicaをクラス1とする
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2, random_state = 0) # 80%のデータを学習データに、20%を検証データにする

lr = LogisticRegression() # ロジスティック回帰モデルのインスタンスを作成
lr.fit(X_train, Y_train) # ロジスティック回帰モデルの重みを学習

上のコードでは省略していますが、ロジスティック回帰のインスタンスを作成する際、いくつかパラメータを設定することができます（下記は、指定できるパラメータの一部）。

パラメータ	説明
penalty	正則化の方法を指定するパラメータ。'l1'か'l2'を指定する。'l1'を指定した場合はL1正則化、'l2'を指定した場合はL2正則化を行う（デフォルト値は'l2'）。L1正則化は多くの説明変数の重みが0になるようにし、特徴選択の手法としても使用される。L2正則化は説明変数の重みが大きくなりすぎるのを防ぎ、過学習を回避するために利用される。
C	正則化の強さを指定するパラメータ。正の値を指定する（デフォルト値は1.0）。Cの値が小さいほど正則化の強さが増す。
random_state	乱数生成器のシードを指定するパラメータ。数値（int型）かRandomStateインスタンスを指定する。モデルの再現性を担保するために、適当な値を指定しておいたほうがよい。

学習により得られた、ロジスティック回帰モデルの切片 $w_0$ はintercept_属性に、説明変数の係数 $w_1$ はcoef_属性に格納されます。学習結果を確認すると、係数は約0.74、切片は約-3.39であることがわかります。

print("coefficient = ", lr.coef_)
print("intercept = ", lr.intercept_)

coefficient =  [[0.73642001]]
intercept =  [-3.39483004]

このモデルに対し、検証データの説明変数の値を引数としてpredictメソッドを実行すると、それぞれの検証データが属するクラスが予測できます。

Y_pred = lr.predict(X_test)
print(Y_pred)

[1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 0 1 0 1 0]

ロジスティック回帰モデルの性能評価

クラス分類の性能評価には、主に以下のような指標が用いられます。下記のコードのように、scikit-learn.metricsライブラリの各種メソッドでこれらの指標の値を算出できます。

混同行列（confusion matrix）
正解率（accuracy）
適合率（precision）
再現率（recall）
F1スコア（F1-score）
AUC（曲線下面積）

混同行列とは、実際のクラスが0,1のデータに対して、クラス0,1に分類されたデータの個数を要素とする行列です。
出力結果の

1行1列目の要素は、実際にクラス0で正しくクラス0に分類されたデータ数（真陰性（TN: True Negative））
1行2列目の要素は、実際にはクラス0だが誤ってクラス1に分類されたデータ数（偽陽性（FP: False Positive））
2行1列目の要素は、実際にはクラス1だが誤ってクラス0に分類されたデータ数（偽陰性（FN: False Negative））
2行2列目の要素は、実際にクラス1で正しくクラス1に分類されたデータ数（真陽性（TP: True Positive））

を表します。真陽性・偽陰性・偽陽性・真陰性の「真・偽」は「分類が当たっているか」を表し、「陽性・陰性」は「クラス1・クラス0」を表します。たとえば、偽陰性は、予測されたクラスは陰性だが分類が間違っていることを表します。

正解率とは、分類したデータの総数のうち、正しく分類されたデータ数の割合です。式で書くと、(TP+TN)/(TP+FN+FP+TN)です。

適合率とは、クラス1に分類されたデータのうち、実際にクラス1であるデータ数の割合です。式で書くと、TP/(TP+FP)です。

再現率とは、実際にクラス1であるデータのうち、クラス1に分類されたデータ数の割合です。式で書くと、TP/(TP+FN)です。

F1スコアは、適合率と再現率の調和平均で定義されます。一般に、適合率と再現率はトレードオフの関係にあるため、F1スコアはこれらのバランスを評価するための指標とみることができます。

from sklearn.metrics import confusion_matrix, accuracy_score, precision_score, recall_score, f1_score

print('confusion matrix = \n', confusion_matrix(y_true=Y_test, y_pred=Y_pred))
print('accuracy = ', accuracy_score(y_true=Y_test, y_pred=Y_pred))
print('precision = ', precision_score(y_true=Y_test, y_pred=Y_pred))
print('recall = ', recall_score(y_true=Y_test, y_pred=Y_pred))
print('f1 score = ', f1_score(y_true=Y_test, y_pred=Y_pred))

confusion matrix = 
 [[ 5  5]
 [ 0 10]]
accuracy =  0.75
precision =  0.6666666666666666
recall =  1.0
f1 score =  0.8

AUCは、ROC (Receiver Operating Characteristic)曲線の下側の面積で、この値が大きいほど分類性能がよいことを示す指標です（最大値は1）。ROC曲線とは、横軸に偽陽性率（FPR: False Positive Rate）を、縦軸に真陽性率（TPR: True Positive Rate）をとった、モデルの分類性能を示した曲線です。クラス分類では、各データがクラス1に属する確率（スコア）を算出し、スコアが閾値を上回っていればデータをクラス1に分類します。したがって、その閾値を変化させると、偽陽性率（実際はクラス0だが誤ってクラス1に分類されるデータの割合）と真陽性率（実際にクラス1で正しくクラス1に分類される確率）が変化します。閾値が極端に大きい場合はどのデータもクラス1に分類されないのでROC曲線は点(0, 0)を通り、閾値が極端に小さい場合はすべてのデータがクラス1に分類されるのでROC曲線は点(1, 1)を通ります。閾値を大きい値から小さい値に変化させていくと、スコアが大きいデータから順にクラス1に分類されるようになります。実際にクラス1のデータのスコアが大きい傾向があれば、真陽性率が偽陽性率よりも早く大きくなり、ROC曲線は点(0, 0)から図の左上を通って点(1, 1)まで描かれます。その結果、ROC曲線の下側の面積であるAUCの値が大きくなります。したがって、AUCは、実際にクラス1のデータに対するスコアが大きめに予測され、実際にクラス0のデータに対するスコアは低めに予測されているかを評価する指標といえます。たとえば、下記のようなコードでROC曲線を描画し、AUCを算出することができます。

from sklearn.metrics import roc_curve, auc

Y_score = lr.predict_proba(X_test)[:, 1] # 検証データがクラス1に属する確率
fpr, tpr, thresholds = roc_curve(y_true=Y_test, y_score=Y_score)

plt.plot(fpr, tpr, label='roc curve (area = %0.3f)' % auc(fpr, tpr))
plt.plot([0, 1], [0, 1], linestyle='--', label='random')
plt.plot([0, 0, 1], [0, 1, 1], linestyle='--', label='ideal')
plt.legend()
plt.xlabel('false positive rate')
plt.ylabel('true positive rate')
plt.show()

また、下記のコードでは、ROC曲線を省略してAUCの値だけを取得しています。

from sklearn.metrics import roc_auc_score
print('auc = ', roc_auc_score(y_true=Y_test, y_score=Y_score))

auc =  0.975

実際にどの指標を利用してクラス分類の性能評価を行うべきかは、クラス分類により何を達成したいか・何を検証したいかという課題設定に依存します。たとえば、実際にクラス1に属するデータが数％しかない状況でクラス1に属するデータを当てたい（例：故障する機器の特定・CVする顧客の特定など）場合、正解率を指標として利用するのは不適当です。その理由は、どんなデータに対してもクラス0と予測してしまえば、99％に近い正解率を得られるためです。

おわりに

この記事では、scikit-learnライブラリでロジスティック回帰モデルを構築し、クラス分類の性能評価をする方法について簡単に触れました。

参考

[第2版]Python機械学習プログラミング達人データサイエンティストによる理論と実装（https://www.amazon.co.jp/dp/B07BF5QZ41/ref=dp-kindle-redirect?_encoding=UTF8&btkr=1 ）