Zenn
Aidemy Tech BlogPublicationへの投稿
📖

教師なしクラスタリングのその先へ~次元削減と定量解析~

Shingo Nakazawa
2022/12/16に公開
Python
scikit-learn
機械学習
データ分析
clustering
tech

これはAidemy Advent Calendar 2022 16日目の記事です。

また、本記事に掲載したコードをまとめたNotebookをGitHubにて公開しています。

はじめに

こんにちは。株式会社アイデミーデータサイエンティストの中沢(@shnakazawa_ja)です。普段は社内のデータ基盤整備、クライアント様のデータ解析、R&Dといった業務に携わっています。

本稿ではPythonを用いたデータのクラスタリングとその先の解析について記載します。
クラスタリングは教師なし学習の代表例として挙げられることが多い一方、その結果の可視化や続く解析、解釈についての情報は少ないように思います。
2022年プロ野球の投手成績を題材に、チュートリアルのクラスタリングから一歩進んだ世界の一例を皆様にご紹介できればと思います。

本稿のメイントピックは「クラスタリング結果の可視化」以降になりますので、背景情報はいらないよという方はそちらから読み始めていただければと思います。

データの準備

データの取得

データはBaseball Reference様よりセ・リーグ、パ・リーグ12球団の2022年投手1軍成績をお借りします。

dataというフォルダを作ってチームごとにcsvとして保存します。(チーム名はBaseball Referenceでの表記にならう)

.
├── 解析用ファイル(.py, .ipynb, etc)
└── data/
    ├── ChibaLotteMarines.csv
    ├── ChunichiDragons.csv
    ├── FukuokaSoftbankHawks.csv
    ├── HanshinTigers.csv
    ├── HiroshimaCarp.csv
    ├── HokkaidoNipponHamFighters.csv
    ├── OrixBuffaloes.csv
    ├── SaitamaSeibuLions.csv
    ├── TohokuRakutenGoldenEagles.csv
    ├── YakultSwallows.csv
    ├── YokohamaBayStars.csv
    └── YomiuriGiants.csv

これをpandasで読み込みましょう。1球団1csvであるため、チーム名に相当するカラムを増やしつつ、全データを結合したDataFrameを作ります。

import pandas as pd
from pathlib import Path

teams = Path('data').rglob('*.csv')
for i, team in enumerate(teams):
    if i == 0:
        df = pd.read_csv(team)
        df['Team'] = team.stem
    else:
        df_tmp = pd.read_csv(team)
        df_tmp['Team'] = team.stem
        df = pd.concat([df, df_tmp], axis = 0)
df

データの前処理

欠損値の処理

データを眺めると以下のような欠損値の存在に気が付きます。

  • 'W-L%'(勝率), 'SO/W'(奪三振/四死球. K/BBと同じ)の値が欠損している投手がいる
  • 'GS'(先発登板数), 'GF'(試合終了投手になった数)は記録が全く無い
  • 中日・岩嵜投手のレコードはほとんどが欠損値となっている

欠損値の扱いも非常に深いテーマではありますが、本稿のメインではないので、ここでは簡単のためにこれらの説明変数/レコードは削除してしまいます。

(余談ですが、'W-L%'も'SO/W'も割り算値なので分母の大小が無視されています。今回これらが欠損値になっている投手のほとんどは登板数が少ない投手です。解析を始める前にそもそもそういった選手の割り算値を、たくさん試合に出ている選手の値と同等に扱うことが適切なのかを考える必要があります。欠損値がなくともカラムごと削除するべきかもしれませんし、割り算値を明示的に与えることがモデルの改善に役立つケースもあります。実践的には、複数の欠損値処理パターンを作り解析、それぞれの出力結果を見て考察、いずれの処理を適用するかを決めたいように思います。(※自分にとって都合の良い結果を恣意的に選ばないようには注意))

本題に戻りまして、'Rk'と'Notes', 'Name-additional'はインデックスや選手を区別するための補足情報等であるので、不要なカラムとして削除してしまいます。

df_tmp = df[df['Name'] != 'Sho Iwasaki']
df_dropped = df_tmp.drop(['Rk','W-L%', 'GS', 'GF', 'SO/W', 'Notes', 'Name-additional'], axis=1)
df_dropped

