「入門はじめての多変量解析」をPythonで写経 Vol.15 ～ ３章「はじめての主成分分析」⑦主成分分析 FOREVER

３章「はじめての主成分分析」

書籍の著者　石村貞夫 先生、石村光資郎 先生

書籍「入門はじめての多変量解析」３章「はじめての主成分分析」の Python写経活動記録 です。 

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。

この記事は、主成分分析の最終話 です。
名残惜しいのでこの書籍からはみ出して、別の書籍の分析事例を参考にして、主成分分析を堪能したいと思います！

ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

今回記事で取り組む書籍

統計検定準１級公式テキストの「ロンドンオリンピック女子７種競技」の主成分分析に取り組みます！

【引用表記】
この記事は統計検定準１級の公式テキスト「統計学実践ワークブック」（以下、「ワークブック」と呼びます）の例題 問 22.1 を引用いたします。
記事に用いるデータは、ワークブックが引用する Wikipedia に掲載されたデータを引用しています。

引用書籍：
日本統計学会公式認定 統計検定準１級対応 統計学実践ワークブック
第１版 第６刷
編者 一般社団法人 日本統計学会
発行所 学術図書出版社

引用データ：
WIKIPEDIA - The Free Eccyclopedia
Athletics at the 2012 Summer Olympics – Women's heptathlon
URL:

主成分分析は続くよ

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# PCA
from sklearn.decomposition import PCA

# 統計
import statsmodels.api as sm

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

分析の概要

今回はいったんテキスト「入門はじめての多変量解析」から離れます。
統計検定準１級のワークブック p.199～ 「2012年ロンドンオリンピック、女子７種競技の得点」の分析例を参考にして、主成分分析を進めます。

７種の各競技の記録はポイント化され、ポイント合計の高い順に順位が決まります。
７つの競技ポイントを変数にした主成分分析で「競技の総合指標」を見つけに行きましょう！

データの概要

■ データの読み込み
Wikipedia の情報を csv ファイル化して読み込みます。
※手作業で csv ファイルを作成しました。

各競技の最高ポイントをハイライト表示します。

### データの読み込み
df = pd.read_csv('./data/2012london.csv')

# データの表示：各競技の最高点をハイライト
df.style.highlight_max(subset=df.columns[4:], color='lightpink')

【実行結果】
列名に７種目が表示されています。
短・中距離走、跳躍、投てきの競技です。
概ね上位者が最高ポイントを獲得している感じです。

■ 可視化で概観を確認
ひとまずポイントの積み上げ棒グラフを描きましょう。
pandas の plot メソッドで描画します。

# 積み上げ棒グラフで可視化

# 全体を（）で括ると途中で改行を入れやすくなる
(df
  .set_index('順位')       # 保険会社列をインデックスに
  .iloc[:, 3:]             # ポイントの列を選択
  .iloc[::-1]              # 行を逆順に
  .plot.barh(              # 横棒グラフの描画
    stacked=True,               # 積み上げ棒グラフ
    alpha=0.7,                  # 透過度
    figsize=(8, 8)              # グラフのサイズ
  )
)

# 凡例をグラフ外に
plt.xlabel('ポイント')
plt.legend(bbox_to_anchor=(1, 1))
plt.grid(axis='x');

【実行結果】
1000ポイントごとの垂直線を補助線にしてみると、100mハードルと走り高跳びの合計は「どんぐりの背比べ」ですが、、砲丸投げ以右でジリジリと差が広がっている感触です。

続いてヒストグラムと散布図行列です。

# ヒストグラムと散布図の描画
sns.pairplot(data=df.iloc[:, list(range(4, 11)) + [0]], height=1.2);

【実行結果】
最下行の順位と各競技ポイントの散布図からは、おおむね、ポイントが大きくなると順位が上がる（縦軸で下に向かう）関係です（当たり前？）
競技間では100mハードルと200m走に正の相関が見られます。

相関行列を確認しましょう。

# 相関係数の表示
df.iloc[:, [0] + list(range(4, 11))].corr().round(3)

【実行結果】

順位と相関が高いのは 走り幅跳び、走り高跳び、100mハードルです。
競技間で相関が高いのは、100mハードルと200m走の 0.626 です。

データの前処理

データを標準化します。
また、変数名等を用意しておきます。

### データの前処理

# 変数の標準化
X = df.iloc[:, 4:]
X_std = (X - X.mean()) / X.std(ddof=1)
display(X_std.head())

