タイタニック号生存者予測 - 主成分分析(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は全主成分ベクトルで一定値存在しているがこれは一様分布に近い分布であるためだと考えられる。つまり、データに特徴がないと言える)など一致した結果であると言えます。
このように主成分分析は特徴量選定に有効です。今回のように、特徴量が少ない場合は次元圧縮の効果も少ないですが、計算コストの削減にも有効です。