Hari S. Sreedeth

Critically ill patients arriving in intensive care rarely fit neat templates. Even within a single ICU, patients differ widely in age, comorbidity burden, physiology and social context. This heterogeneity makes it hard to reason about prognosis, compare like with like, or decide which treatment strategies are most appropriate.

Unsupervised machine learning is an obvious tool for this problem, and a large body of ICU phenotyping work already exists. However, in practice many clustering pipelines are dominated by the blunt signal of multimorbidity: patients with more chronic disease, and who are older, simply get grouped together. The resulting “phenotypes” often mirror gradients in age and comorbidity count rather than revealing distinct, acute physiological patterns that could plausibly change clinical decisions.

The Multimorbidity-Anchored ICU Phenotypes (MAIP) project1 is my attempt to tackle that limitation in a principled, end-to-end way using the SUPPORT-II cohort (over 9,000-patient sample from the Study to Understand Prognoses, Preferences, Outcomes and Risks of Treatment)2. designed MAIP as a full pipeline that (i) explicitly separates chronic disease burden from acute physiology and (ii) carries the resulting clusters all the way through to a human-readable bedside tool.

MMSP (baseline): physiology-only clustering within multimorbidity strata

The MMSP baseline clusters acute physiology features only within each multimorbidity stratum, using the pre-scaled physiology matrix (P_view_scaled.csv). For a given stratum, near-constant physiology columns (within-stratum standard deviation ≤1e-8) were dropped to avoid degenerate PCA components. Dimensionality reduction was then performed with PCA (not FAMD): the number of components was chosen to retain approximately 80% cumulative explained variance, with a minimum of 2 components where feasible and never exceeding the number of retained physiology features.

Clustering was performed on the PCA embedding using PAM k-medoids on Euclidean distances computed in the PCA space. Candidate numbers of clusters were evaluated over K = 2…8. For each K, clustering stability was estimated using bootstrap subsampling (fraction 0.80, B=100) and Adjusted Rand Index (ARI) computed across pairs of bootstrap runs on the overlap of sampled individuals. Internal separation metrics (silhouette, Calinski–Harabasz, Davies–Bouldin) were computed from a representative fit at each K. Final K selection used a stability-first rule: retain K values within 20% of the best stability, then select the best among those by highest silhouette, then highest Calinski–Harabasz, then lowest Davies–Bouldin. The final MMSP partition for the chosen K was obtained using best-of-20 restarts to reduce sensitivity to initialization.

SNF-lite: multi-view similarity fusion within multimorbidity strata

SNF-lite integrates three feature views—comorbidity (C_view.csv), physiology (P_view.csv), and socio-contextual (S_view.csv)—within each multimorbidity stratum. Prior to similarity construction, low-signal variables in the comorbidity and socio-contextual views were removed: constants were dropped, and binary indicators with extreme prevalence (p<0.005 or p>0.995) were excluded. A safeguard check prevents a set of known outcome fields (e.g., death, hospdead, d.time, slos, totmcst) from being present in the C/S views used for fusion.

View-specific patient affinities were then constructed as follows. For comorbidity and socio-contextual views, a Gower-like similarity was computed using: (i) min–max scaling for numeric columns followed by average Manhattan distance, (ii) simple mismatch distance for multi-level categoricals, and (iii) an asymmetric binary treatment for {0,1} indicators using a Jaccard-style distance that ignores 0–0 matches, optionally weighted by prevalence (approximately p(1−p)) to down-weight extremely rare/common flags. For the physiology view, an RBF affinity was computed after standardization; the kernel used a local scaling parameter per individual based on the 7th nearest-neighbour distance (k_local=7) to reduce sensitivity to absolute scale differences across patients.

Fusion used a lightweight SNF routine. Each view affinity was KNN-masked to its top k=25 neighbours (include self), row-normalized to a row-stochastic matrix (diagonal set to 0), and iteratively updated for 10 iterations with mixing parameter α=0.5:

\[ P_v \leftarrow \alpha\,P^{\mathrm{KNN}}_v + (1-\alpha)\,\frac{1}{V-1}\sum_{u \ne v} P_u \]

followed by row-normalization. The final fused affinity was the mean across views and was symmetrized.

Clustering and K selection were performed on the fused affinity via spectral embedding followed by PAM k-medoids on Euclidean distances in the embedding space. Candidate K values were evaluated over K = 2…8, with the same bootstrap ARI stability procedure (B=100, subsample 0.80) and the same internal separation metrics. The same stability-first selection rule (20% tolerance from best stability; then silhouette/CH/DB tie-breaks) was applied. In addition, an eigengap diagnostic was computed from a KNN-sparsified normalized Laplacian (KNN=25) and recorded as a transparency aid. when the eigengap-suggested K achieved a silhouette within 95% of the selected K’s silhouette, the pipeline permitted preferring the eigengap K and refitting labels at that K (as implemented).

Internal clustering performance of the multimorbidity-stratified phenotyping (MMSP) solution and the Similarity Network Fusion (SNF-lite) solution across high, mid and low multimorbidity strata. For each stratum, MMSP metrics are reported at the rule-selected K from phase 1, and SNF-lite metrics are reported at the K chosen by the same stability-first selection rule (search range K = 3–8). Higher stability (ARI), silhouette and Calinski–Harabasz indices, and lower Davies–Bouldin indices, indicate better cluster separation. Across all strata, the SNF-lite solution with K = 3 shows consistently superior internal validity compared with MMSP; combined with prognostic and parsimony analyses, this motivated choosing SNF-lite as the primary phenotyping system.