説明変数のみのDataFrameを作る

選手名、チーム名はクラスタリングをするにあたり不要な情報なので、それらを除いたデータフレームを作成します。

df_features = df_dropped.drop(['Name', 'Team'], axis=1)

正規化

全ての説明変数を同等に扱うため、それぞれの説明変数ごとに最小0/最大1に合わせる正規化を行います。この処理の有無でこの後の次元削減の結果が影響されます。

from sklearn import preprocessing
mm = preprocessing.MinMaxScaler()
df_features_mm = pd.DataFrame(mm.fit_transform(df_features), columns=df_features.columns)
df_features_mm

ここまでで、前処理は完了です。

クラスタリング

次に、ここまでで整備してきたデータをクラスター分けします。
目的によって選択すべき手法は変わりますが、本記事のメインはクラスタリング後の解析であるため、ここでは簡単のためにチュートリアルでのド定番 "k-means"を用い、クラスター数は8個と決めてしまいます。

補足：データの次元が大きい場合はクラスタリングの前に後述する次元削減を行うケースもあります(むしろその方が一般的か)。計算時間がとても長かったり、うまくクラスタリングができないときは検討してみてください。

from sklearn.cluster import KMeans

k = 8 # クラスター数を指示

kmeanModel = KMeans(n_clusters=k, random_state=42)
kmeanModel.fit(df_features_mm)
clusters = kmeanModel.labels_

# クラスターごとに何サンプルあるか
for i in range(k):
    num = list(clusters).count(i)
    print(f'Cluster {i}: n = {num}')
    
>>Cluster 0: n = 24
>>Cluster 1: n = 11
>>Cluster 2: n = 15
>>Cluster 3: n = 35
>>Cluster 4: n = 36
>>Cluster 5: n = 76
>>Cluster 6: n = 61
>>Cluster 7: n = 95

353人の選手を8つのクラスターに分けられました。

クラスタリング結果の可視化

ここからが本記事のメインです！
まず最初に、結果をどう可視化するかという話があります。

そのままプロット

チュートリアルでは適当な説明変数を縦横にとり可視化している場合が多いです。元データの次元数が少ないので、その方法でもそれらしいクラスターが見られます。
今回のデータでも同じ様に行ってみましょう。とりあえず登板数と勝数でプロットしてみましょうか。

import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.family'] = 'Hiragino Sans' # 日本語を表示できるようにする
# Windowsの方は mpl.rcParams['font.family'] = 'MS Gothic' 等

plt.scatter(df_features['G'], df_features['W'],
            c=clusters, cmap='Accent', alpha=0.7) 
plt.title('そのままプロット')
plt.xlabel('登板数')
plt.ylabel('勝数')
plt.colorbar()
plt.show()

よくわからん図ができました。「茶色と紫が登板数が多いのは分かるが、2クラスター間の違いはわからない。青と薄橙も勝ち星が多いのは分かるが……」程度の考察はできますが、これ以上の考察には別の説明変数の組み合わせで可視化したくなります。しかし、全ての組み合わせを1枚1枚見ていくのは現実的ではありません。
では1枚の図でもっとたくさんの情報の表現できないでしょうか？その目的で使える手法の一つが次元削減になります。

次元削減は色々な目的で用いられますが、可視化の文脈では多次元のデータを2次元まで落とすことでデータの特徴を最も表す2軸を抽出、つまり平面にプロットしたときにデータの分離が最も良くなる(と考えられる)2軸を作出することが使用目的になります。

PCA

次元削減の代表として、主成分分析 (principal component analysis: PCA) をまず使ってみましょう。

from sklearn.decomposition import PCA
pca = PCA(random_state=42)
pca.fit(df_features_mm)
score = pd.DataFrame(pca.transform(df_features_mm), index=df_features.index)
score

具体的にどういう変換が行われているのかはググれば良い資料がたくさん出てきます(例)が、可視化の文脈で言うと元のデータが作る多次元空間の中に、データの分散を最も保てる方向に平面を置いて、その平面に対して各レコードを落としたときの数値が第一主成分(上の表で言う一番左のカラム)と第二主成分(左から2番目)のスコアというイメージになります。

