動画や本でデータサイエンスを学んだ後、実際のデータで実践したくなりました。そんな時に出会ったのがKaggleでした。Kaggleではデータサイエンスコンペが開催され、日夜データサイエンティストがしのぎを削っています。
初心者には参加しづらく思えますが、コンペの中には初者のために常時開催されているものもあります。そのひとつが、今回扱うタイタニックのデータセットです。
タイタニックのデータセットには、氏名、年齢、性別、等級、チケット料金等の乗客に関する情報と、乗客の生死の情報(特徴量)が含まれています。これらの情報はどれも理解しやすく、データ量も多くないため、データサイエンス初心者にとってオススメのデータセットと言われています(特徴量同士が相関してしまっている問題はある)。このコンペでは、特徴量から乗客の生死を高精度で予測するモデルを構築することが目的となります。
ライブラリのインポートとデータセットの確認
まずライブラリ一式をインポートします。
#In import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns %matplotlib inline from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA from sklearn.pipeline import make_pipeline from sklearn.model_selection import cross_val_score from sklearn.model_selection import StratifiedKFold from sklearn.linear_model import LogisticRegression from sklearn.tree import DecisionTreeClassifier from sklearn.svm import SVC from sklearn.naive_bayes import GaussianNB from sklearn.ensemble import RandomForestClassifier, VotingClassifier from sklearn.neighbors import KNeighborsClassifier from sklearn.model_selection import GridSearchCV from sklearn.metrics import classification_report, confusion_matrix, roc_curve, roc_auc_score from sklearn.model_selection import learning_curve
面倒なので、Kaggleのトレーニングデータだけを利用することにしました。データフィールドの形式を見たところ、objectがいくつかあります。scikit-learnの多くのモデルは、文字データや欠損データを扱うことができないので、変換や削除する必要があります。
#In
titanic = pd.read_csv('titanic.csv')
titanic.info()
#Out
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 891 entries, 0 to 890
Data columns (total 12 columns):
PassengerId 891 non-null int64
Survived 891 non-null int64
Pclass 891 non-null int64
Name 891 non-null object
Sex 891 non-null object
Age 714 non-null float64
SibSp 891 non-null int64
Parch 891 non-null int64
Ticket 891 non-null object
Fare 891 non-null float64
Cabin 204 non-null object
Embarked 889 non-null object
dtypes: float64(2), int64(5), object(5)
生存に関係しそうな’Age’がたくさん欠損しており、推測値で補完する必要がありそうです。’Cabin’はほとんど欠損しており、この特徴量から生存を推測するのは難しそうです。
#In sns.heatmap(titanic.isnull(), yticklabels=False, cbar=False, cmap='viridis')
文字データや欠損データの補完と削除
まずは、’Cabin’を列ごと削除しました。本当はこうしたデータ操作をする前に、トレーニングデータとテストデータを分けておき、テストデータが完全に未知の状態を用意するのが良いのですが、この先タイタニックに新データが発生するとも思えませんし、めんどくさがりなので許してください。
#In
titanic = titanic.drop('Cabin', axis=1)
次に、’Age’の欠損を補完しなければなりません。方法としては、平均値/中央値/最頻値/ランダム値等、色々あります。しかしながらタイタニックのデータセットでは、生存者の方が性別が明らかな場合が多く、実際に乗客全体の生存率は38%であるのに対し、’Age’が不明な乗客の生存率は29%に止まります。平均値/中央値/最頻値/ランダム値等は、欠損データが規則性なく欠損していることを前提としているため、今回適用するのは良くないかもしれません。データを歪めてしまい、モデルの精度が落ちてしまうと思います。
#In (titanic['Survived'].mean(), titanic[titanic['Age'].isnull()]['Survived'].mean()) #Out (0.3838383838383838, 0.2937853107344633)
そこで、’Name’についている乗客の肩書きから、’Age’の欠損を補完しようと思います。’Name’の構成を見たところ、幸い’,’と’.’でうまく肩書きのみを分離できそうです。
#In pd.DataFrame(titanic.head()['Name']) #Out 0 Braund, Mr. Owen Harris 1 Cumings, Mrs. John Bradley (Florence Briggs Th... 2 Heikkinen, Miss. Laina 3 Futrelle, Mrs. Jacques Heath (Lily May Peel) 4 Allen, Mr. William Henry
肩書きを’Title’として分離後、数え上げました。主要4カテゴリ以外は、マイナーな肩書きであることがわかりました。
#In
#'Name'を','で分離し、分離後の文字データをさらに'.'で分離
titanic['Title'] = titanic['Name'].str.split(', ', expand=True)[1].str.split('.', expand=True)[0]
sns.countplot(y=titanic['Title'])
頻度に応じて’Title’へカテゴリ変数を付与しました。関数がイケてないので、リスト内包表記をdictionaryへ応用しようと思ったのですが、’else: return 5’の部分をうまく表現できず、挫折しました。
#In
#頻度に応じて、'Title'へカテゴリ変数を付与
def title_int(x):
if x == 'Mr':
return 1
elif x == 'Miss':
return 2
elif x == 'Mrs':
return 3
elif x == 'Master':
return 4
else:
return 5
titanic['Title_Int'] = titanic['Title'].apply(lambda x:title_int(x))
‘Title’に応じて、’Age’の欠損を補完しました。’values’の中を見てみると、’Title’相応のそれっぽい年齢で補完できてそうです。また、欠損データがなくなったことにより、図表から黄色の部分がなくなりました。
#In
#キーを作成
keys = [1,2,3,4,5]
#'Age'が欠損していないデータから、'Title'ごとの年齢の平均値を取得
values = np.array(titanic.groupby('Title_Int').mean()['Age'])
#キーと'Title'ごとの年齢の平均値からdictionaryを作成
dic = dict(zip(keys,values))
#'Age'が欠損しているデータへ、dictionaryをマッピング
titanic['Age'] = titanic['Age'].fillna(titanic['Title_Int'].map(dic))
print(values)
[32.36809045, 21.7739726 , 35.89814815, 4.57416667, 42.38461538]
sns.heatmap(titanic.isnull(), yticklabels=False, cbar=False, cmap='viridis')
さらに、文字データを数値データに置換する必要があります。文字データには、’Name’、’Sex’、’Ticket’、’Embarked’があります。’Sex’をカテゴリ変数に置換しました。
#In #'Sex'をカテゴリ変数に置換 titanic['Sex'] = titanic['Sex'].apply(lambda x:0 if x=='male' else 1) titanic['Sex'].unique() #Out array([0, 1])
‘Sex’以外は定量化できなさそうなため、これらは’Cabin’と同様、列ごと削除しました。特徴量が少ないため、特徴量削減(主成分分析/次元削減/決定木・ランダムフォレストによる特徴量の重要度分析等)はやめておきました。ようやくきれいなデータになりました。
#In #データセットから特徴量行列を分離 drop_col = ['PassengerId', 'Name', 'SibSp', 'Parch', 'Ticket', 'Embarked', 'Title', 'Survived'] titanic_features = titanic.drop(drop_col, axis=1) titanic_features.info() #Out <class 'pandas.core.frame.DataFrame'> RangeIndex: 891 entries, 0 to 890 Data columns (total 5 columns): Pclass 891 non-null int64 Sex 891 non-null int64 Age 891 non-null float64 Fare 891 non-null float64 Title_Int 891 non-null int64 dtypes: float64(2), int64(3)
最後に、ターゲットベクトルである’Survived’を分離しました。
#In #データセットからターゲットベクトルを分離 titanic_target = titanic['Survived']
スケーリング
今回使うモデルは、理解しやすいロジスティック回帰にします。このモデルは、特徴量間のスケールに開きがあると、より大きな特徴量に影響され、過学習しやすくなります。例えば、’Age’は0~80の値をとる一方、’Fare’は0~512の値をとるため、’Fare’の影響が大きくなってしまいます。
特徴量間のスケールを揃えるのがスケーリングの目的です。スケーリングにより、各特徴量は平均値0、標準偏差1になるよう変換されます。スケーリング時の注意すべき点として、トレーニングデータで行なったスケーリングをそのままテストデータに適用することが挙げられます。トレーニングデータとテストデータをまとめてスケーリングすると、本来未知であるはずのテストデータを既知としてしまうため、良くないのです。したがって、スケーリングする前にトレーニングデータとテストデータを分離します。
※ターゲットベクトルである’Survived’については、トレーニングデータとテストデータ双方とも、スケーリングする必要はありません。モデルの入力値である右側部分の入力関数にのみスケーリングすべき特徴量が存在するため、仮に単なる数値であるターゲットベクトルをスケーリングしたところで、入力関数の形は変わらないからです。
#In
#スケーリング前にトレーニングデータとテストデータを分離
X_train, X_test, y_train, y_test = train_test_split(titanic_features, titanic_target, test_size=0.25)
#スケーラーを作成
scaler = StandardScaler()
#トレーニングデータをスケーリング
scaled_X_train = scaler.fit_transform(X_train)
#トレーニングデータのスケーリングで、テストデータをスケーリング
scaled_X_test = scaler.transform(X_test)
columns = titanic_features.columns
print('1. scaled_X_train\n{}'.format(pd.DataFrame(scaled_X_train, columns=columns).head()),
'\n\n2. scaled_X_test\n{}'.format(pd.DataFrame(scaled_X_test, columns=columns).head()))
#Out
1. scaled_X_train
Pclass Sex Age Fare Title_Int
0 -1.552449 1.318220 1.824374 0.533983 1.170747
1 0.820631 -0.758599 -0.207585 -0.501537 -0.727455
2 -0.365909 -0.758599 1.673858 -0.388862 -0.727455
3 0.820631 -0.758599 0.196405 -0.490228 -0.727455
4 -1.552449 -0.758599 0.196405 0.385202 -0.727455
2. scaled_X_test
Pclass Sex Age Fare Title_Int
0 0.820631 -0.758599 0.196405 -0.504468 -0.727455
1 -1.552449 1.318220 0.462069 0.932072 1.170747
2 -0.365909 -0.758599 2.050147 -0.411983 3.068950
3 -0.365909 -0.758599 -0.809647 -0.449178 -0.727455
4 0.820631 -0.758599 0.243961 -0.501537 -0.727455
モデルと最適化
１．モデル
ロジスティック回帰は、タイタニックデータセットのようなクラス分類に使われる、わかりやすいモデルです。ロジスティック回帰には右図のシグモイド関数が使われています。まずは入力関数から、シグモイド関数の入力値を求めます。
\[
z=wx+b \verb| ・・・①|
\]
入力値は、’w: 重み’、’x: 特徴量行列’、’b: バイアス(過学習防止)’からなる入力関数の計算結果’z’です。初期の重みは任意の値になります。’z’が何かイメージしにくいのですが、各観測値について、’Age’や’Fare’等の特徴量行列を重みに応じて合算したもので、ターゲットベクトルである’Survived’を予測するのに使われます。この入力値をシグモイド関数へ渡すことで、シグモイド関数は0~1の結果(確率)を返します。閾値を指定しなければ、a(z)が0.5より下の場合は’0: Unsurvived’、0.5より上の場合は’1: Survived’と分類されます。
\[
a(z) = \frac{1}{1 + \mathrm{e}^{-z}} \verb| ・・・②|
\]
２．最適化
初期の重みはテキトーであるため、シグモイド関数も誤ってクラス分類してしまいます。誤ったクラス分類をした際に、入力値を計算している入力関数にペナルティを課し、重みを最適化する必要があります。ペナルティの計算には、残差平方和というコスト関数が使われます。
\[
J(w)=\sum^m_{i=1}\frac{1}{2m}((a(z^{(i)}))-y^{(i)})^2 \verb| ・・・③|
\]
残差平方和は、各観測値について、モデルが予測したクラスと正しいクラスの差異を合計したものです。’a(z): モデルが予測したクラス’を、’y: 正しいクラス’を示しています。例えば、ある観測値について、モデルが100%の自信を持って’0: Unsurvived’だと予測したものの、正しいクラスが’1: Survived’の場合(a(z)=0、y=1の場合)、(0-1)^2=1となります。一方、ある観測値について、モデルが100%の自信を持って’1: Survived’だと予測し、正しいクラスも’1: Survived’の場合(a(z)=1、y=1の場合)、(1-1)^2=0となります。③が最小となるモデルが、良いモデルであると言えるわけです。
③の関数に①②の関数を代入して変形すると、いわゆる対数尤度関数になります。③を最小化することと、④を最小化することは同義です。
\[
J(w)=\sum^m_{i=1} [-y^{(i)}\log(a(z^{(i)}))-(1-y^{(i)})\log(1-a(z^{(i)}))] \verb| ・・・④|
\]
対数尤度関数は、全ての観測値についてのコストの総和であるため少しわかりにくく感じます。ある観測値に限定して考えると、イメージをつかみやすくなります。
右図は、ある観測値について、正しいクラスが’1: Survived’の場合、’a(z): モデルが予測したクラス’を変動させた時の’J(w): コスト’の変化を示しています。負の対数関数のグラフからa(z)=0~1の区間を切り取ったような形で、a(z)=0に近づくにつれコストはプラスへ発散し、a(z)=1に近づくにつれコストは0へ収束します。つまり、コストを下げるには、モデルが正しくクラスを予測する必要があると言うわけです。
\[
J(a(z))=-y\log(a(z))-(1-y)\log(1-a(z))
\]
\[
\begin{eqnarray}J (a(z))={\rm} \left\{\begin{array}{ll} -\log(a(z)) ,&(y = 1)\\ -\log(1-a(z)), &(y = 0) \end{array} \right\} \end{eqnarray}
\]
scikit-learnは賢く、こうしたロジックを1行で表現できます。ロジスティック回帰のモデルを作成し、トレーニングデータを学習させました。
#In #モデルを作成 clf = LogisticRegression() #トレーニングデータでモデルを学習 clf.fit(scaled_X_train, y_train)
予測
学習したモデルへテストデータの特徴量行列を渡し、テストデータのターゲットベクトルである’Survived’のクラスを予測しました。見た目はだいたい合ってそうな感じがします。
#In
#学習したモデルでテストデータのターゲットベクトルのクラスを予測
y_hat = clf.predict(scaled_X_test)
print('1. True Class of Survived\n{}'.format(np.array(y_test)),
'\n\n2. Predicted Class of Survived\n{}'.format(y_hat))
#Out
1. True Class of Survived
[0 1 0 1 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 1
1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0
0 1 1 1 0 1 0 1 0 0 1 1 1 0 0 0 0 0 0 1 0 1 1 0 0 0 1 0 0 0 1 1 0 1 0 1 1
1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 1 1 0 0 1 1 1 0 0 0 1 0
0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 1 1 1 1 1 0 0 0 1 0 0 1 1 0 0 1 1 0 1 1
1 1 1 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 1 0
0]
2. Predicted Class of Survived
[0 1 0 0 0 1 0 0 1 1 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1
0 0 1 0 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 0 1 1 0 0 0 1 1 0 0 1 0 0 0 1 1 0 0
1 0 1 1 0 0 1 1 1 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 1 0 1 0 1 0
0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0 1 0
0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 1 1 1 0 1 1 0 0 0 1 0 1 1 1 0 0 1 1 0 1 1
1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 1 1 1 0
0]
モデルの評価
１．分類レポート
分類レポートは、クラス分類の結果を、精度、再現率、F1値(精度と再現率の調和平均)ごとにまとめたレポートです。基本的にF1値を見れば良く、’0: Unsurvived’に比べて’1: Survived’が低いことから、いびつなモデルになっていることがわかります。’1: Survived’は観測値の数が少ないため、クラス分類を誤った場合、入力関数により大きなペナルティを課すようロジスティック回帰のパラメータ’class_weight=’balanced”を指定すると改善する可能性があります。
#In
#分類レポート
print(classification_report(y_test, y_hat))
#Out
precision recall f1-score support
0 0.85 0.85 0.85 144
1 0.72 0.73 0.73 79
micro avg 0.81 0.81 0.81 223
macro avg 0.79 0.79 0.79 223
weighted avg 0.81 0.81 0.81 223
２．混同行列
混同行列も、クラス分類の結果をまとめたレポートです。混同行列をseabornへ投入することで、それぞれのクラスについてどの程度正しく/誤って分類したかを可視化できます。桁あふれしていないにのなぜか指数表示されるため％で作成しました。
例えば、左上のセルでは、正しいクラスが’0: Unsurvived’で、モデルが’0: Unsurvived’だと正しく予測した割合は、全体の約55％であることを示しています。
#In
#混同行列
matrix = confusion_matrix(y_test, y_hat)
class_names = ['Unsurvived', 'Survived']
dataframe = pd.DataFrame(data=matrix/y_test.count()*100, index=class_names, columns=class_names)
sns.heatmap(dataframe, annot=True, cbar=None, cmap='Blues')
plt.title('Confusion Matrix (%)')
plt.xlabel('Predicted Class')
plt.ylabel('True Class')
plt.tight_layout()
３．受信者動作特性曲線
受信者動作特性曲線(ROC曲線)は、モデルの閾値(確率が何％で’1: Survived’と判断するか)を変動させた時のTrue Positive Rate(真陽性率)とFalse Positive Rate(偽陽性率)の関係を示します(※)。True Positive RateとFalse Positive Rateはトレードオフの関係であるため、True Positive Rateを最大化しつつ、False Positive Rateが最小化されるモデルの閾値がベストです。
以下のグラフでは、True Positive Rateが約80％から伸び悩む一方、False Positive Rateは約20％から急増しています。したがって、分類結果が上記の割合になるようなモデルの閾値がベストであることがわかります。これは、分類レポートや混同行列の結果とも一致しています。
※True Positive Rateは、正しいクラスが’1: Survived’で、モデルが’1: Survived’だと正しく予測した割合です。一方、False Positive Rateは、正しいクラスが’0: Unsurvived’で、モデルが’1: Survived’だと誤って予測した割合です。例えば、一番右上の点は、True Positive RateとFalse Positive Rate双方とも1であることから、モデルの閾値を0にした場合(何でもかんでも’1: Survived’と判断した場合)を示しています。以下の動画でROC曲線がわかりやすく説明されていました。
#In
#ROC曲線
y_hat_proba = clf.predict_proba(scaled_X_test)[:, 1]
false_positive_rate, true_positive_rate, threshold = roc_curve(y_test, y_hat_proba)
plt.plot(false_positive_rate, true_positive_rate, color='blue')
plt.plot([0, 1], ls='--', color='orange')
plt.plot([0, 0], [1, 0], c='.7')
plt.plot([1, 1], c='.7')
plt.title('Receiver Operating Characteristic')
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
４．学習曲線
学習曲線は、トレーニングデータと検証データ(※)モデルの学習状況を示します。次のことがわかりました。
- トレーニングデータの数を増やすにつれ、Training ScoreとCross Validation Scoreの差が収束しているため、過学習は発生していない。
- トレーニングデータの数を増やしても、200件を越したあたりからCross Validation Scoreは向上しないため、追加のデータ収集を行ってもモデル改善には繋がらない。
多くの場合、追加のデータ収集は多大なコストと時間がかかります。例えば、上記の場合のように、追加のデータ収集がモデル改善に繋がらないことがわかっていれば、コストと時間を浪費しなくて済みます。
※交差検証は、データセットをトレーニングデータと、検証データに分離し、トレーニングデータで学習したモデルの検証を行う手法です。1度目の分離でトレーニングデータとしたものは、2度目以降の分離では検証データとし、分割数の回数だけ分離と検証が行われます。持ち回りでトレーニングデータを担当し、学習したモデルで検証データのクラス分類を行い、良し悪しを確認するといった感じです。なお、交差検証にスケール前特徴両行列を使っているのは、交差検証とスケーリングを同時に行うことが難しいからです。パイプラインを使うことで、交差検証にスケーリングを取り込めるものの、検証データもスケーリングされてしまうため良くありません。
#In
#学習曲線
train_sizes, train_scores, test_scores = learning_curve(LogisticRegression(), #ロジスティック回帰
titanic_features, #スケール前特徴量行列
titanic_target, #ターゲットベクトル
cv=5, #交差検証の分割数
scoring='accuracy', #性能評価指標
n_jobs=-1) #全てのコアを利用
train_mean = np.mean(train_scores, axis=1)
train_std = np.std(train_scores, axis=1)
test_mean = np.mean(test_scores, axis=1)
test_std = np.std(test_scores, axis=1)
plt.plot(train_sizes, train_mean, '--', color='blue', label='Training Score')
plt.plot(train_sizes, test_mean, color='orange', label='Cross Validation Score')
plt.fill_between(train_sizes, train_mean - train_std,
train_mean + train_std, color='#DDDDDD')
plt.fill_between(train_sizes, test_mean - test_std,
test_mean + test_std, color='#DDDDDD')
plt.title('Learning Curve')
plt.xlabel('Training Set Size')
plt.ylabel('Accuracy Score')
plt.legend(loc='best')