# 変数名と主成分名のリストの作成
var_names = X.columns
pc_names = [f'PC{i+1}' for i in range(len(var_names))]

【実行結果】
標準化後のデータです。
この記事の主成分分析は標準化後のデータを利用します。

ちょっと寄り道～重回帰分析・決定木分析

主成分分析に突入する前に…
順位に効いてそうな競技を見つけに行きましょう！

📈 重回帰分析・標準化偏相関係数
順位を目的変数、各競技ポイントを説明変数にして重回帰分析を行います。
目的変数・説明変数の両方にデータの標準化を実施して、標準化偏相関係数を取得します。

### 標準化偏回帰係数の確認：目的変数：順位、説明変数：各球技のポイント

# 目的変数の設定、標準化
y = df['順位']
y_std = (y - y.mean()) / y.std(ddof=1)

# 説明変数に定数項を追加
X_std_const = sm.add_constant(X_std)

# 線形回帰の実行
result = sm.OLS(endog=y, exog=X_std_const).fit()

# 結果の表示
result.summary()

【実行結果】
決定係数の高すぎさ、ダービン・ワトソン比（Durbin-Watson）が示すやや残差の自己相関、があるものの、構わず、偏相関係数を見ます。

100mハードル以下の偏相関係数 coef はデータを標準化しているので変数の単位に影響を受けない「標準化偏相関係数」に変身しています。
順位への影響度は絶対値が大きさで見ます。

📈 決定木分析
順位を目的変数、各競技ポイントを説明変数にして決定木分析（回帰木）を行います。
樹形図で順位の決定経路を可視化しましょう。

scikit-learn の DecisionTreeRegressor() でモデルを構築します。

# 回帰木モデルの構築
model = DecisionTreeRegressor()
model.fit(X=X, y=y)

【実行結果】なし

樹形図を描画します。
まずは scikit-learn の標準機能 plot_tree を利用します。

# treeの可視化１
plt.figure(figsize=(10, 7))
plot_tree(model, max_depth=2, feature_names=model.feature_names_in_,
          filled=True, fontsize=10);

【実行結果】
最上位からの一部を切り取っています。

【見方】
一番上の箱は「走り幅跳びが 886.5 ポイント以下」を条件にして、左が Yes、右が No で分岐します。

途中で切れないように、graphviz ライブラリで描画します。

# treeの可視化２
dot_data = export_graphviz(
    model, out_file=None, feature_names=model.feature_names_in_, filled=True)
g = graphviz.Source(dot_data, format='png') 
g

【実行結果】
部分を切り取って掲載しています。
最初に、上位者を分岐した右側を。

右側の末端（最下）に注目します。
value の値が順位です。右側に１～３位が集まっています。

シンプルですね！

続いて下位が集中する左側を見ましょう。
画像をアップにして見て下さい。
27 位、29～30 位が最左に位置します。

決定木系のアルゴリズムはモデル構築時に説明変数（特徴量）の重要度を計算しています。
可視化しましょう。

# 特徴量重要度の可視化

# 並び順の取得
sort_index = model.feature_importances_.argsort()
# 棒グラフの描画
plt.barh(
    model.feature_names_in_[sort_index], model.feature_importances_[sort_index])
# 修飾
plt.xlabel('特徴量重要度', fontsize=12)
plt.grid(axis='x', alpha=0.5)

【実行結果】

条件分岐＝順位決定にとって、走り幅跳びと800m走の重要度が大きいです。

前座で会場が温まったところで真打ちの登場です！

主成分分析

■ 主成分分析の実行
標準化データで、scikit-learn で、サクッと主成分分析を実行します。

### 主成分分析の実行 scikit-learn利用

# PCAインスタンスの生成
pca = PCA()
# PCAを実行して主成分得点を取得
scores = pca.fit_transform(X_std)

【実行結果】なし

■ 固有値・寄与率の確認
固有値と寄与率を確認して、分析に利用する主成分を考えましょう。

# 固有値と寄与率の表示 p.200 表22.4

# データフレームの作成
eigen_vals = pd.DataFrame(
    {'固有値': pca.explained_variance_,
     '寄与率': pca.explained_variance_ratio_,
     '累積寄与率': pca.explained_variance_ratio_.cumsum()},
     index=pc_names)

