クラスタリング（k-means / 階層的 / DBSCAN / GMM）

クラスタリングは、正解ラベルを使わずに、似た症例をグループ化する方法です。脳卒中後の歩行パターン、認知機能の下位項目、嚥下運動の時系列特徴量などを整理する時に役立ちます。一方で、得られた群をそのまま疾患サブタイプと呼ぶには注意が必要です。本稿では、k-means・階層的クラスタリング・DBSCAN・GMM を、臨床研究での読み方を意識して整理します^[6]。

// 01 · LEARN OUTCOMESこの記事で学ぶこと

• k-means・階層的・DBSCAN・GMM の違いを説明できる。
• クラスタ数 k・距離尺度・標準化の重要性を理解できる。
• Silhouette などの指標を、臨床的解釈と分けて読める。
• クラスタを疾患サブタイプと誤読しないための確認点を整理できる。

// 02 · CONCLUSIONまず結論

// 03 · FIGURE図で理解するクラスタリング

k-means は、各クラスタの重心に近い症例をまとめます。丸い群が複数あるデータでは扱いやすい一方で、三日月型や細長い群では苦手です。歩行特徴量のように相関が強い変数では、事前に PCA で構造を見ることもあります。

階層的クラスタリングは、似ている症例や群を順に結合します。デンドログラムを使うと、どこで切ると何群になるかを視覚的に確認できます。N が大きすぎない臨床データでは、説明しやすい選択肢です。DBSCAN は密度の高い部分をクラスタとして扱い、孤立した症例を外れ値として扱える点が特徴で、異常検知 の発想と共通します。GMM は複数のガウス分布が混ざったものとしてデータを見て、各症例に「クラスタ A らしさ 0.7」のような確率を付けられるため、境界が連続的な現象に読みやすいです。

// 04 · CLINICALリハ研究での使いどころ

// 05 · THEORY4 手法の考え方

k-means: 重心に近い症例を集める

k-means は、クラスタ数 k を先に決めます^[1]。各症例を最も近い重心に割り当て、重心を更新する操作を繰り返します。目的は、クラスタ内のばらつきを小さくすることです。

k-means objective

min  Σ_k Σ_{x_i in C_k} || x_i - mu_k ||^2

C_k  : cluster k
mu_k : centroid of cluster k

k-means は、距離に強く依存します。FIM・年齢・歩行速度・加速度特徴量をそのまま混ぜると、単位の大きい変数が結果を支配します。そのため、StandardScaler などで標準化します。

階層的クラスタリング: 近いものから結合する

階層的クラスタリングは、最初は各症例を 1 つの群として扱います。そこから、近い群を順に結合します。Ward 法では、結合後の分散増加が小さい組み合わせを選びます^[2]。N が数百程度のリハデータでは、デンドログラムが説明に使いやすいです。

Ward linkage

At each step:
merge the pair of clusters
that gives the smallest increase
in within-cluster variance.

DBSCAN: 密度の高い場所をクラスタにする

DBSCAN は、半径 eps の中に min_samples 以上の点があるかを見ます^[3]。条件を満たす点をコアポイントと呼び、コアポイントが連結している領域をクラスタにします。どのクラスタにも属さない点は noise として扱われます。

DBSCAN definitions

eps         : neighborhood radius
min_samples : minimum number of samples

core point:
number of neighbors within eps >= min_samples

noise point:
not reachable from any core point

DBSCAN は k を指定しません。その代わり、eps と min_samples に敏感です。eps が小さいと多くが noise になり、大きいと別の群がつながります。結果を図と臨床情報の両方で確認します。

GMM: 混合分布として群を見る

GMM(Gaussian Mixture Model) は^[4]、データが複数のガウス分布から来ていると仮定します。k-means と違い、症例を硬く 1 群に割り当てません。各症例に、各クラスタに属する確率を計算します。

GMM idea

p(x) = Σ_k pi_k * N(x | mu_k, Sigma_k)

pi_k    : mixture weight
mu_k    : mean vector
Sigma_k : covariance matrix

EM algorithm:
E-step: estimate cluster probabilities
M-step: update pi, mu, Sigma

GMM は、群の境界が連続的な場合に便利です。ただし、ガウス分布という仮定を置きます。特徴量の分布が大きく歪む場合は、変換やロバストな前処理を検討します。

評価指標は「正解」ではなく補助線

Silhouette は、同じ群内の近さと、他の群からの遠さを見ます^[5]。Davies-Bouldin は低いほど分離がよい傾向、Calinski-Harabasz は高いほど分離がよい傾向です。既知ラベルがある場合は ARI も使えます。

Common metrics

Silhouette score:
range -1 to 1
larger is often better

Davies-Bouldin index:
smaller is often better

Calinski-Harabasz score:
larger is often better

