4.0 MiB
4.0 MiB
Данные по инсультам¶
Выводим информацию о датасете:
In [352]:
import pandas as pd
df = pd.read_csv("..//..//static//csv//healthcare-dataset-stroke-data.csv")
df
Out[352]:
Атрибуты:
- id – уникальный идентификатор пациента;
- gender – пол пациента: может быть "Male" (мужчина), "Female" (женщина) или "Other" (другой);
- age – возраст пациента (в годах);
- hypertension – наличие гипертонии: 0 – гипертонии нет, 1 – гипертония есть;
- heart_disease – наличие сердечных заболеваний: 0 – заболеваний нет, 1 – заболевание присутствует;
- ever_married – семейный статус пациента: "No" (не состоял в браке) или "Yes" (состоял в браке);
- work_type – тип занятости пациента: "children" (дети), "Govt_job" (государственная служба), "Never_worked" (никогда не работал), "Private" (частная компания) или "Self-employed" (самозанятый);
- Residence_type – место проживания пациента: "Rural" (сельская местность) или "Urban" (город);
- avg_glucose_level – средний уровень глюкозы в крови (в ммоль/л);
- bmi – индекс массы тела пациента;
- smoking_status – статус курения пациента: "formerly smoked" (курил ранее), "never smoked" (никогда не курил), "smokes" (курит), "Unknown" (информация недоступна);
- stroke – факт наличия инсульта: 1 – пациент перенес инсульт, 0 – инсульта не было.
Бизнес-цель: кластеризация пациентов для выявления групп с схожими характеристиками здоровья и рисками инсульта. Что, к примеру, может использоваться для следующего:
- Определение групп пациентов для целенаправленных профилактических мероприятий.
- Оптимизация распределения медицинских ресурсов и создания индивидуализированных программ наблюдения.
Для начала избавимся от пустых значений:
In [353]:
# Количество пустых значений признаков
print(df.isnull().sum())
print()
# Есть ли пустые значения признаков
print(df.isnull().any())
print()
# Процент пустых значений признаков
for i in df.columns:
null_rate = df[i].isnull().sum() / len(df) * 100
if null_rate > 0:
print(f"{i} процент пустых значений: %{null_rate:.2f}")
In [354]:
# Замена значений
df["bmi"] = df["bmi"].fillna(df["bmi"].median())
Визуализация взаимосвязей¶
Для визуализации и выполнения задачи в целом были выбраны столбцы age, avg_glucose_level, bmi, hypertension.
In [355]:
from typing import Any, List
import matplotlib.pyplot as plt
def draw_data_2d(
df: pd.DataFrame,
col1: int,
col2: int,
y: List | None = None,
classes: List | None = None,
subplot: Any | None = None,
):
ax = None
if subplot is None:
_, ax = plt.subplots()
else:
ax = subplot
scatter = ax.scatter(df[df.columns[col1]], df[df.columns[col2]], c=y)
ax.set(xlabel=df.columns[col1], ylabel=df.columns[col2])
if classes is not None:
ax.legend(
scatter.legend_elements()[0], classes, loc="lower right", title="Classes"
)
In [356]:
columns = ['age', 'avg_glucose_level', 'bmi', 'hypertension']
df_reduced = df[columns]
plt.figure(figsize=(16, 12))
draw_data_2d(df_reduced, 0, 1, subplot=plt.subplot(2, 2, 1)) # age vs avg_glucose_level
draw_data_2d(df_reduced, 0, 2, subplot=plt.subplot(2, 2, 2)) # age vs bmi
draw_data_2d(df_reduced, 0, 3, subplot=plt.subplot(2, 2, 3)) # age vs hypertension
draw_data_2d(df_reduced, 1, 2, subplot=plt.subplot(2, 2, 4)) # avg_glucose_level vs bmi
Перед кластеризацией стандартизируем данные:
In [357]:
from sklearn.preprocessing import StandardScaler
columns_to_scale = df_reduced.drop(columns=["hypertension"]).columns
columns_to_keep = ["hypertension"]
scaler = StandardScaler()
data_scaled = scaler.fit_transform(df_reduced[columns_to_scale])
df_scaled = pd.DataFrame(data_scaled, columns=columns_to_scale, index=df_reduced.index)
df_scaled[columns_to_keep] = df_reduced[columns_to_keep]
Иерархическая агломеративная кластеризация¶
Также выведем дендрограмму
In [358]:
import numpy as np
from sklearn import cluster
from scipy.cluster import hierarchy
def run_agglomerative(
df: pd.DataFrame, num_clusters: int | None = 2
) -> cluster.AgglomerativeClustering:
agglomerative = cluster.AgglomerativeClustering(
n_clusters=num_clusters,
compute_distances=True,
)
return agglomerative.fit(df)
def get_linkage_matrix(model: cluster.AgglomerativeClustering) -> np.ndarray:
counts = np.zeros(model.children_.shape[0]) # type: ignore
n_samples = len(model.labels_)
for i, merge in enumerate(model.children_): # type: ignore
current_count = 0
for child_idx in merge:
if child_idx < n_samples:
current_count += 1
else:
current_count += counts[child_idx - n_samples]
counts[i] = current_count
return np.column_stack([model.children_, model.distances_, counts]).astype(float)
def draw_dendrogram(linkage_matrix: np.ndarray):
hierarchy.dendrogram(linkage_matrix, truncate_mode="level", p=3)
plt.xticks(fontsize=10, rotation=45)
plt.tight_layout()
In [359]:
tree = run_agglomerative(df_scaled)
linkage_matrix = get_linkage_matrix(tree)
draw_dendrogram(linkage_matrix)
Попробуем разделить данные на 2 больших кластера, поэтому зададим порог расстояния в 90 единиц.
И визуализируем сами результаты иерархической кластеризации, т.е. распределение кластеров:
In [360]:
result = hierarchy.fcluster(linkage_matrix, 90, criterion="distance")
y_names = ['1', '2']
plt.figure(figsize=(16, 12))
draw_data_2d(df_reduced, 0, 1, result, y_names, plt.subplot(2, 2, 1)) # age vs avg_glucose_level
draw_data_2d(df_reduced, 0, 2, result, y_names, plt.subplot(2, 2, 2)) # age vs bmi
draw_data_2d(df_reduced, 0, 3, result, y_names, plt.subplot(2, 2, 3)) # age vs hypertension
draw_data_2d(df_reduced, 1, 2, result, y_names, plt.subplot(2, 2, 4)) # avg_glucose_level vs bmi
KMeans (неиерархическая четкая кластеризация) для сравнения¶
In [361]:
from typing import Tuple
def print_cluster_result(
df: pd.DataFrame, clusters_num: int, labels: np.ndarray, separator: str = ", "
):
for cluster_id in range(clusters_num):
cluster_indices = np.where(labels == cluster_id)[0]
print(f"Cluster {cluster_id + 1} ({len(cluster_indices)}):")
rules = [str(df.index[idx]) for idx in cluster_indices]
print(separator.join(rules))
print("")
print("--------")
def run_kmeans(
df: pd.DataFrame, num_clusters: int, random_state: int
) -> Tuple[np.ndarray, np.ndarray]:
kmeans = cluster.KMeans(n_clusters=num_clusters, random_state=random_state)
labels = kmeans.fit_predict(df)
return labels, kmeans.cluster_centers_
In [362]:
random_state = 9
labels, centers = run_kmeans(df_scaled, 2, random_state) # также указываем 2 кластера
print_cluster_result(df_scaled, 2, labels)
display(centers)
Также визуализируем результаты:
In [363]:
def draw_cluster_results(
df: pd.DataFrame,
col1: int,
col2: int,
labels: np.ndarray,
cluster_centers: np.ndarray,
subplot: Any | None = None,
):
ax = None
if subplot is None:
ax = plt
else:
ax = subplot
centroids = cluster_centers
u_labels = np.unique(labels)
for i in u_labels:
ax.scatter(
df[labels == i][df.columns[col1]],
df[labels == i][df.columns[col2]],
label=i,
)
ax.scatter(centroids[:, col1], centroids[:, col2], s=80, color="k")
In [364]:
plt.figure(figsize=(16, 12))
draw_cluster_results(df_scaled, 0, 1, labels, centers, plt.subplot(2, 2, 1)) # age vs avg_glucose_level
draw_cluster_results(df_scaled, 0, 2, labels, centers, plt.subplot(2, 2, 2)) # age vs bmi
draw_cluster_results(df_scaled, 0, 3, labels, centers, plt.subplot(2, 2, 3)) # age vs hypertension
draw_cluster_results(df_scaled, 1, 2, labels, centers, plt.subplot(2, 2, 4)) # avg_glucose_level vs bmi
Теперь понизим размерность данных до двух компонент и еще раз осуществим неиерархическую кластеризацию¶
In [365]:
from sklearn.decomposition import PCA
pca_data = PCA(n_components=2).fit_transform(df_scaled)
pca_data
Out[365]:
Визуализация данных после понижения размерности:
In [366]:
plt.figure(figsize=(8, 6))
draw_data_2d(
pd.DataFrame({"Column1": pca_data[:, 0], "Column2": pca_data[:, 1]}),
0,
1
)
Визуализация результатов неиерархической кластеризации для двух кластеров с учетом понижения размерности:
In [367]:
from sklearn.cluster import KMeans
def fit_kmeans(
reduced_data: np.ndarray, num_clusters: int, random_state: int
) -> cluster.KMeans:
kmeans = cluster.KMeans(n_clusters=num_clusters, random_state=random_state)
kmeans.fit(reduced_data)
return kmeans
def draw_clusters(reduced_data: np.ndarray, kmeans: KMeans):
h = 0.02
x_min, x_max = reduced_data[:, 0].min() - 1, reduced_data[:, 0].max() + 1
y_min, y_max = reduced_data[:, 1].min() - 1, reduced_data[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = kmeans.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.figure(1)
plt.clf()
plt.imshow(
Z,
interpolation="nearest",
extent=(xx.min(), xx.max(), yy.min(), yy.max()),
cmap=plt.cm.Paired, # type: ignore
aspect="auto",
origin="lower",
)
plt.plot(reduced_data[:, 0], reduced_data[:, 1], "k.", markersize=2)
centroids = kmeans.cluster_centers_
plt.scatter(
centroids[:, 0],
centroids[:, 1],
marker="x",
s=169,
linewidths=3,
color="w",
zorder=10,
)
plt.title(
"K-means clustering (PCA-reduced data)\n"
"Centroids are marked with white cross"
)
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.xticks(())
plt.yticks(())
In [368]:
kmeans = fit_kmeans(pca_data, 2, random_state)
draw_clusters(pca_data, kmeans)
Анализ оценки инерции для метода локтя (метод оценки суммы квадратов расстояний)¶
In [369]:
import math
def _get_kmeans_range(
df: pd.DataFrame | np.ndarray, random_state: int
) -> Tuple[List, range]:
max_clusters = int(math.sqrt(len(df)))
clusters_range = range(2, max_clusters + 1)
kmeans_per_k = [
cluster.KMeans(n_clusters=k, random_state=random_state).fit(df)
for k in clusters_range
]
return kmeans_per_k, clusters_range
def get_clusters_inertia(df: pd.DataFrame, random_state: int) -> Tuple[List, range]:
kmeans_per_k, clusters_range = _get_kmeans_range(df, random_state)
return [model.inertia_ for model in kmeans_per_k], clusters_range
def _draw_cluster_scores(
data: List,
clusters_range: range,
score_name: str,
title: str,
):
plt.figure(figsize=(8, 5))
plt.plot(clusters_range, data, "bo-")
plt.xlabel("$k$", fontsize=8)
plt.ylabel(score_name, fontsize=8)
plt.title(title)
def draw_elbow_diagram(inertias: List, clusters_range: range):
_draw_cluster_scores(inertias, clusters_range, "Inertia", "The Elbow Diagram")
In [370]:
inertias, clusters_range = get_clusters_inertia(df_scaled, random_state)
display(clusters_range)
display(inertias)
draw_elbow_diagram(inertias, clusters_range)
На графике "Elbow Diagram" (метод локтя) оптимальное количество кластеров определяется точкой, где график начинает "сгибаться", то есть уменьшается прирост качества при добавлении новых кластеров (резкое снижение инерции становится более плавным).
На представленном выше варианте графика видно, что инерция резко падает от 2 до примерно 5 кластеров. После этого снижение инерции становится гораздо менее выраженным. Поэтому в этом случае не будет ошибкой выбрать число от 3 до 5, так как добавление большего количества кластеров уменьшает инерцию незначительно, что может не оправдывать усложнение модели.
Для выбранного же ранее варианта в 2 кластера (в процессе использования алгоритмов) инерция достаточно высокая, поэтому на таком значении, особенно если неизвестны особенности решаемой задачи, лучше не останавливаться.
Выбор количества кластеров на основе коэффициента силуэта¶
In [371]:
from sklearn.metrics import silhouette_score
def get_clusters_silhouette_scores(
df: pd.DataFrame, random_state: int
) -> Tuple[List, range]:
kmeans_per_k, clusters_range = _get_kmeans_range(df, random_state)
return [
float(silhouette_score(df, model.labels_)) for model in kmeans_per_k
], clusters_range
def draw_silhouettes_diagram(silhouette: List, clusters_range: range):
_draw_cluster_scores(
silhouette, clusters_range, "Silhouette score", "The Silhouette score"
)
In [372]:
silhouette_scores, clusters_range = get_clusters_silhouette_scores(df_scaled, random_state)
display(clusters_range)
display(silhouette_scores)
draw_silhouettes_diagram(silhouette_scores, clusters_range)
Коэффициент силуэта рассчитывается с использованием среднего расстояния внутри кластера (а) и среднего расстояния до ближайшего кластера (b) для каждого образца. Лучшее значение — 1, худшее — -1. Значения около 0 указывают на перекрывающиеся кластеры. Отрицательные значения обычно указывают на то, что образец был отнесен к неправильному кластеру.
На графике коэффициента силуэта оптимальное количество кластеров определяется пиком, где значение силуэта максимально, т.к. чем выше значение, тем лучше структура кластеров.
В данном случае из графика и предыдущего вывода списка оценок видно, что максимальное значение коэффициента силуэта наблюдается при 3 или 4 кластерах (около 0.36). Это говорит о том, что при таком количестве кластеров группы имеют наилучшее качество разделения.
Однако, если для задачи требуется большее количество кластеров, можно выбрать другое значение, где коэффициент силуэта все еще достаточно высокий (по сравнению с остальными вариантами). К примеру значения 5 или 7.
Для выбранного ранее варианта разделения на 2 кластера значение коэффициента силуэта равно примерно 0.2908, что указывает на то, что кластеры имеют нечеткую границу, а разделение данных является неоптимальным. Это может быть связано либо с недостаточным количеством кластеров, либо с особенностями самих данных, которые затрудняют их разделение на четко определенные группы.
Пример анализа силуэтов для разбиения от 2 до 12 кластеров¶
In [373]:
from typing import Dict
from sklearn.metrics import silhouette_samples
import matplotlib.cm as cm
def get_clusters_silhouettes(df: np.ndarray, random_state: int) -> Dict:
kmeans_per_k, _ = _get_kmeans_range(df, random_state)
clusters_silhouettes: Dict = {}
for model in kmeans_per_k:
silhouette_value = silhouette_score(df, model.labels_)
sample_silhouette_values = silhouette_samples(df, model.labels_)
clusters_silhouettes[model.n_clusters] = (
silhouette_value,
sample_silhouette_values,
model,
)
return clusters_silhouettes
def _draw_silhouette(
ax: Any,
reduced_data: np.ndarray,
n_clusters: int,
silhouette_avg: float,
sample_silhouette_values: List,
cluster_labels: List,
):
ax.set_xlim([-0.1, 1])
ax.set_ylim([0, len(reduced_data) + (n_clusters + 1) * 10])
y_lower = 10
for i in range(n_clusters):
ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]
ith_cluster_silhouette_values.sort()
size_cluster_i = ith_cluster_silhouette_values.shape[0]
y_upper = y_lower + size_cluster_i
color = cm.nipy_spectral(float(i) / n_clusters) # type: ignore
ax.fill_betweenx(
np.arange(y_lower, y_upper),
0,
ith_cluster_silhouette_values,
facecolor=color,
edgecolor=color,
alpha=0.7,
)
ax.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))
y_lower = y_upper + 10 # 10 for the 0 samples
ax.set_title("The silhouette plot for the various clusters.")
ax.set_xlabel("The silhouette coefficient values")
ax.set_ylabel("Cluster label")
ax.axvline(x=silhouette_avg, color="red", linestyle="--")
ax.set_yticks([])
ax.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
def _draw_cluster_data(
ax: Any,
reduced_data: np.ndarray,
n_clusters: int,
cluster_labels: np.ndarray,
cluster_centers: np.ndarray,
):
colors = cm.nipy_spectral(cluster_labels.astype(float) / n_clusters) # type: ignore
ax.scatter(
reduced_data[:, 0],
reduced_data[:, 1],
marker=".",
s=30,
lw=0,
alpha=0.7,
c=colors,
edgecolor="k",
)
ax.scatter(
cluster_centers[:, 0],
cluster_centers[:, 1],
marker="o",
c="white",
alpha=1,
s=200,
edgecolor="k",
)
for i, c in enumerate(cluster_centers):
ax.scatter(c[0], c[1], marker="$%d$" % i, alpha=1, s=50, edgecolor="k")
ax.set_title("The visualization of the clustered data.")
ax.set_xlabel("Feature space for the 1st feature")
ax.set_ylabel("Feature space for the 2nd feature")
def draw_silhouettes(reduced_data: np.ndarray, silhouettes: Dict):
for key, value in silhouettes.items():
if key > 12:
return
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(18, 7)
n_clusters = key
silhouette_avg = value[0]
sample_silhouette_values = value[1]
cluster_labels = value[2].labels_
cluster_centers = value[2].cluster_centers_
_draw_silhouette(
ax1,
reduced_data,
n_clusters,
silhouette_avg,
sample_silhouette_values,
cluster_labels,
)
_draw_cluster_data(
ax2,
reduced_data,
n_clusters,
cluster_labels,
cluster_centers,
)
plt.suptitle(
"Silhouette analysis for KMeans clustering on sample data with n_clusters = %d"
% n_clusters,
fontsize=14,
fontweight="bold",
)
In [374]:
silhouettes = get_clusters_silhouettes(pca_data, random_state)
draw_silhouettes(pca_data, silhouettes)
