タイタニック号生存者予測- 多変量解析

　本稿では、前回の内容に続き、タイタニック号の生存者予測の議論を行います。前回は、生存状況について一つ一つの特徴量ごとに調べましたが、本データのように特徴量を複数持つ場合は、特徴量間で相関がある可能性が考えられます（疑似相関）。その場合、実際には生存率と直接関係のない特徴量を説明変数としてしまうことで予測に悪い影響を及ぼします（多重共線性：マルチコ現象）。本稿では、特徴量間の相関を考慮し、疑似相関を排除し、適切な特徴量を選択することを目的とします。  

　特徴量間の相関係数

　前回は、特徴量それぞれについて生存率が増加 or 低下したかにより生存率への影響を観察しましたが、pandasでは特徴量間の相関係数をヒートマップで簡単に示すことが可能です。
　まずは、特徴量を評価するために、数字データに変換します。特徴量のラベル付けは、前回の結果を基にしています。

train['Fare'] = pd.cut(train['Fare'], [-1, 10, 60, 100, 1000], labels=False) # 運賃が10＄以下では生存率が低い

train['Sex'] = train['Sex'].replace('male', 0)
train['Sex'] = train['Sex'].replace('female', 1)

train['Embarked'] = train['Embarked'].map({'S':0, 'Q':1, 'C':2})

