「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第６章「機械学習への道」②主成分分析・クラスタリング

第６章「機械学習への道」

書籍の著者　山澤成康 先生

この記事は、書籍「回帰分析から学ぶ計量経済学」第６章「機械学習への道」の Python写経活動 を取り扱います。

今回は教師なし学習の主成分分析とクラスタリングを実践します。
では書籍を開いて回帰分析の旅に出発です🚀

はじめに

書籍「回帰分析から学ぶ計量経済学」のご紹介

このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。

テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。

私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。

回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！

引用表記

この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

第６章 機械学習への道

この記事は第６章の以下の節を取り扱います。

記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる Excel ファイル内のデータをもとにしてCSVファイルを作成し、data フォルダに格納しています。

第６章で用いるライブラリをインポートします。

### インポート

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

## 統計
import scipy.stats as stats

## 機械学習
# irisデータセット
from sklearn.datasets import load_iris
# データの準備
from sklearn.model_selection import train_test_split  # 学習・テストデータの分割
# 交差検証
from sklearn.model_selection import StratifiedKFold   # 層化Kfold交差検証
from sklearn.model_selection import cross_validate    # CVによる学習実行
# 分類アルゴリズム
from sklearn.linear_model import LogisticRegression   # ロジスティック回帰
from sklearn.tree import DecisionTreeClassifier       # 決定木
from sklearn.neighbors import KNeighborsClassifier    # K近傍法
from sklearn.svm import SVC                           # サポートベクターマシン
from sklearn.naive_bayes import GaussianNB            # ナイーブベイズ
from sklearn.ensemble import RandomForestClassifier   # ランダムフォレスト
from sklearn.decomposition import PCA                 # 主成分分析
from sklearn.cluster import KMeans                    # k-means法
from sklearn.cluster import AgglomerativeClustering   # 階層分析
# 評価
from sklearn.metrics import accuracy_score            # 正解率
from sklearn.metrics import classification_report     # 分類評価レポート
from sklearn.metrics import confusion_matrix          # 混同行列
from sklearn.tree import plot_tree                    # 決定木の描画

## 描画
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib

6.6 テスト成績を主成分分析にかける p.199

テキストの学生10人のテストの点数データ（仮想データ）を読み込みます。

### データの読み込み
df_test = pd.read_csv('./data/06_01_testscore.csv', index_col=0,
                  usecols=list(range(7)))
print('df_test.shape:', df_test.shape)
display(df_test)

【実行結果】
６教科のテストの点数です。
化学以降、点数の小数値が気になりますが、仮想なので気にしません。

学生別に教科別得点を折れ線グラフで可視化してみましょう。

### 折れ線グラフを描画

# 教科の並び順を変える
cols = ['数学', '化学', '物理', '日本史', '地理', '国語']

# 学生別の折れ線グラフの描画
sns.lineplot(data=df_test[cols].T)
# 修飾
plt.xlabel('科目')
plt.ylabel('得点')
plt.legend(bbox_to_anchor=(1, 1), title='学生')
plt.grid(lw=0.5)
plt.show()

【実行結果】
数学・化学・物理と日本史・地理・国語はそれぞれ関係がありそうに見えます。

教科ごとの相関係数を見てみましょう。
テキスト p.200 の図表「得点の相関係数行列」に相当します。

### 相関行列の算出
display(df_test.corr().round(3))

【実行結果】
数学・化学・物理の２変数間と、国語・日本史・地理の２変数間に強い正の相関が見られます。

相関行列をヒートマップで可視化してみましょう。

### 相関行列のヒートマップの描画
sns.heatmap(data=df_test.corr().round(3),
            vmin=0.5, vmax=1.15, cmap='Greens',
            annot=True, fmt='.3f', annot_kws={'fontsize': 12});

【実行結果】
強い相関（濃い緑）が一目瞭然です！

それでは主成分分析を実行しましょう。
テキストは「相関行列」を用いることとしているので、標準化したデータを用いて主成分分析を行います。

### 主成分分析の実行