# 結果の表示
eigen_vals.style.format(
    {'固有値': '{:.2f}', '寄与率': '{:.1%}', '累積寄与率': '{:.1%}'})

【実行結果】
累積寄与率 0.7 超え & 固有値の意味ある大きさ 1 以上、とすると第３主成分まで分析するのが良さげです。

寄与率を可視化しましょう。

# 寄与率の可視化

# 寄与率の棒グラフの描画
plt.bar(eigen_vals.index, eigen_vals['寄与率'], label='寄与率')
# 累積寄与率の折れ線グラフの描画
plt.plot(eigen_vals.index, eigen_vals['累積寄与率'], '-o', color='tab:red',
         alpha=0.7, label='累積寄与率')
# 累積寄与率のしきい値の水平点線の描画
plt.axhline(0.7, color='black', linestyle='--', linewidth=0.5,
            label='しきい値（目安）')
# 修飾
plt.xlabel('主成分', fontsize=12)
plt.ylabel('寄与率, 累積寄与率', fontsize=12)
plt.legend();

【実行結果】
第３主成分まで見ていきます！

■ 固有ベクトルの確認
固有ベクトルをサマリーします。

# 固有ベクトル（第1, 2主成分）の表示 p.200 表22.5

# データフレームの作成
eigen_vecs = pd.DataFrame(
    pca.components_.T[:, :3], index=var_names, columns=pc_names[:3])
# 結果の表示
eigen_vecs.round(3)

【実行結果】
第１主成分の正負がワークブックの結果と逆転していますが、本質的には問題ないのでご安心を。

■ 主成分負荷量の確認

# 主成分負荷量（第1, 2主成分）の表示 p.200 表22.6

# 主成分負荷量の算出
# 相関行列ベースの場合、√固有値i × 固有ベクトルij, i:主成分, j:変数
loadings = eigen_vals.loc[['PC1', 'PC2', 'PC3'], '固有値']**(1/2) * eigen_vecs
# 結果の表示
loadings.round(2)

【実行結果】
固有ベクトルで見た考察と同様の傾向があります。

■ バイプロットで分析
第１主成分と第２主成分のバイプロットを可視化します。

# バイプロットの描画

## 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 7), tight_layout=True)

## 30選手の主成分得点の描画
# 主成分得点の散布図の描画（ダミー）
ax.scatter(scores_df['PC1'], scores_df['PC2'], marker='')
# 主成分得点の「順位」のテキストの表示
for x, y, s in zip(scores_df['PC1'], scores_df['PC2'], scores_df.index):
    ax.text(x=x, y=y, s=s, ha='center', va='center', color='royalblue',
            fontsize=12)
# x=0の垂直線、y=0の水平線の描画
ax.axvline(0, color='black', lw=0.5)
ax.axhline(0, color='black', lw=0.5)
# 修飾
ax.set_xlabel('第１主成分', fontsize=12)
ax.set_ylabel('第２主成分', fontsize=12)
ax.set(xlim=(-5, 5), ylim=(-5, 5))

## 7競技の主成分負荷量の描画
# 描画領域の設定：主成分得点と異なるx軸・y軸のスケールで描画するため
twin_ax = ax.twinx().twiny()
# 主成分負荷量の矢印と競技名の描画
for x, y, text in zip(loadings['PC1'], loadings['PC2'], loadings.index):
    # annotateで矢印をプロット
    twin_ax.annotate(
        text='',          # テキスト：なし
        xy=[x, y],        # テキストの位置:[第1主成分負荷量, 第2主成分負荷量]
        xytext=[0, 0],    # 終点の位置：[0, 0]
        arrowprops=dict(  # 矢印の設定
            arrowstyle='->',       # 矢印の形状と向き
            facecolor='tab:red',   # 矢印の表面の色
            edgecolor='tab:red',   # 矢印の輪郭の色
        )
    )
    # textで競技名をプロット
    twin_ax.text(x=x, y=y, s=text, color='tab:red', fontsize=12)
# 修飾
twin_ax.set(xticks=np.arange(-1, 1.1, 0.5), yticks=np.arange(-1, 1.1, 0.5),
            xlim=(-1, 1), ylim=(-1, 1))
plt.show()

【実行結果】
ワークブックと比べて左右が逆転しています。
競技（の主成分負荷量）の方向に注目しましょう。

第１主成分と第３主成分のバイプロットを描画します。

# バイプロットの描画

## 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 7), tight_layout=True)

