Developing the Total Health Profile, a Generalizable Unified Set of Multimorbidity Risk Scores Derived From Machine Learning for Broad Patient Populations: Retrospective Cohort Study

Introduction

Multimorbidity risk scores, which factor in the presence of several chronic conditions, can provide insights into morbidity and mortality [1,2]. In general, the number of co-occurring medical conditions is associated with increased adverse medical outcomes [3-5] and increased use of medical services [6]. This is particularly true for older individuals since the number of co-occurring medical conditions will increase with age [7]. Various approaches to quantifying multimorbidity have been used, including simply counting the number of conditions [8], while more complex tools have also been developed, such as the Charlson Comorbidity Index (CCI) [9] and the Elixhauser Comorbidity Index (ECI) [10]. However, these scores were developed decades ago and are limited to diagnostic information or certain populations, which may limit their utility. A systematic review [11] of multimorbidity scores examined 35 major multimorbidity scores, which could be subclassified by the information they used (eg, prescription data, diagnostic data, self-reported quality of life) and the outcomes they recorded (all-cause mortality, emergency room admissions, and hospital admissions). Patients with multimorbidities are cared for in general practice and by specialists [8,12,13] who use disease-specific risk scores and guidelines. Condition-specific risk scores, such as the Framingham Risk Score [14] for coronary heart disease, can help identify specific interventions to benefit patients and provide actionable information to guide tests and medications. One potential issue with the use of these tools is that sometimes they do not address the overall health care priorities of the patient due to their narrow focus [8,15,16].

Multimorbidity scores tend to only use one type of clinical data, such as diagnoses, prescriptions, or procedures, and rarely integrate them. As a result, they may be missing vital information and relationships in patient information. Although newer methods, such as probabilistic phenotyping [17], may alleviate these concerns, while remaining scalable, these methods are still highly experimental, with a wide variety of methods and little consensus on which ones are most trustable for real-world settings [18]. Using multiple data sets, feature types, and methodological explorations could provide a more comprehensive and robust estimate of multimorbidity risk. Currently, no multimorbidity scores exist that produce granular and overall risk profiles irrespective of age and sex; are derived from a large, representative population of patients; and integrate multiple clinical data sets, including diagnoses, prescriptions, lab results, and procedures using machine learning (ML; building upon previous ML-based strategies and recommendations for multimorbidity analysis by Hassaine et al [19]). Such scores could help health care providers engage in patient-centered care and prescribing, reduce polypharmacy, and guide deprescribing when used together with traditional risk scores and guidelines.

To address this need, we sought to create the total health profile (THP), a set of ML-derived measures of an individual's comprehensive clinical risk. The THP presents clinical risk in five separate models (referred to as “component scores”), producing granular, multimorbid risk scores specific to cardiovascular (“heart score”), respiratory (“lung score”), neuropsychiatric (“neuro score”), renal (“kidney score”), and gastrointestinal (“digestive score”) conditions. These organ systems reflect those involved in the top 10 sources of mortality in the United States [20] and serve to complement existing disease-specific risk scores. We also included, as a sixth score, the total health score (THS), a single view of a patient's overall health across all five of the aforementioned organ systems, which can be compared to existing pure multimorbidity risk scores. Each of these six scores was independently modeled using electronic health record (EHR) data consisting of demographic information, diagnosis codes, lab results, prescriptions, and medical procedural data and required, otherwise, no patient behavior or familial history data. For the unified risk endpoint of all six of the scoring models, we used inpatient (IP) hospital visits. As such, each score’s estimate of clinical risk represents the likelihood of an IP hospital visit over the next 24 months, attributable to the score's clinical category (eg, lung, heart). After training, testing, and calibrating the THS and the five organ system component scores, we analyzed the metrics and generalizability of each score across populations. We also conducted ablation tests of several feature groups to assess their importance in the final set of models. Finally, we discussed the clinical applicability of the THP, limitations of the study, and future work.

Methods

Study Design and Patient Inclusion Criteria

This retrospective cohort study used lab measurements and an administrative claims database of 52 million patients provided by a US health care insurance company. Patients were enrolled in a mixture of commercial, Medicare, Medicaid, and exchange plans. Our study design involved training on retrospective data from a certain time window and assessing performance via a follow-up time window. The retrospective observation period, or the time period in which model features were collected, was defined as January 1, 2016, through December 1, 2017, and the follow-up period, or the time period in which the model labels were collected, as January 1, 2018, through December 31, 2019.

