「入門はじめての多変量解析」をPythonで写経 Vol.11 ～ ３章「はじめての主成分分析」③主成分の解釈、寄与率・累積寄与率、主成分得点

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

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

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

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

この記事は、情報損失量の最小化に基づく主成分分析 に取り組みます。
とりわけ、主成分の解釈、寄与率・累積寄与率、主成分得点 を学びます

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

はじめに

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

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

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

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

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

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

### インポート

# 数値計算
import string                           # 大文字アルファベットを取得
import numpy as np
import pandas as pd

# 描画
import matplotlib.pyplot as plt
from graphviz import Digraph            # パス図（有向グラフ）
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

はじめに

この記事で取り扱う「情報損失量の最小化」に基づく主成分分析は、（個人的感触ですが）メジャーな方法では無いようです。
たとえば、Python の scikit-learn や statsmodels などのライブラリは「分散最大化」に基づいています。
「情報損失量の最小化」と「分散最大化」の結果が一致する部分はありますが、結果が異なる部分もありますし、結果を導くプロセスが異なる面もあります。
ご留意下さいませ。

分析の準備・前回までのダイジェスト

■ データの登録
テキスト p.87 表 3.1.1「介護施設と医療施設（その１）」のデータをお借りします。

### 介護・医療施設データ p.87 表3.1.1

# データの登録
N = 10
data1 = pd.DataFrame(
    {'地域名': list(string.ascii_uppercase[:N]),
     '介護施設': [22, 22, 18, 18, 15, 19, 19, 24, 21, 25],
     '医療施設': [12, 8, 6, 15, 7, 9, 7, 17, 14, 11]},
     index=range(1, N+1))
data1.index.name = 'No.'
data1

【実行結果】
10 の地域における 65 歳以上人口 1 万人あたりの介護施設数、人口 1 万人あたりの医療施設数です。

【補足説明】
地域名設定時に利用した string.ascii_uppercase を補足します。
大文字アルファベットを A から順に取得できる優れものなのです！
大文字アルファベットの先頭 10 文字をリストに格納してみましょう。

### アルファベット大文字を取得
list(string.ascii_uppercase[:10])

【実行結果】

■ その他の設定
その他の準備をまとめて実施します。
地域名をインデックスにするときに利用したい area の作成と、変数だけを取り出したデータフレーム X の作成です。

## 設定と準備

# 変数Xの設定
X = data1.iloc[:, 1:]
display(X.head())

# 地域名
area = data1['地域名']

【実行結果】

■ 固有値・固有ベクトルの算出
情報損失量の最小化に基づく主成分分析の固有値・固有ベクトルを算出します。
前回記事のダイジェストです。

### 固有値・固有ベクトルの取得 ※前回のダイジェスト

# p.100の行列の作成 ※データの分散共分散行列にN-1を掛けて非対角成分の符号を逆転
the_matrix = X.cov() * (N - 1) * np.array([[1, -1], [-1, 1]])

# 上記行列の固有値・固有ベクトルの算出
eig_val1, eig_vec1 = np.linalg.eig(the_matrix)

# 固有ベクトルの変更 ※行の入れ替え、符号の逆転
eig_vec1 = eig_vec1[::-1] * -1

# 結果表示
print('固有値:')
print(eig_val1)
print('固有ベクトル:')
print(eig_vec1)

【実行結果】

■ 第１主成分の数式表現
情報損失量の最小化の場合、最小固有値に対応する固有ベクトルを用いて、第１主成分

が導かれました。

主成分の解釈

※「分散最大化」も同じように主成分の解釈が必要です。

■ 主成分の解釈の必要性
求めた第１主成分は、２つの変数「介護施設」「医療施設」を「総合化」したものと考えられます。
でも総合化で生まれた「第１主成分の意味」はどこにも書かれていないですよね…

主成分分析において、主成分の意味は分析者が解釈して導き出す必要があるのです！

■ 主成分の解釈をサポートする可視化
テキストの「パス図」を用いて、主成分の意味を考えてみましょう。
パス図は「変数と主成分の関係を示す図」として用いられているようです。
描いてみましょう。