# ミドルネーム（敬称）は生存率と関係がありそうだ
train['honorific'] = train['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
train = train.drop(columns = 'Name')
train['honorific'] = train['honorific'].str.replace(' ', '')
train['honorific'] = train['honorific'].map({'col':1, 'Major': 1, 'Sir': 1, 'Master':1, 'Dr':1})
train['honorific'] = train['honorific'].fillna(0)  

train['Familysize'] = train['Parch'] + train['SibSp'] + 1 
train['Familysize'] = train['Familysize'].map({1:1, 2:2, 3:2, 4:2, 5:0 , 6:0, 7:0, 8:0, 11:0})
train['Familysize'] = train['Familysize'].fillna(0)  

# 年齢が10歳以下では生存率が高い、60歳以上では生存率が低い
Age_train = []
for idx, row in train.iterrows():
    if row['Age'] < 10:
        Age_train.append(2)
    elif 10 < row['Age'] < 60:
        Age_train.append(1)
    elif row['Age'] > 60:
        Age_train.append(0)
    else:
        Age_train.append(np.nan)
train['Age'] = Age_train

# チケット番号は数字列の先頭の文字が生存率と影響していそうだ
Ticket_train = []
for idx, row in train.iterrows():
    Ticket_train.append(row['Ticket'].split(' ')[-1])
train['Ticket_num'] = Ticket_train
train['Ticket_num'] = train['Ticket_num'].replace('LINE', 0)
train['Ticket_num1'] = train['Ticket_num'].astype(int)
Ticket_train1 = []
for idx, row in train.iterrows():
    if row['Ticket_num1'] < 3000 or 20000 < row['Ticket_num1'] <300000:
        Ticket_train1.append(1)
    elif 3000 < row['Ticket_num1'] < 10000 or 300000 < row['Ticket_num1'] < 4000000:
        Ticket_train1.append(0)
    elif 10000 < row['Ticket_num1'] < 20000:
        Ticket_train1.append(2)
train['Ticket_num'] = Ticket_train1

sns.set_style(style = 'whitegrid')
sns.set_context('talk')
plt.figure(figsize=(15, 12))
sns.heatmap(train.corr(), annot=True, cmap='RdBu')
plt.show()

　 　Survivedと各特徴量の相関係数を確認すると、性別が最も大きいことがわかります。続いて、乗客クラス、運賃、チケット番号、家族サイズがほぼ同程度の相関であり、相関があると言えそうです。続いて、出港地、年齢、敬称が続きますが、これらの生存率との相関係数は、0.2以下であり、相関は弱いと考えられます。

疑似相関の確認

本題の特徴量の選定を行います。そのために、ヒートマップにて説明変数間で相関係数が大きくなっている特徴量を確認していきます。

PclassとTicket_numberの関係

　ヒートマップから、乗客クラスと運賃および乗客クラスとチケット番号の相関が大きくなっていました。乗客クラスと運賃の相関は直感で想像できそうですが、乗客クラスとチケット番号のように直感では発見できない相関を簡単に見つけ出すことがヒートマップの強みです。
　実際にグラフで2つの特徴量の相関を確認してみると下の結果となります。

sns.catplot(data=train , x='Pclass', hue = 'Ticket_num', kind = 'count', palette = 'gist_heat')
plt.title('Ticket Number by Pclass', fontsize=20)

グラフからTicket_num = 0の乗客の98%がサードクラスであるという顕著な傾向を示しました。乗客クラスという特徴量が曖昧でしたが（ネットでは、ファーストクラスを上流家庭とする例あり）、Pclassはチケットクラスを表していることが示唆されます。つまり、チケット番号を整理することで、Pclassを再現することができるということです。
　Ticket_num = 0かつPclass = 1, 2の乗客を確認してみます。

train[(train['Pclass'] == 1) & (train['Ticket_num'] == 0)]

train[(train['PassengerId'] == 93) | (train['PassengerId'] == 541) | (train['PassengerId'] == 663) | (train['PassengerId'] == 746)]

Pclass = 1の場合、5000番台で新たに区分わけをすることで解決できそうです。

train[(train['Pclass'] == 2) & (train['Ticket_num'] == 0)]

train[(train['PassengerId'] == 309) | (train['PassengerId'] == 433) | (train['PassengerId'] == 875)]

Pclass = 2の場合、3000番台のチケットという条件のみでは、乗客クラスの分類分けはできませんが、先頭のアルファベット列を考慮することで乗客クラスの区分けは可能そうです。
　この結果から、乗客クラスはチケットと強い相関があることから、チケットクラスを表す変数だと判断します。よって、チケット番号と生存率の相関は疑似相関であり、チケットクラス（Pclass）を説明変数として使用することとします。

PclassとFareの関係

これら2つの相関をグラフで確認した結果が下図となります。

sns.catplot(data=train , x='Pclass', hue = 'Fare', kind = 'count', palette = 'ocean_r')
plt.title('FareBand by Pclass', fontsize=20)

　このグラフから、Pclassが小さいほどFareは大きくなる傾向が確認できました。一方で、Pclass = 1かつ FareBand = 1となる乗客が一定数おり、これら2つの特徴量は完全に一致するわけではありませんでした。そこで、２つの特徴量を含んだ新たな特徴量「Fclass = FareBand / Pclass」を導入しました。

train['Fclass'] = (train['Fare'] / train['Pclass']).round(2)
train['Fclass'] = train['Fclass'].map({0:0, 0.33:0, 0.67:1, 0.5:1, 1.0:1, 2.0:2, 3.0:2})
sns.catplot(data=train , x='Fclass', hue = 'Survived', kind = 'count', palette = 'gist_heat')
plt.title('Survived  by Fclass', fontsize=20)

　このような特徴量の変換を行うことで、2つの特徴量の影響を上手く合成することができました。

欠損値補完

　欠損値補完は元データの傾向を変化させてしまう可能性があるため慎重に行う必要があります。簡単な方法は特徴量の平均値を代表値として補完することですが、元の欠損値が多いほどデータの特徴が失われます。本データでは複数特徴量があるため、他の特徴量から欠損値を推測し元のデータ傾向から外れないように補完します。

Ageの欠損値補完

　今回のデータセットにおいて、年齢を推測するヒントとして同乗した家族構成（SibSp、Parch、Familysize ）があります。例えば、家族規模が1である場合、単身での乗船であることから、子どもではないことが推測されます。また、SibSｐ ≠ 0かつ、Parch = 0であれば、配偶者か兄弟が同乗しているため、その相手の年齢に近しいと考えられます。
　しかし、家族単位で年齢情報が欠損していたり、親子関係（どちらが親か）は特定できませんでした。そこで、敬称から年齢の予測を行うことを考えました。実際に未婚女性であればMiss、既婚女性はMrsが使用されるため、敬称は年齢層と関係します。これはヒートマップにおいて、年齢と最も相関が強い変数が敬称であることからも判断できます。 　これを確認するため、年齢が欠損値である乗客の敬称ごとの統計量を計算しました。

honorific_list = train_age['honorific'].values.tolist()
honorific_unique = list(set(honorific_list))
for honorific in honorific_unique:
    age_avg = train[train['honorific'] == honorific]['Age'].mean()
    age_std = train[train['honorific'] == honorific]['Age'].std()
    print('{}：平均値：{} 標準誤差：{}'.format(honorific, age_avg, age_std))

[out]
Mrs：平均値：35.898148148148145 標準誤差：11.433627902196413
Dr：平均値：42.0 標準誤差：12.016655108639842
Miss：平均値：21.773972602739725 標準誤差：12.99029242215268
Mr：平均値：32.368090452261306 標準誤差：12.70879272257399
Master：平均値：4.574166666666667 標準誤差：3.6198716433439615

　この結果から、年齢敬称ごとに年齢差が明確に現れました。特に、Masterは男の子どもを表すようです。従って、敬称ごとの年齢の平均値から標準誤差の範囲で乱数を生成し欠損値に補完しました。

honorific_list = train[train['Age'].isnull()]['honorific'].values.tolist()
age_list = []

for honorific in honorific_list:
    age_avg = train[train['honorific'] == honorific]['Age'].mean()
    age_std = train[train['honorific'] == honorific]['Age'].std()
    age_list.append(random.uniform(age_avg - age_std, age_avg + age_std))
train['Age'][np.isnan(train['Age'])] = age_list
train['Age'] = train['Age'].astype(int)

Embarkedの欠損値補完

　出港地の欠損値は2つと少ないです。そのため、乗客情報を確認し代入値を決定します。
　二人の情報は下記の通りでした。

[train[train['Embarked'].isnull()]

この情報から二人はチケットが同じであることがわかります。但し、家族ではないことがわかります。出港地を推測できそうな特徴量は、チケット番号だと考え、113572近辺のチケットを持つ乗客を調べました。

train[(train['Ticket_num1'] > 113500) & (train['Ticket_num1'] < 113600)]

この結果から、SとCのどちらかだと考えられますが特定には至りません。Amelieという名前からフランス人だと考えCを補完しました。

結果

　

生存者予測

　説明変数を定めたので、実際に予測を行います。予測モデルはランダムフォレストを使用します。まず、sklearnから必要なライブラリをインポートします。

rom sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import precision_score
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix, accuracy_score, precision_score, recall_score, f1_score
from sklearn.model_selection import GridSearchCV

　train_xを説明変数のみとし、train_yを正解ラベルとします。モデルの精度を検証するため、trainデータを訓練用と評価用に分割します。

train_x = train.drop('Survived', axis = 1)
train_y = train['Survived'] 
(train_x, test_x ,train_y, test_y) = train_test_split(train_x, train_y, test_size = 0.2, random_state = 42)

　ランダムフォレストのパラメータを設定しグリッドサーチにて最適なパラメーターを決定します。

parameters = {
    'n_estimators' :[3,5,20],#作成する決定木の数
    'random_state' :[7,42],
    'max_depth' :[3,5,8,10],#決定木の深さ
    'min_samples_leaf': [2,5,10,20],#分岐し終わったノードの最小サンプル数
    'min_samples_split': [2,5,10,20]#決定木が分岐する際に必要なサンプル数
}
#グリッドサーチを使う
clf = GridSearchCV(estimator=RandomForestClassifier(), param_grid=parameters, cv=2, iid=False)
#学習モデルを作成
clf.fit(train_x, train_y)

　訓練データと評価データの精度を確認します。

best_clf = clf.best_estimator_ #ベストパラメーター
print('score: {:.2%}'.format(best_clf.score(train_x, train_y)))
print('score: {:.2%}'.format(best_clf.score(test_x, test_y)))

[out]
score: 83.85%
score: 81.56%

　約8割の精度であり、それなりの精度でした。
　本題のtestデータに対し予測を行います。

pre_x = test.drop('PassengerId', axis = 1) # 説明変数のみにする
pre_y = best_clf.predict(pre_x)
PassengerId = np.array(test["PassengerId"]).astype(int)
my_solution = pd.DataFrame(pre_y, PassengerId, columns = ["Survived"])
my_solution.to_csv("RF_bestclf.csv", index_label = ["PassengerId"])

　初めて出力したExcelをワクワクしながらkaggleに提出したところ結果は0.75358でした。だいぶ下の順位でした。

予測結果考察

　予測結果が芳しくなかった原因を考察します。まず、ヒートマップから、特徴量の相関を確認します。

　ヒートマップからtestデータでは、生存率は性別と強い相関があることがわかります。また、trainデータでFclassという新しい特徴量を作成したものの、testデータでは相関係数が小さくなっていました。つまり、trainデータにおいて、設定した特徴量はtestデータに対し必ずフィットする訳ではないということです。
　ここで、相関が最も強い性別に着目すると下の結果となりました。

　男性は約7%しか生存できておらず、女性は約12%のみが死亡しています。つまり、性別のみから生存者の予測はあらかたできてしまうということです。裏を返すと、性別以外の要因から生存率を判定できる要因を見つけ出さなければ精度改善には至らないということです。そこで、性別ごとに特徴量を再検討していきます。

sns.catplot(data=train, x='Fclass', hue = 'Survived', kind = 'count', col='Sex', palette = 'gist_heat')

　Fclassは男性の場合、生存率が段々と増加する結果となりましたが、女性の場合にはFclass = 3のの乗客に死亡者が集中していました。そのため、女性のFclass区分は再振り分け要と判断します。また、Fclassの元の変数であったPclassとFareの変数も確認します。

　PclassはほとんどFclassと同様の傾向が見られました。Fareでは、男女に分けたことで、Fare = 2,3に大きな差が見られなくなりました。また、Fareはチケット番号が同一である乗客で同料金でした。つまり、同乗している家族の合計金額を表しています。そこで、一人当たりの運賃を表すperFare（Fare/Familysize）という特徴量を導入します。

train['perFare'] = train['Fare'] / train['Familysize']

train['FareBand'] = pd.cut(train['Fare'], [0, 7, 7.5, 7.8, 8, 10, 15, 25, 35, 60, 100, 1000])
sns.catplot(data=train , x='FareBand', kind = 'count', hue = 'Survived', col = 'Sex', palette = 'gist_heat').set_xticklabels(rotation=90)

train['perFareBand'] = pd.cut(train['perFare'], [0, 7, 7.5, 7.8, 8, 10, 15, 25, 35, 60, 100, 1000])
sns.catplot(data=train , x='perFareBand', kind = 'count', hue = 'Survived', col = 'Sex', palette = 'gist_heat').set_xticklabels(rotation=90)

　一人当たりの運賃に変換することで、男性の場合、perFareが25＄以上の乗客の生存率は高いこと、女性では、10＄以上で生存率が増加する傾向を発見することができました。

　続いて、男性と女性で年齢ごとの生存状況は異なる結果となりました。女性の場合、10歳以下の生存率は他年代と比較し高くありませんでした。 　
このグラフからも女性は年齢と生存率の相関が小さいことがわかりました。男性の場合では、10歳以下の生存率が50％を超えるため、その影響を考慮することとします。
　家族サイズについては、男女ともに相関があることがわかります。 　出港地と生存率は相関がないと判断します。

　再予測

　男女ごとに特徴量を設定してきました。そのため、性別ごとに生存者を予測し最終的に個々の結果を統合します。
　性別ごとにデータをわけ、特徴量を変換します。

train_male = train[train['Sex'] == 'male']
train_female = train[train['Sex'] == 'female']
test_male = test[test['Sex'] == 'male']
test_female = test[test['Sex'] == 'female']

dataset = [train_male, train_female, test_male, test_female]
for data in dataset:
    test_Fare = data[(data['Pclass'] == 3) & (data['Age'] > 40)]['Fare'].mean()
    data['Fare'] = data['Fare'].fillna(test_Fare) 
    
    data['Familysize'] = data['Parch'] + data['SibSp'] + 1    
    data['perFare'] = data['Fare'] / data['Familysize']
    
    data['Embarked'] = data['Embarked'].map({'S':0, 'Q':1, 'C':2})
    data['Embarked'] = data['Embarked'].fillna(0) 
    
    data['honorific'] = data['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
    data['honorific'] = data['honorific'].str.replace(' ', '')
    honorific_list = data[data['Age'].isnull()]['honorific'].values.tolist()
    age_list = []
    for honorific in honorific_list:
        age_avg = data[data['honorific'] == honorific]['Age'].mean()
        age_std = data[data['honorific'] == honorific]['Age'].std()
        age_list.append(random.uniform(age_avg - age_std, age_avg + age_std))
    data['Age'][np.isnan(data['Age'])] = age_list
    

dataset_male = [train_male, test_male]
for data_male in dataset_male:
    data_male['Pclass'] = data_male['Pclass'].map({1:1, 2:0, 3:0})
    data_male['honorific'] = data_male['honorific'].map({'col':1, 'Major': 1, 'Sir': 1, 'Master':1, 'Dr':1})
    data_male['honorific'] = data_male['honorific'].fillna(0)  
    data_male['FareBand'] = pd.cut(data_male['perFare'], [-1, 25, 1000], labels = False)
    data_male['AgeBand'] = pd.cut(data_male['Age'], [-5, 10, 100], labels = False)
    
dataset_female = [train_female, test_female]
for data_female in dataset_female:
    data_female['Pclass'] = data_female['Pclass'].map({1:1, 2:1, 3:0})
    data_female['FareBand'] = pd.cut(data_female['perFare'], [-1, 7.3, 10, 1000], labels = False)

dataset_male = [train_male, test_male]
for data_male in dataset_male:
    data_male = data_male.drop(columns = {'Name', 'Sex', 'perFare', 'Age', 'AgeBand', 'Fare', 'SibSp', 'Parch', 'Familysize', 'Ticket', 'Cabin'}, inplace=True)    

dataset_female = [train_female, test_female]
for data_female in dataset_female:
    data_female = data_female.drop(columns = {'Name', 'Sex', 'perFare', 'Age', 'honorific', 'Fare', 'SibSp', 'Parch', 'Familysize', 'Ticket', 'Cabin'}, inplace=True)    

　男女ごとにヒートマップを確認してみます。

plt.figure(figsize=(12, 9))
sns.heatmap(train_male.corr(), annot=True, cmap='RdBu', vmin=-1,vmax=1)
plt.show()

plt.figure(figsize=(12, 9))
sns.heatmap(train_female.corr(), annot=True, cmap='RdBu', vmin=-1,vmax=1)
plt.show()

　これら説明変数に対し、ランダムフォレストを適用します。

train_male_x = train_male.drop(columns = {'Survived', 'PassengerId'}) # 説明変数のみにする
train_male_y = train_male['Survived'] # 正解クラス
train_female_x = train_female.drop(columns = {'Survived', 'PassengerId'}) # 説明変数のみにする
train_female_y = train_female['Survived'] # 正解クラス

# ランダムフォレストのパラメータの候補をいくつか決める
parameters = {
    'n_estimators' :[3,5,20],#作成する決定木の数
    'random_state' :[7,42],
    'max_depth' :[3,5,8,10],#決定木の深さ
    'min_samples_leaf': [2,5,10,20],#分岐し終わったノードの最小サンプル数
    'min_samples_split': [2,5,10,20]#決定木が分岐する際に必要なサンプル数
}

#グリッドサーチを使う
clf_male = GridSearchCV(estimator=RandomForestClassifier(), param_grid=parameters, cv=2, iid=False)
clf_female = GridSearchCV(estimator=RandomForestClassifier(), param_grid=parameters, cv=2, iid=False)

#学習モデルを作成
clf_male.fit(train_male_x, train_male_y)
clf_female.fit(train_female_x, train_female_y)

best_clf_male = clf_male.best_estimator_ 
test_male_x = test_male.drop(columns = {'PassengerId'}) # 説明変数のみにする

test_male_y = best_clf_male.predict(test_male_x)
test_male['Survived'] = test_male_y

best_clf_female = clf_female.best_estimator_ 
test_female_x = test_female.drop(columns = {'PassengerId'}) # 説明変数のみにする
test_female_y = best_clf_female.predict(test_female_x)
test_female['Survived'] = test_female_y

PassengerId_male = np.array(test_male["PassengerId"]).astype(int)
my_solution_male = pd.DataFrame(test_male_y, PassengerId_male, columns = ["Survived"])
PassengerId_female = np.array(test_female["PassengerId"]).astype(int)
my_solution_female = pd.DataFrame(test_female_y, PassengerId_female, columns = ["Survived"])
my_solution = pd.concat([my_solution_male, my_solution_female])
my_solution.to_csv("my_RF_bestclfSex.csv", index_label = ["PassengerId"])

結果は0.78468でした。予測精度は性別で分割することで約３％向上しており有効であることが確認できました。