All patient data were de-identified. Patients selected for inclusion had at least one medical claim in each year of the data collection and follow-up periods and had a known sex, birthdate, and zip code. These inclusion criteria resulted in 14 million patients, from which 1 million patients were randomly selected for analysis using PySpark, resulting in 992,868 patients due to the approximation methods used by PySpark. This patient sample was split into training (n=794,294) and testing (n=198,574) groups corresponding to an 80:20 ratio. Diagnosis codes (International Classification of Diseases, Tenth Revision [ICD-10]), medical procedure codes (Current Procedural Terminology [CPT]), lab data, demographics (social determinants of health [SDoH], patient gender/age), and prescription data (defined by General Product Identifier [GPI] codes) were used from patients who met the selection criteria. Our study, in total, used 88 ICD-10 codes and 30 chronic conditions (derived from ICD-10 codes specified by the ECI), 16 lab types, 764 GPI codes representing 4 GPI prefixes, 14 CPT codes, and 17 demographic markers.

Data Processing

Using the data compiled for the 992,868 patients, we extracted a set of features corresponding to chronic diagnoses, acute diagnoses, IP hospital visits, prescriptions, sociodemographic information, and lab results/physical exam measurements for feature extraction and modeling. A description of all features gathered during the data collection period follows next.

Demographic information was extracted from the United States Census American Community Survey for 2017 at the zip code level. This information included population, household count, and race and ethnicity percentages for that zip code (eg, African American, non-Hispanic White, Hispanic, Asian, Native American), sex percentages per zip code, and economic indicators, including the mean and median income. Demographic data also included the age and sex of the patients. Chronic disease diagnoses were counted as the presence of a chronic disease, while acute diagnoses were counted as the number of those diagnoses in the study period, summed over the component. For instance, 3 atrial fibrillation codes and 2 acute heart failure codes during the 2-year data collection period would have resulted in the number of acute heart diagnoses being 5. Medical procedure features were counted as the count of IP CPT codes that occurred during the data collection period, with otherwise identical score-specific inclusion criteria to the IP hospital labels (discussed in the Model Outcome Labels section). Four groups of prescriptions were included, assigned using the first two digits of the GPI code and indicated by binary presence: antihypertensives, hypoglycemics, lipid-lowering medications, and antithrombotic agents. In all, 16 labs or vitals were included in the study, each one being a numerical feature. If there were multiple results of the same lab data/vitals collected during the data collection period, only the most recent measurement was included. In total, our feature set and labels used the following set of clinical features: diagnoses (88 ICD-10 codes and 30 chronic conditions), labs (16 types), prescriptions (764 GPI codes, representing 4 GPI prefixes), procedures (14 CPT codes), and demographics (15 SDoH and 2 individual patient characteristics).

Except for demographic features, lab values, and vital signs, all input features were filtered on a model-by-model basis to include only score-relevant data (ie, the heart score would be modeled using only physician-curated features related to cardiovascular health). For IP hospital visit features collected during the data collection period, only score-specific IP visit counts were included (ie, the heart score would have as input the number of heart-related IP visitations during the data collection period, not the lung-related, and so on). The set of input features used over all component score models were used as input for the THS model, with an exception for chronic diagnosis features.

Model Outcome Labels

All component score labels were a binary indicator referring to whether a patient had an IP visit within the follow-up period, given that they also had acute or chronic diagnoses within 12 months prior to the IP visit and within 7 days after the IP visit, establishing both a history of that condition and that the IP visit was (likely) related to that condition. These diagnoses would be specific to each component, given by the corresponding ECI comorbidities and ICD-10 codes. For example, a possible positive label for the lung scoring model could be an IP hospital stay CPT code on June 2, 2019, a diagnosis code corresponding to pneumonia 3 months prior to it, and a diagnosis code corresponding to chronic pulmonary disease 2 days after it. The THS label is simply the combination of all the component score labels; if a patient has any positive component score label, the THS label would be positive as well.

Modeling Procedures and Baselines