graphviz ライブラリの 有向グラフ Digraph を利用してパス図を描きます。
テキスト p.102 図 3.4.1「主成分分析のパス図」に相当します。

### 主成分分析のパス図 p.102 図3.4.1

## 設定
# 有向グラフオブジェクトの生成、neatoでnodeの位置調整を実施
g = Digraph(engine='neato')
# nodeの基本属性の設定
g.attr('node', shape='box', fontname='Meiryo UI')

## node：頂点の作成、posで位置固定
g.node('介護施設x1', pos='0, 1!')
g.node('医療施設x2', pos='0, 0!')
g.node('主成分z1', pos='2, 0.5!', shape='oval')

## edge：辺の作成
g.edge('介護施設x1', '主成分z1', label=f'{eig_vec1[0, 0]:.4f}')
g.edge('医療施設x2', '主成分z1', label=f'{eig_vec1[1, 0]:.4f}')

## グラフの表示
g

【実行結果】
右の四角形は２つの変数、右の楕円は主成分。
矢印は変数から主成分に向けられていて、係数 $${a_1, a_2}$$（固有ベクトル）が添えられています。

■ 第１主成分の解釈
テキストにならって図を読み解きましょう。

この関係性から第１主成分を解釈すると
「介護施設・医療施設の両方の大きさ」⇒「地域の介護・医療の充実度」
と考えられます！

テキストは「福祉の充実度」と解釈しています。
この解釈でパス図を上書きしてみましょう。
テキスト p.103 図 3.4.4「主成分に名前を付ける！」に相当します。

### 主成分に名前を付ける！ p.102 図3.4.4

## 設定
# 有向グラフオブジェクトの生成、neatoでnodeの位置調整を実施
g = Digraph(engine='neato')
# nodeの基本属性の設定
g.attr('node', shape='box', fontname='Meiryo UI')

## node：頂点の作成、posで位置固定
g.node('介護施設x1', pos='0, 1!')
g.node('医療施設x2', pos='0, 0!')
g.node('福祉の充実度', pos='2, 0.5!', shape='oval', 
       fillcolor='lavender', style='filled')

## edge：辺の作成
g.edge('介護施設x1', '福祉の充実度', label=f'{eig_vec1[0, 0]:.4f}')
g.edge('医療施設x2', '福祉の充実度', label=f'{eig_vec1[1, 0]:.4f}')

## グラフの表示
g

【実行結果】

前回記事の主成分の直線を添えたデータの散布図で、主成分の解釈を深堀りしましょう。

【解釈】
赤い直線が第１主成分 $${z_1}$$ 、真ん中の赤い点が平均（中心）です。
第１主成分 $${z_1}$$ が「地域の介護・医療の充実度」と解釈すると…

• 中心より右上の方向に位置する地域（青い点）は「介護・医療が充実している地域」と考えられます。

• 中心より左下の方向に位置する地域は「介護・医療が充実しているとはいいにくい地域」と考えられます。

寄与率・累積寄与率

※「分散最大化」と計算方法が異なります。

■ 寄与率の概要
「寄与率」は各主成分がデータ全体をどの程度説明しているかを示す指標です。
下の図で確認しましょう。

各主成分は「新しい情報量」でデータを説明します。
「情報損失量最小化」に基づく主成分分析の「固有値 $${\lambda}$$ 」は「情報損失量の二乗和」を示します。
※前回記事で「固有値が情報損失量の二乗和 $${U}$$ から導かれたこと」を思い出しましょう。

ということで、テキストの公式をお借りします。

📊 寄与率の定義（情報損失量最小化の場合）

🔢 寄与率の計算
分析データの第１主成分の寄与率を計算しましょう。

### 寄与率の算出 p.105

# 2変数の平均値の算出
x1_bar, x2_bar = X.mean().values

# 元の情報量の平方和の算出
sum_squared_info_form = sum([(x1 - x1_bar)**2 for x1 in X['介護施設']]) \
                      + sum([(x2 - x2_bar)**2 for x2 in X['医療施設']])

