29 · Unüberwachtes Lernen und Phänotypisierung von Patient:innen

import pandas as pd
import numpy as np
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

labs   = pd.read_csv("https://schradern.github.io/data-science-coach/data/labor.csv")

NUMERIC = ["alter", "sofa_score", "crp_mg_l", "bmi",
           "leukozyten_g_l", "kreatinin_mg_dl", "laktat_mmol_l"]

df = pd.read_csv("https://schradern.github.io/data-science-coach/data/kohorte.csv").merge(labs, on="patient_id", how="left")
X_raw = df[NUMERIC].values

prep = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("scale",  StandardScaler()),
])
X = prep.fit_transform(X_raw)

library(tidyverse); library(recipes)

labs   <- read_csv("https://schradern.github.io/data-science-coach/data/labor.csv", show_col_types = FALSE)
df <- read_csv("https://schradern.github.io/data-science-coach/data/kohorte.csv", show_col_types = FALSE) |> left_join(labs, by = "patient_id")
numeric_cols <- c("alter","sofa_score","crp_mg_l","bmi",
                  "leukozyten_g_l","kreatinin_mg_dl","laktat_mmol_l")

rec <- recipe(~ ., data = df[numeric_cols]) |>
  step_impute_median(all_predictors()) |>
  step_normalize(all_predictors())
X_scaled <- rec |> prep() |> bake(new_data = NULL) |> as.matrix()

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

inertias, silhouettes = [], []
K_RANGE = range(2, 9)

for k in K_RANGE:
    km = KMeans(n_clusters=k, random_state=SEED, n_init=10)
    labels = km.fit_predict(X)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X, labels))

# In dieser Kohorte ist der Silhouettenkoeffizient bei k=2 am höchsten (0,130),
# knapp vor k=7 (0,123) und k=5 (0,122); k=3 liegt mit 0,115 sogar darunter.
# Alle Werte liegen zwischen 0,10 und 0,13 -> schwache, überlappende
# Clusterstruktur, kein k trennt die Kohorte sauber. Trotzdem wird unten
# bewusst mit k=3 weitergearbeitet (drei Schweregrad-Stufen), eine klinisch
# begründete, keine rein statistisch optimale Wahl (siehe Stolperstein).
best_k = 3
km_final = KMeans(n_clusters=best_k, random_state=SEED, n_init=10).fit(X)
df["cluster"] = km_final.labels_

library(factoextra)
fviz_nbclust(X_scaled, kmeans, method = "silhouette")
km <- kmeans(X_scaled, centers = 3, nstart = 10, iter.max = 100)
df$cluster <- km$cluster

from scipy.cluster.hierarchy import linkage, dendrogram, fcluster
import matplotlib.pyplot as plt

Z = linkage(X, method="ward")
# Cut at k=3 to match k-means
labels_hier = fcluster(Z, t=3, criterion="maxclust")

dist_mat <- dist(X_scaled, method = "euclidean")
hc <- hclust(dist_mat, method = "ward.D2")
plot(hc, labels = FALSE, main = "Dendrogramm (Ward)")
abline(h = 15, col = "red", lty = 2)   # Schnittlinie bei k=3
df$cluster_hier <- cutree(hc, k = 3)

from sklearn.decomposition import PCA

# PCA: linear, preserves global variance
pca = PCA(n_components=2, random_state=SEED)
X_pca = pca.fit_transform(X)
explained = pca.explained_variance_ratio_
print(f"PC1+PC2 explain {explained.sum():.1%} of variance")

# UMAP: non-linear, preserves local structure
try:
    import umap
    reducer = umap.UMAP(n_components=2, random_state=SEED)
    X_2d = reducer.fit_transform(X)
    method_label = "UMAP"
except (ImportError, TypeError) as exc:
    # ImportError: umap-learn not installed. TypeError: an installed
    # umap-learn can be incompatible with the installed scikit-learn version
    # (a real version-skew failure, not just a missing package).
    from sklearn.manifold import TSNE
    X_2d = TSNE(n_components=2, random_state=SEED, perplexity=30).fit_transform(X)
    method_label = "t-SNE (UMAP nicht verfügbar/inkompatibel)"
    print(f"Note: UMAP failed ({exc}), using {method_label}")

library(umap)  # install.packages("umap") falls nicht vorhanden
pca_res  <- prcomp(X_scaled, scale. = FALSE)
X_pca_r  <- pca_res$x[, 1:2]
# UMAP (Fallback: PCA nutzen, falls Paket fehlt)
tryCatch({
  umap_res <- umap(X_scaled)
  X_umap   <- umap_res$layout
}, error = function(e) {
  message("umap nicht verfügbar, nutze PCA-2D als Fallback")
  X_umap <<- X_pca_r
})

import pandas as pd

profile = df.groupby("cluster")[NUMERIC + ["verweildauer_tage", "verstorben_30d"]].mean().round(2)
print(profile.T)

# Qualitative labels — assigned AFTER reading the profile, by ranking
# clusters on mean SOFA (illness severity). k-Means cluster IDs (0, 1, 2, …)
# are arbitrary and NOT guaranteed to come out in severity order (they don't,
# for this cohort/seed: the highest-SOFA cluster is ID 1, not ID 2), so a
# fixed {1: "mild", 2: "moderate", 3: "severe"}-style mapping keyed on raw
# IDs would silently mislabel the phenotypes.
severity_order = profile["sofa_score"].sort_values().index.tolist()
phenotype_names = dict(zip(severity_order, ["Leichter Verlauf", "Mittelschwer", "Hochrisiko"]))
df["phenotype"] = df["cluster"].map(phenotype_names)

profile <- df |>
  group_by(cluster) |>
  summarise(across(c(alter, sofa_score, crp_mg_l, laktat_mmol_l, verweildauer_tage,
                     verstorben_30d), mean, na.rm = TRUE))
print(profile)