A novel method for handling pre-existing conditions in multivariate prediction model development for COVID-19 death in the Department of Veterans Affairs.

The Coronavirus Disease 2019 (COVID-2019) has changed lives worldwide. In the USA, COVID-19 is estimated to have 17.2–163.4 times higher mortality rate compared to seasonal influenza; the difference in numbers is based on whether deaths due to pneumonia among those infected with influenza were included [ 1 , 2 ]. In the Department of Veterans Affairs (VA), COVID-19 has been associated with 4.0-, 2.9-, and 1.3 times higher risk of death, mechanical ventilation, and admission to intensive care, respectively [ 3 ]. Several mathematical models have been proposed for predicting death from COVID-19 [ 4–8 ]. Since they are very helpful for patient management and resource allocation, their use has proliferated. One of the most important components of these models is the set of pre-existing conditions. Many diseases increase the mortality rate because they diminish the host response to infection, cause end-organ dysfunction that is further compromised by COVID-19, or severely limit life expectancy and functional status. One method for handling these conditions is to gather them under broad categories (e.g. “malignancy”) and derive a regression coefficient for each grouping [ 6 , 7 ]. Another approach is to use a comorbidity scoring system such as the Charlson comorbidity index (CCI) or Elixhauser comorbidity index [ 4 , 5 , 8–13 ]. With this approach, the groupings are assigned weights, and the weights are added to generate a summary score. The summary score is then used as a covariate in the model. CCI and Elixhauser have been extensively validated and applied to many conditions [ 14 ]. Moreover, mathematical modeling has shown that they provide all that is necessary to handle confounding from the underlying diagnoses [ 15 ]. Having said this, their development was based on solely inpatient data [ 16 , 17 ]. Also, neither method handles protective conditions. Lastly, rare diseases are not well represented. As a result, a person with a common disease of limited impact could be given a poorer prognosis than another with a rare but lethal condition. A more sophisticated approach is to do a systematic survey of all pre-existing conditions, determine which have an impact on outcomes, and generate a predicted probability of death that represents the aggregate risk posed by those that are statistically significant. Like a propensity score, this predicted probability of death (PDeathDx) can be used as a variable in the model. Fortunately, advanced computer methods have made it possible to generate these estimates. The purpose of this study is to describe a novel approach to handling comorbidities in COVID-19 prediction models and compare the performance of PDeathDx with conventional comorbidity measures in distinguishing fatal from nonfatal cases.