ARI:
agreement with external labels

しかし、指標が高い群が臨床的に意味を持つとは限りません。年齢や施設差だけで分かれている可能性もあります。クラスタごとの FIM・歩行速度・疾患重症度・施設構成を必ず確認します。

手法選択は「データの形」と「研究目的」から考える

どのクラスタリング手法を使うかは、解析前にある程度の方針を置きます。歩行速度や TUG・BBS などの連続変数で群が丸く分かれると考えるなら k-means が候補です。重症例から軽症例まで連続的に広がる場合は、GMM の方が読みやすいことがあります。

症例数が N=100〜300 程度で、研究会や論文の図として説明しやすさを重視するなら、階層的クラスタリングが向きます。デンドログラムはどの症例同士が近いかを追いやすいからです。ただし、症例数が大きくなると図が読みにくくなります。

歩行センサーや筋電図では、計測エラーや体動アーチファクトが混ざります。その場合、DBSCAN は外れ値候補を noise として扱えます。ただし、noise が臨床的異常を意味するとは限りません。機器エラー・装着位置・測定手順の違いも確認します。

重要なのは、最初から 1 つの手法だけに固定しないことです。k-means・階層的・GMM で大きく似た構造が見えるなら、結果の安定性を説明しやすくなります。逆に、手法ごとに全く違う群が出る場合は、そのデータに強いクラスタ構造がない可能性もあります。

クラスタ名は最後に付ける

「低機能群」「認知低下群」「歩行不安定群」のような名前は、読者にとって分かりやすいです。しかし名前を先に置くと、解析結果をその名前に合わせて読んでしまいます。まずは Cluster 1・Cluster 2 のように中立的に示します。

その後、各群の特徴を表にします。FIM motor・FIM cognition・歩行速度・TUG・BBS・MMSE・年齢・発症からの日数・施設を並べます。その特徴が一貫している場合に限り、Discussion で臨床的な呼び方を提案します。

さらに、クラスタ名には幅を持たせます。「歩行不安定群」と断定するより、「歩行不安定性が高い可能性のある群」と書く方が安全です。教師なし学習では目的変数がないため、群の意味は解析後に解釈するものだからです。

// 06 · IMPLEMENTATION · PYTHONPython 実装例

以下は scikit-learn を用いた教育用の仮想データ例です^[7]。実データでは、個人情報保護・倫理審査・施設ルール・研究計画書・利用規約を確認します。クラスタリングでも、前処理を全データに fit してから解析すると、評価や再現性の見積もりが楽観的になることがあります（09·02 データリーケージ 参照）。

import numpy as np
import pandas as pd

from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN
from sklearn.mixture import GaussianMixture
from sklearn.metrics import silhouette_score, davies_bouldin_score
from sklearn.metrics import calinski_harabasz_score
from sklearn.model_selection import train_test_split

RANDOM_STATE = 42

# Educational example only.
# Confirm ethics, privacy, and institutional rules.
df = pd.read_csv("example_rehab_dataset.csv")

feature_cols = [
    "age", "fim_motor_admission", "fim_cognition_admission",
    "gait_speed", "tug_seconds", "bbs_score",
    "step_length_asymmetry", "cadence",
    "accel_variability", "grip_strength"
]

X = df[feature_cols].copy()

num_pipe = Pipeline(steps=[
    ("imputer", SimpleImputer(strategy="median")),
    ("scaler", StandardScaler())
])

preprocess = ColumnTransformer(
    transformers=[("num", num_pipe, feature_cols)]
)

X_dev, X_holdout = train_test_split(
    X, test_size=0.25, random_state=RANDOM_STATE
)

# Fit preprocessing on the development set only.
X_dev_z = preprocess.fit_transform(X_dev)
X_holdout_z = preprocess.transform(X_holdout)

def cluster_scores(X_array, labels):
    labels = np.asarray(labels)
    valid = labels != -1
    labels_valid = labels[valid]
    X_valid = X_array[valid]
    n_clusters = len(set(labels_valid))
    if n_clusters < 2:
        return np.nan, np.nan, np.nan, n_clusters
    return (
        silhouette_score(X_valid, labels_valid),
        davies_bouldin_score(X_valid, labels_valid),
        calinski_harabasz_score(X_valid, labels_valid),
        n_clusters
    )

records = []

for k in range(2, 7):
    km = KMeans(n_clusters=k, n_init=20, random_state=RANDOM_STATE)
    labels = km.fit_predict(X_dev_z)
    sil, db, ch, nc = cluster_scores(X_dev_z, labels)
    records.append(["kmeans", f"k={k}", nc, sil, db, ch, np.nan])