## 30選手の主成分得点の描画
# 主成分得点の散布図の描画（ダミー）
ax.scatter(scores_df['PC1'], scores_df['PC3'], marker='')
# 主成分得点の「順位」のテキストの表示
for x, y, s in zip(scores_df['PC1'], scores_df['PC3'], scores_df.index):
    ax.text(x=x, y=y, s=s, ha='center', va='center', color='royalblue',
            fontsize=12)
# x=0の垂直線、y=0の水平線の描画
ax.axvline(0, color='black', lw=0.5)
ax.axhline(0, color='black', lw=0.5)
# 修飾
ax.set_xlabel('第１主成分', fontsize=12)
ax.set_ylabel('第３主成分', fontsize=12)
ax.set(xlim=(-5, 5), ylim=(-5, 5))

## 7競技の主成分負荷量の描画
# 描画領域の設定：主成分得点と異なるx軸・y軸のスケールで描画するため
twin_ax = ax.twinx().twiny()
# 主成分負荷量の矢印と競技名の描画
for x, y, text in zip(loadings['PC1'], loadings['PC3'], loadings.index):
    # annotateで矢印をプロット
    twin_ax.annotate(
        text='',          # テキスト：なし
        xy=[x, y],        # テキストの位置:[第1主成分負荷量, 第2主成分負荷量]
        xytext=[0, 0],    # 終点の位置：[0, 0]
        arrowprops=dict(  # 矢印の設定
            arrowstyle='->',       # 矢印の形状と向き
            facecolor='tab:red',   # 矢印の表面の色
            edgecolor='tab:red',   # 矢印の輪郭の色
        )
    )
    # textで競技名をプロット
    twin_ax.text(x=x, y=y, s=text, color='tab:red', fontsize=12)
# 修飾
twin_ax.set(xticks=np.arange(-1, 1.1, 0.5), yticks=np.arange(-1, 1.1, 0.5),
            xlim=(-1, 1), ylim=(-1, 1))
plt.show()

【実行結果】

第２主成分と第３主成分のバイプロットを可視化します。

# バイプロットの描画

## 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 7), tight_layout=True)

## 30選手の主成分得点の描画
# 主成分得点の散布図の描画（ダミー）
ax.scatter(scores_df['PC2'], scores_df['PC3'], marker='')
# 主成分得点の「順位」のテキストの表示
for x, y, s in zip(scores_df['PC2'], scores_df['PC3'], scores_df.index):
    ax.text(x=x, y=y, s=s, ha='center', va='center', color='royalblue',
            fontsize=12)
# x=0の垂直線、y=0の水平線の描画
ax.axvline(0, color='black', lw=0.5)
ax.axhline(0, color='black', lw=0.5)
# 修飾
ax.set_xlabel('第２主成分', fontsize=12)
ax.set_ylabel('第３主成分', fontsize=12)
ax.set(xlim=(-5, 5), ylim=(-5, 5))

## 7競技の主成分負荷量の描画
# 描画領域の設定：主成分得点と異なるx軸・y軸のスケールで描画するため
twin_ax = ax.twinx().twiny()
# 主成分負荷量の矢印と競技名の描画
for x, y, text in zip(loadings['PC2'], loadings['PC3'], loadings.index):
    # annotateで矢印をプロット
    twin_ax.annotate(
        text='',          # テキスト：なし
        xy=[x, y],        # テキストの位置:[第1主成分負荷量, 第2主成分負荷量]
        xytext=[0, 0],    # 終点の位置：[0, 0]
        arrowprops=dict(  # 矢印の設定
            arrowstyle='->',       # 矢印の形状と向き
            facecolor='tab:red',   # 矢印の表面の色
            edgecolor='tab:red',   # 矢印の輪郭の色
        )
    )
    # textで競技名をプロット
    twin_ax.text(x=x, y=y, s=text, color='tab:red', fontsize=12)
# 修飾
twin_ax.set(xticks=np.arange(-1, 1.1, 0.5), yticks=np.arange(-1, 1.1, 0.5),
            xlim=(-1, 1), ylim=(-1, 1))
plt.show()

【実行結果】
１位・２位は原点に近めであり、順位と主成分の関係性を捉えにくいかもしれません。

いろんな手法で分析できました。
とても楽しかったです！！！

記事の最後はChatGPTが締めくくります。
今回はステージの演奏者に耳をすまして。

📘 ChatGPTのひとこと：

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。