第一主成分 (principal component 1: PC1), 第二主成分 (PC2)をそれぞれ横・縦軸にとってプロットします。

plt.scatter(score.iloc[:,0], score.iloc[:,1],
            c=clusters, cmap='Accent', alpha=0.7) 
plt.title('PCA plot')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.colorbar()
plt.show()

先程よりもクラスターの分離が良くなりました。

それぞれのPCに元の説明変数がどのように寄与しているかを調べれば、この状態でもそれぞれのクラスターがどういう性質を持つか十分に議論ができそうです。

しかし、まだクラスター2 (薄橙) と4 (青)はほぼ重なっていたり、可視化の観点ではもっと工夫ができそうです。

t-SNE

次に紹介するのは非線形の次元削減手法のt-SNEです。(t-Distributed Stochastic Neighbor Embedding. ティスニーと読みます)

非線形であるがゆえに、元の説明変数間の線形の関係を保ちたかったり、値を戻す必要がある場合は使えませんが、データ間の関係を低い次元に落とし表現することはPCAより上手にできるケースが多いです。

scikit-learnに含まれているので簡単に実行できます。

from sklearn.manifold import TSNE
embedding = TSNE(random_state=42).fit_transform(df_features_mm)
plt.scatter(embedding[:, 0], embedding[:, 1],
    c=clusters, cmap='Accent', alpha=0.7)
plt.title('t-SNE plot')
plt.colorbar()
plt.show()

PCAよりもクラスターの分離がハッキリとしました。特に2と4、1と6の分離が顕著です。

UMAP

もう一つ、UMAP (ユーマップ) という手法も使ってみましょう。

こちらもt-SNE同様、非線形の手法ですが、t-SNEよりも速く、似たデータはより凝集するのが特徴とされています。(が、そう単純なものではなく、色々と議論が続いています。参考1、参考2)

こちらのライブラリを使うと簡単に実行できます。ライブラリ名がumapではなくumap-learnである点に注意してください。(umapはPython2向けの別パッケージ)

import umap
mapper = umap.UMAP(random_state=42)
embedding = mapper.fit_transform(df_features_mm)
plt.scatter(embedding[:, 0], embedding[:, 1],
    c=clusters, cmap='Accent', alpha=0.7)
plt.title('UMAP plot')
plt.colorbar()
plt.show()

今回のデータにおいてはt-SNEよりも綺麗にクラスターの表現ができているように見えます。

(補足)t-SNEとUMAPのパラメータ設定

ちなみに、今回はパラメータチューニングを一切行っていませんが、いずれの手法もその設定が結果に劇的な影響を与えます。
パラメータを変えるとどのような振る舞いを見せるのかをインタラクティブに体験できるサイトがそれぞれありますのでご紹介しておきます。

後者のサイトではデータの解釈に関しても重要な示唆が提示されており、特に以下の2文は頭においておくとよいでしょう。

  1. Cluster sizes in a UMAP plot mean nothing (UMAPのプロットにおけるクラスタサイズには意味がない)
  2. Distances between clusters might not mean anything (クラスタ間の距離は意味をなさないかもしれない)

前述の参考1、参考2も合わせ、非線形次元削減の結果の解釈性は今後も議論が続きそうです。

クラスターの特徴の同定

ここまででクラスタリングが終わり、投手を8つのタイプに分類できました。そして可視化を行い、クラスターがある程度綺麗に分かれる形で2次元に表現できました。
次のステップとしては、それぞれのクラスターがどういう特徴を持つのかが当然知りたくなります。

説明変数の大小の可視化

最初の1手として、先程作成した2次元プロット上でそれぞれの説明変数の大小を色で表現します。今回はベースとしてUMAPのプロットを採用します。
ここで作るプロットと、クラスターで色分けしたグラフとを並べて見ることで、それぞれのクラスターの特徴の定性的な考察ができます。

fig = plt.figure(figsize=(15,15))
cols = df_features.columns
for i, col in enumerate(cols):
    ax = fig.add_subplot(5, 5, i+1, title=col)
    ax.scatter(embedding[:, 0], embedding[:, 1],
        c=df_features[col], cmap='plasma', alpha=0.8)
