PR曲線とROC曲線を理解する（後編）

こんばんは🐼

今回は、3週間ほど前に書いた「PR曲線とROC曲線を理解する（前編）」という記事の続編です。

ontheblink.hatenablog.com

前回はPR曲線・ROC曲線について理解するために欠かせない混同行列や、適合率・再現率といった指標について、その意味と共になぜ必要なのかを説明しました。

今回の目的としては、

• PR曲線・ROC曲線とは何か
• ビジネス上どのように役立てることができるのか

ということについて、説明・考察していきたいと思います。

PR曲線・ROC曲線を描く目的

PR曲線・ROC曲線について説明する前に、前回の記事のおさらいも含めてPR曲線(ROC曲線)をなぜ描くのかについてお話します。

それは、「作成した機械学習をビジネス上どのように役立てるか」というところに繋がってきます。

• スパムメールを除去するタスクにおいて、必要なメールを削除してしまうよりも、受信トレイにスパムメールが残ってしまう方がマシ。
• ガンの陽性を発見するタスクにおいて、陽性の患者を陰性と診断してしまうよりも、陰性の患者を陽性と診断してしまう方がマシ。

といったように、現実では同じ予測の誤りでも、その誤りによって生じるコストやリスクが異なるため、精度(accuracy)だけではモデルの良し悪しの判断が付かないケースは往々にしてあります。

この時、「ビジネス上どちらの誤りをどの程度優先して減らすか」ということを含めてモデル評価するために、PR曲線・ROC曲線を描画します。

どちらの誤りも減らすということは根本的なモデルの性能を向上させる他ないので、基本的には片方の誤りを減らすと片方の誤りは増えるというトレードオフの関係にあります。

サンプルデータセット

本記事では、不均衡データに対するモデルの性能評価について簡潔な例を示すために、手頃に手に入るskicit-learnのトイデータセットのラベルを不均衡にして使用しようと思います。

import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

# digitsデータセットの読み込み
digits = load_digits()
y = digits.target == 8

X_train, X_test, y_train, y_test = train_test_split(
    digits.data, y, random_state=42, stratify=y)

数字が「8である(1:陽性)」「8でない(0:陰性)」の2種類のラベルを付けることで、不均衡データに対する2値分類タスクの例として使います。

np.bincount(y)

全サンプルのうち陽性クラスが174サンプル(9.7%)、陰性クラスが1623サンプル(90.3%)であり、不均衡データとなっていることがわかります。

モデリング

モデルの性能を比較するために、様々なモデルで訓練を行い、テストデータに対する認識精度(accuracy)を出力します。

ダミー分類器

はじめに問答無用で多数派クラスと予測するダミー分類器を作ります。

from sklearn.dummy import DummyClassifier
# strategy = 'most_frequent'で常に多数派クラスを予測
dummy_clf = DummyClassifier(strategy='most_frequent').fit(X_train, y_train)
pred_dummy = dummy_clf.predict(X_test)
print('test accuracy:', dummy_clf.score(X_test, y_test))

認識精度は90.2%となります。内容を知らない状態で精度が90.2%ということだけを聞くと「性能がいいモデルだ」と感じる人が多いかもしれませんが、実情は多数派クラスを出力しているだけです。多数派クラスが99.9%のデータであれば認識精度は99.9%となります。

この簡単な例だけでも認識精度のみを使ってモデルを評価することの危険さが解ると思います。

ロジスティック回帰

from sklearn.linear_model import LogisticRegression

param_grid = [{'C': [0.001, 0.01, 0.1, 1, 10]}]

grid_search = GridSearchCV(LogisticRegression(random_state=42, solver='lbfgs', max_iter=10000), param_grid, cv=5)
grid_search.fit(X_train, y_train)
print(f"Best parameters: {grid_search.best_params_}")
print(f"Best validation score {grid_search.best_score_}")

logreg_clf = grid_search.best_estimator_
pred_logreg = logreg_clf.predict(X_test)
print('test accuracy:', logreg_clf.score(X_test, y_test))

グリッドサーチという手法を使って、適切なハイパーパラメータを探索していますがここではPR曲線・ROC曲線の理解が主テーマなので、パラメータの探索は最低限に行っています。（以降のモデルについても同様） 正則化パラメータC: 0.01のロジスティック回帰分類器モデルはテストデータに対して96.4%の認識精度が得られるようです。

決定木

from sklearn.tree import DecisionTreeClassifier

param_grid = [{'max_depth': [1, 3, 5, 7, 9]}]

grid_search = GridSearchCV(DecisionTreeClassifier(random_state=42), param_grid, cv=5)
grid_search.fit(X_train, y_train)
print(f"Best parameters: {grid_search.best_params_}")
print(f"Best validation score {grid_search.best_score_}")

