Python scikit-learnのランダムフォレストで受診予約のNo-Showを予測する

Kaggleのデータセットを使って、ランダムフォレストで受診予約のNo-Showを予測します。

データセットのロード

今回はKaggleで公開されているMedical Appointment No Showsを使っていきます。 このデータは、受診予約で1レコードとなっており、患者の情報(年齢・性別やかかっている病気など)や予約どおりに現れたか・現れなかったか、などが15の含まれています。 全体では30万レコードになります。

このデータをcsvとしてダウンロードし、DataFrameに読み込ませます。

# ライブラリのインポート
import pandas as pd
import matplotlib.pyplot as pyplot
% matplotlib inline
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import f1_score
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import make_scorer
import seaborn

# CSVファイルからDataFrameへロード
original_df = pd.read_csv('No-show-Issue-Comma-300k.csv')
original_df.head(6)

以下のデータが含まれています。 typoが目立ちます。

  • 患者の情報
    • 年齢(Age)、性別(Gender) 、ハンディキャップ有無(Handcap)、奨学金有無(Scholarship)
    • アルコール依存(Alcoolism)、喫煙(Smokes)、高血圧(HyperTension)などの、生活習慣に関わる情報
    • 糖尿病(Diabetes)、結核(Tuberculosis)などの、病気に関わる情報
  • 予約の情報
    • 予約登録日時(AppointmentRegistration)、予約日(ApointmentData)、予約日の曜日(DayOfTheWeek)、SMSリマインド?(Sms_Reminder)、予約日と予約登録日までの日数(AwaitingTime)
    • 今回の目的変数となる予約に現れたかどうかの情報はStatusに入っており、現れた場合は'Show-Up'、現れなかった場合は'No-Show'

Show-Upが209,269件、No-Showが90,731件で、概ね7:3の比率です。

original_df['Status'].value_counts()
#Show-Up    209269
#No-Show     90731
#Name: Status, dtype: int64

各データの分布を確認します。 Ageがマイナスになっていたり、2値と思われていたHandcapやSms_Reminderが2以上の値を含んでいたり、若干のデータ不備があるようです。

original_df.hist(figsize=(12, 12))

念のため、月ごとに予約数とNo-Show数に経時的な変化が無いか確認しましたが、特に見られませんでした。

# 経時的な変化は無いか? → 無さそう(12月だけちょっと多いか)
appointments = pd.DataFrame(index=original_df.index)
appointments['AppointmentDate'] = pd.to_datetime(original_df['ApointmentData']).apply(lambda a: '{}/{:02}'.format(a.year, a.month))
appointments['NoShow'] = original_df['Status'].apply(lambda d: 1 if d == 'No-Show' else 0)
temp = pd.DataFrame(appointments.groupby('AppointmentDate').count())
temp['All'] = appointments.groupby('AppointmentDate').count()
temp['NoShow'] = appointments.groupby('AppointmentDate').sum()
temp.reset_index(inplace=True)
plt.plot(temp.index, temp['All'])
plt.plot(temp.index, temp['NoShow'])
plt.show()

f:id:ohke:20170730103451p:plain

目的変数と説明変数の抽出

ロードしたDataFrameから目的変数・説明変数を抽出します。

まずは元データをほとんどそのままコピーして作ります(日時データのAppointmentRegistrationとApointmentDataを除いています)。 Genderの2値化(Mの場合に1とする)、DayOffTheWeekのone-hot-encoding、および、typoの修正も行います。

  • 例えば、DayOffTheWeekがMondayの場合は、AppointmentMondayが1になり、それ以外は0となるように列を分解します(one-hot-encoding)
features_df = pd.DataFrame()

# 目的変数の抽出(No-Showなら1)
features_df['Outcome'] = original_df['Status'].apply(lambda s: 1 if s == 'No-Show' else 0)

# 元データを説明変数に追加(typoも同時に修正する)
features_df['Age'] = original_df['Age']
features_df['Male'] = original_df['Gender'].apply(lambda g: 1 if g == 'M' else 0) # 2値変数化
features_df['Diabetes'] = original_df['Diabetes']
features_df['Alcoholism'] = original_df['Alcoolism']
features_df['HiperTension'] = original_df['HiperTension']
features_df['Handicap'] = original_df['Handcap']
features_df['Scholarship'] = original_df['Scholarship']
features_df['Smokes'] = original_df['Smokes']
features_df['SmsReminder'] = original_df['Sms_Reminder']
features_df['Tuberculosis'] = original_df['Tuberculosis']
features_df['AwaitingTime'] = original_df['AwaitingTime']