fig.tight_layout()
plt.show()

こうしてみると

  • 左の島は勝ち (W) 負け (L) イニング数 (IP) が多い
  • 下に行くに従って登板数 (G) が増える
  • 右下の方は登板数や申告敬遠 (IBB) が多い
  • 右下の独立した島は突出してセーブ数 (SV) が多い
  • 右上は防御率 (ERA) 与四球率 (BB9) が高い

などといったことが分かります。

定量評価とVolcano plot

定性的な評価と合わせ、定量的な評価も行います。加えて、評価結果をVolcano plotという方法で可視化してみましょう。
Volcano plotは

  • 対象クラスターにおける説明変数の平均値が、他のクラスターの平均値と比較して何倍あるか(何分の1か)
  • その差が統計的仮説検定で有意差がつくか

を同時に可視化する手法です。一般に横軸が数値の比率のlog2、縦軸がp値の-log10を取ることが多いです。

ここであの悪名高いp値が登場します。かつ、典型的な多重検定であるので、検定結果に対してBenjamini-Hochberg (BH) 法などでp値の補正が必要となります。(補正された結果をq値と呼ぶこともあるようです)

(p値に関して色々と思うところはあるのですが、議論しだすと一人でアドベントカレンダーを埋めることになってしまうので、ここでは割愛します)

パッケージを使ってもいいのですが、馴染みのない表現方法かもしれませんので中身を知るために書いてみます。

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

df_features_with_clusters = df_features.copy()
df_features_with_clusters['cluster'] = clusters
cols = df_features.columns

# 閾値を定める。ここでは補正後のp値 (q値) が0.05以下かつ数値の比 (Fold Change:fc)が2倍以上あるか、を閾値とします。
q_threshold = 0.05
fc_threshold = 2

fig = plt.figure(figsize=(10,10))
# クラスターごとに評価する
for i in range(k):
    p_values = []
    fcs = []
    # 変数 = プロット上の1点
    for col in cols:
        # 検定
        group_1 = df_features_with_clusters[df_features_with_clusters['cluster'] == i][col]
        group_2 = df_features_with_clusters[df_features_with_clusters['cluster'] != i][col]
        p_value = stats.ttest_ind(group_1, group_2, equal_var=False)[1]
        p_values.append(p_value)
        
        # Fold change. 平均での比較が不適切であればここをmedian等に変える
        fc = group_1.mean()/group_2.mean()
        fcs.append(fc)

    # p-valueの補正
    q_values = multipletests(p_values, method='fdr_bh')[1]

    # 閾値を超えたものは色を変える
    colors = []
    for col, q_value, fc in zip(cols, q_values, fcs):
        # 対象がその他の2倍大きいときはオレンジ
        if q_value < q_threshold and fc > fc_threshold:
            colors.append('orange')
        # その他が対象の2倍大きいときは水色
        elif q_value < q_threshold and fc < 1/fc_threshold:
            colors.append('skyblue')
        # 大きな違いがない場合は灰色
        else:
            colors.append('gray')

    ax = fig.add_subplot(3,3,i+1)
    ax.scatter(np.log2(fcs), -np.log10(q_values),
    c=colors)
    
    # 図をきれいに見せるためのあれこれ。好みの世界
    max_val = max(abs(np.nanmin(np.log2(fcs)[np.log2(fcs) != -np.inf])), max(np.log2(fcs)))
    ax.set_xlim([-max_val-1, max_val+1]) # -infがあるので。-inf = そのクラスターでは全員が0
    ax.set_ylim(ax.get_ylim())
    # 閾値に点線をつける
    ax.hlines([-np.log10(q_threshold)], -max_val-1, max_val+1, 'gray', 'dashed', linewidth=0.5, alpha=0.5)
    ax.vlines([np.log2(fc_threshold), np.log2(1/fc_threshold)], ax.get_ylim()[0], ax.get_ylim()[1], 'gray', 'dashed', linewidth=0.5, alpha=0.5)

    # ラベルとアノテーション
    ax.set_title(f'Cluster {i}')
    ax.set_xlabel('logFC')
    ax.set_ylabel('-log10(q-values)')
    for j, label in enumerate(cols):
        ax.annotate(label, (np.log2(fcs)[j], -np.log10(q_values)[j]), size=9)