# 設定
n_components = 6  # 主成分数
pc_names = [f'PC{i}' for i in range(1, n_components+1)]

# データの標準化 ※scipyのzscoreを使用
df_test_ss = pd.DataFrame(stats.zscore(df_test, ddof=1),
                          columns=df_test.columns, index=df_test.index)

# pcaの実行
pca = PCA(n_components=n_components)
result = pca.fit_transform(df_test_ss)

# 固有値、固有ベクトル、主成分得点、主成分負荷量のデータフレーム化
eig_val_df = pd.DataFrame({'固有値': pca.explained_variance_,
                           '寄与率': pca.explained_variance_ratio_,
                           '累積寄与率': np.cumsum(pca.explained_variance_ratio_)},
                           index=pc_names)
eig_vec_df = pd.DataFrame(pca.components_.T,
                          columns=pc_names, index=df_test_ss.columns)
pc_score_df = pd.DataFrame(result, columns=pc_names, index=df_test_ss.index)
loadings_df = eig_vec_df * eig_val_df['固有値'].apply(np.sqrt).T

# 主成分分析の結果の表示
print('主成分分析')
print('n =', pca.n_samples_, '\n')
print('相関行列の固有値解析')
display(eig_val_df.round(4))
print('固有ベクトル')
display(eig_vec_df.round(3))
print('主成分得点')
display(pc_score_df.T[:2].round(2))   # PC1～PC2を表示
print('主成分負荷量')
display(loadings_df.round(3))

【実行結果】
固有ベクトルの正負の符号の反転が見られますが、テキストの「gretl」の分析結果と同一の結果を得られました。
第２主成分までの累積寄与率が 99.6 ％という驚異的な寄与率になりました。

（ちなみに情報）
上記の固有値と固有ベクトルがデータの相関行列の固有値解析の結果と一致することを確認しましょう。

### データの相関行列の固有値解析

# 指数表記を禁止
np.set_printoptions(suppress=True)

# データの相関行列の固有値解析
eig_vals, eig_vecs = np.linalg.eig(df_test_ss.corr())
# 固有値の大きい順位並び替えるときのインデックス
sort_idx = np.argsort(eig_vals)[::-1]
# 固有値、固有ベクトルの表示
print('固有値:')
print(eig_vals[sort_idx].round(4))
print('固有ベクトル:')
print(eig_vecs[:, sort_idx].round(3))

# 指数表記をデフォルト値に戻す
np.set_printoptions(suppress=False)

【実行結果】
若干正負符号が異なりますが、絶対値は一致しました。

学生の主成分得点と教科の主成分負荷量をプロットします。
バイプロットみたいなものです（矢印は無いですが。。。）

### 主成分負荷量と主成分得点の描画

## 設定と準備 Listsの初期化
list1, list2 = [], []

## 科目の描画
# 主成分負荷量の散布図の描画（科目）
plt.scatter(loadings_df['PC1'], loadings_df['PC2'])
# 文系・理系の分離
for i in range(len(eig_vec_df)):
    if loadings_df.iloc[i, 1] > 0:
        list1.append(loadings_df.index[i])
        xy1 = loadings_df.iloc[i, :2]
    else:
        list2.append(loadings_df.index[i])
        xy2 = loadings_df.iloc[i, :2]
# 科目の表示
for text, xy in zip([','.join(list1), ','.join(list2)], [xy1+0.1, xy2+0.1]):
    plt.annotate(text=text, xy=xy)

## 学生の描画
# 主成分得点の散布図の描画（学生）
plt.scatter(pc_score_df['PC1'], pc_score_df['PC2'])
# 学生の表示
for i in range(len(pc_score_df)):
    plt.annotate(text=pc_score_df.index[i], xy=pc_score_df.iloc[i, :2]+0.04)

# x軸=0、y軸=0の描画
plt.axhline(0, color='black', lw=0.5)
plt.axvline(0, color='black', lw=0.5)
# 修飾
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('主成分負荷量（科目）と主成分得点（学生）のプロット')
plt.show()