tree_clf = grid_search.best_estimator_
pred_tree = tree_clf.predict(X_test)
print('test accuracy:', tree_clf.score(X_test, y_test))

決定木ではロジスティック回帰よりも認識精度が低く94%となります。しかしここでも注意が必要なのは、常に多数派クラスを出力するダミー分類器ですら90.2%の認識精度を出せるのだから、これだけではまだ「ロジスティック回帰よりも決定木が悪いモデルである」と結論付けることは出来ないということです。

不均衡データでは一般に多数派クラスに寄せて予測すると精度が出やすいので、少数派クラスをどの程度正しく予測できているかどうかは精度からは詳しく判定できないことに注意が必要です。

ランダムフォレスト

from sklearn.ensemble import RandomForestClassifier

forest_clf = RandomForestClassifier(n_estimators=400, random_state=42)
forest_clf.fit(X_train, y_train)

pred_forest = forest_clf.predict(X_test)
print('test accuracy:', forest_clf.score(X_test, y_test))

ランダムフォレストの認識精度は96.9%となりました。続いて2値分類タスクのモデルを正しく評価するための有効手段である混同行列と各種性能指標を表示します。

混同行列・各種性能指標の表示

PR曲線は前回の記事で説明した適合率・再現率のトレードオフの関係を、同じくROC曲線は後述する真陽性率と偽陽性率の関係を視覚的に表したものです。

そのためまずは可視化のもととなる混同行列と各種性能指標の表示を行います。これらはsklearn.metricsモジュールのconfusion_matrix関数および、classification_report関数で実現できます。

ダミー分類器

from sklearn.metrics import confusion_matrix, classification_report
import matplotlib.pyplot as plt
import seaborn as sns

targets = ['not 8', '8']

matrix = confusion_matrix(y_test, pred_dummy)
# 混同行列の可視化
sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets)

# 各種性能指標の表示
print(classification_report(y_test, pred_dummy,
                            target_names=targets))

当然ながら全てを「8以外(not 8)」と予測したので、「8である(8)」と予測したサンプルの数は0になっています。

classification_reportでは「8以外」を陽性としたときの適合率・再現率・F値と、「8」を陽性としたときの適合率・再現率・F値の双方が出力されます。

これにより、少数派クラスについて全く予測ができていない（再現率がゼロである）ことが分かります。現実では少数派クラスが見つけたいもの（陽性）であることがほとんどです。

ロジスティック回帰

他のモデルについても同様に、混同行列と各種指標を表示します。

matrix = confusion_matrix(y_test, pred_logreg)
sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets)

print(classification_report(y_test, pred_logreg,
                            target_names=targets))

ロジスティック回帰では、「8」を陽性としたときの再現率は0.68まで増加しており、68%の「8」を正しく予測できていることがわかります。また、適合率が0.94と高いことから「8」であると予測したものが、実際には「8以外」であることは少ないといえます。

十分な精度が得られているかは別で検討が必要ですが、幾分意味のあるモデルとなっているといえるでしょう。

決定木

matrix = confusion_matrix(y_test, pred_tree)
sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets)

print(classification_report(y_test, pred_tree,
                            target_names=targets))

決定木の「8」を陽性としたときの再現率はロジスティック回帰とほぼ同等ですが、適合率は大きく下がっており、全体としてロジスティック回帰より性能が悪いことが分かります。

ランダムフォレスト

matrix = confusion_matrix(y_test, pred_forest)
sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets)

print(classification_report(y_test, pred_forest,
                            target_names=targets))

ランダムフォレスト分の性能指標はロジスティック回帰とよく似ています。「8」を陽性としたときの適合率が100%であることから、より曖昧な「8」についても「8」と予測することで、適合率と再現率のトレードオフを調整することができます。（そのような用途でPR曲線が便利だと考えられます）

PR曲線・ROC曲線の描画

いよいよPR曲線・ROC曲線の描画です。PR曲線・ROC曲線の解釈については、前編の記事で説明した指標（classification report関数で出力される指標）についての理解が出来ていれば難しいことはありません。

探しているクラス（今回の例でいえば「8」）を陽性とした上で、

• PR曲線 : X軸を再現率(recall)、Y軸を適合率(precision)としてあらゆる閾値についてプロットしたグラフ
• ROC曲線 : X軸を真陽性率(=再現率)、Y軸を偽陽性率としてあらゆる閾値についてプロットしたグラフ

この「あらゆる閾値について」という所がポイントで、各種の分類アルゴリズムは分類の際に、サンプルが陽性に属する確率というものを出力します。