fig.suptitle('Volcano plots')
fig.tight_layout()
plt.show()

  • オレンジの点 = 他のクラスターと比較して、対象クラスターで有意に大きな数値を持つ
  • 水色の点 = 他のクラスターと比較して、対象クラスターで有意に小さな数値を持つ

を意味します。今回は点が少ないのですが、何百何千も説明変数があるデータをこの手法で可視化すると噴火する火山のように見えることからVolcano plotと呼ばれます。
さて、ここから読み取れる特徴は

  • クラスター0は登板数 (G) を始め多くの指標が他より少ないが、防御率 (ERA) や四死球率 (BB9) が高い
  • クラスター1はセーブ数 (SV) が多い
  • クラスター2は完封 (SHO), 完投 (CG) が多い
  • クラスター3は負け数 (L) が多く、セーブが少ない
  • クラスター4はクラスター2と似た傾向で、イニング数 (IP) や勝ち (W) 負けが多い。一方で、完投・完封は飛び抜けて多くはない
  • クラスター5は目立って多い指標は無く平均的。少ない指標は完投・完封・セーブ・勝ち数・負け数
  • クラスター6は登板数と申告敬遠 (IBB) が多い
  • クラスター7は登板数を始め多くの指標が他より少ない

といったところでしょうか。

考察とアノテーション

ここまでの結果にドメイン知識、今回であれば野球の知識を導入して考察していきます。
まずそれぞれのクラスターについてですが、上の結果から

  • クラスター0: 1軍レベルにまだ満たない投手
  • クラスター1: 守護神
  • クラスター2: エース
  • クラスター3: ローテ入りはしていない先発投手
  • クラスター4: 先発ローテ
  • クラスター5: 中継ぎ
  • クラスター6: 頼れる中継ぎエース
  • クラスター7: 1軍定着を狙う投手

といったような考察ができます。申告敬遠が中継ぎエース(と思われるクラスター)で多いのは面白いですね。終盤の大事なところを担うからこそ、歩かせることも多いという事実を反映した結果になっているように思います。

UMAPの図の上にアノテーションを付けてみます。

cluster_centers = mapper.transform(kmeanModel.cluster_centers_)
cluster_names = ['1軍挑戦', '守護神', 'エース', '谷間の投手', '先発ローテ', '中継ぎ', '中継ぎエース', '1軍定着狙う']

plt.scatter(embedding[:, 0], embedding[:, 1],
    c=clusters, cmap='Accent', alpha=0.7)
plt.title('UMAP plot')
for i, label in enumerate(cluster_names):
    plt.annotate(label, cluster_centers[i], size=8)
plt.colorbar()
plt.show()

もう一つ、別の見せ方でも出してみましょう。

df_features_mm_clusters = df_features_mm.copy()
df_features_mm_clusters['cluster'] = clusters
x = []
y = []
targets = []
colors = []
for i, col in enumerate(cols):
    for j, cluster_name in enumerate(cluster_names):
        target_value = df_features_mm_clusters[df_features_mm_clusters['cluster']==j][col].mean()
        x.append(j)
        y.append(i)
        targets.append(np.exp(1+target_value*4)) # ここは図を見ながら調整
plt.scatter(x, y, s=targets, c=targets, cmap='plasma')
plt.xticks(list(range(k)), cluster_names, fontsize=8)
plt.yticks(list(range(len(cols))), cols, fontsize=8)
plt.show()


円が大きいほど、そのクラスターでその説明変数が大きいことを意味します。

こうしてみると、クラスターの名付けは特徴を捉えられていそうです。
実務的にはこの図を見てからクラスター名をつける場合も多いかもしれませんね。「クラスター1はSVタイプだ」「クラスター2はSHOタイプだ」といったような。
また、今回は全ての説明変数を表示していますが、説明変数が多い場合には何かしらの基準で選んだ(例えばVolcano plotで当たりをつけた)変数だけを選んで縦軸に表すことで、その変数がそれぞれのクラスターでどういう違いを持っているかを示すために使える表現法でもあります。

