はじめに

以前の記事 ¹ で seaborn の penguins のデータセットを使って重回帰分析を行いました。そこで今回は同じデータセットに対して、K-平均法を使用して、このデータをグループ化し、ペンギンに関するパターン分析を行います。

K-平均法とは

実際のコードに入る前にK-平均法とその評価手法（エルボー法とシルエット分析）について軽く触れておきます。K-平均法は教師なしの分割アルゴリズムです。データ群をクラスタリング（分類）する、つまりラベルのないデータをグループまたはクラスターに整理するために使用されます。 $k$ はモデル内の重心の数、つまりクラスターの数を表します。例えば、今回であればペンギンの種類が3種類と分かっているのであれば、$k = 3$としても良いでしょう。

エルボー法とシルエット分析

前段で「ペンギンの種類が3種類と分かっているのであれば、$k = 3$ 」と書きましたが、本当に（データセットを統計的に分析して）クラスターの数は3でよいのでしょうか？実は「アデリーペンギンとヒゲペンギンが同一種族（とデータ上みなすことが出来て）として、ジェンツーペンギンが完全に別種族する方がデータ統計上、良いのではないか？」ということがあります。K-平均法は教師なし学習なので、モデルを比較するためのラベル付きデータはありません。そこで、クラスタリングのモデルを評価し、どのグルーピングが意味を持つかを決定するために、エルボー法とシルエット・スコアを使用します。

エルボー法

エルボー法はクラスター分析においてデータセット内のクラスターの数を決定する際に使用される手法です。異なるモデルの慣性モーメントを視覚的に比較するために折れ線グラフを使用します。K-平均法では、異なるkの値の比較をして行われます。
例：

この場合ですと5, 6, 7あたりが肘のように見えるので、$k = 5, 6, 7$が候補となります。

シルエット分析 (Silhouette Analysis)

エルボー法は主観的であり、信頼性が低いと考えられています。そこでシルエット分析です。シルエット分析とは、シルエットスコア（Silhouette Score）を比較する手法です。シルエットスコアはオブジェクトが自身のクラスター（凝集度）と他のクラスター（分離度）とでどの程度類似しているかを示す尺度です。そして、シルエットスコアは、モデル内のすべての観測値のシルエット係数の平均として定義されます。シルエット係数は以下の式で定義されます：

$a = $インスタンスと同じクラスタ内の他のインスタンスとの平均距離
$b = $インスタンスから、最も近い他のクラスタ（つまり、インスタンスが割り当てられているクラスタを除く）の各インスタンスまでの平均距離
$max(a,b) = $ a, bのいずれか大きい方の値。

Silhouette Coefficient= \frac{b − a}{max(a,b)}

シルエット係数は、-1 から +1 の範囲で値を取り、その値が+1に近いほど、その点は自分のクラスタ内の他の点に近く、他のクラスタ内の点からよく離れていることを意味します。シルエットスコアは、すべての観測値のシルエット係数の平均なので、シルエットスコアが大きいほど任意のクラスタ内のポイントが互いに近く、クラスタ自体が互いにより分離しています。つまり、よくクラスタリングされていると言えます。

参照

ここではサラッと触りだけ触れましたが、詳細まで定義されていないと落ち着かない人はちゃんとした学術論文 ² を読むことをおススメします。

インポートとデータの読み込み

パッケージのインポート

K-means, silhouette_score, StandardScaler に必要なライブラリをインポートしていきます。

# データ分析に必要な標準ライブラリ群
import numpy as np
import pandas as pd
import seaborn as sns

# モデルと評価に必要なライブラリ群
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import OneHotEncoder

# ビジュアル化に必要なライブラリ群
import matplotlib.pyplot as plt
import seaborn as sns

# データの確認
pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 1000)
print(penguins.head(n = 10))
penguins = sns.load_dataset("penguins", cache=False)

出力例