All scores were calculated using a gradient-boosted tree classifier, with default hyperparameters, using the Scikit-Learn Python 3.6 package (version 0.24.1). Using demographics, diagnoses, lab values, procedures, and prescription data as input and IP visits as binary labels, separate models were trained for each score and subsequently calibrated using an isotonic regression with 3-fold cross-validation over the training set. Discriminative results from the models were obtained using the optimal threshold point of the training set (given by the threshold that yielded the smallest difference between the true-positive rate and the false-positive-rate) and applied to the testing set. All missing values were mean-imputed, and all input features for each model were mean-normalized using the training data.

We had multiple baselines: a logistic regression model with default hyperparameters using the statsmodel package (version 0.12.0) with otherwise identical feature sets, and a comparison of the performance of the THS to commonly used scores of a similar nature, specifically the CCI and the ECI, in predicting the hospitalization endpoint. We also conducted multiple ablation tests on the feature groups: a set of gradient-boosted tree classifiers, all with procedural data, but having only one set of either lab, SDoH, prescription, or diagnosis information. For the baseline gradient-boosted comparison with the combined feature model for patient subgroups, CI calculations were generated using 100 bootstrap iterations of 10% of the given demographic. The patient subgroups analyzed were patients with two or more of any comorbidity and one or more prescriptions of hypertensive, hyperglycemic, lipid-lowering, or antithrombotic medications.

Radial Plots

Radial plots were generated using three patients who were closest to each of the centroids of a fitted, randomly initialized K-means model, with a K value of 3. The K-means algorithm used the Scikit-Learn Python 3.6 package (version 0.24.1), and the plots themselves were generated using Plotly.

Model Discrimination and Generalization/Sensitivity Analyses

Models were assessed on three levels: discriminative performance, calibration, and generalizability in performance across different demographics. To assess the discriminative performance of each model in the THP, we calculated the area under the receiver operating curve (AUROC), sensitivity, and specificity using Scikit-Learn for the testing set of 198,574 patients. We also plotted the AUROC for all scores on the testing set (Figure 1). All CIs for the discriminative metrics were generated using 500 bootstrap samples of 20,000 from the testing data set. We selected the AUROC as our primary metric because it represents a comprehensive measure of the true-positive-rate and false-positive-rate trade-off without needing an optimal threshold point. Since our outcomes exhibited strong class imbalance, which may have led to overly optimistic AUROC values, we used sensitivity and specificity as secondary model measures. To assess calibration performance, we created calibration plots using Scikit-Learn, graphing predicted probability versus positive label percentage across 10 uniform probability bins. We assessed calibration, as measured by calibration plots, as the primary measurement of clinical utility as it gives a clear idea of how these scores can be used to identify sick patients, avoid alarm fatigue, and be interpreted as a probabilistic likelihood. To assess the generalization performance of each model, we studied how the performance and scores of the models vary across age, zip code income, and sex categories. We plotted how the AUROC varies among age groups (decade age groupings), median income groups (low, medium and high), and binary gender (male or female). Additionally, we computed statistical significance Z-tests for AUROC pairwise differences between all groups within each category. Due to the lower sample sizes of the groups, CI calculations were generated using 100 bootstrap samples of 10% of the given demographic of the testing group. The THS and component scores were then analyzed by plotting the distribution of scores as a function of age and disease burden (measured by the presence of pre-existing comorbidities during the observation period). Specifically, we looked at the distributions of the THS and the component scores among various age groups for patients with zero comorbidities found during the data collection period and patients with at least one ECI comorbidity related to the given component found during the data collection period.

‎View this figure

Physician-Guided Feature Selection and Curation

To select features to incorporate into the THS, a physician-guided curation method was incorporated, which involved selecting conditions, reviewing clinical practice guidelines for important conditions, and identifying clinical measures, tests, and pharmacological interventions in those guidelines. No statistical feature selection techniques were used, as those offer improved accuracy mainly in cases with relatively small training data sets or models that are sensitive to unsuspected feature correlations [21].

An overview of the manually guided feature selection process is described below:

1. Selection of disease categories/subscores: The top causes of death across the United States were reviewed from the 2018 mortality statistics from the National Center for Health Statistics, and five main categories that the causes of death could be classified into were identified: cardiac (heart), respiratory (lung), neuropsychiatric (neuro), gastrointestinal (digestive), and renal (kidney).
2. We then obtained the leading causes of IP conditions for 2011-2013 using the Agency for Health care Research and Quality (AHRQ) Healthcare Cost and Utilization Project (HCUP) database. We cross-referenced the top 30 codes for each year with the 5 categories we developed in step 1. We did not include codes corresponding to obstetrical conditions, complications related to birth and delivery, multisystem malignancy, and musculoskeletal conditions.
3. Using the AHRQ HCUP codes, we extracted the ICD-10 codes and selected additional medically related conditions (eg, selecting the ICD-10 codes for ischemic stroke codes, in addition to ICD-10 codes for hemorrhagic stroke).
4. To obtain corresponding prescription, lab, and procedure codes, we then reviewed clinical practice guidelines for the identified conditions (eg, stroke, chronic obstructive pulmonary disease) from the United States Preventive Services Task Force [18,22]. These guidelines were then reviewed and lab data and procedures corresponding to diagnosis and management were extracted by a physician. These were then manually mapped to the corresponding drug (GPI codes), procedure (CPT codes), and lab (Logical Observation Identifiers Names and Codes) codes.

MI-CLAIM Checklist

This work meets the Minimum Information about Clinical Artificial Intelligence Modeling (MI-CLAIM) requirements for sharing design, data/optimization, model performance, model examination, and level of reproducibility [23].

Role of the Funding Source

The funding source collected the raw data.

Ethics Approval

The analysis presented here is not to be characterized as human subject research. We are presenting the results of an analysis conducted for a health plan’s health care operations in accordance with the Health Insurance Portability and Accountability Act of 1996 (HIPAA). Only aggregated results of the business analysis are provided, and no individually identifiable information (protected health information or otherwise) was used in the development of this presentation.

Results

Patient Cohort

In all, 992,868 patients matched the inclusion criteria (Table 1). The majority were female (n=560,165, 56.4%), with 432,703 (43.6%) males, concordant with the 2019 census results [20]; the median age (41 years) was slightly higher than the national average; and the number of comorbidities tended to increase with age, consistent with previous findings [7]. The mean patient age was 39.07 years (95% CI 39.02-39.12), the mean number of comorbidities was 1.71 (95% CI 1.71-1.72), and the percentage of patients with IP visits was 1.65% (females 1.7%, males 1.7%). IP visits were also positively correlated with age for each organ system (Spearman correlation=0.314), which is concordant with previous studies [21]. The sample-weighted summary of zip-code-level demographics had a median income of $66,829, and the racial averages by zip code were 73.8% White, 5.9% Asian, and 11.9% African American, closely matching census averages (Table 1).

Overall Model Performance

All models outperformed the logistic regression baseline and were well specified with AUROCs of 0.83 for the THS, 0.89 for heart, 0.86 for lung, 0.84 for neuro, 0.90 for kidney, and 0.83 for digestive functions (Figure 1). All six models were well calibrated. Additional metrics (sensitivity and specificity) can be found in Table 2.

One benefit of the THP is that it is personalized to the patient to allow for nuanced interpretation based on the affected organ system, in addition to robust predictive performance. Figure 2 demonstrates an illustrative radial plot example of three patients who were around the same age (50-60 years old) and had the same rough THS (>0.8). Unlike grouped scoring systems, the THP enables the clinician to understand the personalized drivers for that score, thereby enabling clinical decisions that are specific to the individual patient. The score of patient 1 was driven primarily by neuro and heart issues, while the score of patient 2 was affected by kidney, neuro, and heart diseases, and the score of patient 3 was mostly affected by heart, lung, digestive, and kidney maladies.

‎View this figure

Important Model Features

We obtained the most important features of the THS and of each component model. For the THS, the biggest drivers of the model were the use of prior IP hospital visits (0.41) and age (0.15), followed by an acute heart-related diagnosis (0.11), uncomplicated hypertension (0.06), and acute neurological conditions (0.02); see Figure 3. The key features for each component model were directly relevant elements. For example, the most important features of the kidney model included a diagnosis of renal failure, as well as age, last recorded measurement of estimated glomerular filtration rate (eGFR) and hemoglobin A_1c (HbA_1c), and any acute kidney-related diagnosis. Generally, age and a prior history of hospitalization for issues relating to the organ system in question were the most important features across most component models.