# 予約日の曜日をone-hot-encoding
d = pd.get_dummies(original_df['DayOfTheWeek'])
features_df['AppointmentMonday'] = d['Monday']
features_df['AppointmentTuesday'] = d['Tuesday']
features_df['AppointmentWednesday'] = d['Wednesday']
features_df['AppointmentThursday'] = d['Thursday']
features_df['AppointmentFriday'] = d['Friday']
features_df['AppointmentSaturday'] = d['Saturday']
features_df['AppointmentSunday'] = d['Sunday']

各列の相関係数を求めます。

相関係数pandas.DataFrame.corr()で取得できるので、seaborn.heatmap()で可視化します。

  • 最も相関が高いのは年齢(Age)で、年齢が高いほうがNo-Showは少ないようです(負)が、それ以外に有力な説明変数は無さそうです
  • 今回の分析とは関係ないですが、年齢と高血圧・糖尿病、喫煙とアルコール依存症も相関が高いことが伺えます(このあたりは通説と概ねマッチしてますね)
pyplot.figure(figsize=(15,15))
seaborn.heatmap(features_df.corr(), annot=True)

ランダムフォレストでの予測(1回目)

この時点の特徴量を使って、ランダムフォレストで予測させてみます。

ランダムフォレストは、複数の決定木で多数決することで分類する方法です。 各決定木はランダムに復元抽出されたサンプルと、同じくランダムに選択された特徴量を使って学習を行うため、それぞれ異なった決定木となります。

scikit-learnではRandomForestClassifierが提供されており、決定木の個数や各決定木で使う特徴量の個数、各決定木の深さなどを指定できますが、まずはデフォルトで学習・テストさせてみます。 ランダムフォレストでは、データの正規化・標準化を考える必要が無いので、今回のように連続値変数(AgeやAwaitingTime)と2値変数(その他)が混在しているケースでも、簡単に試すことができる利点があります。

結果の分布に偏りがあるため、今回はスコアだけではなく、F値も見ます(No-Showは全体の3割なので、全て0と予測すると、それだけでスコアは0.7になってしまうためです)。 F1値は適合率(あるラベルに分類されたデータの内、正しく分類された検証データの割合)と再現率(あるラベルに分類されるべき検証データの内、そのラベルに分類された割合)の調和平均で、scikit-learnではsklearn.metrics.f1_scoreを使います。

結果を見てみると、学習データではスコア0.80だったのに対して、テストデータでは0.65に留まっており、過学習の傾向があります。 混合行列でも、TP(右下の値、正しくNo-Showと予測した数)が4797で、全体の21%程度しか予測できていません。

# 説明変数と目的変数の分離
X = features_df.ix[:, 'Age':]
y = features_df['Outcome']

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

# ランダムフォレストの作成
forest = RandomForestClassifier(min_samples_leaf=3, random_state=0)
forest.fit(X_train, y_train)

# 評価
print('Train score: {}'.format(forest.score(X_train, y_train)))
print('Test score: {}'.format(forest.score(X_test, y_test)))
print('Confusion matrix:\n{}'.format(confusion_matrix(y_test, forest.predict(X_test))))
print('f1 score: {:.3f}'.format(f1_score(y_test, forest.predict(X_test))))
# Train score: 0.804
# Test score: 0.647
# Confusion matrix:
# [[43714  8633]
#  [17856  4797]]
# f1 score: 0.266

説明変数の追加

以下の仮説に基いて、説明変数を8つ追加してみます。 予約登録した曜日・時間帯に予約登録した人の方が、そうでない人よりもNo-Showは低いだろうという推測に基づいています。

  • 予約登録した曜日(RegistrationMonday〜RegistrationSunday, 7変数)
    • 予約登録した曜日が1(それ以外は0)
  • 予約登録した時間帯(RegistrationWorktime)
    • 営業時間帯(9時-17時と仮定)に予約登録していると1(それ以外は0)
regs = pd.to_datetime(original_df['AppointmentRegistration'])
features_df['RegistrationMonday'] = regs.apply(lambda r: 1 if r.date().weekday() == 0 else 0)
features_df['RegistrationTuesday'] = regs.apply(lambda r: 1 if r.date().weekday() == 1 else 0)
features_df['RegistrationWednesday'] = regs.apply(lambda r: 1 if r.date().weekday() == 2 else 0)
features_df['RegistrationThursday'] = regs.apply(lambda r: 1 if r.date().weekday() == 3 else 0)
features_df['RegistrationFriday'] = regs.apply(lambda r: 1 if r.date().weekday() == 4 else 0)
features_df['RegistrationSaturday'] = regs.apply(lambda r: 1 if r.date().weekday() == 5 else 0)
features_df['RegistrationSunday'] = regs.apply(lambda r: 1 if r.date().weekday() == 6 else 0)

