「入門はじめての多変量解析」をPythonで写経 Vol.14 ～ ３章「はじめての主成分分析」⑥多重共線性の影響を検討

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

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

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

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

この記事は、主成分分析が多重共線性の影響を受けるかどうか を確かめるものです。
多重共線性は重回帰分析に悪影響を及ぼす事象でした。

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

はじめに

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

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

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

この記事は３章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

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

### インポート

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

# PCA
from sklearn.decomposition import PCA
from statsmodels.multivariate.pca import PCA as sm_PCA

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

主成分分析が多重共線性の影響を受けるか？ p.118～

■ 多重共線性
多重共線性は変数間の相関係数が高い場合に生じる事象です。
以前、重回帰分析が多重共線性の影響を受けることを検討しました。
テキストは、説明変数に多重共線性が存在する場合、「分散共分散行列の逆行列が存在しないため重回帰分析ができない」ことを問題視しました。

テキストは、主成分分析においても「分散共分散行列の逆行列が存在しないことが影響するか」の確認に進みます。

１つめのデータで主成分分析

■ データの準備
テキスト p.118 表 3.11.1「多重共線性のあるデータ（その１）」をお借りします。

### 多重共線性のあるデータ（その１） p.118 表3.11.1

# データの登録
data4 = pd.DataFrame(
    {'x1': [2, 1, 2, 3, 4, 4, 1, 3, 2, 1],
     'x2': [2, 1, 5, 9, 1, 4, 1, 9, 2, 7],
     'x3': [4, 2, 5, 8, 7, 8, 2, 8, 4, 4],
     'x4': [2, 1, 3, 5, 3, 4, 1, 5, 2, 3]},
     index=range(1, 11))
data4.index.name = 'No.'

# 結果の表示
data4

【実行結果】
４つの変数、10 個のデータです。

相関係数を確認します。

# 相関係数の表示
data4.corr().round(3)

【実行結果】
$${x_3, x_4}$$、$${x_1, x_3}$$ の相関係数が大きいようです。

データの分散共分散行列に逆行列が存在しないことを確認します。
分散共分散行列のランクが変数の数 $${4}$$ 未満の場合、逆行列が存在しません。

# 分散共分散行列のランク
cov4_rank = np.linalg.matrix_rank(data4.cov(ddof=1))
print(f'分散共分散行列のランクは {cov4_rank} です')

【実行結果】
ランク $${2 < 4}$$ ですので、逆行列が存在しません。

■ ヘルパー関数の定義
固有値計算と主成分負荷量計算の関数を作ります。

### ヘルパー関数

# 固有値データフレーム作成関数
def make_variance_df(eig_vals):
    return pd.DataFrame({'固有値': eig_vals,
                         '寄与率': eig_vals / sum(eig_vals),
                         '累積寄与率': np.cumsum(eig_vals) / sum(eig_vals)},
                        index=[f'PC{i+1}' for i in range(len(eig_vals))])

# 主成分負荷量データフレーム作成関数 ※想定：データが標準化されている
def make_loadings_df(eig_vals, eig_vecs, df, use_statsmodels=False):
    # statsmodelsの場合、N-1でスケーリングする
    scale_val = len(df) - 1 if use_statsmodels else 1
    # eig_vecs.shape=(変数数, 主成分数)
    return pd.DataFrame(
        eig_vecs * (eig_vals / scale_val)**(1/2), index=df.columns,
        columns=[f'PC{i+1}' for i in range(eig_vecs.shape[1])])

scikit-learn、statsmodels、分散共分散行列の固有値・固有ベクトルの３つのケースを確認します。
各ライブラリが多重共線性の影響でエラーになったら困ります…

作成する指標は固有値と主成分負荷量です。
主成分負荷量は統計解析ソフト SPSS の成分行列に相当します。

🖲️scikit-learn

### sklearnのPCA p.119 表3.11.2

# データの標準化
X4_std = (data4 - data4.mean()) / data4.std(ddof=1)

# PCAインスタンスの生成
n_components = 4
pca4 = PCA(n_components=n_components)
# PCAの実行、主成分得点の算出
score4 = pca4.fit_transform(X4_std)

# 固有値・寄与率の表示
display(make_variance_df(pca4.explained_variance_).round(4))

# 主成分負荷量の表示
display(make_loadings_df(
    pca4.explained_variance_, pca4.components_.T, X4_std).round(3))

【実行結果】
無事、算出できました。
下の表は主成分負荷量です。

第１主成分と第２主成分が求まりました。
累積寄与率は最初の２つの主成分で合計１になっており、２つの主成分でデータのすべてを説明できています！

「主成分の数は分散共分散行列のランクまで」となるようです。
ChatGPTに訊いてみました。

🖲️statsmodels

### statsmodelsのPCA  p.119 表3.11.2

# PCAの実行
pca4_sm = sm_PCA(X4_std, standardize=False)