【実行結果】
第１主成分は総合的な学力（マイナス値：学力が高い、プラス値：学力が高くない）、第２主成分は文系の強さ（プラス値）か理系の強さ（マイナス値）かの指標になりそうです。

6.7 クラスター分析 p.203

■ データの読み込み
あやめデータセットを使います。

### データの読み込み

# データセットの取得
iris = load_iris()

# データフレーム化
df = pd.DataFrame(np.hstack([iris.data, iris.target.reshape(-1, 1)]),
                  columns=iris.feature_names + ['class'])
df['class'] = df['class'].astype(int)

# アヤメの名前の取得
iris_names = iris.target_names

# 結果の表示
print('df.shape: ', df.shape)
display(df.head())
print('iris_names: ', iris_names)

【実行結果】

■ 非階層分析 p203
k-means法であやめデータセットを３つのクラスタに分けます。
クラスタリングは「教師なし学習」ですので、正解値（class）を使用せず、似たデータ点でクラスタ（グループ）を推定します。

### k-means法によるクラスタリングの実行

# クラスタリングの実行
result = KMeans(n_clusters=3, random_state=123).fit(df.iloc[:, :-1])
result.labels_

【実行結果】

クラスタを可視化しましょう。
実績値（正解値）も分かるようにしています。

### クラスタラベルの描画
df1 = df.copy()
df1['cluster'] = result.labels_
sns.scatterplot(data=df1, x='petal length (cm)', y='petal width (cm)', 
                hue='cluster', palette=['tab:green', 'tab:blue', 'tab:red'],
                style='class', markers=['D', 'o', 's'])
plt.show()

【実行結果】
クラスタは色（青、赤、緑）で表現しています。
いい感じにクラスタを識別できていると思います！

■ 階層分析 p.204
scikit-learn の凝集的クラスタリング AgglomerativeClustering を利用します。
クラスタ間の距離はウォード法（デフォルト値）、データ間の距離はユークリッド距離（デフォルト値）です。

### 凝集的クラスタリングによるクラスタリングの実行

# クラスタリングの実行
result = AgglomerativeClustering(n_clusters=3).fit(df.iloc[:, :-1])
result.labels_

【実行結果】

クラスタを可視化しましょう。

### クラスタラベルの描画
df2 = df.copy()
df2['cluster'] = result.labels_
sns.scatterplot(data=df2, x='petal length (cm)', y='petal width (cm)', 
                hue='cluster', palette=['tab:green', 'tab:blue', 'tab:red'],
                style='class', markers=['D', 'o', 's'])
plt.show()

【実行結果】
いい感じにクラスタリングできています！

デンドログラムを描画しましょう。
階層的クラスタリングはデンドログラムでクラスタ形成の様子を可視化できてよきです！

公式サイトのデンドログラム描画関数をお借りします。

### デンドログラム描画関数
# https://scikit-learn.org/stable/auto_examples/cluster/plot_agglomerative_dendrogram.html

from scipy.cluster.hierarchy import dendrogram

def plot_dendrogram(model, **kwargs):
    # Create linkage matrix and then plot the dendrogram

    # create the counts of samples under each node
    counts = np.zeros(model.children_.shape[0])
    n_samples = len(model.labels_)
    for i, merge in enumerate(model.children_):
        current_count = 0
        for child_idx in merge:
            if child_idx < n_samples:
                current_count += 1  # leaf node
            else:
                current_count += counts[child_idx - n_samples]
        counts[i] = current_count

    linkage_matrix = np.column_stack(
        [model.children_, model.distances_, counts]
    ).astype(float)

    # Plot the corresponding dendrogram
    dendrogram(linkage_matrix, **kwargs)

再度モデリングを行って、デンドログラムを描画します。

### デンドログラムの描画

# 凝集的クラスタリングの実行
model = AgglomerativeClustering(distance_threshold=0, n_clusters=None)
model.fit(df.iloc[:, :-1])

# デンドログラムの描画
plt.figure(figsize=(15, 7))
plot_dendrogram(model, color_threshold=10)

【実行結果】
150のデータ点を横並びにするのは（小さくなりすぎて）見にくいですかね。。。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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