for k in range(2, 7):
    ward = AgglomerativeClustering(n_clusters=k, linkage="ward")
    labels = ward.fit_predict(X_dev_z)
    sil, db, ch, nc = cluster_scores(X_dev_z, labels)
    records.append(["ward", f"k={k}", nc, sil, db, ch, np.nan])

for eps in [0.5, 0.8, 1.1]:
    dbs = DBSCAN(eps=eps, min_samples=10)
    labels = dbs.fit_predict(X_dev_z)
    sil, db, ch, nc = cluster_scores(X_dev_z, labels)
    records.append(["dbscan", f"eps={eps}", nc, sil, db, ch, np.nan])

for k in range(2, 7):
    gmm = GaussianMixture(
        n_components=k, covariance_type="full",
        random_state=RANDOM_STATE
    )
    labels = gmm.fit_predict(X_dev_z)
    sil, db, ch, nc = cluster_scores(X_dev_z, labels)
    records.append(["gmm", f"k={k}", nc, sil, db, ch, gmm.bic(X_dev_z)])

summary = pd.DataFrame(records, columns=[
    "method", "parameter", "n_clusters",
    "silhouette", "davies_bouldin",
    "calinski_harabasz", "bic"
])
print(summary.sort_values(["method", "silhouette"], ascending=[True, False]))

# PCA is used only for visualization.
pca = PCA(n_components=2, random_state=RANDOM_STATE)
X_dev_pca = pca.fit_transform(X_dev_z)

final_model = KMeans(n_clusters=3, n_init=20, random_state=RANDOM_STATE)
dev_cluster = final_model.fit_predict(X_dev_z)
holdout_cluster = final_model.predict(X_holdout_z)

profile = X_dev.copy()
profile["cluster"] = dev_cluster
print(profile.groupby("cluster").median(numeric_only=True))

# Review clusters by FIM, gait speed, BBS, TUG, age, and site.
# Do not name clusters before checking clinical plausibility.

この実装では、まず開発用データだけで前処理を fit しています。holdout は、最終的なクラスタ割り当ての再現性確認に使います。k-means は predict できますが、階層的クラスタリングや DBSCAN は新規データへの割り当てに追加工夫が必要です。

研究計画段階では、クラスタリングを行う理由も明確にします。「予後が異なる群を探す」のか、「運動パターンを整理する」のか、「外れ値を見つける」のかで、手法と評価が変わります。目的が曖昧なまま複数手法を試すと、後から都合のよい図だけを選んだように見えます。解析前に主解析と感度分析を分けておくと、説明の透明性が高まり、再解析時にも有用です。

論文用には、クラスタごとの要約表を作ります。年齢・性別・疾患重症度・FIM・歩行速度・BBS・TUG・施設・測定機器を群ごとに並べ、そこで初めて、群名を付けるかどうかを検討します。PCA や UMAP を図示に使う場合も、「見た目が分かれている」だけで判断しません。可視化パラメータ・乱数 seed・前処理・欠測補完の方法を Methods に書きます（03·14 次元削減 も参照）。

感度分析も重要です。k-means では、k を 2〜6 で変えた時に主要な群の解釈がどの程度変わるかを見ます。DBSCAN では eps と min_samples を少し変え、GMM では BIC だけでなく各群の症例数も確認します。症例数が極端に少ない群は、外れ値の集まりかもしれません。さらに bootstrap で症例を再抽出し、クラスタ構造が安定するかを見る方法もあります。同じような群が繰り返し出るなら結果は説明しやすく、再抽出のたびに群が大きく変わる場合は、「サブタイプ発見」よりも「探索的な傾向」として書く方が自然です。

// 07 · MYTHSよくある誤解

誤解: クラスタリングで得られた群はサブタイプ

クラスタは、入力した特徴量と手法で見えた群です。臨床サブタイプと呼ぶには、別データでの再現性・臨床転帰との関係・専門家による解釈が必要です。

誤解: k-means なら N が小さくても使いやすい

N が小さいと、初期値や少数の外れ値に結果が揺れます。また、球状クラスタという前提が合わないと、見かけ上きれいでも意味が不安定になります。

誤解: Silhouette が最大の k が正解

Silhouette は補助指標です。k が 2 で最大でも、臨床的には 3 群の方が説明しやすい場合があります。指標と臨床解釈を分けて記載します。

誤解: DBSCAN は k を決めなくてよいので楽

DBSCAN は eps と min_samples に敏感です。スケールや密度が変わると結果も変わります。パラメータを変えた感度分析が重要です。

// 08 · WRITING論文 Methods に書くこと

クラスタリング研究では、モデル名だけでは不十分です。どの特徴量を使ったか・どの時点のデータか・欠測値をどう扱ったか・標準化したかを明記します。入院時 FIM と退院時 FIM を混ぜると、時間情報が混ざり解釈が崩れます。

