オーバーサンプリングの問題点 - smote VS GICaPS

今回は、不均衡データでのクラス分類の方法および問題点について投稿します。

不均衡データと主要な分析手法

不均衡データとは、クラス分類したい目的変数に大きな偏りがあるデータを指します。例としてよく挙げられるのが、ガンの予測です。100人がガン検診を受診した場合、実際にガンである患者が5人のみだったとします。この結果を予測する場合、全患者が陰性であると予測すると正解率は95%と高確率である一方、陽性患者を誰一人当てることができていません。マイノリティ側(この例では、陽性)の予測精度に関心がある場合、この結果は好ましくありません。また、機械学習モデルを用いる上でも、マジョリティ側の予測精度を向上させたほうが、全体の精度が向上するため、マイノリティ側の予測精度がおざなりになるという問題があります。

そこで、不均衡データを均衡データに変換してしまおうという考えが提案されています。主に、アンダーサンプリングとオーバーサンプリングという手法です。
アンダーサンプリングとは、マイノリティ側のデータ数と一致するようにマジョリティ側のデータを削減する方法です。上記のガンの陽性・陰性の予測例では、マイノリティ側(陽性)は5人ですので、マジョリティ側(陰性)も5人に合わせ全体を10人としデータを再構築します。アンダーサンプリングの問題は、マイノリティ側のデータ数が非常に少ない場合、全体のデータ数が激減してしまうことです。
反対に、オーバーサンプリングはマジョリティ側のデータ数にマイノリティ側のデータ数を合わせる方法です。つまり、マジョリティ側(陰性)は95人ですので、マイノリティ側(陽性)を95人になるようにし、全体で190のデータ数を構築します。アンダーサンプリングでは、データを削減するという方法だったため、全データが実際のデータでしたが、オーバーサンプリングの場合、データを増やすため擬似データを使用することになります。不均衡の傾向が顕著な場合、擬似データの割合が大きくなってしまいます。 どちらの手法にも課題はありますが、疑似データをより正確に生成することが可能となれば、オーバーサンプリングは今後不均衡データを扱う際に主流となる可能性が高いのではないかと考えています。そのため、本記事では、現状のオーバーサンプリングの手法の考察と最新の研究内容について紹介したいと思います。

smote : Synthetic Minority Over-sampling Technique

smoteとは、現在オーバーサンプリングを行う上で主要な手法です。というのも、pythonのモジュールにあるため、簡単に実装することができるためです。私も最近になって、初めて不均衡データというものを知り実装してみましたが非常に簡単でした。一方、モデルの予測精度は向上しませんでした。ビッグデータを扱う上では、一つ一つのデータというのをじっくり見る機会というのは少ないと思います。そのため、簡単に実装できるモジュールを安易に利用するのではなく、その原理と問題点を理解することが求められます。

smoteの原理

smoteの原理は公式のページを参照して頂くことが最も間違いありません。今回は、公式ページを参照して簡単に紹介します。

f:id:datascience_findings:20210419212842p:plain
smoteの原理

newサンプルを生成する仕組みは上図に示すように非常に単純です。2つのマイノリティデータを選択し、2点間の線上に乱数λ(0<λ<1)を振るだけです。このアルゴリズムをマイノリティデータがマジョリティデータと一致するまで繰り返すのです。

imbalanced-learn.org

smoteの実装

実際に簡単な例を作成し実装します。必要なモジュールをimportします。

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
sns.set_context('talk')
import random
from imblearn.over_sampling import SMOTE # smoteのモジュールとなります

不均衡データは、クラス1が300個、クラス0が10個となるように作成しました。また、説明変数は2つ用意しました。

# クラス0のデータを作成
class0_num = 10
class0_a0 = [random.randint(200, 500) for i in range(class0_num)]
class0_a1 = [random.randint(200, 500) for i in range(class0_num)]
class0 = [0 for i in range(class0_num)]
# クラス1のデータを作成
class1_num = 300
class1_a0 = [random.randint(0, 500) for i in range(class1_num)]
class1_a1 = [random.randint(0, 500) for i in range(class1_num)]
class1 = [1 for i in range(class1_num)]
# クラス0と1を結合
class1_a0.extend(class0_a0)
class1_a1.extend(class0_a1)
class1.extend(class0)

作成したデータは下図のようになります。

lt.scatter(class1_a0, class1_a1, color='olive')
plt.scatter(class0_a0, class0_a1, color='tomato')
plt.show()

f:id:datascience_findings:20210419214628p:plain
赤のデータがクラス0のマイノリティを表しています。今回のデータでは、10:300なので、クラス0のデータを290個増やします。

smote = SMOTE()
origin = pd.DataFrame(data = [class1_a0, class1_a1]).T
origin_smote, class_smote = smote.fit_resample(origin, class1) # smoteの引数は説明変数、目的変数となります

上記コードを実行するだけで、origin_smote, class_smoteとして均衡データが作成されます。実際に散布図で確認します。

index_0 = [i for i, x in enumerate(class_smote) if x == 0]
index_1 = [i for i, x in enumerate(class_smote) if x == 0]

plt.scatter(class1_a0, class1_a1, color='olive')
plt.scatter(origin_smote.loc[index_0][0], origin_smote.loc[index_0][1], color= 'navy')
plt.scatter(origin.loc[class1_num:][0], origin[class1_num:][1], color='crimson')
plt.title('all data')
plt.show()

plt.scatter(origin_smote.loc[index_0][0], origin_smote.loc[index_0][1], color= 'navy')
plt.scatter(origin.loc[class1_num:][0], origin[class1_num:][1], color='crimson')
plt.title('class0 data')
plt.show()
f:id:datascience_findings:20210419220315p:plainf:id:datascience_findings:20210419220321p:plain

青で示されるプロットがsmoteにより生成された擬似データとなります。all dataで生成データの傾向を観察すると、クラス0のデータの範囲内で分散しているように見えます。しかし、クラス0のみで観察した場合、少しいびつなデータとなっていると感じるのではないでしょうか。というのもドーナッツ型のようなデータ補完がされているためです。2点の赤いプロットを結ぶ線上にしかデータは生成されないため、originalデータの位置よって、いびつなデータとなる可能性があります。
同じ条件で乱数を振り、複数回データを生成した結果を載せます。

f:id:datascience_findings:20210419221803p:plainf:id:datascience_findings:20210419221809p:plain
f:id:datascience_findings:20210419222013p:plainf:id:datascience_findings:20210419221956p:plain
f:id:datascience_findings:20210419222302p:plainf:id:datascience_findings:20210419222317p:plain

何度か繰り返しましたがどの場合においても、いびつな傾向が感じられました。これは、優先的に決定される2点の赤プロットが近しい2点であるためです。そのため、対角線上にデータを生成するのではなく、外枠にデータが生成される傾向となります。

smoteの問題点

データの外枠から補完されることの問題点は、説明変数間の相関が変化(低下)してしまうことにあります。実際にdata3でのオリジナルデータとsmoteデータの相関係数を計算すると大きく低下していることとがわかります。

print('original dataの相関係数')
print(origin.loc[class1_num:].corr())
print()
print('smote dataの相関係数')
print(origin_smote.loc[class1_num:].corr())
[out]
origin dataの相関係数
         0        1
0  1.00000 -0.25939
1 -0.25939  1.00000

smote dataの相関係数
          0         1
0  1.000000 -0.090797
1 -0.090797  1.000000

このように、一見上手くデータを補完することができているように見えたとしてもマイノリティ側のデータ傾向は大きく変わっている可能性があります。

説明変数の次元数と相関係数の関係

上記までの例では、smoteを簡単に理解するため説明変数は2次元としていました。しかし、実際に扱うデータは高次元であるため、その影響を確認します。
data3に説明変数と一つ追加し3次元とします。

# a2という説明変数を追加
class0_a2 = [random.randint(200, 500) for i in range(class0_num)]
class1_a2 = [random.randint(0, 500) for i in range(class1_num)]
class1_a2.extend(class0_a2)
origin_3 = pd.DataFrame(data = [class1_a0, class1_a1, class1_a2]).T
# smoteを実行
origin3_smote, class_smote = smote.fit_resample(origin_3, class1) 
# 相関行列を出力
print('original dataの相関係数')
print(origin_3.loc[class1_num:].corr().round(3))
print()
print('smote dataの相関係数')
print(origin3_smote.loc[class1_num:].corr().round(3))
[out]
original dataの相関係数
       0      1      2
0  1.000 -0.259 -0.080
1 -0.259  1.000 -0.008
2 -0.080 -0.008  1.000

smote dataの相関係数
       0      1      2
0  1.000 -0.413 -0.049
1 -0.413  1.000  0.145
2 -0.049  0.145  1.000

2次元データでは説明変数a0, a1の相関係数はsmoteにより低下していましたが、3次元データの場合大きく増加しました。
実際に両データを比較します。

f:id:datascience_findings:20210420182353p:plainf:id:datascience_findings:20210420182400p:plain

左が説明変数が2次元である場合、右が説明変数が3次元である場合のsmoteの実行結果です。2つの図からわかるように、たしかに3次元の場合がより外枠だけでなく、中心も満遍なく補完されています。これは3次元となったことで、smoteにより結ばれる2点の選択肢が増えたためです。
その傾向を確かめるため、説明変数を10次元とし検証します。

# +7次元分のデータを作成
origin_p = origin_3
for p in range(7):
    class0_ap = [random.randint(200, 500) for i in range(class0_num)]
    class1_ap = [random.randint(0, 500) for i in range(class1_num)]
    class1_ap.extend(class0_ap)
    class1_ap_df = pd.DataFrame({p+3:class1_ap})
    origin_p = pd.concat([origin_p, class1_ap_df], axis = 1)
originp_smote, class_smote = smote.fit_resample(origin_p, class1)
print('original dataの相関係数')
print(origin_p.loc[class1_num:].corr().round(3))
print()
print('smote dataの相関係数')
print(originp_smote.loc[class1_num:].corr().round(3))
[out]
original dataの相関係数
       0      1      2      3      4      5      6      7      8      9
0  1.000 -0.259 -0.080  0.265 -0.071  0.233 -0.123 -0.529  0.030  0.187
1 -0.259  1.000 -0.008  0.387  0.379  0.300 -0.210  0.095 -0.318 -0.273
2 -0.080 -0.008  1.000 -0.382 -0.370  0.558 -0.271 -0.508 -0.501  0.493
3  0.265  0.387 -0.382  1.000  0.302  0.247  0.175 -0.032  0.209 -0.108
4 -0.071  0.379 -0.370  0.302  1.000 -0.424  0.715  0.036 -0.011 -0.677
5  0.233  0.300  0.558  0.247 -0.424  1.000 -0.491 -0.600 -0.599  0.385
6 -0.123 -0.210 -0.271  0.175  0.715 -0.491  1.000  0.066  0.039 -0.340
7 -0.529  0.095 -0.508 -0.032  0.036 -0.600  0.066  1.000  0.391 -0.183
8  0.030 -0.318 -0.501  0.209 -0.011 -0.599  0.039  0.391  1.000 -0.143
9  0.187 -0.273  0.493 -0.108 -0.677  0.385 -0.340 -0.183 -0.143  1.000

smote dataの相関係数
       0      1      2      3      4      5      6      7      8      9
0  1.000 -0.375  0.026  0.136 -0.160  0.179 -0.135 -0.569  0.105  0.233
1 -0.375  1.000 -0.013  0.406  0.329  0.333 -0.196  0.123 -0.368 -0.263
2  0.026 -0.013  1.000 -0.415 -0.510  0.633 -0.428 -0.561 -0.517  0.619
3  0.136  0.406 -0.415  1.000  0.355  0.156  0.229  0.059  0.218 -0.212
4 -0.160  0.329 -0.510  0.355  1.000 -0.502  0.775  0.182  0.033 -0.754
5  0.179  0.333  0.633  0.156 -0.502  1.000 -0.585 -0.620 -0.619  0.444
6 -0.135 -0.196 -0.428  0.229  0.775 -0.585  1.000  0.193  0.101 -0.475
7 -0.569  0.123 -0.561  0.059  0.182 -0.620  0.193  1.000  0.409 -0.263
8  0.105 -0.368 -0.517  0.218  0.033 -0.619  0.101  0.409  1.000 -0.149
9  0.233 -0.263  0.619 -0.212 -0.754  0.444 -0.475 -0.263 -0.149  1.000