species     island  bill_length_mm  bill_depth_mm  flipper_length_mm   body_mass_g     sex  
0  Adelie  Torgersen            39.1           18.7              181.0        3750.0    Male  
1  Adelie  Torgersen            39.5           17.4              186.0        3800.0  Female  
2  Adelie  Torgersen            40.3           18.0              195.0        3250.0  Female  
3  Adelie  Torgersen             NaN            NaN                NaN           NaN     NaN  
4  Adelie  Torgersen            36.7           19.3              193.0        3450.0  Female  
5  Adelie  Torgersen            39.3           20.6              190.0        3650.0    Male  
6  Adelie  Torgersen            38.9           17.8              181.0        3625.0  Female  
7  Adelie  Torgersen            39.2           19.6              195.0        4675.0    Male  
8  Adelie  Torgersen            34.1           18.1              193.0        3475.0     NaN  
9  Adelie  Torgersen            42.0           20.2              190.0        4250.0     NaN

探索的データ解析

データ探索

データセットに含まれるカテゴリ変数の種類の数を確認しておきます。

print(penguins['species'].unique())

print(penguins['island'].unique())

print(penguins['sex'].unique())

出力例

# species
['Adelie' 'Chinstrap' 'Gentoo']

# island
['Torgersen' 'Biscoe' 'Dream']

# sex
['Male' 'Female' nan]

MALE と Male のような大文字、小文字の違いによるカテゴリ変数の登録は無さそうです。ですが、nanとあるように、欠損値がありそうです。

補足

もし、大文字、小文字が混在しているようであれば

penguins_subset['sex'] = penguins_subset['sex'].str.upper()

のようにして変換を実施。

欠損値の確認

K-平均法では欠損値はないことが前提となります。よって、データレコードに欠損値がないか確認します。

penguins['species'].value_counts(dropna = False)

出力例

species               0
island                0
bill_length_mm        2
bill_depth_mm         2
flipper_length_mm     2
body_mass_g           2
sex                  11
dtype: int64

欠損値のあるレコードがあるので、それを削除します。

penguins['species'].value_counts(dropna = False)

出力例

species              0
island               0
bill_length_mm       0
bill_depth_mm        0
flipper_length_mm    0
body_mass_g          0
sex                  0
dtype: int64

欠損値がなくなったことを確認できました。

データのエンコード

K-means では、クラスタリングに数値列が必要です。カテゴリ列sexを数値に変換します。species列はクラスタリングアルゴリズムの特徴として使用されていないため、変換する必要はありません。

# 「性別」列をカテゴリから数値に変換します。
penguins_subset['sex'] = OneHotEncoder(drop='first').fit_transform(penguins_subset[['sex']]).toarray()

データセットからカテゴリ列islandを削除します。同じ種のペンギンが性別に基づいて異なる身体的特徴を示すかどうかを確認しようとしているので、islandは不要と判断します。

penguins_subset = penguins_subset.drop(['island'], axis=1)

それと、species列は数値ではないことに留意します。また、species列を削除は行いません。これは後でクラスターを理解するのに有用と考えられるからです。

特徴をスケーリングする

K-means では、類似性の尺度として観測間の距離を使用するため、モデリングの前にデータをスケーリングすることが重要です。 scikit-learn の StandardScaler 関数などのサードパーティツールを使用します。 StandardScaler は、各ポイント $X_i$ を、その特徴の平均観測値を減算し、標準偏差で割ることでスケーリングします:

x_{scaled} = \frac{x_i - mean(x)}{σ}

これにより、すべての変数の平均が 0、分散/標準偏差が 1 になります。

注: 種の列は特徴ではないため、スケーリングする必要はありません。

まず、「種」列を除くすべての特徴を DataFrame X にコピーします。

X = penguins_subset.drop(['species'], axis=1)

StandardScaler を使用して X の特徴をスケーリングし、スケーリングされたデータを新しい変数 X_scaled に割り当てます。

X_scaled = StandardScaler().fit_transform(X)

データモデリング

K-平均法