As of 30 September 2021, there were 347 220 COVID-19 patients in the CSDR. Of these, 339 772 (97.9%) had at least one pre-existing condition, forming the basis of this study. The mean age at the time of diagnosis was 58.6 ± 16.7 years, 84.1% were male, 22.9% were members of a racial minority, 9.0% were Hispanic, 0.7% were on supplemental oxygen, 11.8% were current smokers, and 9.1% had been fully vaccinated at least 14 days before testing positive. Presumed to have the delta variant, 21.5% acquired their infections after 1 July 2021. Overall, 18 120 patients (5.33%) died within 60 days of their NAAT. For the study cohort, 82 578 233 ICD-9 codes had been entered into the electronic record before the VA converted to ICD-10 codes in 2015. Of these, 81 671 483 (98.9%) were successfully converted to an ICD-10 equivalent using the CMS crosswalk. Additionally, 78 976 269 ICD-10 codes were entered at least 14 days before testing positive. After consolidation, 29 162 710 separate diagnoses were given to patients representing 41 341 ICD-10 codes. The sample size was insufficient to test each ICD-10 code for its prognostic significance. Therefore, using ICD hierarchy, individual codes were aggregated to 1890 category conditions which were assigned to the group for the first time on 19 184 437 occasions. Of the 1890 category conditions, 425 involved ≥100 subjects and had a lower boundary for the CI ≥ 1.50 or upper boundary ≤ 0.80 ( Table 1 ). One diagnosis [Z11 (screening for other viral diseases)] was given to 120 308 subjects and had a high RR for death; it was removed because it may have been used for COVID-19 testing before there was a suitable ICD-10 code. Preliminary regressions indicated that an additional 25 provided a perfect prediction of the outcome when present A56 (Other sexually transmitted chlamydial diseases), A60 (Anogenital herpesviral [herpes simplex] infections), A63 (Other predominantly sexually transmitted diseases, not elsewhere classified), A64 (Unspecified sexually transmitted disease), A74 (Other diseases caused by Chlamydiae), F12 (Cannabis-related disorders), F15 (Vestibular Assessment), F90 (Attention-deficit hyperactivity disorders), G43 (Migraine), J03 (Acute tonsillitis), K58 (Irritable bowel syndrome), L70 (Acne), N76 (Other inflammation of vagina and vulva), N83 (Noninflammatory disorders of ovary, fallopian tube and broad ligament), N84 (Polyp of female genital tract), N87 (Dysplasia of cervix uteri), N93 (Other abnormal uterine and vaginal bleeding), R51 (Headache), R87 (Abnormal findings in specimens from female genital organs), R92 (Abnormal and inconclusive findings on diagnostic imaging of breast), X50 (Overexertion and strenuous or repetitive movements), Z30 (Encounter for contraceptive management), Z32 (Encounter for pregnancy test and childbirth and childcare instruction), Z33 (Pregnant state), and Z60 (Problems related to social environment). Another 20 were affected by collinearity. A diagnostic grid was therefore assembled where each patient was assigned a value for each of the remaining 378 category conditions. If the diagnosis was present, (RR − 1) was assigned as its value and 0 otherwise. Stepwise logistic regression showed that 153 were statistically significant, independent predictors of death ( Table 2 ). PDeathDx was derived for each patient; AUROC = 0.811 ± 0.002. Single variable logistic models were constructed for age at diagnosis, Charl2Yrs, CharlEver, Elix2Yrs, and ElixEver, and AUROCs determined for their predicted probabilities. PDeathDx was less powerful than age as a discriminator (AUROC = 0.812 ± 0.001; P  < 0.001) but was superior to the Charl2Yr (AUROC = 0.727 ± 0.002; P  < 0.001), CharlEver (AUROC = 0.753 ± 0.002; P  < 0.001), Elix2Yr (AUROC = 0.694 ± 0.002; P  < 0.001); and ElixEver (AUROC = 0.731 ± 0.002; P  < 0.001). As can be seen from Figure 2 , PDeathDx has the best performance each of the models among patients who have lower probability of death, which is relevant to the majority of patients. Comparison of AUROC ICD-10 category conditions predictive of death on univariate analysis Multivariate model for COVID-19 death based upon pre-existing conditions Ranked in descending order by contribution to the logit function. Table 1 ranks the category conditions by their RR on univariate analysis. Degenerative neurologic disease and severe functional disability are represented in multiple category conditions within the 20 highest risk conditions while only one solid tumor (malignant neoplasm of other and unspecified parts of biliary tract) appears. In terms of reduced risk of death, several functional diagnoses, gynecologic disorders, and sexually transmitted diseases are present. Table 2 shows the ICD-10 codes that were independently predictive of death rank ordered by their contribution to the logit when present. Many high-risk conditions found upon univariate analysis became less so, to the point of some becoming protective, when adjusted for the effects of others. Hypertension was the most important independent risk factor for death and represented a greater threat than coronary artery disease. Degenerative neurologic diseases were prominently represented at the top of the list, while malignancies comprised the bulk of high-risk conditions.