では最後に確認として、実際にクラスターに含まれる選手を見てみましょう。

df_with_clusters = df_dropped.copy()
df_with_clusters['cluster'] = clusters
df_target_cluster = df_with_clusters[df_with_clusters['cluster'] == 2]
df_target_cluster

クラスター2を取り出すと、たしかにエース格の投手が選ばれていますね。

さらにその先へ

教師なし学習の典型例であるクラスタリングの、その後に続く可視化、解析、そして解釈の一例をご紹介してきました。
データ集めから解釈まで通して書いた結果こんなに長い記事になってしまい恐縮です…ここまで読んでいただきありがとうございました。

今回ご紹介したところからさらに先の展開もいくつか提案してみます。

まず、こうした解析を通して、こういうデータが追加で欲しいという着想を得ることができます。
例えば、今回のクラスタリングの中で「中継ぎ」と名付けたクラスはかなり大雑把な括りとなっており、もう少し細かく分類したいなと感じます。
お気づきかもしれませんが、実は今回のデータにはホールドポイントが含まれていません。すなわち、優れた中継ぎ投手を簡単に差別化する指標が存在しないんですね。もしホールドポイントの記録もあれば、中継ぎエースも守護神や先発のように独立した島ができる可能性は高いだろうと思いますし、中継ぎの中でも中継ぎエースとまではいわないが接戦も任されるタイプ、敗戦処理、などとより細かい区分けができたかもしれません。
また、先発登板数のデータがあれば、今回「中継ぎ」にまとめられた中でも「先発もできるタイプ」「中継ぎ特化タイプ」を分けられるかもしれません。(ピンクの点が一部先発の島に入っているのは、こういう区分けの存在を示唆しているのだと思います)。
このように、現在手持ちのデータで一旦できあがったクラスター分類を見た上で、直感的に分かれそうだけどうまく分けられていないクラスターを分類できそうな説明変数を考察することは一つの有意義な展開かと思います。

また、データを３次元で見てみると新しいデータの切り口が得られることもあります。
例えば、今回の2次元のUMAPプロットに年齢の軸を加えれば「未来のエース・守護神候補」が見つかるかもしれません。

数年分のデータを合わせてクラスタリング・プロットするのも面白そうですね。その上で同じ選手がどのような変遷を辿っているかを追うことで、「典型的な選手のキャリアパターン」をいくつか見つけることができ、そこから　「この若手選手はパターン2に該当するから、来年はこの辺りに位置する成績を残しそう」みたいな予測ができるかもしれません。

そして最も強力なアプリケーションになると個人的に考えているのが疑似的な時系列変化の解析です。
今回のUMAPプロットでは右上から左下に向かって低い成績から高い成績へと、成績の連続的な変化が表現されていましたね(t-SNEでは左上から右下)。さらに先発タイプか抑えタイプかで枝分かれしていました。
これが、例えば購買履歴データの可視化だとすると、右上は最近会員登録したばかりの人で、そこからヘビーユーザーになると左下に進み、購入している商品の種類によって左か下かに分かれていく、というようにユーザー間の差異を擬似的な時系列変化として捉えることができるケースもあります。


(左: UMAP, 右: t-SNEのプロット上に疑似的な変化の方向を示す矢印を引いた)

こうしたデータの捉え方ができると 「このクラスターは今後このクラスターへと変遷していくだろう」「ということは、このクラスターに属する人は、今後こういうものを買う可能性が高い」「そこにリーチできる施策を打とう！」 といったかたちで実際のアクションに繋げていくこともできるかもしれません。

本稿では教師なしクラスタリングを一歩進めた応用の一例を紹介してきました。クラスタリングに連なる手法は単なるデータの分類に留まらず、データ間の連続的な性質の変遷、ともすれば擬似的な時系列関係をも表現できうる高いポテンシャルを持っています。
ぜひ皆様のお手元のデータでも本稿の解析を再現し、データに対する理解を深め楽しんでください！

関連記事

修正履歴

  • 2023-11-28
    • 誤字修正
  • 2023-12-04
    • 関連記事リンク追加
Aidemy Tech Blog
Aidemy Tech BlogPublication

株式会社アイデミーの技術ブログです。

Discussion