次に、K-平均法を当てはめて、K のさまざまな値に対する慣性モーメントを評価します。データにクラスターがいくつ存在するかがわからない可能性があるため、まず K-平均法を当てはめて、 $k$ のさまざまな値に対する慣性値を調べます。これを行うには、 num_clusters と x_vals (X_scaled) を受け取り、各 $k$ 値の慣性モーメントのリストを返す kmeans_inertia という関数を作成します。

関数内で K 平均法を使用する場合は、random_state を設定します。ここでは 42 としました。こうすることで、結果の再現性を担保することができるようになります。

# K-平均法を適合し、k のさまざまな値に対する慣性を評価します。
num_clusters = [i for i in range(2,11)]

def kmeans_inertia(num_clusters, x_vals):
    
    inertia = []
    for num in num_clusters:
        kms = KMeans(n_clusters=num, random_state=42)
        kms.fit(x_vals)
        inertia.append(kms.inertia_)
    
    return inertia

kmeans_inertia 関数を使用して、k = 2 から 10 までの慣性のリストを取得します。

inertia = kmeans_inertia(num_clusters, X_scaled)
print(inertia)

出力例

[885.622414365225, 578.8284278107235, 477.22956735281934, 284.5464837898288, 218.0515629783023, 201.51585567445818, 196.58908645819542, 180.65751845628574, 170.15703589755495]

エルボー法

次に、 num_clusters と inertia の関係を示す折れ線グラフを作成します。この関係を視覚化するには、seaborn または matplotlib を使用します。

# 折れ線グラフの作成
plot = sns.lineplot(x=num_clusters, y=inertia, marker = 'o')
plot.set_xlabel("Number of clusters");
plot.set_ylabel("Inertia");

プロットのエルボー（肘）は通常、曲線の中で最も急な角度または曲がりを持つ部分です。線が急激に曲がっている場所、つまり $k = 6$ であるように見えます。

シルエットスコア

前段のエルボー（肘）が $k = 6$ であることを判断するためにシルエットスコアを調べます。そこで、シルエットスコアを評価する関数、各 $k$ 値のスコアのリストを返す関数を作成します。引数は整数値（クラスター数）のリストとそのデータの配列になります。

def kmeans_sil(num_clusters, x_vals):
    sil_score = []
    for num in num_clusters:
        kms = KMeans(n_clusters=num, random_state=42)
        kms.fit(x_vals)
        sil_score.append(silhouette_score(x_vals, kms.labels_))

    return sil_score

sil_score = kmeans_sil(num_clusters, X_scaled)
print(sil_score)

出力例

[0.44398088353055243, 0.45101024097188375, 0.4489699212061028, 0.5199985748608681, 0.5223086008347771, 0.47386350642293157, 0.4715443426463867, 0.4160561489496056, 0.4183063433691049]

$k = 6$ の時、シルエットスコアが最大値の0.5223086008347771を取るので、クラスター数は6とするのが良いようです。
注：配列の戻り値の1番目は $k = 2$ の時のシルエットスコア。

念のためにグラフ化してみます。

plot = sns.lineplot(x=num_clusters, y=sil_score, marker = 'o')
plot.set_xlabel("# of clusters");
plot.set_ylabel("Silhouette Score");

plt.show()

最適なk 値

最適な k 値を決定するために、データセットを6つのクラスターモデルにあてはめます。その上で適合モデルの一意のラベルを出力し、結果を確認します。

kmeans6 = KMeans(n_clusters=6, random_state=42)
kmeans6.fit(X_scaled)

print('Unique labels:', np.unique(kmeans6.labels_))

出力例

Unique labels: [0 1 2 3 4 5]

次に、penguins_subset にクラスターの割り当てを示す新しい列クラスターを作成します。これにより各クラスターのラベルの意味を理解し、クラスタリングが意味をなすかどうかを判断することができるようになります。
注：penguins_subset を使用するのは、スケーリングされていないデータを解釈する方が簡単な場合ことが多いためです。

# 新しいカラム `cluster` の追加
penguins_subset['cluster'] = cluster_6.labels_
penguins_subset.head()