Cases were defined as VA patients who had at least one positive nucleic acid amplification test (NAAT) at the VA or elsewhere, identified in the VA’s COVID-19 Shared Data Resource (CSDR). We defined death within 60 days of the NAAT as the primary outcome. The outcome was retrieved from the post-index conditions file of the CSDR, which assigns a 1 to those who died and 0 otherwise. Likewise, the 2-year CCI score (Charl2Yrs), lifetime CCI score (CharlEver), 2-year Elixhauser score (Elix2Yrs), and lifetime Elixhauser score (ElixEver) were retrieved from the CSDR for each patient. Pre-existing conditions were identified by reviewing all diagnoses entered into the electronic medical record for outpatient visits, on patient problem lists, or at the time of hospital discharge. “Pre-existing” refers to entries made at least 14 days prior to the diagnosis of COVID-19. This precaution excludes any entries that may have been made during the pre-symptomatic phases of COVID-19. International Classification of Diseases, Ninth Revision (ICD-9) codes were converted to ICD-10 using a crosswalk provided by the Centers for Medicare & Medicaid Services (CMS). A “category condition” was defined as all characters preceding the decimal point for ICD-10 codes or the ICD-9 equivalent. A patient was considered to have or not have each category condition prior to the COVID-19 diagnosis. A computer program was used to identify all patients with a given condition who died or survived as well as all patients without the condition who died or survived. For a visual depiction of workflow, see Figure 1 . The software used these cell frequencies to derive the relative risk (RR) of death associated with the condition and its confidence interval (CI). CIs were adjusted for multiple comparisons by the Bonferroni procedure [ 18 ]. A category condition was considered to have a significant effect on the outcome if there were at least 100 cases and if the lower limit for the CI was ≥1.5 or the upper limit was ≤0.80, denoting risk or protective conditions, respectively. Diagnosis grid overview. Cat: category condition; Pt: patient; Dth: death; RelRisk: relative risk. Another software program was used to create a diagnostic grid in which each patient was assigned a value for all significant category conditions. Since a relative risk of 1 has no effect the scale for relative risk was centered on 1, then that value was entered for the corresponding category condition, if present. Those without the condition were treated as if they had a condition with no prognostic significance and assigned a value of 0. The diagnostic grid was exported to a statistical program (StataMP 17, StataCorp LP, College Station, TX, USA). Stepwise logistic regression was used to identify those that were independently predictive of death. This procedure adjusted the contribution of each condition for the presence of the others, assigned a predicted probability of death to each patient (PDeathDx), and used PDeathDx to generate an area under a receiver operating characteristic curve (AUROC). The product of coefficient × (RR − 1) is the contribution of each condition to the logit function if present. It was therefore used to rank order the importance of the category conditions. Differences in categorical variables were tested by chi-square analysis; differences in continuous variables were analyzed by the unpaired Student’s t -test or rank sum test. Separate logistic models were used to examine the effect of age, Charl2Yrs, CharlEver, Elix2Yrs, and ElixEver as single predictors of death; the AUROC for each was compared to the AUROC for PDeathDx. This study followed the Transparent Reporting of a multivariate prediction model for Individual Prognosis or Diagnosis guidelines specific to development. The New Mexico VA Health Care System Institutional Review Board approved this study, granting a waiver of informed consent.