この確率は各分類器のpredict_probaメソッドで出力できます。（但し一部の分類器はpredict_probaメソッドを持たず、decision_functionメソッドで予測の不確実性を評価します）

デフォルトでは陽性の確率が50%以上のサンプルを陽性と予測しますが、この閾値は変更可能です。

閾値を増やすと、より陽性の確率が高いもののみを陽性と判定するため、適合率が高まり、再現率が下がる傾向になります。逆に閾値を減らすと、陽性かどうか曖昧なものまで陽性と判定するため、適合率が下がり、再現率が高まる傾向になります。

PR曲線はこの閾値の変更による適合率と再現率のトレードオフの関係を可視化したものです。

ROC曲線についてもトレードオフの関係を表しています。真陽性率は再現率と同じですが、偽陽性率については「本当は陰性のサンプルの中で陽性と予測された割合」を表しています。

偽陽性率は誤りの割合なので低いほど良い指標です。初見ではやや混乱するので定期的に確認して覚えていきましょう。

偽陽性率は、1 - (陰性を陽性と見た時の再現率)と言い換えることも出来ます。（陰性を陽性と見た時の再現率、つまり本当に陰性のサンプルの中で陰性と予測された割合は「特異度」と呼びます）

そのため自分はROC曲線が何を表しているか混乱した時は、陽性の再現率と陰性の再現率のトレードオフと覚えています。（厳密に言うと再現率は「実際に陽性のサンプルのうち何%を正しく予測できたのか」なので、再現率と特異度のトレードオフですが、特異度という言葉はあまり使われないので）

長くなりましたが、ロジスティック回帰とランダムフォレストについてそれぞれの曲線を描画してみましょう。

sns.set();
from sklearn.metrics import precision_recall_curve, roc_curve

fig, axes = plt.subplots(1, 2, figsize=(12, 6))

# PR曲線の描画 (ロジスティック回帰, ランダムフォレスト)
precision_logreg, recall_logreg, thresholds_logreg = precision_recall_curve(y_test, logreg_clf.predict_proba(X_test)[:, 1])
precision_forest, recall_forest, thresholds_forest = precision_recall_curve(y_test, forest_clf.predict_proba(X_test)[:, 1])
default_threshold_logreg = np.argmin(np.abs(thresholds_logreg - 0.5))
default_threshold_forest = np.argmin(np.abs(thresholds_forest - 0.5))
axes[0].plot(recall_logreg, precision_logreg, label='LogisticRegression')
axes[0].plot(recall_forest, precision_forest, label='RandomForestClassifier', linestyle='--')
axes[0].plot(recall_logreg[default_threshold_logreg], precision_logreg[default_threshold_logreg],
             'x', c='k', markersize=10, label='default threshold(50%)')
axes[0].plot(recall_forest[default_threshold_forest], precision_forest[default_threshold_forest],
             'x', c='k', markersize=10)
axes[0]
axes[0].legend()
axes[0].set_title('PR Curve', fontsize=14)
axes[0].set_xlabel('recall')
axes[0].set_ylabel('precision')

# ROC曲線の描画 (ロジスティック回帰, ランダムフォレスト)
fpr_logreg, tpr_logreg, thresholds_logreg = roc_curve(y_test, logreg_clf.predict_proba(X_test)[:, 1])
fpr_forest, tpr_forest, thresholds_forest = roc_curve(y_test, forest_clf.predict_proba(X_test)[:, 1])
default_threshold_logreg = np.argmin(np.abs(thresholds_logreg - 0.5))
default_threshold_forest = np.argmin(np.abs(thresholds_forest - 0.5))
axes[1].plot(fpr_logreg, tpr_logreg, label='LogisticRegression')
axes[1].plot(fpr_forest, tpr_forest, label='RandomForestClassifier', linestyle='--')
axes[1].plot(fpr_logreg[default_threshold_logreg], tpr_logreg[default_threshold_logreg],
             'x', c='k', markersize=10, label='default threshold(50%)')
axes[1].plot(fpr_forest[default_threshold_forest], tpr_forest[default_threshold_forest],
             'x', c='k', markersize=10)
axes[1].legend()
axes[1].set_title('ROC Curve', fontsize=14)
axes[1].set_xlabel('TPR(True Positive Rate)')
axes[1].set_ylabel('FPR(False Positive Rate)')

plt.show()

X印でデフォルトの閾値の際の値を表しています。PR曲線は右上に広がる（適合率も再現率も高い）ほど、ROC曲線は左上に広がる（真陽性率が高く、偽陽性率が低い）ほど良いモデルであると言えます。