相関行列から、10次元でのa0とa1のsmote dataの相関係数は3次元だった場合と比較し低下していますが、2次元smote dataと比較すると相関係数は増加しています。この結果から、説明変数の次元を増やした場合、相関係数の振れ幅がoriginal dataの相関係数に収束していくと考えられます。これは、下図に示すように次元が増えることでマイノリティクラスの領域内が次第に埋められていくためだと考えます。

f:id:datascience_findings:20210420212451p:plainf:id:datascience_findings:20210420212500p:plain

オーバーサンプリングの研究動向

上記までのようにsmoteはシンプルなアルゴリズムである反面、補完データに偏りやマイノリティクラス内での不均衡が生じ得ます。そこで、「A Method for Handling Multi-class Imbalanced Data byGeometry based Information Sampling and Class Prioritized Synthetic Data Generation (GICaPS)」の論文で、データ生成がより自然(均等)となるようなアルゴリズムについて述べられていますので紹介します。

arxiv.org

GICaPS (Geometry based Information Sampling and Class Prioritized Synthetic Data Generation)

 GICaPSとは、データの幾何学情報(分布)に基づきサンプリングを行う手法です。アンダーサンプリングとオーバーサンプリングを組み合わせた手法が論文内で提案されています。しかし、今回はオーバーサンプリングに着目しそのアルゴリズムを理解します。本論文でsmoteの問題点として、クラス間の干渉の影響が考慮されていないということを挙げています。これは、smoteのアルゴリズムからもわかるようにマイノリティ側のデータ配置のみしか考慮していません。つまり、マジョリティとマイノリティの境界と思われる周辺に擬似データが生成されると機械学習モデルの精度を低下させる危険があります。  簡単な例を下に示します。

f:id:datascience_findings:20210420223524p:plainf:id:datascience_findings:20210420224051p:plain

 この例では、クラス1は0〜300の範囲、クラス0では200〜500としデータを生成しました。そのため、originalデータから確認できるようにクラス0と1ではデータ傾向は大きく異なります。しかし、smoteにより擬似データを生成するとその境界は考慮されていないため、クラスが重なる領域にoriginalでは1個のみでしたが、35個のデータが生成されてしまいました。これを回避する手法として、no man’s landという概念を導入しています。

no man’s landを用いたオーバーサンプリング

no man’s landとは、翻訳すると誰のものでもない土地を意味します。今回の場合では、class0, 1どちらにも属さない領域となります。no man’s landの領域の決定方法は論文内で下図のように説明されています。 f:id:datascience_findings:20210421210335p:plain
図中の(a)が最もシンプルな例です。k近傍法で切り取られた空間Vで、マイノリティクラスの2点(Xm, Xv)を結んだ際に、マジョリティクラスQが存在する場合、長方形(高さ:線分に最も近いQの線分との垂直距離、幅:線分間にQが存在する幅)がNo man's landを表します。
論文で紹介されているNo man's landを求めるより一般化した方法を紹介します。下図では、マジョリティクラスの任意の点を線分ab上に投射した際の長さを求める方法が述べられています。 f:id:datascience_findings:20210421230455p:plain 実装し、線分abとno man's landのlを求めます。

# a, b点を決定し、ベクトルabを求める
a = np.array([[class0_a0[1]], [class0_a1[1]]])
b = np.array([[class0_a0[0]], [class0_a1[0]]])
ab = b - a
# t1, t2を決定し、at1, at2ベクトルを求める
t1 = np.array([[class1_a0[157]], [class1_a1[157]]])
t2 = np.array([[class1_a0[46]], [class1_a1[46]]])
at1 = t1 - a
at2 = t2 - a
# 論文の式から、pt1, pt2を求める
pt1 = (np.dot(ab.T, at1) * ab) / np.dot(ab.T, ab)
pt2 = (np.dot(ab.T, at2) * ab) / np.dot(ab.T, ab)
apt1 = pt1 + a
apt2 = pt2 + a

# 散布図を作成
plt.scatter(class1_a0, class1_a1, color='olive', label = 'class1')
plt.scatter(origin.loc[class1_num:][0], origin[class1_num:][1], color='crimson', label = 'class0')
plt.title('smote data')
plt.legend(loc="upper left", bbox_to_anchor=(1.02, 1.0,), borderaxespad=0)

plt.scatter(a[0], a[1], color = 'fuchsia')
plt.text(a[0], a[1], 'a')
plt.scatter(b[0], b[1], color = 'fuchsia')
plt.text(b[0], b[1], 'b')
plt.scatter(t1[0], t1[1], color = 'lime')
plt.text(t1[0], t1[1], 't1')
plt.scatter(t2[0], t2[1], color = 'lime')
plt.text(t2[0], t2[1], 't2')

f:id:datascience_findings:20210512225514p:plain

plt.scatter(a[0], a[1], color = 'fuchsia')
plt.text(a[0], a[1], 'a')
plt.scatter(b[0], b[1], color = 'fuchsia')
plt.text(b[0], b[1], 'b')
plt.scatter(t1[0], t1[1], color = 'lime')
plt.text(t1[0], t1[1], 't1')
plt.scatter(t2[0], t2[1], color = 'lime')
plt.text(t2[0], t2[1], 't2')
plt.annotate('', xy = (apt2[0], apt2[1]), xytext = (a[0], a[1]), arrowprops = dict(width = 0.5, headwidth = 7, color = 'black'))
plt.annotate('', xy = (apt1[0], apt1[1]), xytext = (a[0], a[1]), arrowprops = dict(width = 0.5, headwidth = 7, color = 'black'))
plt.text(apt1[0], apt1[1], 'pt1')
plt.text(apt2[0]-50, apt2[1], 'pt2')
plt.xlim([0, 500])
plt.ylim([0, 500])
plt.title('pt1 and pt2')

f:id:datascience_findings:20210512225648p:plain
上記図から、pt1およびpt2を求めることができました。しかし、今回選択したt2は線分ab上で垂直に交わりません。no man' landの対象となり得るのは線分ab上で垂直に交わるt1のようなプロットです。これを満たす条件として、下記2つを考えます。
1.線分aptの長さが線分abより短い
2.線分aptと線分abが同方向を向く
また、点aに対しマジョリティクラス全てを点tとみなすわけではありません。論文の通り、Vmという空間を定義し、Vm内のプロットに対しno man's landを決定していきます。上記例では、t2を含めVmを定義した場合、データ数が膨大となりオーバーサンプリング能力と計算量が見合わなくなってしまいます。そのため、KNN(K-近傍法)を用いて適切な範囲のVmを設定する必要があります。
これらを踏まえた上でGICAPsを実装します。

GICaPSの実装