features_df['RegistrationWorktime'] = regs.apply(lambda r: 1 if r.hour >= 9 and r.hour < 17 else 0)

もう一度、相関係数を表示させてみます。 追加した変数も決定的とは言えなさそうです。

重要度の測定と説明変数の削減

再度、randomForestClassifierで学習・テストさせてみますと、全体的に若干の改善が見られます。

X = features_df.ix[:, 'Age':]
y = features_df['Outcome']

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

forest = RandomForestClassifier(random_state=0)
forest.fit(X_train, y_train)
print('Train score: {:.3f}'.format(forest.score(X_train, y_train)))
print('Test score: {:.3f}'.format(forest.score(X_test, y_test)))
print('Confusion matrix:\n{}'.format(confusion_matrix(y_test, forest.predict(X_test))))
print('f1 score: {:.3f}'.format(f1_score(y_test, forest.predict(X_test))))
# Train score: 0.829
# Test score: 0.641
# Confusion matrix:
# [[42923  9424]
#  [17528  5125]]
# f1 score: 0.276

RandomForestClassifierは、feature_importances_というプロパティを持っており、各説明変数の重要度を持っています。 棒グラフで可視化してみると重要度が高い方からAge、AwaitingTime、Male、・・・と順に続いています。

values, names = zip(*sorted(zip(forest.feature_importances_, X.columns)))

pyplot.figure(figsize=(12,12))
pyplot.barh(range(len(names)), values, align='center')
pyplot.yticks(range(len(names)), names)

下5つの変数についてはほぼ影響を与えないようなので、削除します。

features_df.drop(['RegistrationSunday', 'AppointmentSunday', 'Tuberculosis', 'RegistrationSaturday', 'AppointmentSaturday'], axis=1, inplace=True)

ハイパーパラメータの探索

最後にグリッドサーチでハイパーパラメータを探索します。

  • n_estimators: 決定木の数(デフォルトは10)を設定しますが、基本的に多ければ多いほど良くなる(=計算時間とのトレードオフ)なので、今回は10に固定します
  • max_features: 各決定木で分類に使用する説明変数の数で、今回は4パターン(1, ‘auto’=全ての説明変数の数の2乗根、None=全ての説明変数の数と同じ)にします
  • max_depth: 各決定木の深さを表し、深ければ深いほど複雑な分岐になります(=過学習を起こしやすい)
  • min_samples_leaf: 決定木の葉に分類されるサンプル数を決めるパラメータで、3パターン(1, 2, 4)で試します

今回は、スコア方法をF1にします。 GridSearchCVのscoringオプションで変更できます。 またcv=4として、4グループずつに分けて交差検証します。

# ハイパーパラメータ
forest_grid_param = {
    'n_estimators': [100],
    'max_features': [1, 'auto', None],
    'max_depth': [1, 5, 10, None],
    'min_samples_leaf': [1, 2, 4,]
}

# スコア方法をF1に設定
f1_scoring = make_scorer(f1_score,  pos_label=1)

# グリッドサーチで学習
forest_grid_search = GridSearchCV(RandomForestClassifier(random_state=0, n_jobs=-1), forest_grid_param, scoring=f1_scoring, cv=4)
forest_grid_search.fit(X_train, y_train)

# 結果
print('Best parameters: {}'.format(forest_grid_search.best_params_))
print('Best score: {:.3f}'.format(forest_grid_search.best_score_))
# Best parameters: {'max_depth': None, 'max_features': 1, 'min_samples_leaf': 1, 'n_estimators': 10}
# Best score: 0.276

得られたベストパラメータで学習・テストさせてみましたが、デフォルトとほぼ変わらない結果となりました。 これ以上の精度を目指すなら、説明変数を増やす必要があるかもしれません。

X = features_df.ix[:, 'Age':]
y = features_df['Outcome']

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

best_params = forest_grid_search.best_params_

forest = RandomForestClassifier(random_state=0, n_jobs=-1, 
                                max_depth=best_params['max_depth'], 
                                max_features=best_params['max_features'], 
                                min_samples_leaf=best_params['min_samples_leaf'],
                                n_estimators=best_params['n_estimators'])