デフォルトの閾値ではあまり差がなかったロジスティック回帰とランダムフォレストですが、ランダムフォレストにはより良い閾値がありその閾値の元でロジスティック回帰よりも良い結果をもたらすことがPR曲線・ROC曲線の描画により初めて分かりました。

AUC

ROC曲線は左上に広がるほど良いモデルであることから、モデルの有効な指標としてROC曲線の下の領域の面積(0から1)が使われ、この指標はAUCスコアと呼ばれます。 AUCスコアはroc_auc_score関数により算出できます。

from sklearn.metrics import roc_auc_score

print('AUC score(logreg):', roc_auc_score(y_test, logreg_clf.predict_proba(X_test)[:, 1]))
print('AUC score(forest):', roc_auc_score(y_test, forest_clf.predict_proba(X_test)[:, 1]))

AUCスコアにより、閾値を考慮したモデルの性能を捉えられていることがわかります。

閾値の選択

ランダムフォレストにはより良い閾値の選択の余地がありそうだと分かりました。PR曲線から適合率を下げずに0.80程度まで再現率を増やすことができると読み取れます。

いくつかの閾値について、PR曲線にプロットを追加してみましょう。

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)

precision_forest, recall_forest, thresholds_forest = precision_recall_curve(y_test, forest_clf.predict_proba(X_test)[:, 1])
threshold_125_forest = np.argmin(np.abs(thresholds_forest - 0.125))
threshold_22_forest = np.argmin(np.abs(thresholds_forest - 0.22))
threshold_36_forest = np.argmin(np.abs(thresholds_forest - 0.36))
default_threshold_forest = np.argmin(np.abs(thresholds_forest - 0.5))
ax.plot(recall_forest, precision_forest, label='RandomForestClassifier')
ax.plot(recall_forest[default_threshold_forest], precision_forest[default_threshold_forest],
        'x', c='k', markersize=10, label='threshold 50%(default)')
ax.plot(recall_forest[threshold_36_forest], precision_forest[threshold_36_forest],
        'o', c='k', markersize=10, fillstyle="none", label='threshold 36%')
ax.plot(recall_forest[threshold_22_forest], precision_forest[threshold_22_forest],
        '^', c='k', markersize=10, fillstyle="none", label='threshold 22%')
ax.plot(recall_forest[threshold_10_forest], precision_forest[threshold_125_forest],
        's', c='k', markersize=10, fillstyle="none", label='threshold 12.5%(min)')
ax.legend()
ax.set_ylim([0, 1.1])
ax.set_title('PR Curve', fontsize=14)
ax.set_xlabel('recall')
ax.set_ylabel('precision')

plt.show()

PR曲線の閾値の最小は必ずしもゼロでないことに注意が必要です。

thresholds_forest.min()

適合率と再現率のトレードオフの選択は、冒頭に述べた通りビジネス上の影響を考慮して決定する必要があります。

上記のプロットからは、「8」を必ず正しく予測することが重要なタスクでは閾値0.36を、「8」の1つ残らず予測することが必要な(予測漏れがない)タスクでは閾値0.125を、バランスをとって閾値0.22といった閾値が考えられます。

例: 閾値を0.36に設定した場合の混同行列と分類レポートの表示

これはテストデータで閾値を決定し、テストデータで評価を行っている悪い例ですが、閾値の変更により指標が変わる様子を示すために記載しておきます。

pred_forest_36 = forest_clf.predict_proba(X_test)[:, 1] >= 0.36

matrix = confusion_matrix(y_test, pred_forest_36)
sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets)

print(classification_report(y_test, pred_forest_36,
                            target_names=targets))

適合率は1.00のまま、再現率が0.68から0.82に向上していることがわかります。適切な閾値の設定により精度が向上する一例と言えそうです。

テストデータの使用に関する注意点

モデルの評価の際には、基本的に「テストデータを使って性能評価を行うのは最後の1回だけ」を意識しておかないといけません。

テストデータの予測結果を基に描画したPR曲線・ROC曲線からモデルの閾値を決定すれば、それはテストデータに合わせた閾値となってしまうため楽観的な精度が得られます。

実際にはテストデータとは別に検証データを用意し、検証データに対するPR曲線・ROC曲線を基に閾値の選択を行いましょう。

まとめ

書きたいこと書いてると記事が長くなってしまいましたね。

PR曲線やROC曲線は2値分類タスクにおいて、ビジネス上の誤判定の影響を考慮する際には、不均衡データを扱う場合に限らず使用されます。

閾値の変更を行う際は、ハイパーパラメータの決定を行う時と同様に検証データを使いましょう。初心者のうちは、テストデータの情報を使ってまやかしの精度を叩き出してしまわないよう常に気を配らないとやらかします。（自分もよくやらかします）

では、いってきます🐼