出力例

  species  bill_length_mm  bill_depth_mm  flipper_length_mm  body_mass_g  sex cluster  
0  Adelie            39.1           18.7              181.0       3750.0  1.0       5  
1  Adelie            39.5           17.4              186.0       3800.0  0.0       2  
2  Adelie            40.3           18.0              195.0       3250.0  0.0       2  
3  Adelie            36.7           19.3              193.0       3450.0  0.0       2  
4  Adelie            39.3           20.6              190.0       3650.0  1.0       5

groupby を使用して、cluster が species によって区別できるかどうかを確認します。

print(penguins_subset.groupby(by=['cluster', 'species']).size())

出力例

cluster  species  
0        Chinstrap    32
1        Gentoo       58
2        Adelie       73
         Chinstrap     2
3        Adelie        2
         Chinstrap    34
4        Gentoo       61
5        Adelie       71
dtype: int64

上記の通り、groupby の結果は、各clusterをspeciesで区別できることを示しています。しかしながら、同じspeciesでありながら、群が分かれてしまっていることが確認できます。その前に、一度これらの結果を視覚化して、分かりやすくしておきます。

penguins_subset.groupby(by=['cluster', 'species']).size().plot.bar(title='Clusters differentiated by species',
                                                                   figsize=(6, 5),
                                                                   ylabel='Size',
                                                                   xlabel='(Cluster, Species)');

plt.show()

やはり同じspeciesでありながら別群になっていることがグラフからも読み取れます。そこで、更にgroupby を使用して、各clusterがspeciesとsexで区別できるかどうかを確認します。

print(penguins_subset.groupby(by=['cluster','species', 'sex']).size().sort_values(ascending = False))

出力例

cluster  species    sex
2        Adelie     0.0    73
5        Adelie     1.0    71
4        Gentoo     1.0    61
1        Gentoo     0.0    58
3        Chinstrap  1.0    34
0        Chinstrap  0.0    32
2        Chinstrap  0.0     2
3        Adelie     1.0     2
dtype: int64

やはりspeciesとsexによって群が分かれていることが確認できそうです。念のためにグラフ化します。

penguins_subset.groupby(by=['cluster','species','sex']).size().unstack(level = 'species', fill_value=0).plot.bar(title='Clusters differentiated by species and sex',
                                                                                                                      figsize=(6, 5),
                                                                                                                      ylabel='Size',
                                                                                                                      xlabel='(Cluster, Sex)')
plt.legend(bbox_to_anchor=(1.3, 1.0))

plt.show()

cluster 1 と 3 は1種、単一の性別ではありませんでしたが、groupbyのアルゴリズムによって大体speciesとsexで区別された群が生成されることが示されています。

6つの群があることは理にかなっています。なぜなら、この調査データには、3つの種とそれぞれに性別の違いがあることが分かっているためです。 (男女の性別 × 3つの異なる種 = 6 つの群)。

K-平均法の3Dクラスタリングの可視化

最後の分類された6つの群がどれぐらいの距離、どのように群化されているかを視覚的に分かるようにしておきます。ここでは3Dクラスタリングを用いて可視化します。

fig = plt.figure(figsize=(12,10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X[:, 0], X[:, 1], X[:, 4], c=labels)
ax.scatter(centers[:, 0], centers[:, 1], centers[:, 4], marker='x', s=200, c='red', alpha=1)

ax.legend(labels=['K-means','cluster centers','kd'])

plt.title("Kmeans clustering")
plt.show()

結論

上記より、K-平均のモデルを使用することで、事前にパターンを知らなくてもデータ内の隠れた構造を明らかにすることが出来ることが分かりました。めでたしめでたし。

以上。

Pythonを使って重回帰分析を行う - Qiita https://qiita.com/tatsu_sekine/items/4c50a0b0e0257761cc2f ↩
Silhouettes: A graphical aid to the interpretation and validation of cluster analysis - ScienceDirect https://www.sciencedirect.com/science/article/pii/0377042787901257 ↩