査読では、クラスタ数の決定根拠を問われやすいです。Elbow・Silhouette・BIC・デンドログラムのどれを使ったか、なぜ最終的にその k を採用したかを説明します。「最も見やすかったから」では弱くなります。

もう 1 つの論点は、クラスタの臨床的意味です。たとえば歩行特徴量で 3 群が得られても、年齢だけで分かれている可能性があります。施設や測定機器で分かれている可能性もあります。その場合、病態のサブタイプとは言いにくいです。

Results では、クラスタリング指標だけでなく、群ごとの FIM・歩行速度・TUG・BBS・MMSE・疾患重症度を表にします。Discussion では、探索的解析として扱うこと、因果効果や治療反応性を示すものではないことを明記します。

予後予測を目的とする場合は、クラスタを説明変数として使う段階と、予測モデルを評価する段階を分けます。クラスタ作成を全データで先に行い、そのラベルを使って交差検証すると、情報漏れに近い設計になることがあります。データリーケージ の観点で確認します。

// 09 · CHECKLIST実装前チェックリスト

• 01クラスタリングの目的が「探索」か「予測前処理」かを分けた
• 02FIM・歩行速度・TUG などの測定時点を揃えた
• 03StandardScaler を Pipeline 内で使った
• 04k・eps・min_samples・covariance_type を記録した
• 05複数指標で k の妥当性を見た
• 06クラスタごとの臨床背景表を作った
• 07施設差や測定機器差で分かれていないか確認した
• 08別データ・holdout・bootstrap で再現性を見た

// 10 · QUIZ理解度チェック

1. Q1k-means を使う前に特に重要な前処理はどれですか。
   • 目的変数を二値化する
   • 特徴量を標準化する
   • すべての特徴量をカテゴリ化する
   • 欠測値をそのまま残す
   SHOW ANSWER

   B. k-means は距離に依存するため、標準化が重要です。
2. Q2DBSCAN の特徴として適切なのはどれですか。
   • クラスタ数 k を必ず指定する
   • 外れ値を noise として扱える
   • すべての点を必ずクラスタに入れる
   • ガウス分布を仮定する
   SHOW ANSWER

   B. DBSCAN は密度ベースで、孤立した点を noise として扱えます。
3. Q3GMM の説明として適切なのはどれですか。
   • 各症例を確率的にクラスタへ割り当てる
   • デンドログラムを必ず使う
   • eps を主なパラメータにする
   • 距離尺度を使わないため標準化は無関係
   SHOW ANSWER

   A. GMM は各クラスタに属する確率を計算できます。
4. Q4クラスタリング結果を論文で示す時の注意として適切なのはどれですか。
   • Silhouette が高ければ臨床サブタイプと断定できる
   • クラスタ名は解析前に決めてよい
   • 群ごとの臨床背景を確認する
   • 施設差は教師なし学習では問題にならない
   SHOW ANSWER

   C. クラスタの意味は、臨床背景や再現性と合わせて確認します。

// 11 · FAQよくある質問

クラスタリングで得られた群はサブタイプですか？

そのままではサブタイプの証明になりません。クラスタは、入力した特徴量・距離尺度・手法・パラメータで「見えた群」です。臨床サブタイプとして提示するには、別データでの再現性、臨床転帰との関係、専門家による解釈が必要です。

k-means・階層的・DBSCAN・GMM はどう使い分けますか？

球状クラスタが複数想定されるなら k-means、N が数百で説明しやすさを重視するなら階層的（Ward + デンドログラム）、外れ値・密度の濃淡を扱いたいなら DBSCAN、群の境界が連続的で確率的解釈をしたいなら GMM が候補です。複数手法で結果が概ね一致する場合に解釈の安定性が増します。

クラスタ数 k はどう決めますか？

Elbow（クラスタ内ばらつきの減少）・Silhouette・Davies-Bouldin・Calinski-Harabasz・GMM では BIC などの指標を併用します。ただし、数値だけで決めず、臨床的に説明しやすい k かどうか、別データで再現するかも合わせて判断します。

// 12 · RELATED関連記事

// REF参考文献

1. MacQueen JB. Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability 1967;1:281-297.
2. Ward JH Jr. Hierarchical grouping to optimize an objective function. Journal of the American Statistical Association 1963;58(301):236-244.
3. Ester M, Kriegel HP, Sander J, Xu X. A density-based algorithm for discovering clusters in large spatial databases with noise. Proceedings of the Second International Conference on Knowledge Discovery and Data Mining 1996:226-231.
4. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society Series B 1977;39(1):1-38.
5. Rousseeuw PJ. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics 1987;20:53-65.
6. Hastie T, Tibshirani R, Friedman J. The Elements of Statistical Learning. 2nd ed. Springer; 2009.
7. Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 2011;12:2825-2830.