forest.fit(X_train, y_train)
print('Train score: {:.3f}'.format(forest.score(X_train, y_train)))
print('Test score: {:.3f}'.format(forest.score(X_test, y_test)))
print('Confusion matrix:\n{}'.format(confusion_matrix(y_test, forest.predict(X_test))))
print('f1 score: {:.3f}'.format(f1_score(y_test, forest.predict(X_test))))
# Train score: 0.828
# Test score: 0.638
# Confusion matrix:
# [[42551  9796]
#  [17368  5285]]
# f1 score: 0.280

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)

Python 回帰木でセッション数を予測するモデルを作成する

前回の投稿では線形回帰を使ってセッション数を予測しましたが、今回は回帰木を使ってみます。

Python GoogleAnalyticsのデータを使って線形回帰でセッション数を予測するモデルを作る - け日記

回帰木による学習・テスト

前回の投稿では、本ブログの1日あたりのセッション数(来訪回数)を予測するために、2つの説明変数(4/1を0として何日目かを表すnth、休日を表すholiday)を2次の多項式に拡張して、線形回帰で学習・テストさせました。

今回使う回帰木は、説明変数に閾値を設けて学習サンプルを分割していき(多くの実装は二分割)、リーフでの目的変数の平均値をそのリーフに分類されたサンプルの予測値とするものです。

散布図を見てみますと、概ね4つのグループ(四角枠)に分けられそうです。 回帰木で学習すると、例えばholidayが0ならばAグループ、0ならばBグループ、Aグループの内nthが40以下ならA-1グループ、41以上ならA-2グループ、・・・というように、今ある説明変数だけでもざっくりと分割できると予想されます。

scikit-learnを使った回帰木の学習では、sklearn.tree.DecisionTreeRegressorを使います。

  • 学習データとテストデータを3:1で分割しています
  • 回帰木の深さが深いほど過学習を起こしやすくなるため、最大深さ3としています
    • 4グループに分割できればいいので、二分木であれば深さ2までで良さそうですが、おそらく最初にholidayを使って平日(上部3グループ)と休日(下部1グループ)に分けられるので、平日の3グループをさらに分割するためには深さ2では足りないと予想されたため、3にしました
  • 決定木は閾値による分類のため、説明変数の正規化・標準化を必要ありません
from sklearn.tree import DecisionTreeRegressor

# features_dfの取得については前回の記事を参照
X = features_df[['nth', 'holiday']]
y = features_df['sessions']
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

# 最大深さ3の回帰木を作成
regressor = DecisionTreeRegressor(max_depth=3)

# 学習・テスト
regressor.fit(X_train, y_train)
print('Train score: {:.3f}'.format(regressor.score(X_train, y_train)))
print('Test score: {:.3f}'.format(regressor.score(X_test, y_test)))
# Train score: 0.949
# Test score: 0.862

# プロット
plt.scatter(X['nth'], y)
plt.plot(X['nth'], pipeline.predict(X))
plt.show()

結果として決定係数R2は、学習データで0.949、テストデータで0.862となり、前回の線形回帰を使ったモデル(学習データで0.916、テストデータで0.832)よりも少し改善されました。 また、プロットした図でも、概ね最初に予想した4グループに分割されていることが確認できます。

ただし、回帰木では長期的な増加傾向を反映できていないので、7月以降のセッション数をこれを使って予測すると、おそらく線形回帰よりも性能は悪くなると予想されます。

f:id:ohke:20170720081901p:plain

回帰木を描画

最後に、学習で獲得した回帰木を描画し、どういった条件でサンプルを分割させているのかを見てみます。 回帰木の描画ではdotファイルを出力し、graphvizで画像ファイル(png)に変換します。

Graphviz | Graphviz - Graph Visualization Software

まずgraphvizをインストールしておきます。

$ brew install graphviz

Pythonからはsklearn.tree.export_graphvizをインポートして、export_graphvizでdotファイルを出力します。

from sklearn.tree import export_graphviz

export_graphviz(regressor, out_file='tree.dot', feature_names=['nth', 'holiday'])

最後に、graphvizのdotコマンドでpngファイルへ出力します。

$ dot -Tpng tree.dot -o tree.png

以下のような回帰木が図として出力されます。 valueがそのグループの平均値で、mseが平均二乗誤差です。

当初の見立て通り、最初にholidayを閾値として平日と休日のグループに分け、以降はnthの値を閾値として分割していることがわかります。 例えば、4/3(nth=3, holiday=0)のセッション数は、一番左下のノードに分類されて74.3と予測されます。

f:id:ohke:20170720082249p:plain