ABSTRACT
Objective To derive a predicted probability of death (PDeathDx) based upon complete sets of ICD-10 codes assigned to patients prior to their diagnosis of COVID-19. PDeathDx is intended for use as a summary metric for pre-existing conditions in multivariate models for COVID-19 death.
Methods Cases were identified through the COVID-19 Shared Data Resource (CSDR) of the Department of Veterans Affairs. The diagnosis required at least one positive nucleic acid amplification test (NAAT). The primary outcome was death within 60 days of the first positive test. We retrieved all diagnoses entered into the electronic medical record for visits, on problem lists, and at the time of hospital discharge if they were at least 14 days prior to the NAAT. ICD-9 codes were converted to ICD-10 equivalents using a crosswalk provided by the Centers for Medicare/Medicaid Services. ICD-10 codes were converted to their category diagnoses defined as all columns to the left of the decimal point. Each patient was considered to have or not have each category diagnosis prior to the NAAT. A computer program calculated the number of cases for each category diagnosis, the relative risk (RR) of death, and its confidence interval (CI) using a Bonferroni adjustment for multiple comparisons. RRs were re-centered by subtracting 1 so that high-risk conditions had a positive value while protective conditions had a negative one. Diagnoses found to be significant were entered into a logistic model for death in a stepwise fashion. Each patient was assigned (RR-1) to each category diagnosis if they had the condition or 0 otherwise. The resulting model was used to derive PDeathDx for each patient and the area under its receiver operating characteristic (ROC) curve calculated. Single variable logistic models were also derived for age at diagnosis, the Charlson 2-year (Charl2Yr) and lifetime (CharlEver) scores, and the Elixhauser 2-year (Elix2Yrs) and lifetime (ElixEver) scores. Stata was used to compare the ROCs for PDeathDx and each of the other metrics.
Results On September 30, 2021 there were 347,220 COVID-19 patients in the CSDR. 18,120 patients (5.33%) died within 60 days of their diagnosis. After consolidating ICD-9 and ICD-10 codes, 29,162,710 separate diagnoses were given to the subjects representing 41,341 ICD-10 codes. This set was reduced to 1,890 category diagnoses assigned to the group for the first time on 19,184,437 occasions. Of the 1,890 category diagnoses, 425 involved >= 100 subjects and had a lower boundary for the CI >= 1.50 (a high-risk condition) or upper boundary <= 0.80 (a protective condition). Stepwise logistic regression showed that 153 were statistically significant, independent predictors of death. PDeathDx was slightly less powerful than age as a discriminator (ROC = 0.811 +/- 0.002 vs 0.812 +/- 0.001, respectively; P < 0.001) but was superior to the Charl2Yr (ROC = 0.727 +/- 0.002; P < 0.001), CharlEver (ROC = 0.753 +/- 0.002; P <= 0.001), Elix2Yr (ROC = 0.694 +/- 0.002; P < 0.001); and ElixEver (ROC = 0.731 +/- 0.002; P < 0.001). Univariate analysis and multivariate modeling showed that many of the most high-risk conditions are under-represented or not included in the Charlson Index. These include hypertension, dementia, degenerative neurologic disease, or diagnoses associated with severe physical disability.
Conclusions Our method for handling pre-existing conditions in multivariate analysis has many advantages over conventional comorbidity indices. The approach can be applied to any condition or outcome, can use any categorical predictors including medications, creates its own condition weights, handles rare as well as protective conditions, and returns actionable information to providers. The latter include the specific ICD-10 groups, their contribution to the risk, and their rank order of importance. Finally, PDeathDx is equivalent to age as a discriminator of outcomes and outperforms 4 other comorbidity scores. If validated by others, this approach provides an alternative and more robust approach to handling comorbidities in multivariate models.
INTRODUCTION
Many mathematical models have been proposed for predicting death from COVID-19 infection (1-5). They have gained popularity because they are very useful for managing patients and allocating scarce resources. 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 (“malignancy”) and derive a regression coefficient for each grouping (3, 4). Another approach is to use a comorbidity scoring system such as the Charlson Comorbidity Index (CCI) or Elixhauser (ELIX) score (1, 2, 5). 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 ELIX have been extensively validated and applied to many conditions (6). Moreover, mathematical modeling has shown that they provide all that is necessary to handle confounding from the underlying diagnoses (7). However, there are still problems with either approach:
CCI was developed to predict mortality one year after hospitalization. It is less clear that the same conditions are the most influential for outcomes mediated by an immune response such as breakthrough infections.
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.
There is often little evidence that members of each group have the same prognostic significance.
Neither method handles conditions that are protective.
These approaches do not return actionable information to clinicians. While they may become aware of their patients’ risk, providers are not informed of the sources of such risk. As a result, they cannot take steps to mitigate that risk through further studies or treatments.
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 non-fatal cases.
METHODS
Cases were identified through VA’s COVID-19 Shared Data Resource (CSDR). Membership in this registry requires at least one positive nucleic acid amplification test (NAAT). The primary outcome was death within 60 days of the diagnosis. 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 ELIX score (Elix2Yrs), and lifetime ELIX score (ElixEver) were retrieved from 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. ICD-9 codes were converted to ICD-10 using a crosswalk provided by the Centers for Medicare/Medicaid Services. A “category diagnosis” 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 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. The software used these cell frequencies to derive the relative risk of death associated with the condition (RR) and its confidence interval (CI). CI were adjusted for multiple comparisons by the Bonferroni method. A category diagnosis 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. Another software program was used to create a diagnostic grid in which each patient was assigned a value for all significant category diagnoses. Recall that a condition with a relative risk of 1 has no effect on the outcome. The scale for relative risk was therefore centered on 1 by subtracting 1 from the relative risk and entering that value 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 (Stata). 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 (ROC area). 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 diagnoses. 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. Stata provides a function that compares the ROC areas for pairs of these variables.
RESULTS
On September 30, 2021, there were 347,220 COVID-19 patients in VA’s COVID-19 Shared Data Resource. 339,772 (or 97.9%) had at least one pre-existing condition and form the basis for this report. 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; 95.8% were veterans; 0.7% were on supplemental oxygen; and 11.8% were current smokers. 9.1% had been fully vaccinated at least 14 days prior to the COVID-19 diagnosis. 21.5% acquired their infections after July 1, 2021 and were presumed to have the delta variant. Overall, 18,120 patients (5.33%) died within 60 days of their diagnosis.
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 (or 98.9%) were successfully converted to an ICD-10 equivalent using the CMS crosswalk. In addition, 78,976,269 ICD-10 codes were entered at least 14 days prior to the COVID-19 diagnosis. After consolidation, 29,162,710 separate diagnoses were given to the subjects representing 41,341 ICD-10 codes. The sample size was insufficient to test each ICD-10 code for its prognostic significance. The individual codes were therefore reduced to 1,890 category diagnoses which were assigned to the group for the first time on 19,184,437 occasions.
Of the 1,890 category diagnoses, 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 another 26 provided a perfect prediction of the outcome when present, while 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 diagnoses. 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 and its ROC area determined. This analysis showed that PDeathDx provided excellent discrimination between those who died and those who survived (ROC area = 0.811 ± 0.002). Single variable logistic models were also constructed for age at diagnosis, Charl2Yrs, CharlEver, Elix2Yrs, and ElixEver, and ROC areas determined for their predicted probabilities. PDeathDx was less powerful than age as a discriminator (ROC = 0.812 ± 0.001; P < 0.001) but was superior to the Charl2Yr (ROC = 0.727 ± 0.002; P < 0.001), CharlEver (ROC = 0.753 ± 0.002; P <= 0.001), Elix2Yr (ROC = 0.694 ± 0.002; P < 0.001); and ElixEver (ROC = 0.731 ± 0.002; P < 0.001).
* Most Common ICD-10 Code and Description for Each Category Diagnosis
*Most Common ICD-10 Code and Description for Each Category Diagnosis
**Ranked in Descending Order by Contribution to the Logit Function
Tables 1 and 2 show that this approach provides much more clinical information than the indices. Table 1 ranks the category diagnoses by their RR on univariate analysis. Degenerative neurologic disease and severe functional disability are prominently represented within the top 20 most influential conditions. This observation is significant because they are not given much weight in comorbidity indices. On the other hand, only one solid tumor (carcinoma of the ampulla of Vater) appears within the top 20. On the other end of the scale, several functional diagnoses, gynecologic disorders, and sexually transmitted diseases reduced the risk of death. 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 became 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.
DISCUSSION
We propose a novel method for handling pre-existing conditions in COVID prediction models that has several theoretical advantages over conventional comorbidity indices. Table 3 compares the attributes of the standard comorbidity indices and PDeathDx. Unlike a standard index, our approach creates a different model for every disease. The outcome can be changed by re-defining a target group. Univariate statistical analysis is used to assign condition weights. The method can be applied to any type of categorical predictor including medications and procedures. It handles rare as well as protective conditions. In addition to identifying high-risk codes, the output includes a probability for each pattern of conditions. It can therefore detect circumstances where disease combinations are problematic. Finally, the model can be used to return actional 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.
However, PDeathDx requires programming at a high level, takes more time, and uses more computing resources. The first 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 second 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 large number of rows.
PDeathDx provided powerful discrimination between COVID-19 patients who died and survived, out-performed the other comorbidity indices, had an ROC area equivalent to that generated by age, and remains the second most powerful predictor of all those that we have examined thus far. The most influential category conditions included entries not usually included in comorbidity indices. This observation suggests that it is risky to pick comorbidities for analysis without a systematic review of all those experienced by the cohort. Moreover, 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 for 3 reasons. First, it was intended to represent pre-existing conditions in an overarching model containing multiple domains. One domain was comprised of demographic characteristics including age at diagnosis. Next, we did not want age to displace the diagnoses 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 modelled at all, 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.
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. Finally, to develop the most efficient model, it is reasonable to start with a comprehensive list of diagnoses and work backwards. That strategy allows one to determine if time-limited sampling compromises the ability of PDeathDx to discriminate between survivors and non-survivors. On the other hand, starting with a fixed time frame for diagnoses is arbitrary and provides no information about the benefits of more remote data. Fortunately, our computing resources were more than sufficient to handle all pre-existing conditions.
The major limitation of our approach is that it handles all pre-existing diagnoses – not just the most recent ones. Thus, 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. Of course, our conclusions are limited to patients with characteristics like the veteran population. Further studies should be done on other populations and disease states before the method should be widely applied. If validated by others, our method could provide a more robust alternative to comorbidity scores for handling pre-existing conditions in multivariate models.
Data Availability
All data produced in the present study are available upon reasonable request to the authors.
