Python scikit-learnのSVMで糖尿病データセットを分類する

ピマ・インディアンの糖尿病のデータセットを使って、SVMで糖尿病か否かを分類します。 合わせて、欠測値の補完や、SVMに欠かせないデータの標準化、グリッドサーチによるハイパーパラメータの探索も行っています。

ピマ・インディアンの糖尿病データセット

今回はピマ・インディアンの糖尿病データセットを使います。 糖尿病のデータセットとしては、scikit-learn付属のもの(sklearn.datasets.load_diabetes)も書籍等でよく見かけますが、それとは別でUCI(カリフォルニア大学アーバイン校)から提供されており、実際にサンプルされたデータです。

ピマ・インディアンの女性の症例768サンプル(内糖尿病と診断されたのは268サンプル)で、各サンプルは8つの計測値(年齢、妊娠回数、BMIなど)と1つの診断結果(糖尿病か否か)の値を持ちます。

UCI Machine Learning Repository: Pima Indians Diabetes Data Set

今回の機械学習とは全く関係ないのですが、ピマ・インディアンについて少し調べてみました。

アメリカアリゾナ州に定住しているピマ・インディアンは、1900年代後半から食生活の変容によって肥満と2型糖尿病が急増しました。 その後の研究で、エネルギー代謝のコントロール分子であるβ3アドレナリン受容体が発見され、ピマ・インディアンのケースはこのβ3アドレナリン受容体の遺伝子変異の発現によるものとされています(ちなみに日本人も同様の変異を持ちやすいようです)。

おそらくこうした背景もあって、まとまったデータセットとして採取・提供されいているのでしょう。

データセットのロード

データセットCSVファイルとして提供されていますので、pandas.read_csvでDataFrameにロードします。

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
from sklearn.metrics import confusion_matrix
from sklearn.preprocessing import Imputer
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV

diabetes_df = pd.read_csv(
    'http://archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data', 
    header=None)

ロードすると9列×768行のDataFrameに展開されます。 各列の意味は以下の通りで、8つの説明変数(0〜7)と1つの目的変数(8)になっています。

  • 0 Number of times pregnant: 妊娠回数
  • 1 Plasma glucose concentration a 2 hours in an oral glucose tolerance test: 血糖濃度(経口ブドウ糖負荷試験後2時間の値)
  • 2 Diastolic blood pressure (mm Hg): 最低血圧(mm/Hg)
  • 3 Triceps skin fold thickness (mm): 上腕三頭筋皮下脂肪の厚さ(mm)
  • 4 2-Hour serum insulin (mu U/ml): 血清インスリン濃度(経口ブドウ糖負荷試験後2時間の値)
  • 5 Body mass index (weight in kg/(height in m)2): BMI
  • 6 Diabetes pedigree function: 糖尿病血統要因
  • 7 Age (years): 年齢
  • 8 Class variable (0 or 1): 糖尿病が陽性ならば1(目的変数)

説明変数(X)と目的変数(y)に分離しておきます。

X = diabetes_df.iloc[:, :8]
y = diabetes_df.iloc[:, 8:].values.flatten() # 1次元に展開

print('X shape: {}, y shape: {}'.format(X.shape, y.shape))
# X shape: (768, 8), y shape: (768,)

SVMで学習・テスト

学習アルゴリズムとして、今回はサポートベクタマシン(SVM)を使います。

SVMは、各データがプロットされた空間(この場合は8次元空間)に超平面の境界(決定境界)を設けて陰性のデータと陽性のデータを分離しますが、このとき各データと決定境界のマージンが最大となるようにする学習アルゴリズムです。 こちらの方のスライドがとてもわかり易いです。

https://www.slideshare.net/mknh1122/svm-13623887

まずは何も手を加えていない状態のデータを使って、SVMで学習・テストさせてみます。

結果はやはり良くなく、分類精度は学習データで100%だったのに対してテストデータでは68%となり、典型的な過学習に陥ってました。 混同行列でテストデータの分類結果を見ると、テストサンプル292を全て0(陰性)で分類していることがわかります。

# 学習データとテストデータの分離
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

# 学習とテスト
svc = SVC()
svc.fit(X_train, y_train)

print('Train score: {:.3f}'.format(svc.score(X_train, y_train)))
print('Test score: {:.3f}'.format(svc.score(X_test, y_test)))
# Train score: 1.000
# Test score: 0.677
print('Confusion matrix:\n{}'.format(confusion_matrix(y_test, svc.predict(X_test))))
# Confusion matrix:
# [[130   0]
# [ 62   0]]

欠測値の補完と標準化

まずはデータをpandas.tools.plotting.scatter_matrixで、各変数毎の偏りや変数間の相関性などを一覧にしてみます。

pd.tools.plotting.scatter_matrix(diabetes_df, alpha=0.3, figsize=(16,16))
plt.show()