We present a novel method for handling pre-existing conditions in COVID-19 prediction models, currently available to VA clinicians that has several theoretical advantages over conventional comorbidity indices. (Individuals interested in the model should contact the corresponding author for more information.) Tables 1 and 2 show that this approach provides more clinical information than Charlson or Elixhauser comorbidity indices, handling rare and protective conditions. Table 3 compares the attributes of these indices with PDeathDx. Unlike conventional indices, our methodology can create a different model for every disease state and outcome, flexibilities that have been advocated in comorbidity indices since 1993 [ 19 ]. Whereas the Charlson and Elixhauser indices have static conditions in their models, PDeathDx’s methodology starts with all documented diagnostic codes, permitting flexibility of conditions to be included in the model for each disease state and outcome based on their corresponding RRs and CIs. Following this line of thought, probabilities, analogous to condition weights in the conventional comorbidity indices, will change across diseases and outcomes of interest. Furthermore, since a part of the output also includes a probability for each pattern of conditions, it can detect circumstances where disease combinations are problematic. Finally, the model can be used to return actionable information to providers. Instead of total or condition scores, this information includes the predicted probability of the outcome, the specific diagnoses contributing to the prediction, and their rank order of importance. Such information may prompt the provider to explore a mechanism of injury or an intervention to mitigate the risk. Comparison of methods for handling comorbidities PDeathDx provided powerful discrimination between COVID-19 patients who died and survived; outperforming the other comorbidity indices, it is the second most powerful predictor of those that we have examined thus far. Models using conventional comorbidity indices often assign little weight or usually do not include some of the highest risk conditions (e.g. degenerative neurological diseases); the same is true of conditions associated with COVID-19 severity (e.g. pemphigoid). Moreover, our model revealed poor prognoses were assigned to those with dementia, degenerative neurological diseases, and severe disabilities—conditions not usually associated with cardiorespiratory injury or impaired immunity. These disorders are often associated with a poor quality of life that dictates a conservative approach to treatment. If this decision is common for COVID-19 patients, then death may be more indicative of the patient’s baseline condition than severity of illness. In that event, this outcome may not be suitable for studies of interventions. Finally, multivariate analysis showed that some high-risk conditions became protective when adjusted for the effect of other components of the model. It would have been a mistake to treat the effect of these conditions as independent and additive. We did not include age in our model for PDeathDx. First, it was intended to represent pre-existing conditions in an overarching model containing multiple domains. One domain comprised demographic characteristics including age at diagnosis. More importantly, we did not want age to displace the conditions highly correlated with age in the model. To understand the mechanisms that lead to a fatal outcome, explanatory variables should take precedence over disease markers. If age is modeled, it should be for a residual effect once causal factors have been included. Finally, age-based models do not support decision making at the point of care. For example, a physician might feel comfortable withdrawing care if the patient had Alzheimer’s disease, but not simply because the patient was old. Clinicians make recommendations based upon an assessment of the underlying conditions. It is the patient’s prerogative to decide what is appropriate based upon age and quality of life. Still, in light of a reviewer’s comment, we added age to the model. With age alone having better performance than PDeathDx, we tested its inclusion to the original model. As expected, its addition led to better performance: AUROC = 0.848; P  < 0.001. We identified two studies using clinical information (i.e. laboratory results, vital signs) with a higher AUROC, each limited to the inpatient setting and comprising a few hundred patients [ 20 , 21 ]. Thus, we anticipate the next step in our approach; including demographics, laboratory values, and vital signs; would be the inpatient counterpart to this model, with a higher AUROC than found in this study. By using the current model in the outpatient setting, expeditiously identifying higher risk patients could allow time for preventive measures (e.g. conversations about potentially better management of pre-existing high-risk conditions, individualized messaging about masking or vaccination) to occur before infection or re-infection. We included all diagnoses in the medical record because some conditions have an effect over the patient’s lifetime even if they have not been recently active. Examples include intravenous drug use or sickle cell trait. We also wanted to test our computing resources to see if they were capable of handling problems of this size. (Fortunately, our computing resources were more than sufficient.) Finally, to develop the most efficient model, it is reasonable to start with a comprehensive list of diagnoses and work backward. That strategy allows one to determine if time-limited sampling compromises the ability of PDeathDx to discriminate between survivors and nonsurvivors. Alternatively, starting with a fixed time frame for diagnoses is arbitrary and provides no information about the benefits of more remote data. However, we still note this to be the major limitation of our approach: A person with chronic renal failure (CRF) who undergoes a transplant and regains normal renal function will still be included in the analysis of CRF. Second, PDeathDx requires programming at a high level, takes more time, and uses more computing resources. One hurdle is the systematic search for candidate diagnoses among millions of entries in the electronic record and performing screening statistical tests on thousands of root diagnoses. The other is creating a diagnostic array in which each patient is given a score for every significant diagnosis experienced by the cohort. The table must then be transferred to a statistical program capable of handling hundreds of predictors and a large number of rows. Third, the results of this analysis are dependent upon the accuracy of the documentation of pre-existing conditions, testing positive, and mortality. Documentation exists for 99% capture rate for deaths occurring internal and external to the VA [ 22 ]. Fourth, our conclusions are limited to patients with characteristics like the veteran population, which comprised older, majority male, patients. However, we note a study found similarities with the Medicare population [ 23 ]. With an AUROC = 0.811 and better performance than conventional comorbidity indices, this method shows promise. Allocating better estimates of patient risk which are informed from others’ experiences within the same healthcare system provides a complementary approach to those already in existence to support the learning health system. The approach needs to be replicated in different healthcare systems as differing care across regions and countries, unique patient characteristics in each system, and emerging variants over time may reveal different risk and protective conditions (and thus corresponding probabilities). Further studies should be done on other populations and conditions before the method should be widely applied. If validated, our method could provide a more robust alternative to comorbidity scores for handling pre-existing conditions in multivariate models.