‎View this figure

Generalizability Across Subgroups

We examined differences in performance for the THS as well as each component score for population subgroups based on sex, median neighborhood income, and age (Figures 4A, 4B, and 4C, respectively). Statistical AUROC comparisons [24] with Bonferroni corrections revealed that there were no significant differences in model performance on the basis of sex or income. There were occasional statistical differences in performance on the basis of age, primarily related to the 80-90-year-old population, which had significantly fewer data points than any of the other age groups (see Table 1).

‎View this figure

Ablation Tests

We conducted ablation tests of several feature groups to assess the need for them in the final set of models. All ablation test models used procedural data as input, as we assumed access to this information is a given due to the outcome prediction being procedural as well. The combined feature set outperformed the lab, SDoH, and prescription-only models. The combined feature set also outperformed the diagnosis-only model, though not significantly. This statistically nonsignificant outperformance was similarly observed in multiple patient subgroups focused on at-risk patients.

Performance Comparisons to Charlson and Elixhauser Risk Scores

In addition to a model baseline (logistic regression), we also compared the performance of the THS to commonly used scores of a similar nature, specifically the CCI and the ECI. Both the CCI and the ECI generate a risk score based on different weights associated with certain diseases based on ICD-10 codes. The THS (AUROC=0.82) outperformed both the CCI (AUROC=0.74) and the ECI (AUROC=0.65); see Table 3. We further compared performance across subgroups. Of the baselines, the CCI is a consistently better predictor than the ECI. Across subgroups, the THS also consistently outperformed the score baselines for every age, sex, and income bracket. The improvement was perhaps most pronounced for ages between 20 and 50, as well as for individuals living in lower-income neighborhoods.

Discussion

Principal Findings

There is a continued need for an updated clinical score that profiles patients based on multimorbidities that are equitable across populations and nuanced enough to facilitate precision medicine. To facilitate clinical decision making across patient populations, we created an automated, generalizable, integrated, multimorbidity risk profile across several clinical domains. The THP is composed of cardiovascular, respiratory, neuropsychiatric, renal, and gastrointestinal clinical risk subprofiles, as well as a sixth score, the THS, representing the overall combinatorial risk of the five organ-specific scores. We followed ML best practices to train six integrated models on large-scale EHR data with the forecasted probability of a risk endpoint, organ-specific IP hospital visits, over a 2-year window as the target. We chose IP hospital stays as our risk endpoint because reductions in overall health, whether due to multiple health conditions [6,25] or aging [26], are associated with increased hospital visits [27,28]. The primary contribution of this work goes beyond the models themselves by matching clinical knowledge to data that are available at scale, across a diverse cohort of patients. In our experiments, we found that the profile demonstrates high performance in terms of the AUROC on the aggregate held-out testing set. Importantly, there was consistent performance of this scoring system across sexes, diverse patient ages, and income levels. The THS model and each of the component models learned to generate predictions by focusing on appropriate, clinically relevant patient features. The THP is personalized based on individual organ system risk drivers, and visualizations, such as radar plots, can be used to facilitate explainability and encourage confidence of clinical decision making, providing meaningful feature importance. The THS outperformed relevant baselines, specifically the other commonly used multimorbidity scoring systems CCI and ECI, for every age, sex, and income bracket. Finally, we conducted multiple ablation tests, while retaining procedure features, to determine the relative contribution of feature groups to the THP. In this experiment, we found that while the combined feature set predictive performance outperformed the prescription, lab, and SDoH ablations, it was largely similar to the diagnosis ablation. However, we hypothesized that we would find larger differences in performance among at-risk populations and found, in patients with multiple comorbidities and on certain prescriptions, minor but consistent increases in predictive power using the combined feature set versus the diagnosis feature set, implying that risk prediction is improved on more complex patients, given more complex data. As more features, including more labs, diagnosis, and prescriptions, are added to the THP, future work will more closely examine which demographics benefit from it.