Prognostic performance was evaluated using 5-fold cross-validated Harrell’s C-index for 6-month mortality within each multimorbidity stratum. The base model included age, sex, number of comorbidities (num.co) and APS, and was compared against versions augmented with MMSP or SNF-lite phenotypes. Across low (N = 4,183; 2,765 deaths), mid (N = 3,822; 2,648 deaths) and high multimorbidity (N = 1,100; 788 deaths), SNF-lite clusters consistently delivered larger gains in cross-validated C-index than MMSP clusters, with the strongest uplift in the low- and high-burden groups.

Five-fold cross-validated Harrell’s C-index for 6-month mortality by multimorbidity stratum. The base model includes age, sex, number of comorbidities (num.co) and the Acute Physiology Score (APS). MMSP and SNF-lite phenotypes are added as categorical predictors; we report mean C-index and mean ΔC versus the base model across folds.

MMSP and SNF-lite differ not only in how they discover structure, but also in how compact and interpretable the resulting phenotypes are. Below, I summarise the number and balance of clusters and the complexity of the surrogate decision trees in each multimorbidity stratum.

SNF-lite recovered three clinically coherent phenotypes within each multimorbidity stratum. In the high-burden cohort, clusters separate patients with diabetic multi-organ failure from chronic cardio–pulmonary and cardiometabolic cancer profiles. Mid-multimorbidity phenotypes span older COPD-dominant chronic respiratory disease, multi-organ failure with malignancy and cirrhosis, and compensated heart failure with solid tumours. Even among low-multimorbidity patients, SNF-lite distinguishes low-acuity solid tumour cases, neurologic catastrophe/coma, and acute multi-organ failure in otherwise low-comorbidity patients. Together, these nine phenotypes map onto familiar ICU patterns while making the spectrum of chronic burden and acute severity explicit.

Within each multimorbidity stratum, I estimated adjusted associations between SNF-lite phenotype membership and 365-day all-cause mortality using Cox proportional hazards models, with cluster as the exposure and adjustment for age (standardised), sex, multimorbidity count (standardised), and acute physiology severity (APS, standardised). Cluster membership remained independently prognostic after adjustment. In the high multimorbidity cohort (N=1100), both High_MM_1 and High_MM_2 had lower mortality hazards than the reference High_MM_0 (High_MM_1 vs High_MM_0: HR 0.57 [0.39–0.82]; High_MM_2 vs High_MM_0: HR 0.51 [0.37–0.70]). In the mid stratum (N=3822), Mid_MM_1 showed the highest adjusted risk relative to Mid_MM_0 (HR 1.87 [1.56–2.23]), with a smaller but still significant elevation for Mid_MM_2 (HR 1.30 [1.11–1.53]). In the low stratum (N=4183), the coma-dominant phenotype Low_MM_1 had markedly higher adjusted risk than Low_MM_0 (HR 2.36 [1.93–2.89]), while Low_MM_2 had lower adjusted hazard than Low_MM_0 (HR 0.59 [0.49–0.73]), consistent with long-term risk being driven by underlying disease composition in addition to acute severity. APS was the strongest and most consistent covariate across strata (HR ~1.4–1.6 per 1 SD increase).

Proportional hazards diagnostics based on Schoenfeld residuals indicated departures from strict PH assumptions even after truncation to 365 days; hazard ratios are therefore interpreted as average effects over the 1-year horizon and are presented alongside unadjusted KM curves as supporting evidence.

In addition to the adjusted Cox models reported above, I ran supporting outcome analyses to contextualise these findings, including unadjusted Kaplan–Meier curves with logrank tests, and non-parametric comparisons of ICU length of stay and total medical costs across phenotypes. These supplementary outputs are provided in the Supplementary section (figures and tables), alongside cluster profile summaries used for phenotype interpretation.

Beyond mortality, I assessed whether the identified phenotypes corresponded to distinct patterns of healthcare resource consumption. Using Kruskal-Wallis tests, I compared distributions of hospital length of stay (LOS) and total medical costs across phenotypes within each stratum to quantify the operational impact of these clinical subtypes.

High Stratum: Significant heterogeneity was observed in both LOS ($P=0.015$) and costs ($P=0.006$). The high-risk reference phenotype High_MM_0 was associated with the highest resource burden (Median LOS: 9.0 days; Median Cost: $23,412), while the lower-risk phenotypes High_MM_1 and High_MM_2 showed reduced utilization (Median LOS: 7.0 days for both).

Mid Stratum: Resource utilization varied markedly ($P < 0.001$). Consistent with its high mortality risk, Mid_MM_1 accrued the highest median costs ($22,042) and longest stay (9.0 days), nearly double the costs of the lower-risk reference Mid_MM_0 ($12,189).

Low Stratum: The starkest contrast mirrored the survival findings ($P < 0.001$). The high-risk Low_MM_1 phenotype had a median LOS of 8.0 days and costs of $21,114, whereas the protective Low_MM_2 phenotype demonstrated the lowest resource utilization in the entire cohort (Median LOS: 4.0 days; Median Cost: $9,207).