# 情報損失量の平方和（第１主成分の固有値）の算出
sum_squared_info_loss1 = eig_val1[0]

# 主成分z1の寄与率
(sum_squared_info_form - sum_squared_info_loss1) / sum_squared_info_form

【実行結果】
寄与率は $${76.6\%}$$ です。

■ 分散最大化による主成分分析の寄与率・累積寄与率
少し、分散最大化に寄り道します。

「分散最大化」の場合の固有値は「各主成分がデータのバラツキをどれだけ説明しているかを表す指標」です。
「主成分が元のデータをどれだけ説明しているかを表す指標」との説明もよく見かけます。
このような固有値を使って寄与率を求めます。
ざっくり寄与率は「固有値 ÷ 固有値の合計」であり、元のデータ全体のうち主成分が説明できている割合なのです。

📊 寄与率・累積寄与率の定義（分散最大化の場合）
説明変数の数＝主成分の最大数＝$${p}$$ とします。
固有値 $${\lambda}$$ は大きな値の順に並び替え済みとします。

分母は全固有値の合計です。
累積寄与率の分子は第 $${i}$$ 主成分までの固有値の累計です。

累積寄与率の大きさは新しい情報量＝主成分得点の分散だそうです。
データの要約や次元削減の際には、少ない新しい変数＝主成分を選びたいでしょう。
そんなときに累積寄与率が活躍します。

例えばデータの 70% が説明できていればよい、という風にしきい値を決めておきます。
この場合、累積寄与率が 70% に達するまでの主成分を選択すればよいことになります。

サクッと分散最大化の場合の寄与率を算出しましょう。

### 分散最大化の寄与率の算出

# 分散最大化の固有値の算出 ※データの分散共分散行列の固有値
eig_vec1_bunsan, _ = np.linalg.eig(X.cov())

# 固有値を降順でソート
eig_vec1_bunsan = eig_vec1_bunsan[eig_vec1_bunsan.argsort()[::-1]]

# 寄与率の算出 ※各主成分の固有値 / 固有値の合計
kiyo_rate = eig_vec1_bunsan / eig_vec1_bunsan.sum()

# 結果の表示
print('固有値　　　　　　:', eig_vec1_bunsan)
print('第１主成分の寄与率:', kiyo_rate[0])
print('第２主成分の寄与率:', kiyo_rate[1])

【実行結果】

分散最大化の場合の第１主成分の寄与率は、情報損失量最小化の寄与率と（ほぼほぼ）同じになりました！
70% をしきい値にする場合、第１主成分を選択すれば良いことになります！

主成分得点

※「分散最大化」も同じ計算式で主成分得点を求めます。

■ 主成分得点の概要
下の図の「新しい情報量」が主成分得点です。
データの中心（平均）と主成分軸上のデータ点（赤い点）の長さが主成分得点です。

📊 主成分得点の公式
テキストの計算式をお借りします。
第１主成分の主成分得点は、変数の平均 $${\bar{x}_1, \bar{x}_2}$$、第１主成分の係数（固有ベクトル） $${a_1, a_2}$$ を用いて、

となります。

🔢 主成分得点の計算と読み解き
分析データの第１主成分得点を Python で計算します。

### 主成分得点の算出 p.106

# 主成分得点の算出 ※固有ベクトル1 × 中心化したx1 + 固有ベクトル2 × 中心化したx2
score = (X - X.mean()) @ eig_vec1[:, 0]

# データフレーム化
score = score.rename('主成分得点').to_frame()
score.index = area

# 結果の表示
score

【実行結果】
各地区の第１主成分の主成分得点を計算しました。

【読み解き】
第１主成分の解釈は「地区の介護・医療の充実度」でした。
また、主成分得点の原点は０です。
つまり、この主成分得点が正の値の場合には、介護・医療が充実している度合いを示すと言えるでしょう。
反対に、この主成分得点が負の値の場合には、介護・医療が充実していない度合いを示すと言えそうです。

絶対値の大きな地域を見てみると…
地域 H はとても充実して、地域 E、C の充実度は低い、と考えられます。

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

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