The THP’s multimorbidity scores can be distinguished from traditional multimorbidity scores in three ways: First, they are derived from a comprehensive set of diagnostic, prescription, lab, and medical procedure data. This is in contrast to other multimorbidity scores that use only one set of information, such as diagnoses (as is the case with the ECI [10] and the CCI [9]) or prescriptions (such as Rx-Risk [29]). Second, these scores were derived from a large and diverse cohort of 794,294 patients with medical data spanning decades. Third, the THP was calculated from patients of both sexes and from across the age spectrum (3-90 years), rather than focusing on mostly geriatric populations as with traditional multimorbidity. To the best of our knowledge, this is the first time that ML was used to integrate multiple types of physician-curated clinical information from a large, diverse population and produce a multimorbidity score that can help guide patient care irrespective of sex and age.

As part of the overall multimorbidity score in the THP, we calculated robust, organ-system-specific scores that provide a more granular picture of health. We believe that these organ-specific multimorbidity scores can complement existing condition-specific scores in clinical practice by providing additional validation for treatment decisions for cardiovascular, respiratory, neuropsychiatric, gastrointestinal, and renal domains. Along these lines, we note that these disease-specific scores often use patient-reported outcomes as part of their input [14], with some even using them exclusively [3]. Although EHR software systems may have health-based modules to automatically compute such scores at the population level, these self-reported data are frequently unavailable or unreliable [30], making it difficult to scale these scores to the population level with a high degree of efficacy. Although the THP cannot be directly compared against these alternative risk scores, as they typically focus on diagnoses versus emergency events, the fact that the THP consistently achieved relatively high AUROCs is nevertheless promising with regard to its ability to complement these more specific risk scores. Specifically, it says something well established in multimorbidity scores but understudied in more specific risk scores: foregoing patient input (which typically contains useful information) entirely, in exchange for more scalable data, can still lead to strong results. Moreover, these alternate risk scores are also typically hyperspecific, limiting their clinical utility to a subset of patients—likely due to them being built on similarly restricted cohorts (eg the American Heart Association pooled cohort equations for atherosclerotic cardiovascular disease derived from cohorts exclusively in the age range of 40-79 years). As our approach has no constraints upon individual patients’ age or sex, and are built using a similarly diverse cohort, risk profiles that are applicable to a far larger population can be easily derived. Of course, we assessed generalizability only among three well-known dimensions (age, sex, and income), and there are far more subtle biases that have been observed even among established risk scores, such as the CHADS2VASC stroke score underestimating risk in patients with chronic renal disease [31]. Further study will be needed to fully examine these sorts of biases in our proposed risk models, but even in this case, the scalability of our approach will only make this research simpler to perform.

Limitations

Data-related limitations of this study include unmeasured variables and incomplete observations. Regarding the former, in this study, we did not include lifestyle behavioral data, such as nutrition, smoking, and physical activity. Although reporting of these factors is known to be inconsistent and unreliable [32], especially in healthy populations (which typically lack recent EHR/claims medical history), they play an important role in clinical outcomes. We believe this would be most addressable through the collection of passive data from wearable sensors, which future work will include. On a similar note, although we were able to use aggregate statistics for race and economic status based on zip-code-derived census data, we were unable to track them at an individual level. Though this form of zip code aggregation has been shown to be useful in clinical risk assessment [33], individual SDoH data could increase the precision and accuracy of THP multimorbidity scores. Future studies of the THP will examine the impact of longer observation and follow-up windows on strategies for clinical intervention. Finally, we note the unreliability of claims data at large, as they are typically produced with financial incentives that are not necessarily aligned with patient care, though they are still often used for risk assessment problems [34,35].

Conclusion

In summary, we combined practical clinical knowledge with modern ML on large-scale data to produce THP multimorbidity scores to aid in decision making across generalizable patient populations. We believe that the THP will allow for more targeted prioritization of care-gap closure, the assessment of comprehensive risk profiles for a greater number of patients, and facilitation of better physician-patient interactions and joint decision making via feature explainability. Although prospective studies will be required to measure the utility of this approach, our intention is that the THS may be used as a preliminary risk stratifier to rapidly prioritize patients for care from a population health management perspective [36]. Once a patient is engaged with a care provider, the organ-specific scores can be used to guide, and explain, individualized clinical interventions based on existing best practices. This would provide the foundation for an integrated continuum between population health and personalized medicine. Finally, we also note the promise that the THP has for clinical research at large, reflecting the rare opportunity to study holistic clinical risk at an extreme scale, potentially unveiling clinically valuable insights.