GICaPSを実装するため、まずオーバーサンプリングの対象となるデータ(origin)を用意します。上記までの例で使用したデータを今回も用います。
実装したプログラムが下記となります。
まず、k=5のKNNを用い全マイノリティデータに空間Vmを定義しました。次に、k=5のデータの中に含まれるマイノリティクラスのデータを除外しました。(no man's landに関係するのはマジョリティデータのみであるためです) そして残ったデータで線分ab上で垂直に交わるデータがtとなります。
no man's landのlは、tのうち線分atの距離が最大となるtmaxと最低となるtminを距離マトリックスから抽出し線分ptminと線分ptmaxの差とします。 この作業を全てのマイノリティクラスのデータに対し適用することで、全てのマイノリティクラス間(線分ab上)のlを求めることができます。線分abに対する線分lの大きさの比がその線分ab間でオーバーサンプリングするデータ数の割合となります。 このようにしてGICaPSでは、クラス間のデータ位置を考慮しオーバーサンプリングすることが可能となります。

H = 300 # オーバーサンプリングするデータ数
p = 1 # no man's landの範囲調整パラメータ
# クラスごとのインデックスを取得
index_0 = [i for i, label in enumerate(origin['class_label']) if label == 0]
index_1 = [i for i, label in enumerate(origin['class_label']) if label == 1]

# class0のデータを抽出
origin_0 = origin[origin['class_label'] == 0]
origin_0np = np.array(origin_0.drop('class_label', axis = 1))
# class1のデータを抽出
origin_1 = origin[origin['class_label'] == 1]
origin_1np = np.array(origin_1.drop('class_label', axis = 1))
# オーバーサンプリングの対象となるデータ
origin_list = np.array(origin.drop('class_label', axis = 1))
# 全てのデータ間の距離マトリックスを作成
distance_list = distance.pdist(origin_list)
distance_matrix = squareform(distance_list)


# knn空間(k=5)を作成
knn_list = [] # class0でknn5個を抽出したもの
S_list = []    
# knnの最大距離、最小距離の要素を抽出
for i in range(len(index_0)): # i  = a 少数クラスのデータ数分繰り返す
    idx_6 = np.argpartition(distance_matrix[index_0[i]], 6)[:6] # 行で距離が短い下位6つの要素のインデックスを取得
    idx_5 = idx_6[~(idx_6 == i)] # matrixの対角線は0なので除外する
    knn_list.append(idx_5)
    # knn内のクラス1のみを抽出
    knn_list_majority = [] # 
    for knn_list_row in knn_list:
        knn_majority = [index for i, index in enumerate(knn_list_row) if index not in index_0] # 少数派を除外
        knn_list_majority.append(knn_majority)
    Sa_list = []
    for j in range(len(index_0)): # j = b
        a = origin_0np[i].reshape(2, 1)
        b = origin_0np[j].reshape(2, 1)
        ab = b - a
        ab_length = np.linalg.norm(ab)
        # tmax, tmin を求める
        distance_row = distance_matrix[index_0[i]]
        knn_distance_row = distance_row[knn_list_majority[i]]
        if len(knn_distance_row) > 0: # knnにクラス1が含まれている場合
            idx_list = np.argpartition(knn_distance_row , kth = len(knn_distance_row)-1)
            ab_element = []   # 線分ab上にあるtを抽出 
            for idx in idx_list:
                knn_idx = knn_list_majority[i][idx] 
                t = origin_list[knn_idx].reshape(2, 1)
                at = t - a
                at_internal = np.dot(at.T, ab)
                at_length = np.linalg.norm(at)
                theata = np.arccos(at_internal/(at_length * ab_length))
                
                if (0 < theata < math.pi /2)  and at_length / ab_length < 1: # 線分atと線分abのなす角が0以上90度以下(abと同じ向き)かつ絶対値が1以下(abの線分内)の場合
                    ab_element.append(knn_idx)
   
            if len(distance_matrix[i][ab_element]) > 1: # knnにクラス1が2つ以上含まれている場合 (1つしかないとl=0となりnomanslandなくなる)
                idx_max = np.argpartition(distance_matrix[i][ab_element], kth = len(distance_matrix[i][ab_element])-1)[-1] # knnの範囲内で距離が最も大きいmajority_indexを取得
                idx_min = np.argpartition(distance_matrix[i][ab_element], kth = len(distance_matrix[i][ab_element])-1)[0]
                knn_idx_max = ab_element[idx_max]
                knn_idx_min = ab_element[idx_min]
                tmax = origin_list[knn_idx_max].reshape(2, 1)
                tmin = origin_list[knn_idx_min].reshape(2, 1)
                atmax = tmax - a
                atmin = tmin -a
                ptmax = (np.dot(ab.T, atmax) * ab) / np.dot(ab.T, ab) 
                ptmin = (np.dot(ab.T, atmin) * ab) / np.dot(ab.T, ab)
                aptmax = ptmax + a
                aptmin = ptmin + a
                l = np.linalg.norm(aptmax - aptmin) # 線分ab上のno man's land の区間
                S = p * (ab_length - l) # 線分ab上のオーバーサンプリング対象区間の大きさ
                Sa_list.append(S)  # aに対する全てのbでのSをSa_listに追加

            else:
                S = p * ab_length  #tmax, tminが存在しない場合は線分ab自体がSとなる
                Sa_list.append(S) # Sa_listは任意のaに対するSを格納
        else:
            S = p * ab_length 
            Sa_list.append(S)
    S_list.append(Sa_list) # S_lsitは全てのa, bのSを格納

Sa_sum_list = [] # ΣVmSを算出
for i in S_list:
    Sa_sum_list.append(sum(i))

Nm_list = [] # 任意のaでいくつオーバーサンプリングするか算出
S_sum = sum(Sa_sum_list) # ΣxmΣVmSを算出
for Sa_sum in Sa_sum_list:
    Nm = H * Sa_sum / S_sum
    Nm_list.append(Nm)

Nv_list = [] # 任意の線分abでいくつオーバーサンプリングするか計算
for i in range(10):
    Nv_list_a = []
    for j in range(10):
        Nv = (Nm_list[i] * S_list[i][j]) / Sa_sum_list[i]
        Nv_list_a.append(Nv)
    Nv_list.append(Nv_list_a)
    
X_new = [] # 新しいデータを生成
for i in range(10):
    for j in range(10):
        Nv = Nv_list[i][j]
        N = round(Nv)
        a = np.array([[class0_a0[i]], [class0_a1[i]]])
        b = np.array([[class0_a0[j]], [class0_a1[j]]])
        ab = b - a
        for y in range(N):
            r = np.random.random(a.shape[0])
            rm = np.array([[r[0]], [r[1]]])
            X = a + y * (ab / Nv) + rm
            X_new.append(X)

GICaPSとsmoteの比較

マイノリティデータとマジョリティデータのオリジナルデータを下図に示します。 f:id:datascience_findings:20210515004825p:plain

オリジナルデータに対しsmoteとGICaPSでそれぞれオーバーサンプリングした結果が下図になります。

f:id:datascience_findings:20210515004809p:plainf:id:datascience_findings:20210515004844p:plain

GICaPSとsmoteを重ねた図が下図になります。 f:id:datascience_findings:20210515010736p:plain 上図より、GICaPSでは、クラス間のデータ位置が考慮されているため、マジョリティクラスが存在する範囲のオーバーサンプリング数はsmoteより少ないことがわかります。また、smoteでは、マイノリティクラスデータを結ぶ線分上かつマイノリティデータ間でKNNが適用されるため、オーバーサンプリングした結果が不均一となっていました。しかし、GICaPSでは、全てのマイノリティデータ間でオーバーサンプリングを行うため不均一性が和らいでいることが確認できます。
以上より、GICaPSはオーバーサンプリング手法として有効であると考えます。しかし、予測に疑似データを使用する危険性を無くせたわけではありません。オリジナルデータおよびオーバーサンプリング後のデータを確認することが必要だと考えます。

Random Forest(Classifier)

今回の記事では、分類器の一つである決定木を応用したランダムフォレストについて学んだ知見をまとめます。

ランダムフォレストのアルゴリズム

ランダムフォレストとは決定木の一つの手法です。決定木とは前回の交差検証の記事で示したようにノードのサンプルの不純度を最小とするようにサンプルを分類する手法でした。しかし、決定木では一つの木のみしか作成されず、分類する基準は下図に示すように当該ルートのサンプルによって一意に定められてしまうため、未知データに対する汎化性が乏しくなるという問題があります。(左側ノードのperFare≦26.144の分類基準は既知データにフィットしすぎている可能性がある)

f:id:datascience_findings:20210320094018p:plain
決定木による分類

そこで、決定木を発展させたランダムフォレストという手法により、決定木を複数作成し、決定木間で多数決を取ることで汎化性を向上させた分類基準を得ることが可能となりました。また、ランダムフォレストの決定木はブートストラップによる再標本化により作成されます。

ブートストラップ

ブートストラップとは再標本化の一つの手法であり、人為的に母集団と似た標本を作成するという考え方です。具体的には、母集団からランダムに要素を抽出、またその際、要素の重複を許容します。このようにして、任意サイズの標本を作成することで母集団と似た傾向を持つがそれぞれ異なる標本群を作成することが可能となります。(ランダムフォレストの場合、標本サイズは母集団サイズの0.62倍と決まっています)

ランダムフォレストの実装

今回もいつもと同様タイタニック号のデータセットを用いランダムフォレストを実装します。 まずは、特徴量エンジニアリングをします。

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

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


train['Parch'] = train['Parch'].astype(int)
train['SibSp'] = train['SibSp'].astype(int)
train['Familysize'] = train['Parch'] + train['SibSp'] + 1
test['Parch'] = test['Parch'].astype(int)
test['SibSp'] = test['SibSp'].astype(int)
test['Familysize'] = test['Parch'] + train['SibSp'] + 1

test_Fare = train[(train['Pclass'] == 3) & (train['Age'] > 40)]['Fare'].mean()
test['Fare'] = test['Fare'].fillna(test_Fare) 

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


train['honorific'] = train['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
train['honorific'] = train['honorific'].str.replace(' ', '')
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)

test['honorific'] = test['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
test['honorific'] = test['honorific'].str.replace(' ', '')
honorific_list = test[test['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))
test['Age'][np.isnan(test['Age'])] = age_list
test['Age'] = test['Age'].fillna(28)
test['Age'] = test['Age'].astype(int)

説明変数と目的変数を定義し、訓練データを訓練データと検証データに分割します。

trainE = train.loc[:, ['perFare', 'Pclass', 'Embarked', 'Sex', 'Familysize']]
trainO = train['Survived']
# データの分割
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(trainE, trainO, test_size=0.25, random_state=42)

ホールド・アウトにより分割した訓練データに対しランダムフォレストを適用します。そして、学習した決定木を用いて、オリジナルデータ(母集団)、学習に用いた訓練データ、検証データ(未知データ)に対し予測精度を算出しました。

# ランダムフォレストの作成
from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(random_state = 42)
clf = rf.fit(x_train, y_train)
print('n_estimators : 10')
print('score(オリジナルデータ): {:.3%}'.format(clf.score(trainE, trainO)))
print('score(訓練データ0.75): {:.3%}'.format(clf.score(x_train, y_train)))
print('score(検証データ0.25): {:.3%}'.format(clf.score(x_test, y_test)))
[out]
n_estimators : 10
score(オリジナルデータ): 89.675%
score(訓練データ0.75): 92.665%
score(検証データ0.25): 80.717%

ここで、前回作成した決定木と予測精度の比較を行います。

# 決定木の作成
my_tree = tree.DecisionTreeClassifier()
clf = my_tree.fit(x_train, y_train)

print('score(オリジナルデータ): {:.3%}'.format(clf.score(trainE, trainO)))
print('score(訓練データ0.75): {:.3%}'.format(clf.score(x_train, y_train)))
print('score(検証データ0.25): {:.3%}'.format(clf.score(x_test, y_test)))
[out]
score(オリジナルデータ): 89.562%
score(訓練データ0.75): 92.665%
score(検証データ0.25): 80.269%

ランダムフォレストと決定木ではほとんど予測精度に差違を確認できませんでした。要因の一つとしてランダムフォレストのデフォルトの設定では、決定木が10個しか作成されないため、差が生じなかったのではないかと考えられます。

そこで、作成する決定木の数と検証データの予測精度の関係を算出しました。

rfclf = []
for i in range(5, 100, 5):
    rf = RandomForestClassifier(n_estimators = i, random_state = 42)
    clf = rf.fit(x_train, y_train)
    rfclf.append(clf.score(x_test, y_test))

dx = 0.1
xmin, xmax = 0, 100
x = np.arange(xmin, xmax+dx, dx)
plt.plot(range(5,100,5), rfclf)
plt.hlines(0.80269, xmin, xmax, "blue", linestyles='dashed') 
plt.title('Relationship between n_estimators and prediction accuracy')
plt.xlabel('n_estimators')
plt.ylabel('Verification Score')

f:id:datascience_findings:20210322221051p:plain
図中の破線は単純な決定木の予測精度を表しています。この結果から、たしかに初期段階(n_estimators <25)では、作成する決定木の数により予測精度が向上する傾向が確認できます。しかし、それ以降数を増加させても予測精度は低下しました。これは母集団に類似した標本群に対し過学習してしまったためだと考えられます。そのため、ランダムフォレストで作成する決定木の数は十分に検討する必要があります。
実際にn_estimators = 20のときのランダムフォレストを可視化してみました。このようにそれぞれの決定木の形状変化が見て取れるため、ランダムフォレストを適用する意味を視覚的にも認識できます。 f:id:datascience_findings:20210322214757g:plain

feature_importances

ランダムフォレストの機能の一つに、feature_importancesがあります。これは、説明変数がどの程度目的変数に影響しているかを示す特徴量重要度です。feature_importancesはルートノードが分岐する際の不純度の変化量を全ノードに対し計算し総和で除すことで任意の特徴量の重要度を算出しています。
実際には、不純度変化量は下記の式で計算され f:id:datascience_findings:20210322225149p:plain
feature_importancesは下記式で計算されます。 f:id:datascience_findings:20210322225212p:plain
実際にpythonで実装します。

rf = RandomForestClassifier(n_estimators  = 25, random_state = 45)
clf25 = rf.fit(x_train, y_train)
features = trainE.columns
importances = clf.feature_importances_
indices = np.argsort(importances) # importancesを並び替える


plt.barh(range(len(indices)), importances[indices], align='center')
plt.yticks(range(len(indices)), features[indices])
plt.show()

f:id:datascience_findings:20210322225921p:plain
この結果から、最も重要度が高い特徴量はperFareであることがわかりました。しかし、実際にはランダムフォレストでの特徴量重要度の計算は分岐に使用される回数が多い特徴量ほど不純度変化量が増加してしまいます。そして、perFareという量的データ(Familysizeも量的データだがその範囲は小さいため除く)は分岐の基準となりやすいです。反対に、質的データであるSexはトップの分岐にて0か1で完全に分類されるため決定木での登場は一度のみしかないという不利な点があります。そのため、feature_importancesは特徴量選定の一つのメソッドですが、変数の種類により重要度に差が生じることは留意する必要があります。

交差検証

今回は機械学習モデルの作成時に必須な考え方である交差検証について学んだことをまとめます。

学習モデルの汎化性

交差検証の役割は、学習モデルの汎化性を向上させることです。汎化性とは、学習に用いていないデータ(未知データ)を予測する能力です。一般に、学習モデルの汎化性を検証するために、訓練(train)データと検証(validation)データと試験(test)データを用意します。検証データを用意しない場合、未知データに対する予測精度を確かめることができません。機械学習では過学習という訓練データにフィットし過ぎてしまう危険があり、その場合、訓練データとテストデータで予測精度に大きな乖離が生じます。 そこで、汎化性を検証するため検証データ(学習に用いず、正解ラベルがあるデータ)を用意します。しかし、単純な方法(訓練データの3割を検証データとする等)により検証データを用意しただけでは、ランダムなデータ選択であってもデータのばらつきにより十分に汎化性を評価することはできません。交差検証とは、訓練データおよび検証データの分割を反復し検証データの違いによる予測精度のばらつきを評価することであり、より正確に学習モデルの汎化性を検証することができます。
検証データを用いた汎化性の評価手法は様々あるため、主要なものを今回取り上げます。

ホールド・アウト検証(train_test_split)

sklearnライブラリのモジュールの一つであるtrain_test_splitは、データ(訓練用データ)をtrain data とtest dataに分割する機能を持ちます。実際に適用します。今回もkaggleのタイタニック号生存者予測のデータを用います。まず、特徴量を生成します。

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

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


train['Parch'] = train['Parch'].astype(int)
train['SibSp'] = train['SibSp'].astype(int)
train['Familysize'] = train['Parch'] + train['SibSp'] + 1
test['Parch'] = test['Parch'].astype(int)
test['SibSp'] = test['SibSp'].astype(int)
test['Familysize'] = test['Parch'] + train['SibSp'] + 1

test_Fare = train[(train['Pclass'] == 3) & (train['Age'] > 40)]['Fare'].mean()
test['Fare'] = test['Fare'].fillna(test_Fare) 

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


train['honorific'] = train['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
train['honorific'] = train['honorific'].str.replace(' ', '')
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)

test['honorific'] = test['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
test['honorific'] = test['honorific'].str.replace(' ', '')
honorific_list = test[test['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))
test['Age'][np.isnan(test['Age'])] = age_list
test['Age'] = test['Age'].fillna(28)
test['Age'] = test['Age'].astype(int)

学習に使用する説明変数をtrainEとし、正解ラベルをtrainOとしました。

trainE = train.loc[:, ['perFare', 'Pclass', 'Embarked', 'Sex', 'Familysize']]
trainO = train['Survived']

train_test_splitにより、訓練データと検証データに分割します。分割サイズはデフォルトの3:1の割合とします。

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(trainE, trainO, random_state=42)

データの中身(もともとの訓練データ、訓練データを分割した訓練データ、検証データ)を確認します。

print('元の訓練データ')
print(trainE.describe().round(2))
print()
print('訓練データ')
print(x_train.describe().round(2))
print()
print('検証データ')
print(x_test.describe().round(2))

f:id:datascience_findings:20210316225239p:plain
3つのデータの特徴量ごとの平均・標準偏差に大きな差違は見られませんでした。つまり、予測精度に差は生じないのではないかと考えられます。そこで、分割後の訓練データに対し、決定木を学習させ、そのモデルを3つのデータに適用しました。

# 決定木の作成
my_tree = tree.DecisionTreeClassifier()
clf = my_tree.fit(x_train, y_train)

# 正解率
print('score(オリジナルデータ): {:.2%}'.format(clf.score(trainE, trainO)))
print('score(訓練データ): {:.2%}'.format(clf.score(x_train, y_train)))
print('score(検証データ): {:.2%}'.format(clf.score(x_test, y_test)))
[out]
score(オリジナルデータ): 89.67%
score(訓練データ): 92.66%
score(検証データ): 80.72%

学習に用いた訓練データおよび、訓練データが75%含まれるオリジナルデータはいずれも約90%の正解率でした。一方、未知データである検証データでは約80%と大幅に正解率が低下し事前の予測とは異なる結果となりました。
つまり、データフレームごとの特徴量の統計量の差異を確認しただけでは予測精度を推定することはできないということです。

予測精度と分割サイズの影響

訓練データのサイズを0.75としていましたが、0.5、0.9と変更し予測精度を比較しました。

x_train, x_test, y_train, y_test = train_test_split(trainE, trainO, test_size=0.5, random_state=42)
clf = my_tree.fit(x_train, y_train)
print('score(オリジナルデータ): {:.2%}'.format(clf.score(trainE, trainO)))
print('score(訓練データo.5): {:.2%}'.format(clf.score(x_train, y_train)))
print('score(検証データ0.5): {:.2%}'.format(clf.score(x_test, y_test)))
print()
x_train, x_test, y_train, y_test = train_test_split(trainE, trainO, test_size=0.1, random_state=42)
clf = my_tree.fit(x_train, y_train)
print('score(オリジナルデータ): {:.2%}'.format(clf.score(trainE, trainO)))
print('score(訓練データo.9): {:.2%}'.format(clf.score(x_train, y_train)))
print('score(検証データ0.1): {:.2%}'.format(clf.score(x_test, y_test)))
[out]
score(オリジナルデータ): 85.75%
score(訓練データ0.5): 94.16%
score(検証データ0.5): 77.35%

score(オリジナルデータ): 91.36%
score(訓練データ0.9): 92.26%
score(検証データ0.1): 83.33%

この結果から、train dataサイズが0.5と小さくなるとオリジナルデータと検証データの予測精度が低下し、train dataサイズが0.9と大きくなると予測精度は向上しました。これは、訓練時により多くのデータで学習することで精度が向上するためです。しかし、検証データの割合が小さいことからこの検証結果を真の汎化性と判断してよいか疑問があります。(実際にこの学習モデルでtest dataを予測しkaggleに提出した結果、正解率は80%に届かなかったことから過度な汎化性を示していると考えられます)

交差検証

上記のホールド・アウト検証(一度の検証)では、正確な汎化性を評価することができませんでした。そこで、検証を繰り返す交差検証という手法が多くの場合用いられます。 主要なK-分割交差検証と層化k分割交差検証を取り上げます。

K-分割交差検証

K-分割交差検証とは、元データをk個の訓練データと検証データに分割しk個全てのデータセットに対し検証を行う手法です。これにより、データの組み合わせにより予測精度に差が生じるかがわかります。(差が生じる場合は、データのばらつきが大きく、汎化性に乏しいと考えられます)
まず、実際にデータを分割します。

from sklearn.model_selection import KFold

kf3 = KFold(n_splits=3, shuffle=True, random_state = 42)
print('分割数:3')
for kf_train, kf_test in kf3.split(train):
    print('訓練データサイズ:{} 検証データサイズ:{}'.format(len(kf_train), len(kf_test)))
print()    
kf5 = KFold(n_splits=5, shuffle=True, random_state = 42)
print('分割数:5')
for kf_train, kf_test in kf5.split(train):
    print('訓練データサイズ:{} 検証データサイズ:{}'.format(len(kf_train), len(kf_test)))
[out]
分割数:3
訓練データサイズ:594 検証データサイズ:297
訓練データサイズ:594 検証データサイズ:297
訓練データサイズ:594 検証データサイズ:297

分割数:5
訓練データサイズ:712 検証データサイズ:179
訓練データサイズ:713 検証データサイズ:178
訓練データサイズ:713 検証データサイズ:178
訓練データサイズ:713 検証データサイズ:178
訓練データサイズ:713 検証データサイズ:178

このように設定した分割数kに応じて、1/kの大きさの検証データがランダム(shuffle = Trueの場合)に作成されます。また、検証データを入れ替えるためk回データ分割が行われていることがわかります。
実際に分割されたデータの中身を確認します。

kf3 = KFold(n_splits=89, shuffle=True, random_state = 42)

print('分割数:89')
for i, (kf_train, kf_test) in enumerate (kf3.split(train)):
    print('{}番目の検証データ:{}'.format((i+1), kf_test))
    if i == 9: 
        break
[out]
分割数:89
1番目の検証データ:[ 39 136 137 208 290 300 333 439 709 720 840]
2番目の検証データ:[110 244 294 344 485 621 653 696 853 886]
3番目の検証データ:[ 30  63 192 396 447 538 673 682 819 877]
4番目の検証データ:[ 23 120 141 198 204 235 620 739 793 842]
5番目の検証データ:[ 67  86 210 350 362 448 477 659 790 837]
6番目の検証データ:[ 44  70 211 280 299 360 408 585 802 820]
7番目の検証データ:[ 72 168 196 312 422 426 446 532 591 772]
8番目の検証データ:[254 266 311 357 539 605 767 833 835 889]
9番目の検証データ:[ 66 165 174 215 250 309 319 493 778 822]
10番目の検証データ:[ 76 281 321 327 338 388 541 545 625 712]

このように分割された中身には要素番号のみが格納されています。また、検証データが重ならないように抽出されていることがわかります。(89分割のため最初の10個のみ表示)

続いて交差検証による予測精度のばらつきおよび予測の過程(決定木の可視化)を行います。

# 可視化に必要なモジュールをインポート
from sklearn.model_selection import GridSearchCV
import pydotplus as pdp

kf3 = KFold(n_splits=3, shuffle=True, random_state = 42)
print('分割数:3')
for i, (kf_train, kf_test) in enumerate (kf3.split(train)):
    # 訓練データ
    x_train = trainE.loc[kf_train]
    y_train = trainO.loc[kf_train]
    # 検証データ
    x_test = trainE.loc[kf_test]
    y_test = trainO.loc[kf_test]
    # 決定木を学習し予測
    clf = my_tree.fit(x_train, y_train)
    print('{}番目のscore:{:.2%}'.format((i+1), clf.score(x_test, y_test)))

    # 学習した決定木を可視化
     dot_data = tree.export_graphviz(clf, 
                                    out_file=None, # ファイルは介さずにGraphvizにdot言語データを渡すのでNone
                                    filled=True, # Trueにすると、分岐の際にどちらのノードに多く分類されたのか色で示してくれる
                                    rounded=True, # Trueにすると、ノードの角を丸く描画する
                                    feature_names=trainE.columns[:5], # 分類の基準となった特徴量を示す
                                    class_names=['Survived', 'Dead'], # 生存したか否かを示す
                                    )
    graph = pdp.graph_from_dot_data(dot_data)
    file_name = "./tree_visualization{}.png".format(i+1)
    graph.write_png(file_name)
[out]
分割数:3
1番目のscore:79.12%
2番目のscore:76.43%
3番目のscore:80.47%

<一度目の分割>決定木1 f:id:datascience_findings:20210320092917p:plain <1度目の分割(トップ切り取り)> f:id:datascience_findings:20210320094018p:plain <2度目の分割>決定木2 f:id:datascience_findings:20210320105209p:plain <3度目の分割>決定木3 f:id:datascience_findings:20210320105222p:plain

Kfoldにより3度分割した決定木を示しました。決定木の要素であるノードの一番上には分類する特徴量の基準が示されます。基準はジニ不純度により決定します。不純度が低い場合ジニ不純度は0に、高い場合には1に近づきます。不純度は下記の式で表されます。 f:id:datascience_findings:20210320094951p:plain
ここで、tは不純度を求めたいノード、cは目的変数(Survived or Dead)のクラス数、Nはデータのサンプル数、niはクラスごとに属するデータ数を表します。 つまり、トップのノードにおける不純度は下記の式で表されます。

# 不純度
gini = 1-((372/594)**2 + (222/594)**2)
gini
[out]
0.46811549841852873

そして、訓練データを最初にSex≧0.5(男性か女性か)で分類することで、2層目のノードではジニ不純度が0.305(女性)、0.4(男性)となり元の不純度より低くなっていることがわかります。決定木では、ジニ不純度が最も小さくなるような条件を選定し学習します。今回は、Kfoldにより3度データを分割していますがどの場合でも最初の分類では、Sexとなりました。つまり、SexはSurvivedと最も相関が高い特徴量であるということを示しており、これまでの記事で示した相関係数とも一致します。
また、2層目のノードで分類する特徴量は女性ではperFare、男性ではPclassと異なりました。特に男性のPclas≦2.5での分類は決定木の収束が最も早くSurvivedに強く寄与する特徴量であることがわかりました。 多変量解析の記事で示したように、説明変数をデータ全体に対し決定するのではなく、性別で分類した後に決定することが予測精度の向上に繋がった結果とも一致します。

決定木と予測精度の関係

また、決定木と予測精度の関係を考察します。予測精度は決定木2が決定木1,3と比較し低い結果となりました。これは決定木の幅に関係すると考えます。決定木の幅が広いということは、決定木の深さが増してもノードが収束していかないということを意味します。実際に決定木2では深さ8でノードが40ありましたが決定木3ではノードは28でした。つまり、決定木2の訓練データは全体的に不純度が高く分割されたデータであると解釈できます。 このように、分割するデータセットによって、予測精度に差が生じることがわかりました。モデルを学習する際には、交差検証により汎化性を検証することが重要であることが実際に確かめられました。

層化k分割交差検証

続いて、層化k分割交差検証を適用します。これは、分割するデータセットの目的変数の出現比率を一定とし分割する手法です。タイタニック号の全乗客の生存率は約38%でした。

np.round(train['Survived'].mean(),2)

[out]
0.38

先程の単純なK個への分割ではデータセットごとの生存率に差が生じてしまいます。仮に分割した訓練データ全てが生存者となった場合にはモデルが学習を行う必要がなくなります。そこで、層化k分割交差検証を適用することで、正解ラベルが揃っている場合の予測精度のばらつきを確認することができます。ばらつきが小さくなれば、選定した説明変数によってデータ全体をざっくり網羅できていると考えられます。 説明変数は先程までと同様perFare、Pclass'、Embarked、Sex、Familysizeの場合(CASE1)とEmbarkedのみとした場合(CASE2)で比較しました。

skf = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
score = []
print('CASE1')
for i, (skf_train, skf_test) in enumerate (skf.split(trainE, trainO)):
    # 訓練データ
    x_train = trainE.loc[skf_train]
    y_train = trainO.loc[skf_train]
    # 検証データ
    x_test = trainE.loc[skf_test]
    y_test = trainO.loc[skf_test]
    # 決定木を学習し予測
    clf = my_tree.fit(x_train, y_train)
    score.append(clf.score(x_test, y_test))

    print('{}番目のscore:{:.2%}'.format((i+1), clf.score(x_test, y_test)))
    print('{}番目の訓練データの生存者の比率:{:.2%}'.format((i+1), y_train.mean()))   
    print('{}番目の検証データの生存者の比率:{:.2%}'.format((i+1), y_test.mean()))       
    print()
print('予測精度の標準偏差:{:.2}'.format(np.std(score)))
CASE1
1番目のscore:78.11%
1番目の訓練データの生存者の比率:38.38%
1番目の検証データの生存者の比率:38.38%

2番目のscore:78.79%
2番目の訓練データの生存者の比率:38.38%
2番目の検証データの生存者の比率:38.38%

3番目のscore:79.12%
3番目の訓練データの生存者の比率:38.38%
3番目の検証データの生存者の比率:38.38%

予測精度の標準偏差:0.42
trainE2 = train.loc[:, ['Embarked']]
trainO2 = train['Survived']
score = []
skf = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
print('CASE2)
for i, (skf_train, skf_test) in enumerate (skf.split(trainE2, trainO2)):
    # 訓練データ
    x_train = trainE2.loc[skf_train]
    y_train = trainO2.loc[skf_train]
    # 検証データ
    x_test = trainE2.loc[skf_test]
    y_test = trainO2.loc[skf_test]
    # 決定木を学習し予測
    clf = my_tree.fit(x_train, y_train)
    score.append(clf.score(x_test, y_test))
    
    print('{}番目のscore:{:.2%}'.format((i+1), clf.score(x_test, y_test)))
    print('{}番目の訓練データの生存者の比率:{:.2%}'.format((i+1), y_train.mean()))   
    print('{}番目の検証データの生存者の比率:{:.2%}'.format((i+1), y_test.mean()))       
    print()
print('予測精度の標準偏差:{:.2}'.format(np.std(score)))
CASE2
1番目のscore:64.98%
1番目の訓練データの生存者の比率:38.38%
1番目の検証データの生存者の比率:38.38%

2番目のscore:61.95%
2番目の訓練データの生存者の比率:38.38%
2番目の検証データの生存者の比率:38.38%

3番目のscore:64.65%
3番目の訓練データの生存者の比率:38.38%
3番目の検証データの生存者の比率:38.38%

予測精度の標準偏差:1.36

この結果から、CASE2の場合は正解率が劣ることはもちろん、標準偏差もCASE1と比較し約3倍となりました。例え、目的変数が含まれる割合を等しく分割したとしても、説明変数が不十分であれば分割されるデータセットによって精度に差が生じるということだと思われます。反対に、予測精度のばらつきを十分に小さくなるまでは説明変数を再考する余地があるのではないかと考えます。

タイタニック号生存者予測 - 主成分分析(PCA)

本稿では、多変量解析で用いられる主成分分析(PCA)という手法を実際に適用し効果を確認することを目的とします。PCAとは、端的には多次元データが持つ情報を主成分ベクトルに変換する手法です。  PCAを用いるメリットは、元データ(特徴量)の情報要約することができるということです。これにより、重要な特徴量の選定、多次元データの次元圧縮を可能とします。 今回PCAを適用するデータは前回までと同様kaggleのタイタニック号生存者予測のデータを用います。

2変数でのPCA

データを主成分ベクトルに変換するためには、主成分を理解する必要があります。
そこで、簡略化のために2変数(AgeとFare)のみを使用します。まず、2変数の関係を散布図より確認します。

df = train.loc[:, ['Age', 'Fare']]
sns.regplot(x=df.Age, y=df.Fare)

f:id:datascience_findings:20210313195124p:plain
この散布図から横軸と縦軸のスケールが一致していないことがわかります。データ間の関係を語る上では、標準化という作業は基本的には欠かせません。例えば、異なる試験間で100点満点なのか200点満点なのかというようにデータスケールが異なれば1点の重みが変わります。標準化は、データの値と平均との差を取り、標準偏差で除すことで変換することができます。

# 標準化
dfs= df.iloc[:, :].apply(lambda x: (x-x.mean())/x.std(), axis=0)
sns.regplot(x=dfs.Age, y=dfs.Fare)

f:id:datascience_findings:20210313200847p:plain

続いて、この標準化データに対しPCAを適用し主成分ベクトルを得ていきます。
(第一)主成分ベクトルとは元データの情報を最も表現可能な新たな軸を意味しています。情報量が多いとは分散(ばらつき)が大きいことと同義です。つまり、PCAとは元データに対し分散が大きくなる主成分軸を探していけばよいのです。

これをわかりやすくグラフで表現しました。 主成分ベクトルは PC1 = w1x×Age+w1y×Fare ※w1x:第一主成分ベクトルのAgeの割合、w1y:第一主成分ベクトルのFareの割合
で表されます。 簡単な例とするため、AgeとFareは2つのデータのみでPCAを適用します。

# 2点のデータを抽出
x = (dfs['Age'][110:112]).round(2).reset_index(drop = True)
y = (dfs['Fare'][110:112]).round(2).reset_index(drop = True)

分散は、データの平均とのズレの2乗和をデータ数で除すことで求められます。今回はデータ数が2つなので下記のように求められます。 f:id:datascience_findings:20210620005934p:plain
また、w1xとw1yには下記の制約があります。 f:id:datascience_findings:20210620010401p:plain
これを踏まえてPC1を求めた結果が下のグラフです。

x_ave =x.mean()
y_ave =y.mean()
n = 2 # データ数
m = 30 # w1xの刻み数
w1x = 0
V_list = []
ims = []
for i in range(m):
    w1x = w1x + 1 /m
    w1y = (1 - w1x ** 2) ** (1/2)
    # 分散を計算
    V = ((x[0] * w1x + y[0] * w1y - (x_ave * w1x + y_ave * w1y)) ** 2 + \
           (x[1] * w1x + y[1] * w1y - (x_ave * w1x + y_ave * w1y)) ** 2 ) / n
    V_list.append(V)
    # データと主成分軸との直交点を計算
    x0 = (x[0] * w1x + y[0] * w1y) * w1x / (w1x ** 2 + w1y ** 2)
    y0 = w1y * x0 / w1x
    x1 = (x[1] * w1x + y[1] * w1y) * w1x / (w1x ** 2 + w1y ** 2)
    y1 = w1y * x1 / w1x
    
    plt.figure(figsize=(10, 10))
    plt.plot([-1.5, 1.5], [0, 0], 'black') # x軸
    plt.plot([0, 0], [-1.5, 1.5], 'black') # y軸
    plt.plot([-w1x,w1x], [- w1y, w1y]) # 主成分軸
    plt.plot([x0, x1], [y0, y1]) # 主成分軸上のデータ範囲
    plt.scatter(x, y) # データプロット
    plt.scatter(w1x, w1y) # 主成分ベクトルの重みをプロット
    draw_circle = plt.Circle((0, 0), 1, fill=False, ec='silver') # 半径1の円
    plt.gcf().gca().add_artist(draw_circle)
    plt.text(x[0] + 0.1, y[0] + 0.1, "({},{})".format(x[0],y[0]))
    plt.text(x[1], y[1], "({},{})".format(x[1],y[1]))
    plt.text(w1x, w1y, "({},{})".format(np.round(w1x, 2),np.round(w1y, 2)))
    plt.text(w1x / 3, w1y /3, "V = {}".format(V.round(3))) # 分散を表示
    plt.xlabel('Age')
    plt.ylabel('Fare')
    plt.savefig('img{}.png'.format(i))
    plt.show()

f:id:datascience_findings:20210620010847g:plain
上記グラフから主成分軸の回転に伴う分散の変化が見てとれます。図中のオレンジ線はデータと主成分軸が直交する点の長さを示しています。分散が大きいということはデータの広がりが大きいということがオレンジ線から読み取れます。
また、オレンジ点は(w1x, w1y)を表しており、半径1の円周上に乗っていることがわかります。上記でw1xおよびw1yに制約を設けたと記載しましたが、この理由も図から理解することができます。w1xおよびw1yは主成分軸の角度を表す値なので主成分軸(青線/オレンジ線)上であればどの値をとっても問題ありません。しかし、これら値が一意に定まりません。そこで、w1xおよびw1yというベクトルに大きさの制約を設けることで一意に定めることができるのです。(図の線上かつ円周上の制約)

上記例では、w1xを刻ませることで主成分軸を見つけることができましたが、数学的には分散共分散行列からラグランジュの未定乗数法を使用して求めます。ラグランジュの未定乗数法とは、等式の制約を利用して関数を最大化する値を求める方法です。PCAの導出としては理解すべき内容ですが、PCAのイメージは上記グラフであるということを抑えておけば数式を理解できなくてもPCAを使用することができます。

続いて、 上記例では、2点のみしか使用していませんでしたが、全データに対しPCAを適用します。
python の場合、sklearnにPCAモジュールがありますので、こちらを用います。

from sklearn.decomposition import PCA

#主成分分析の実行
pca = PCA()
pca.fit(dfs)
# データを主成分空間に写像
feature = pca.transform(dfs)

df_PC = pd.DataFrame(feature, columns=["PC{}".format(x + 1) for x in range(len(dfs.columns))])
sns.regplot(x=df_PC.PC1, y=df_PC.PC2)

f:id:datascience_findings:20210313214143p:plain
第一主成分、第二主成分に分解してプロットした結果が上図になります。このグラフは、元のグラフと形状は大きく変わっていないように見えます。というのも、2変数の主成分分析とは、元の軸を回転させただけに過ぎないためです。じゃあ、意味ないではないかと思われるかもしれませんが、重要なのは、今回作成したPC1(第一主成分)はAgeとFare両方の情報を持つ軸であるということです。今回の場合、PC1の寄与率は下に示すように0.54とそこまで大きくありませんが、仮に寄与率が0.9(90%)あればPC1という特徴量だけで十分に元データを語ることができます。

 # 寄与率
pd.DataFrame(pca.explained_variance_ratio_, index=["PC{}".format(x + 1) for x in range(len(dfs.columns))],columns=['Contribution rate'])
[out]
    Contribution rate
PC1 0.548031
PC2 0.451969

このようにPCAが強力なツールであるのは、主成分に複数の特徴量のデータを持たせることで新たな特徴量を生成することができるということです。また、それと同時に、寄与率が低い主成分を取り除いたとしても元データに影響しないため、その分特徴量の次元を圧縮することができます。

n次元データでのPCA

 本題のタイタニック号のデータに対してPCAを適用します。基本的に2変数の場合と同様です。 まずは、タイタニック号データを数値データに変換します。

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

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

train['Parch'] = train['Parch'].astype(int)
train['SibSp'] = train['SibSp'].astype(int)
train['Familysize'] = train['Parch'] + train['SibSp'] + 1

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

train['honorific'] = train['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
train['honorific'] = train['honorific'].str.replace(' ', '')
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)

今回PCAに用いる特徴量を抽出し新しいデータフレームを作成します。

df2 = train.loc[:, ['Age', 'perFare', 'Pclass', 'Embarked', 'Sex', 'Familysize']]

データを標準化し、PCAを実行します。

df2s= df2.iloc[:, :].apply(lambda x: (x-x.mean())/x.std(), axis=0)

pca = PCA()
pca.fit(df2s)

feature = pca.transform(df2s)

主成分ごとの寄与率および累積寄与率は下図表のようになりました。

寄与率
pd.DataFrame(pca.explained_variance_ratio_, index=["PC{}".format(x + 1) for x in range(len(df2s.columns))],columns=['Contribution rate'])
[out]
    Contribution rate
PC1 0.303899
PC2 0.227552
PC3 0.153157
PC4 0.124604
PC5 0.114656
PC6 0.076133
# 累積寄与率を図示する
import matplotlib.ticker as ticker
plt.gca().get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
plt.plot([0] + list( np.cumsum(pca.explained_variance_ratio_)), '-o')
plt.xlabel('Number of principal components')
plt.ylabel('Cumulative contribution rate')
plt.show()

f:id:datascience_findings:20210313233316p:plain
ちなみに主成分ベクトルは用いた特徴量の数だけ作成されます。(今回は6個)これは、特徴量の数だけ軸が存在するため、その数だけ固有値が存在するからです。また、主成分ベクトルは分散が大きい(情報の損失が少ない)方向から決定していきます。そのため、第一主成分が最も寄与率が大きく、第六主成分の寄与率が最も小さくなります。寄与率が小さいということは元データに対し影響が小さいということなので、除くことができます。一般には累積寄与率が8〜9割を越えれば十分に元データを説明できるみたいです。今回の場合、累積寄与率が8割となる第四主成分までを考慮することとします。
主成分ベクトルと元の特徴量との関係は固有ベクトルから読み取ることができます。

# PCA の固有ベクトル
Eigenvector = (pd.DataFrame(pca.components_, columns=df2s.columns[0:], index=["PC{}".format(x + 1) for x in range(len(df2s.columns))])).round(2)
Eigenvector
   Age perFare Pclass Embarked Sex Familysize
PC1 0.40   0.56   -0.59  0.33    0.10  -0.25
PC2 -0.44 0.19    -0.11  0.24    0.63    0.56
PC3 0.32   0.02   -0.34  -0.80   0.14  0.35
PC4 -0.37  0.42   -0.15  -0.04   -0.72   0.36
PC5 -0.52  0.38   0.05   -0.44   0.21 -0.58
PC6 0.37    0.57  0.71 -0.05   0.08    0.17

考察

第一主成分が含む主な特徴量は、Pclass、perFareであることがわかります。第二主成分ベクトルはSexとFamilysize。第三主成分ベクトルはEmbarked、第四主成分ベクトルはSexが支配的です。また、Ageは全ての主成分ベクトルで約0.3~0.5 の範囲で強い特徴は見られませんでした。これら結果から、前回までのタイタニック号の分析でわかっていたSurvivedとSexの相関が最も強いこと(Sexは第二、第四主成分ベクトルの主変数であり、合計するとこれら主成分ベクトルの寄与率は0.37であり最も強い影響を持つと言える)、一方AgeはSurvivedとの相関が小さいこと(Ageは全主成分ベクトルで一定値存在しているがこれは一様分布に近い分布であるためだと考えられる。つまり、データに特徴がないと言える)など一致した結果であると言えます。
このように主成分分析は特徴量選定に有効です。今回のように、特徴量が少ない場合は次元圧縮の効果も少ないですが、計算コストの削減にも有効です。

独立性の検定 - 特徴量選定

特徴量選定の方法として、前回までは目的変数と説明変数の相関係数に着目してきました。しかし、相関係数は説明変数として妥当だと判断する基準値は持ちません。

そこで、今回は統計的手法である独立性の検定を用いたいと思います。独立性の検定とは、異なる事象の発生に関係性があるかを判定する検定のことです。
つまり、特徴量選定を行う上で、目的変数と説明変数に対し独立性の検定を適用することで、帰無仮説(変数間は独立である)が採択されれば、特徴量として相応しくないと判断し、帰無仮説が棄却されれば、特徴量に選定すべきと判断可能ではないかと考えました。

独立性の検定の適用

今回のデータも前回までと同様、kaggleのタイタニック号生存者予測を用います。
独立性の検定の適用は、相関係数が大きかったSurvivedとSexおよび相関係数が小さいSurvivedとEmbarkedに適用しその結果を比較します。

SurvivedとSex

pythonで独立性の検定を実装していきます。 まずは、検定に必要なクロス集計表を作成します。具体的には、性別の独立性を検定するので、女性かつ生存した人数・女性かつ死亡した人数・男性かつ生存した人数・男性かつ死亡した人数を求めます。

a = (train[train['Survived'] == 1]['Sex'] ==1).sum()
b = (train[train['Survived'] == 1]['Sex'] ==0).sum()
c = (train[train['Survived'] == 0]['Sex'] ==1).sum()
d = (train[train['Survived'] == 0]['Sex'] ==0).sum()

df = pd.DataFrame([[a,b],   
                   [c,d]])
df
[out]
   0   1
0  233 109
1  81  468

上記の値が実測値(実測度数)となります。
独立性の検定では理論度数(独立と仮定したときの値)も必要です。つまり、女性(314人)と男性(577人)それぞれの生存率は全乗客(891人)の生存率と一致するということです。理論度数のクロス集計表は以下のようになります。

#全乗客の生存率
survived = (train['Survived'] == 1).sum()
all = len(train)
rate = survived/all

aa = (a+c)*rate
cc = (a+c)*(1-rate)
bb = (b+d)*rate
dd = (b+d)*(1-rate)

dff = pd.DataFrame([[aa,bb],   
                   [cc,dd]])
dff.round(2)
[out]
      0    1
0  120.53 221.47
1  193.47 355.53

独立性の検定では検定量χ^ 2値が使用されます。χ^ 2値の求め方は\frac{(実測度数-理論度数)^ 2}{理論度数}  を各要素計算した和で表されます。

((a-aa)**2)/aa +((b-bb)**2)/bb +((c-cc)**2)/cc + ((d-dd)**2)/dd
[out]
263.05057407065567

となりました。χ^ 2検定での今回の自由度は1です。しかし、χ^ 2=236は付表に記載されていないのでp値はわかりませんでした。
そこで、ここまでは地道に計算してきましたが、scipyの統計モジュール機能を用いて計算します。

import scipy.stats as st
x2, p, dof, expected = st.chi2_contingency(df,correction=False)

print(f'p値    = {p :.60f}')
print(f'カイ2乗値 = {x2:.2f}')
print(f'自由度   = {dof}')
print(expected)
[out]
p値    = 0.000000000000000000000000000000000000000000000000000000000037
カイ2乗値 = 263.05
自由度   = 1
[[120.52525253 221.47474747]
 [193.47474747 355.52525253]]

計算した結果、p値は10の−60乗ほどと物凄く小さな値となりました。一般に独立性の検定での有意水準は1か5%ほどであるので物凄く小さいです。つまり、この結果から、生存率と性別の関係の実測度数は独立と仮定した場合に対し、10の-60乗ほどの確率でしか起こり得ない事象であることから、帰無仮説を棄却しこの2つの変数には関係があると結論付けられます。
この結果は、相関係数が0.5ほどあったことからも容易に想像できましたが、続いて相関係数が0.2ほどのSurvivedとEmbarkedに関して独立性を確認します。

SurvivedとEmbarked

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

o = (train[train['Survived'] == 0]['Embarked'] ==0).sum()
p = (train[train['Survived'] == 0]['Embarked'] ==1).sum()
q = (train[train['Survived'] == 0]['Embarked'] ==2).sum()
r = (train[train['Survived'] == 1]['Embarked'] ==0).sum()
s = (train[train['Survived'] == 1]['Embarked'] ==1).sum()
t = (train[train['Survived'] == 1]['Embarked'] ==2).sum()

df_embarked = pd.DataFrame([[o,r],   
                            [p,s],
                            [q,t]])
df_embarked
[out]
    0  1
0  427    217
1  47 30
2  75 95

Embarkedを上記のように数値データに変換し、クロス集計表を作成しました。
χ^ 2検定を適用した結果を下に示します。

x2, p, dof, e = st.chi2_contingency(df_embarked,correction=False)
print(f'p値    = {p :.10f}')
print(f'カイ2乗値 = {x2:.2f}')
print(f'自由度   = {dof}')
[out]
p値    = 0.0000008294
カイ2乗値 = 28.01
自由度   = 2

SurvivedとEmbarkedに対し適用したp値は約8×10^-7乗であり、SurvivedとSexと比較し大きくなっていることが確認できました。しかし、それでもなお独立性の検定で一般に定められる有意水準1%と比較すると非常に小さい値であり、EmbarkedはSurvivedに対し、独立ではないと判定することができそうです。

結論

 独立性の検定を相関係数が0.5であるSurvivedとSex、相関係数が0.2であるSurvivedとEmbarkedに対しそれぞれ適用した結果、有意水準を1%とすると検定統計量χ^ 2はどちらのp値も1%には程遠い小さい値となりました。この結果から、SexおよびEmbarkedはSurvivedに対しどちらも独立ではないと判断できます。
また、独立性の検定は主に、変数が独立であるはずのものが、独立ではないのではないかと疑わしい状況において適用されるものです。例えば、A店とB店でドーナツの質量が異なるのではないかなどです。つまり、今回のように相関係数が0.2もある変数間で独立性の検定を適用すれば、当然非常に小さい値になるということです。
 よって、当初の想定ははずれ独立性の検定を特徴量選定に適用するのはあまり有益ではないという結果となりました。しかし、相関係数と独立性の検定でのp値の関係性を掴むことができたことは大きな学びです。

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

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

 特徴量間の相関係数

 前回は、特徴量それぞれについて生存率が増加 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()

f:id:datascience_findings:20210227222250p:plain    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)

f:id:datascience_findings:20210228102439p:plain
グラフから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)]

f:id:datascience_findings:20210228104301p:plain
f:id:datascience_findings:20210228104252p:plain
Pclass = 1の場合、5000番台で新たに区分わけをすることで解決できそうです。

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

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

f:id:datascience_findings:20210228104752p:plain
f:id:datascience_findings:20210228104757p:plain
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)

f:id:datascience_findings:20210228135847p:plain
 このグラフから、Pclassが小さいほどFareは大きくなる傾向が確認できました。一方で、Pclass = 1かつ FareBand = 1となる乗客が一定数おり、これら2つの特徴量は完全に一致するわけではありませんでした。そこで、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)

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

欠損値補完

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

Ageの欠損値補完

 今回のデータセットにおいて、年齢を推測するヒントとして同乗した家族構成(SibSp、Parch、Familysize )があります。例えば、家族規模が1である場合、単身での乗船であることから、子どもではないことが推測されます。また、SibSp ≠ 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()]

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

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

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

結果

 f:id:datascience_findings:20210228222616p:plain

生存者予測

 説明変数を定めたので、実際に予測を行います。予測モデルはランダムフォレストを使用します。まず、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でした。だいぶ下の順位でした。

予測結果考察

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

f:id:datascience_findings:20210304213543p:plainf:id:datascience_findings:20210304213547p:plain

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

f:id:datascience_findings:20210304214901p:plainf:id:datascience_findings:20210304214906p:plain

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

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

f:id:datascience_findings:20210304220433p:plain
 Fclassは男性の場合、生存率が段々と増加する結果となりましたが、女性の場合にはFclass = 3のの乗客に死亡者が集中していました。そのため、女性のFclass区分は再振り分け要と判断します。また、Fclassの元の変数であったPclassとFareの変数も確認します。 f:id:datascience_findings:20210304221015p:plain
f:id:datascience_findings:20210304221021p:plain
 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)

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

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

 再予測

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

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()

f:id:datascience_findings:20210308210541p:plainf:id:datascience_findings:20210308210545p:plain


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

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でした。予測精度は性別で分割することで約3%向上しており有効であることが確認できました。

タイタニック号生存者予測 - 特徴量の扱い

 

 

特徴量の扱い方について学んだことを本稿にまとめることを目的とします。 目的変数と説明変数をひとつひとつ確認することで考察していきます。

今回は、kaggleのコンペティションに使用されたタイタニック号の生存者予測のデータを用います。

www.kaggle.com

提供されるデータにはtrain dataおよびtest dataの2つがあります。train dataは予測モデルを構築するためのデータ、test dataには予測するためのデータ(生存情報なし)が格納されています。 今回のタイタニックコンペティションでは、生存した乗客を予測することが目的となります。

まず、データの中身を確認します。train dataを表に示します。横軸に乗客の情報(特徴量)が並んでいます。縦軸が乗客数であり、891人がであることがわかります。

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
train = pd.read_csv("/Users/***/Desktop/kaggle/train.csv")
train

f:id:datascience_findings:20210221105947p:plain
Train Data of Tytanic

それぞれの特徴量は、 PassengerId:乗客に付与された番号。kaggle側にて付与したものだと考えられますが、現段階では、生存率との相関関係は不明です。
Survived:生存結果であり、1が生存者、0が死亡者を表します。この特徴量を予測することが今回の目的です。
Pclass:乗客のクラス(ランク)を表します。1,2,3という番号がありそうです。
Name:乗客の名前を表します。
Sex:乗客の性別を表します。
Age:乗客の年齢を表します。
SibSp:乗客の兄弟、配偶者の数を表します。
Parch:乗客の両親、子どもの数を表します。
Ticket:乗客のチケット番号を表します。
Fare:乗客の運賃を表します。
Cabin:乗客の部屋番号を表します。
Embarked:乗客が乗船した出港地を表します。

現段階において、どの特徴量が生存率に影響するかは不明です。本稿では、これら特徴量と生存率の関係を考察していきます。

SurvivedとPclassの関係

Pclassは乗客のランクを表す特徴量です。ランクが高い乗客であるほど、優先的に救助された可能性が考えられるため、生存率との相関が強いと考えられます。図1にPclassの分布を示します。

sns.catplot(data=train , x='Pclass', kind = 'count', ax = 'White', palette = 'Blues')

f:id:datascience_findings:20210221154450p:plain
図1 Pclassの分布
図1から、Pclass=3の人数が乗客の半数以上を占めることがわかります。Pclass=1, 2(上流階級)は人数が少ないことがわかりました。 次に、図2にPclassごとの生存状況を表2にPclassごとの生存率を示します。(Survived : 0 = 死亡者、1 = 生存者)

sns.catplot(data=train , x='Pclass', hue = 'Survived', kind = 'count', palette = 'gist_heat')
plt.show()
rate = train.groupby('Pclass')['Survived'].apply(lambda d: d.value_counts(normalize=True)*100)
rate = pd.DataFrame(rate.reset_index())
rate = rate.loc[rate['level_1'] == 1]
rate.drop('level_1', axis =1).reset_index(drop = True).round(2)
plt.show()

f:id:datascience_findings:20210221161401p:plain
Pclassごとの生存者数と死亡者数
f:id:datascience_findings:20210221165348p:plain
表2 Pclassと生存率の関係
図2と表2から、Pclassが小さくなる(階級が高い)ほど生存率が増加することがわかります。よって、Pclassは生存率に寄与する特徴量であると判断できます。

SurvivedとSexの関係

 性別の事前に想定される影響は、優先的に女性が救助される、もしくは、体力に勝る男性の生存率が高いが考えられます。図3に男性と女性の総数、図4に性別ごとの生存状況を示します。表3に性別ごとの生存率を示します。

sns.catplot(data=train , x='Sex', kind = 'count', ax = 'White', palette = 'Blues')
plt.show()
sns.catplot(data=train , x='Sex', hue = 'Survived', kind = 'count', palette = 'BuPu')
plt.show()

f:id:datascience_findings:20210221171035p:plainf:id:datascience_findings:20210221170958p:plain
図3 性別ごとの総数 図4 性別ごとの生存者数と死亡者数
f:id:datascience_findings:20210221173401p:plain
表3 性別と生存率の関係
図3から女性は男性と比較し乗員数は少ないことがわかります。また、図4と表3から、女性の生存率は男性と比較し非常に高いことから、性別もまた生存率を予測する特徴量であると判断できます。

SurvivedとEmvarkedの関係

 Emvarked(出港地)が生存率にどのように影響するか想定することは難しいです。このような相関を持つか判断し難い特徴量の扱いには慎重になるべきだと考えられます。  まず、unique関数を使用し、出港地を確認します。

train['Embarked'].unique()

f:id:datascience_findings:20210221175043p:plain ここで、SはSouthampton(イギリス)、CはCherbourg(フランス)、QはQueenstown(アイルランド)を表します。また、nanは欠損値であり出港地が不明な乗客の存在を表します。 図5に出港地ごとの生存状況、表4に出港地ごとの生存率を示します。

sns.catplot(data=train , x='Embarked', hue = 'Survived', kind = 'count', palette = 'gist_heat')
f:id:datascience_findings:20210221180552p:plainf:id:datascience_findings:20210221180557p:plain

図からC地点の乗客は生存率が50%を超える一方、他2つの地域は約35%の生存率と大きく異なることがわかりました。しかし、タイタニック号の航路についてwikipediaを確認した結果、3場所ともにヨーロッパと土地的差は見られませんでした。

f:id:datascience_findings:20210613222632p:plain
タイタニック号航路
Sの生存率が最も低い考えられる一つの要因としては、出発地であることから乗客だけでなく船員が多く乗っていた可能性はありますが、Cが特別生存率が高い要因はわかりませんでした。  なぜ、Cの乗客の生存率が高いか推定できない場合には、他の特徴量が関係している(疑似相関)の可能性が考えられます。 ここまでの2つの特徴量(PclassとSex)とEmbarkedに疑似相関があるか確認を行った結果を図6、7に示します。

sns.barplot(data=train , x='Sex', y='Survived', hue = 'Embarked', palette = 'Dark2')
sns.barplot(data=train , x='Pclass', y='Survived', hue = 'Embarked', palette = 'Dark2')
f:id:datascience_findings:20210221194402p:plainf:id:datascience_findings:20210221194448p:plain

図6、7では、barpltを用いて縦軸に生存率を示しています。グラフに表示されるエラーバーはグラフ頂点から95%の可能性でデータが含まれる範囲を表します。性別ごとに出港地の生存状況を比較した結果、性別に因らずC地点の生存率が高いことがわかりました。つまり、出港地と性別に疑似相関は認められません。  次に、階級ごとに出港地との生存率を比較しました。その結果、Q地点のPclass = 1,2において、エラーバーが0から1となり信頼できる生存率が取得できませんでした。この原因として、n数(Q地点でのPclass = 1,2の乗客数)が少ないことが考えられます。その総数をカウントした結果を下図に示します。

sns.catplot(data=test , x='Embarked', hue = 'Pclass', kind = 'count', palette = 'gist_earth')
plt.title('Passengers of Embarked by Pclass', fontsize=20)

f:id:datascience_findings:20210221195951p:plain
このグラフから、事前の推察の通り、n数が小さいことが確認できました。 この結果も含め、考察を行うと、Pclass = 1,2では、エラーバーが出港地に因らず重なることから、出港地と生存率に相関関係はみられません。しかし、Pclass = 3においては、S地点の乗客は生存率が他の出港地と比較し明らかに低いことがわかりました。

SurvivedとCabinの関係

 train dataの概要から客室の情報に欠損値が確認できます。また、アルファベットと数字の組み合わせで表されることがわかります。Embarkedの調査と同様、unique関数を使用し、特徴量にどのうような値が格納されているか確認します。

cabin = train['Cabin'].unique()
cabin = cabin.astype('str')
print(np.sort(cabin))

print(train['Cabin'].isnull().sum())
f:id:datascience_findings:20210222184338p:plainf:id:datascience_findings:20210222185017p:plain

この結果から、先頭のアルファベットはAからGまでであることがわかります。続く数字は1桁から3桁の数字であることが確認できました。また、欠損値は687個(77%)と非常に多いことがわかります。
先頭のアルファベットが生存率に影響を下図、表に示します。

train['Cabin'] = train['Cabin'].str[0]
sns.catplot(data=train , x='Cabin', hue = 'Survived', kind = 'count')

rate = train.groupby('Cabin')['Survived'].apply(lambda d: d.value_counts(normalize=True)*100)
rate = pd.DataFrame(rate.reset_index())
rate = rate.loc[rate['level_1'] == 1]
rate.drop('level_1', axis =1).reset_index(drop = True).round(2)
f:id:datascience_findings:20210222204418p:plainf:id:datascience_findings:20210222204423p:plain

これら図表から、客室の存在により、生存率が向上していることがわかります。また、B、D、Eから始まる客室は生存率が75%程度と非常に高いことがわかります。また、客室の中で、Aの生存率は50%以下でした。タイタニック号の客室は下図のように上層から下層に向け、AからGに推移することがネット情報から確認できています。救命ボートは、デッキに搭載されていることから、Aの客室は他の客室と比較し、避難が容易であったと考えられますが、実際の結果は異なります。
 これは、欠損値が多いためであると考えられます。欠損値は687個(全データは891個)あり、大半のデータが反映されておりません。これは客室のユニークな部屋数は147室であり、各アルファベットの層の客室番号に欠番が見られることからも推測できます。客室は生存率に影響し得る特徴量であるため、欠損値を補完し予測に用いることが望ましいですが、今回のように欠損値が大半を占める場合は不適切な相関を得る可能性があるため、特徴量として相応しくないと判断できます。

SurvivedとAgeの関係

 年齢の分布と欠損値の数を確認します。

sns.catplot(data=train, x='Survived', y='Age', kind = 'swarm', palette = 'gist_heat')
plt.title('Survived Rate by Age', fontsize=20)

train['Age'].isnull().sum()
out : 177

f:id:datascience_findings:20210224093334p:plain
図から、死亡者と生存者で分布の形は大きく異ならないようにみられることから、年齢と生存率に強い相関はないと考えれます。また、欠損値は177個と少なくありませんので、扱いに注意が必要となります。 実際に年代ごとの生存状況のグラフと生存率の表を確認します。

train['AgeBand'] = pd.cut(train['Age'], 8, labels = False)
train['AgeBand'] = train['AgeBand'] * 10
sns.catplot(data=train , x='AgeBand', hue = 'Survived', kind = 'count', palette = 'gist_heat')
plt.title('Survived Rate by AgeBand', fontsize=20)

rate = train.groupby('AgeBand')['Survived'].apply(lambda d: d.value_counts(normalize=True)*100)
rate = pd.DataFrame(rate.reset_index())
rate = rate.loc[rate['level_1'] == 1]
rate.drop('level_1', axis =1).reset_index(drop = True).round(2)
f:id:datascience_findings:20210224100959p:plainf:id:datascience_findings:20210224101141p:plain

表から、0年代(0~9歳)では、生存率が高くなっており、優先的に救助された可能性が考えられます。反対に高齢層(60~80歳)は生存率が低く、極寒環境で体力が持たなかった可能性が考えられます。その他の年代では、生存率は約40%程度であり、891人全ての乗客の生存率と大きな相違はありません。

SurvivedとSibSpおよびParchの関係

 ここでは、同乗した兄弟・配偶者の人数と親・子の人数をまとめて家族サイズとして評価します。この2つの特徴量は、兄弟か配偶者、親か子という重要とも考えられる変数が混在してしまっており、混在する情報をなくすため、家族サイズ(SibSp + Parch + 1)を新しいパラメータとして定義します。ここで+1は、自身の数です。

train['Familysize'] = train['Parch'] + train['SibSp'] + 1
sns.catplot(data=train , x='Familysize', hue = 'Survived', kind = 'count', palette = 'gist_heat')
plt.title('Survived Rate by Familysize', fontsize=20)

rate = train.groupby('Familysize')['Survived'].apply(lambda d: d.value_counts(normalize=True)*100)
rate = pd.DataFrame(rate.reset_index())
rate = rate.loc[rate['level_1'] == 1]
rate.drop('level_1', axis =1).reset_index(drop = True).round(2)
f:id:datascience_findings:20210224113843p:plainf:id:datascience_findings:20210224113850p:plain

これら図表から、家族サイズが2,3,4人においては、生存率が高くなっており、その他家族規模では生存率は38%より低くなっています。しかし、家族規模4人以上では、人数が少ない(全数の5%以下)ため、n数として不足していると考えられます。また家族規模が11人であれば、11人の倍数のPassengerIdが得られるはずですが、実際には下の結果が得られました。

print('家族規模が4人である乗客の人数:{}人'.format((train['Familysize'] == 4).sum()))
print('家族規模が5人である乗客の人数:{}人'.format((train['Familysize'] == 5).sum()))
print('家族規模が6人である乗客の人数:{}人'.format((train['Familysize'] == 6).sum()))
print('家族規模が7人である乗客の人数:{}人'.format((train['Familysize'] == 7).sum()))
print('家族規模が8人である乗客の人数:{}人'.format((train['Familysize'] == 8).sum()))
print('家族規模が11人である乗客の人数:{}人'.format((train['Familysize'] == 11).sum()))

f:id:datascience_findings:20210224115249p:plain

この結果から、これらクラスに最低13個の欠損値が含まれます。仮に家族規模が11人である欠損値の最低4人が全員生きていると仮定すると生存率は50%を超えます。そのため、現段階では、家族規模により大きく生存率が変化するとは断定できません。

SurvivedとFareの関係

 運賃ごとの生存状況を下図に示します。

sns.catplot(data=train, x='Survived', y='Fare', kind = 'swarm', palette = 'Dark2')
plt.title('Survived Rate by Fare', fontsize=20)

f:id:datascience_findings:20210224133750p:plain
上図より、100$以上の乗客では、生存率が高そうであることが伺えます。しかし、大半の乗客が100$以下に密集しているいることから、この区間を詳しく見る必要があります。
下記のように区分わけをします。

train['FareBand1'] = pd.cut(train['Fare'], [0, 25, 50, 75, 100, 1000])
sns.catplot(data=train , x='FareBand1', kind = 'count', palette = 'Blues_r').set_xticklabels(rotation=90)
plt.title('Survived Rate by FareBand1 ', fontsize=20)

train['FareBand2'] = pd.cut(train['Fare'], [0, 7, 7.5, 7.8, 8, 10, 15, 25, 35, 60, 100, 1000])
sns.catplot(data=train , x='FareBand2', kind = 'count', palette = 'Blues_r').set_xticklabels(rotation=90)
plt.title('Survived Rate by FareBand2 ', fontsize=20)

fare = train['Fare'].unique()
np.sort(fare)
f:id:datascience_findings:20210224140004p:plainf:id:datascience_findings:20210224140009p:plain

f:id:datascience_findings:20210224140747p:plain Fareband1では0〜100$の区間を4分割しましたが、25$以下の乗客数が多数でした。そこでFareband2のグラフに示すようにさらに細分化しました。このグラフから、運賃には特定の範囲(7〜10$)にの分布が高くなっています。 人数を大まかに分けられたFareBand2について、生存状況を調べると下の図表になりました。

train['FareBand2'] = pd.cut(train['Fare'], [0, 7, 7.5, 7.8, 8, 10, 15, 25, 35, 60, 100, 1000])
sns.catplot(data=train , x='FareBand2', kind = 'count', hue = 'Survived', palette = 'gist_heat').set_xticklabels(rotation=90)
plt.title('Survived Rate by FareBand2 ', fontsize=20)

rate = train.groupby('FareBand2')['Survived'].apply(lambda d: d.value_counts(normalize=True)*100)
rate = pd.DataFrame(rate.reset_index())
rate = rate.loc[rate['level_1'] == 1]
rate.drop('level_1', axis =1).reset_index(drop = True).round(2)
f:id:datascience_findings:20210224163235p:plainf:id:datascience_findings:20210224163228p:plain

この図表から、10$以下の乗客は生存率が38%を大きく下回っており、10$以上で運賃に比例して生存率が増加する傾向が見られました。よって、運賃は生存率の特徴量として妥当であると考えられます。

SurvivedとTicketの関係

 ユニーク関数でチケットの中身を確認しますが、数が多いため、ここでは、200個を表示します。

train['Ticket'].unique()[:200]

f:id:datascience_findings:20210224165649p:plain
結果から、チケットの表記は数字の列および先頭にアルファベットが付与される場合がありそうです。そのため、アルファベットと数字の列に分け考えます。
 まず、数字列から考察を行います。

Ticket_num = []
for idx, row in train.iterrows():
    Ticket_num.append(row['Ticket'].split(' ')[-1])
train['Ticket_num'] = Ticket_num
train['Ticket_num'] = train['Ticket_num'].replace('LINE', 0)
train['Ticket_num1'] = train['Ticket_num'].astype(int)
sns.catplot(data=train, x='Survived', y='Ticket_num1', kind = 'swarm', palette = 'gist_heat') 

sns.catplot(data=train, x='Survived', y='Ticket_num1', kind = 'swarm', palette = 'gist_heat') 
plt.title('Survived Rate by Ticket number1', fontsize=20)
plt.yscale('log')
plt.ylim(1,1E7)
f:id:datascience_findings:20210224220718p:plainf:id:datascience_findings:20210224221119p:plain

左のグラフでは、数字列が桁数により大きく異るため、分布の形がわかりません。そのため、右のグラフでは縦軸を対数表示とすることで、桁数の変化に対しても分布を観察することができています。このグラフから、数字列は主に四桁以上であることが伺えます。また、10の5乗以上の区間では、数列のバンドが形成されており、先頭の数字のみで分けられる可能性が考えられます。

train['Ticket_num'] = pd.cut(train['Ticket_num1'], [0, 3000, 10000, 20000, 300000, 4000000])
sns.catplot(data=train , x='Ticket_num', hue = 'Survived', kind = 'count', palette = 'gist_heat').set_xticklabels(rotation=60)
plt.title('Survived Rate by NumberBand', fontsize=20)

rate = train.groupby('Ticket_num')['Survived'].apply(lambda d: d.value_counts(normalize=True)*100)
rate = pd.DataFrame(rate.reset_index())
rate = rate.loc[rate['level_1'] == 1]
rate.drop('level_1', axis =1).reset_index(drop = True).round(2)
f:id:datascience_findings:20210227190927p:plainf:id:datascience_findings:20210227190931p:plain

この図表から、先頭が3であるチケットの生存率は基準生存率を下回ることが確認できます。また、1万台のチケットは生存率が高くなっていました。

2つ前のグラフ中ののチケット番号が100以下(外れ値と思われる)となる乗客の情報を確認します。

train[train['Ticket_num1'] < 100]

f:id:datascience_findings:20210224222406p:plain
このデータフレームから100以下にプロットされたチケットは「S.O./P.P. 3」特殊なチケットでした。また、「LINE」という数字列があり、このチケットはFare = 0であることがわかりました。ここから、乗客ではなく、船員という新たな特徴量を考える必要性も考えられます。
そこで、運賃が0$である乗客(船員)を確認します。

train[train['Fare'] == 0]

f:id:datascience_findings:20210224223935p:plain
船員であれば、出発地(S)から乗船すること、同乗した家族はいないことが推測できますが、全員この条件を満たしていました。この条件の場合、船員の生存者は一人のみであり、有意な特徴量であると考えられますが、これだけでは船員として少ないため他の条件も考慮する必要がありそうです。

続いて、チケットの先頭に付くアルファベットを考察します。

Ticket_str = []
for idx, row in train.iterrows():
    Ticket_str.append(row['Ticket'].split(' ')[0])
train['Ticket_str'] = Ticket_str
train['Ticket_str'] = train['Ticket_str'].where(train['Ticket_str'] != train['Ticket_num'])
train['Ticket_str'].unique()

sns.catplot(data=train , x='Ticket_str', hue = 'Survived', kind = 'count', size=4, aspect=3, palette = 'gist_heat').set_xticklabels(rotation=90)
plt.title('Survived Rate by Ticket strings', fontsize=20)
f:id:datascience_findings:20210225215158p:plainf:id:datascience_findings:20210225215203p:plain

 ユニーク関数により種類を確認した結果、これら文字列に傾向は見られませんでしたが、STONやCAなど含まれるものが複数あることから、グループ分けはできそうです。また、「SC/Paris」から地名(出港地)の意味を持つ可能性が考えられます。
 グラフの結果を見ると、各アルファベット列の頻度は非常に少ないことがわかります。頻度が最も多く、生存率が高い「PC」を持つデータを分析します。

trainPC = train[train['Ticket'].str.contains('PC')]
trainPC.describe()

f:id:datascience_findings:20210225220406p:plain
PCを含む場合、全てPclass = 1であり、運賃の平均値は全乗客の平均の約4倍でした。このことから、PCと生存率は疑似相関であることが考えられます。その他アルファベット列は、n数が少ないこと、アルファベット列の有無により一概に生存率が向上する傾向が見られないことから特徴量としては不十分であると考えます。

SurvivedとNameの関係

 外国人の名前は大きくファミリーネーム、ミドルネーム、ファーストネームに分けられます。ファミリーネームを分析することで家族を特定できそうです。また、trainデータを見ると、ミドルネームにMr.やMrs.があり、一つの特徴量となりそうです。今回のデータセットにおいては、カンマ表記がされているため、先頭がファミリーネームとなります。
まず、ファミリーネームから確認します。

train['family_name'] = train['Name'].str.split(',', expand = True)[0]
len(train['family_name'].unique())
out : 667

ユニークなファミリーネームの数は667個あり、ファミリーネームから生存率の考察を行うことは難しいと判断します。
続いて、ミドルネーム(敬称)を考察します。

train['honorific'] = train['Name'].str.split(',', expand = True)[1].str.split('.', expand = True)[0]
sns.catplot(data=train , x='honorific', hue = 'Survived', kind = 'count', size=4, aspect=2, palette = 'gist_heat').set_xticklabels(rotation=90)
plt.title('Survived Rate by Honorific', fontsize=20)

f:id:datascience_findings:20210225224309p:plainf:id:datascience_findings:20210225224539p:plain

この結果から、Mrの生存率16%は、男性のみの生存率19%と大差なく、Mrs, Missの生存率(79%、70%)もまた女性のみの生存率74%と大差がないことを確認しました。
これらを除いた乗客の情報を確認します。

train.query('honorific not in ("Mr", "Mrs", "Miss")').describe()

この結果から、全乗客とMr,Mrs,Missを除いた乗客の生存率は後者が高くなりました。2つのデータ間での大きな違いは年齢と運賃です。2つのデータ間に平均値の差は見られませんが、中央値に違いが見られます。年齢は後者データでは9歳であり、幼子が多いことが確認できます。これは、「Age」の特徴量を扱った際の傾向(10歳以下では生存率が増加)と一致します。また、運賃の中央値は後者データでは29$であり、全乗客の運賃の中央値の約2倍でした。これもまた、「Fare」の特徴量を扱った際に見られた傾向(運賃が高いほど生存率が高い)と一致しました。
 後者データの「Age」には欠損値が5つあります。欠損値補完の際にこの結果は有用であると考えられます。

 まとめ

 各特徴量について、それぞれ生存率を確認し、その結果から、下記の傾向が得られました。 1.全乗客の生存率  891名の乗客の生存率は約38%であり、この値を基準とし、特徴量の相関を確認します。 2.SurvivedとPclassの関係  Pclassの値が小さくなる(階級が上がる)ほど生存率が増加することを確認しました。Pclass = 1,2では基準生存率以上であり、Pclass = 3は基準値以下であり、生存率の説明変数として妥当だと考えます。
3.SurvivedとSexの関係
 女性の生存率は約74%と、基準値を大幅に上回りました。反対に、男性の生存率は約19%であり、基準値を下回りました。この結果から、性別は生存率の説明変数として妥当だと考えます。
4.SurvivedとEmvarkedの関係
 出港地がC(フランス)の場合、生存率は基準値を上回りました。Q(アイルランド)の生存率はほぼ基準値と一致し、S(イギリス:スタート地点)の生存率は基準値を下回りました。結果から、出港地は生存率の説明変数として妥当だと考えられますが、その根拠を実現象から想定できません。よって、これは疑似相関の可能性が考えられます。また、出港地は欠損値を2つ含んでおり、補完する必要があります。
5.SurvivedとCabinの関係
 客室の特徴量は先頭のアルファベットと続く数字列の組み合わせであり、先頭のアルファベットは船体のフロアを表していることをタイタニック号の図面から確認しました。図面から、救命ボートが近い上層フロア(A>F)に近づくほど、生存率が増加すると考えましたが、その傾向は見られませんでした。理由の一つとして、Cabinは欠損値を687個(約78%)含んでおり、正確な傾向を反映していないためだと考えます。また、傾向が見られないことから欠損値補完も不能であり、客室は生存率の説明変数として妥当ではないと判断します。
6.SurvivedとAgeの関係
 年代ごとに分けた結果、10歳以下は基準生存率を上回り、60歳以上の生存率は基準値を下回りました。10〜60歳では、生存率はほぼ基準値と一致しました。しかし、欠損値は177個あり、適切に補完した上で、生存率との相関を議論する必要があります。
7.SurvivedとSibSpおよびParchの関係
 SibSpとParchは同乗した人数を表すことから一つの特徴量(Familysize = SibSp + Parch + 1)に変換しました。この結果から、家族規模が2,3,4人の場合、生存率が基準値を上回りました。単身または5人以上の家族の場合、生存率を下回る傾向が見られました。この結果から、家族規模と生存率に相関が認められそうですが、出港地の場合と同様直接要因を想定できません。よって、家族規模については、より深堀り(家族単位で生存するか家族の特定人物(子どもや女性)が生存するか)検討する必要があると判断します。
8.SurvivedとFareの関係
 運賃が10$以下の乗客が約37%を占め、生存率は基準値を下回りました。10$以上の乗客は、運賃が上昇するにつれ、生存率が増加する傾向が確認できました。よって、運賃は生存率の説明変数として妥当であると考えます。
 また、運賃が0である乗客を確認したところ、 9.SurvivedとTicketの関係
 チケットは数字列と先頭の文字列(891枚の内230枚)で表されます。数字列は桁数および最初の数字で区分わけできそうであることを確認しました。数字列が2桁以下の乗客を確認した結果、運賃が0である乗客のがいることから、船員という新たな特徴量を検討する必要性が示唆されました。
 文字列がPCである場合、生存率は基準値を上回りました。これら乗客の情報を確認すると、Pclassが全員1であり、運賃が高いという傾向を確認しました。しかし、文字列ごとのn数は少なく、文字列の持つ意味が不明であり、文字列ごとの生存率に見られる傾向は異なることから説明変数として妥当でないと判断します。
10.SurvivedとNameの関係
 名前のファミリーネームと敬称に着目しました。ファミリーネームは667個の種類があり、この変数のみで考察は不能だと判断しました。敬称は、Mr、Mrs、Missは頻度が多いものの、それぞれ男性、女性の生存率からの変化は小さく特徴量として不適切であると考えられます。これらを除いた場合、年齢の中央値は9と大幅に小さくなりました。年齢の欠損値補完の参考にできそうです。

課題

 今回、特徴量ごとに生存率との関係を評価しました。しかし実際には、第三の特徴量を介して擬似的に相関がある場合も考えられます。それらを考慮し、適切な変数に変換することを次回の課題とします。
 また、予測を行うためには欠損値を補完する必要があります。最も簡単な補完処理は全数の平均値を代入することですが、欠損値が多いほど単純な補完では、実現象からずれが生じてしまいます。特に、今回のデータのように多数特徴量がある場合には、それら特徴量から補完する条件を絞り込むことも可能です。次回以降では、特徴量間の相関から欠損値を補完する方法について検討も行います。