# 固有値・寄与率の表示
res = make_variance_df(pca4_sm.eigenvals.values)
display(res.round(4))

# 主成分負荷量の表示
display(make_loadings_df(
    res['固有値'].values, pca4_sm.eigenvecs.values, X4_std, True).round(3))

【実行結果】
無事、算出できました。
固有値の値は scikit-learn と異なります。

🖲️分散共分散行列の固有値・固有ベクトル

ワーニング非表示のおまじないを唱えます。

# 追加インポート
import warnings
# ワーニング非表示
warnings.filterwarnings('ignore')

計算します。

### 分散共分散行列の固有値・固有ベクトル・主成分負荷量   p.119 表3.11.2

# 分散共分散行列の固有値・固有ベクトルの算出
eig_vals4_cov, eig_vecs4_cov = np.linalg.eig(X4_std.cov())
sort_idx = eig_vals4_cov.argsort()[::-1]

# 固有値・寄与率の表示
display(make_variance_df(eig_vals4_cov[sort_idx]).round(4))

# 主成分負荷量の表示
display(make_loadings_df(
    eig_vals4_cov[sort_idx], eig_vecs4_cov[:, sort_idx], X4_std).round(3))

【実行結果】
無事、算出できました。

２つめのデータで主成分分析

テキストに用意された２つめのデータの主成分分析に移ります。
結論的には、１つめのデータと同様に、主成分分析を実行できます。

■ データの準備
テキスト p.120 表 3.11.3「多重共線性のあるデータ（その２）」をお借りします。

### 多重共線性のあるデータ（その２） p.120 表3.11.3

# データの登録
data5 = pd.DataFrame(
    {'x1': [2, 1, 2, 3, 4, 4, 1, 3, 2, 1],
     'x2': [2, 1, 5, 9, 1, 4, 1, 9, 2, 7],
     'x3': [4, 2, 5, 8, 7, 8, 2, 8, 4, 4],
     'x4': [2, 1, 3, 5, 3, 4, 1, 5, 2, 3],
     'x5': [4, 2, 4 ,6, 8, 8, 2, 6, 4, 2]},
     index=range(1, 11))
data5.index.name = 'No.'

# 結果の表示
data5

【実行結果】
５つの変数、10 個のデータです。

相関係数を確認します。

# 相関係数の表示
data5.corr().round(3)

【実行結果】
相関係数の高い変数だらけです…

データの分散共分散行列に逆行列が存在しないことを確認します。
分散共分散行列のランクが変数の数 $${5}$$ 未満の場合、逆行列が存在しません。

# 分散共分散行列のランク
cov5_rank = np.linalg.matrix_rank(data5.cov(ddof=1))
print(f'分散共分散行列のランクは {cov5_rank} です')

【実行結果】
ランク $${2 < 5}$$ ですので、逆行列が存在しません。
第２主成分まで算出可能となるでしょう。

🖲️scikit-learn

### sklearnのPCA p.121 表3.11.4

# データの標準化
X5_std = (data5 - data5.mean()) / data5.std(ddof=1)

# PCAインスタンスの生成
n_components = 5
pca5 = PCA(n_components=n_components)
# PCAの実行、主成分得点の算出
score5 = pca5.fit_transform(X5_std)

# 固有値・寄与率の表示
display(make_variance_df(pca5.explained_variance_).round(4))

# 主成分負荷量の表示
display(make_loadings_df(
            pca5.explained_variance_, pca5.components_.T, X5_std)
        .round(3))

【実行結果】
無事、算出できました。
第１主成分と第２主成分が求まりました。

🖲️statsmodels

### statsmodelsのPCA  p.121 表3.11.4

# PCAの実行
pca5_sm = sm_PCA(X5_std, standardize=False)

# 固有値・寄与率の表示
res = make_variance_df(pca5_sm.eigenvals.values)
display(res.round(4))

# 主成分負荷量の表示
display(
    make_loadings_df(
        res['固有値'].values, pca5_sm.eigenvecs.values, X5_std, True)
    .round(3))

【実行結果】
無事、算出できました。
固有値の値は scikit-learn と異なります。

🖲️分散共分散行列の固有値・固有ベクトル

### 分散共分散行列の固有値・固有ベクトル・主成分負荷量   p.121 表3.11.4

# 分散共分散行列の固有値・固有ベクトルの算出
eig_vals5_cov, eig_vecs5_cov = np.linalg.eig(X5_std.cov())
sort_idx = eig_vals5_cov.argsort()[::-1]

# 固有値・寄与率の表示
display(make_variance_df(eig_vals5_cov[sort_idx]).round(4))

# 主成分負荷量の表示
display(
    make_loadings_df(
        eig_vals5_cov[sort_idx], eig_vecs5_cov[:, sort_idx], X5_std)
    .round(3))

【実行結果】
無事、算出できました。
scikit-learn と比べて、正負の符号が逆転していますが、値は一致しています。

記事の最後は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の教科書です。
よかったらぜひ、お試しくださいませ。

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