はじめに
以前の記事 1 で 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に近いほど、その点は自分のクラスタ内の他の点に近く、他のクラスタ内の点からよく離れていることを意味します。シルエットスコアは、すべての観測値のシルエット係数の平均なので、シルエットスコアが大きいほど任意のクラスタ内のポイントが互いに近く、クラスタ自体が互いにより分離しています。つまり、よくクラスタリングされていると言えます。
参照
ここではサラッと触りだけ触れましたが、詳細まで定義されていないと落ち着かない人はちゃんとした学術論文 2 を読むことをおススメします。
インポートとデータの読み込み
パッケージのインポート
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 ↩