いくつかのことが言えるかと思います。

  • 説明変数1〜5に値が0となるデータが含まれている
    • いずれの変数も0が計測されるとは考えられないため、これらはおそらく欠測値
    • Data Set DescriptionでもMissing Attribute Values: Yesと説明されています
  • 目的変数に大きな影響を与える入力変数は無さそう(8行目または8列目の分布図)
    • 説明変数1〜6では陽性の方が陰性よりも少し高いか、という程度の差

欠測値を含むデータについて除外することもできますが、そうなると大半のデータを失ってしまいますので、別の値で補完する方法を考えます。 双峰性こそ見られませんが、説明変数4(血漿中インスリン濃度)のように正規分布(釣り鐘状)になっていない値もありますので、平均値ではなく、中央値を使うこととします。

欠測値の補完にはsklearn.preprocessing.Imputerを使います。 missing_valueで保管される値を0とし、strategyでmedian(中央値)を指定しています。

med_imp = Imputer(missing_values=0, strategy='median', axis=0)
med_imp.fit(X.iloc[:, 1:6])
X.iloc[:, 1:6] = med_imp.transform(X.iloc[:, 1:6])

また、SVMは特徴量のスケールが揃っていないと上手く学習できないアルゴリズムなので、正規化(全ての値を0〜1に収まるように変換)または標準化(平均0、分散1となるように変換)が必要です。 説明変数4(血漿中インスリン濃度)や6(糖尿病血統要因)では外れ値とも言えるような値も見られますので、外れ値の影響を受けやすい正規化ではなく、標準化を選択します。

値の標準化には、sklearn.preprocessing.StandardScalerを使います。 他の学習アルゴリズムなどと同様、学習データのみで学習(fit)させてから、変換(transform)する必要があります。

# 学習データとテストデータの分離
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

# 標準化
std_scl = StandardScaler()
std_scl.fit(X_train)
X_train = std_scl.transform(X_train)
X_test = std_scl.transform(X_test)

欠測値の補完と標準化を終えたところで、再度SVMで学習・テストさせてみますと、テストデータで分類精度76%となり、62ケースの内30ケースは陽性と判断できていました。

svc = SVC()
svc.fit(X_train, y_train)

print('Train score: {:.3f}'.format(svc.score(X_train, y_train)))
print('Test score: {:.3f}'.format(svc.score(X_test, y_test)))
# Train score: 0.825
# Test score: 0.755
print('Confusion matrix:\n{}'.format(confusion_matrix(y_test, svc.predict(X_test))))
# Confusion matrix:
# [[115  15]
#  [ 32  30]]

ハイパーパラメータの選択

SVMには、インスタンス生成時に指定できるオプションがいくつかありますが、そのうち学習に大きく影響を与えるのがCとgammaで、これらはハイパーパラメータとして問題に応じて調整する必要があります。

  • Cは、正則化のためのパラメータで、値を大きくすれば誤分類によるペナルティは大きくなります(=正則化が弱い複雑な決定境界になります)
  • gammaは、説明変数の高次元化に使われるRBFカーネル(ガウスカーネル)のパラメータで、値を大きくすると決定境界はより複雑になります

forループなどでパラメータを1つずつ変えて最高スコアとなるパラメータを見つけ出すこともできますが、データの分離(ハイパーパラメータの選択では、学習データとテストデータに加えて、ハイパーパラメータの評価用のデータも別に用意する必要がある)や交差検証も合わせて行ってくれるsklearn.model_selection.GridSearchCVを使うほうが手軽です。

今回は、Cとgammaでそれぞれ0.001, 0.01, 0.1, 1, 10, 100の6つの値をとるように設定(合計6×6=36パターン)します。 また10分割の交差検証で評価するようにしています。

結果として、C=1およびgamma=0.01の場合でテストデータの分類精度は78%となり、陽性62ケース中33ケースを正しく分類できるように改善されたことがわかります。

svc_param_grid = {
    'C': [0.001, 0.01, 0.1, 1, 10, 100],
    'gamma': [0.001, 0.01, 0.1, 1, 10, 100]
}

svc_grid_search = GridSearchCV(SVC(), svc_param_grid, cv=10)
svc_grid_search.fit(X_train, y_train)

print('Train score: {:.3f}'.format(svc_grid_search.score(X_train, y_train)))
print('Test score: {:.3f}'.format(svc_grid_search.score(X_test, y_test)))
# Train score: 0.764
# Test score: 0.781
print('Confusion matrix:\n{}'.format(confusion_matrix(y_test, svc_grid_search.predict(X_test))))
# Confusion matrix:
# [[117  13]
#  [ 29  33]]
print('Best parameters: {}'.format(svc_grid_search.best_params_))
print('Best estimator: {}'.format(svc_grid_search.best_estimator_))
# Best parameters: {'C': 1, 'gamma': 0.01}
# Best estimator: SVC(C=1, cache_size=200, class_weight=None, coef0=0.0,
#   decision_function_shape=None, degree=3, gamma=0.01, kernel='rbf',
#   max_iter=-1, probability=False, random_state=None, shrinking=True,
#   tol=0.001, verbose=False)