Approximate reciprocal relationship between two cause-specific hazard ratios in COVID-19 data with mutually exclusive events

1 Introduction

Survival analysis studies the longitudinal event data. Regression in survival analysis investigates whether a factor contributes to the hazard (rate of risks) of the event under study. The hazard ratio (HR) is the ratio of two hazards, one with the factor taking the at-risk value (e.g. smoking) and another without (e.g. not-smoking). Since hazard and HR describes the instantaneous risk, or rate of event occurrence (e.g. death from a specific disease), it is a very different concept from the life-long risk of having that event [1]. As a result, regression in survival analysis (e.g. Cox regression) is different from the static case-versus-control regression analysis (e.g. logistic regression). Take the following two statements as example: “smoking makes lung cancer patients die faster” and “smoking makes a person more likely to die from lung caner than a non-smoker”; the first would be a conclusion from a survival analysis, whereas the second from a case-control type of analysis.

The COVID-19 pandemic since 2020 [2–4] provides a unique longitudinal event data. First of all, a COVID-19 patient admitted to a hospital sees his/her outcome relatively quickly: either the patient survived or does not in a matter of days. As a result, there are very few right-censored data where the outcome is still unknown at the time of data collection. Of course, there exist chronic or long COVID-19 survivors who are not completely cured [5–9], but they are unlikely to die from COVID-19 in the future.

The second feature of COVID-19 longitudinal event data is that the two events, death and release from the hospital, are not the traditional “competing risk” events [10, 11]. Although not strictly defined as such, competing risks events are often two unfavorable events with one occurring before the occurrence of another. In COVID-19 data, the event of being released from hospital is a favorite event, and a description of them in principle should not be connected to words like “risk” or “hazard”. A higher HR for the event of releasing from hospital implies a faster recover, thus a factor that contributes to this higher HR actually provide protection. A better name to describe death and release event pairs could just be “competing events” without using the word “risk”.

Also, in the COVID-19 data, the event of death and the event of released from hospital are mutually exclusive. Although organ transplant and death can be mutually exclusive when the organ involved is (e.g.) heart, transplant of many other types of organ is an event that preceed, and may have an impact on, death. Even in the case of heart transplant, survival after the operation is not completely guaranteed. Regardless of these detail, factors affecting transplant timing are basically external, whereas those affecting mortality without a transplant are mostly internal. Our COVID-19 death/release mutually exclusive event pair does not have a correspondence in death/organ-transplant pairs.

As first proposed in Kalbfleisch and Prentice [12] (see also, e.g. [10, 13]), given the event time (T, K) where T is the time to event (or to censoring), K = 1, 2 for two event types (we may understand that k = 0 means right-censored), the cause-specific hazard function can be defined as:

where K for event type (1 or 2), and conditioning T ≥ t means the function is defined before the occurring of either event, and x is a covariate. Equation (1) is the per unit time probability that only event-K occurs. When K = 1 events are considered (e.g. death), K = 2 events (e.g. release) are converted to right-censored (which is essentially equivalent to K = 0), and h K = 1 = d e a t h c s is specified; similarly when K = 2 events are the focused events (e.g. release), death events are right censored, resulting in h K = 2 = r e l c s [12]. On the surface, converting the other not-considered event as right-cencored sample may seem to violate a basic assumption of survival analysis, i.e. the non-informative censoring. But this definition of hazard is better suited to studying of the effect of covariate [10, 14] and it is very convenient to apply as standard programs can be used.

When studying the impact of the effect of a covariate x, we assume an exponential contribution from x on the baseline hazard in Eq. (1) (Cox proportional hazard model), i.e.

We study the two cause-specific hazard ratios (HRs), one for time-to-death (treating release event as right-censored) and another for time-to-release (then treating death event as right-censored):

where x is the independent binary factor (for continuous factor, the definition is similar, with two hazards evaluated at x-levels differing by 1 unit), and the time dependency is supposedly canceled. The question we ask: what is the relationship between the two cause-specific hazard ratios, csHR _k=death and csHR _k=rel? Intuitively, if a factor value leads to faster death, the same factor may lead to slow recovery, and vice versa. In other words, a larger csHR _k=death may imply a smaller csHR _k=rel. Our working hypothesis is that csHR _k=death ≈ 1/csHR _k=rel, which is also hypothesized to be approximately equal to the odds-ratio from a logistic regression analysis.

If the event time is not long, as is the case for COVID-19 patients, at certain fixed time all samples will have one of the two events (death and release) occurred. Odds ratio (OR) for binary variable x is defined as:

where N _death(x = 1) is the number of deceased patients with x = 1, N _rel(x = 1) is the number of released patients with x = 1, etc. The subscript “death” can be dropped because it is understood that we are interested in how the presence of a x value contributes to the death. For continuous variable x, x = 0 and x = 1 is replaced by a unit change of the variable value x = x ₀ and x = x ₀ + 1, averaged over all x ₀. OR is not defined in the survival analysis, but is intrinsically related to the HR. A possible relationship between OR and HR is another hypothesis.

We will use a survival data of n = 450 COVID-19 test our working hypotheses. Previously we asked the question on whether the two survival analyses and one logistic regression can all identify a risk factor [15]. Here we are asking more quantitative questions concerning these analyses.

2 Data

The COVID-19 patients dataset (n = 450) used in this study was collected from the Tokat State Hospital. Electronic medical records, including patient demographics, clinical manifestation, comorbidities, laboratory tests results were collected. According to the severity and outcomes of the patients, they were divided into two groups: deceased/nonsurviving, released/survived severe groups. Clinical and laboratory data were comparable among the two groups. Ministry of Health permission and the ethics committee of Tokat GaziOsmanPasa University Ethics committee permission was obtained with the number 83116987/360 on March 4, 2021.

Besides age and gender, we selected these 18 laboratory testing measurements at the time of hospital admission in the survival analysis n value is the number of samples with measurement value): (from the complete blood count panel) white blood cell (WBC) count (n = 432), neutrophil (NEU) count (n = 383), lymphocyte (LYM) count (n = 382), hemoglobin (HGB) (n = 432), platelet (PLT) (n = 432), mean corpuscle volume (MCV) (n = 432), mean platelet volume (MPV) (n = 430); (from metabolic panel) glucose (n = 449), alanine adinotransferase (ALT) (n = 435), aspartate aminotrasferase (AST) (n = 429), blood urea nitrogen (urea or BUN) (n = 436), creatine (n = 436), calcium (n = 393), potassium (n = 433), sodium (n = 435); (others) ferritin-1 (FER1) (n = 407), d-dimer (n = 390), and lactate dehydrogenase (LDH) (n = 316). The NLR (neutrophil/lymphocyte ratio) and PLR (platelet/lymphocyte ratio) are derived quantities. There are other test measurements, either only available for fewer number of patients, or due to other reasons, that are not included in the analysis.

3 Methods

The single-variable Cox regression in Eq. (2) can be extended to multivariate context, and the cause-specific log hazard ratio is defined as:

The left-hand side of Eq. (5) contains time t whereas the right-hand side does not, which is the proportional hazard assumption, that the time-dependent part is canceled. i is the index for factors (the index i on the left hand size of Eq. (5) is implied). From Eq. (5), it can be seen that csH R death = e β K = 1 , i , csH R release = e β K = 2 , i , i.e. HR is an exponential function of the regression coefficient β. In practice, cause-specific Cox regression is carried out by labeling event = 2 as right-censored event = 0 when focusing on event = 1, and vice versa when focusing on event = 2.

The R (https://www.r-project.org/) statistical package survive [16] is used in the analysis. The coxph is used for Cox regression, cox.zph is used for testing proportional hazard assumption, based on Schoenfeld residual method [17].

4 Results

Summary statistics: Of the n = 450 samples, 353 (78%) are released from the hospital, 97 (22%) were deceased, with a deceased/released ratio of 27.5%. There is no right-censored patients in this collection. A 57.6% of patients are male, and this proportion is slightly lower for the released group (56.1%), as compared to the deceased group (62.9%). The mean/median age of all patients is 59.7/66, and the mean/median age for the released group is much younger (55.8/60) than that for the deceased group (74.1/75).

Table 1 provides a summary statistics of all 18 factors, for all patients, released group, and deceased group. Both mean and median value are listed. The following variables do not seem to follow a normal distribution ferritin-1, d-dimer, glucose, ALT, AST, urea, creatine, LDH, NEU, and perhaps LYM and PLT. These can either be due to a long-tail in the distribution (there are few samples with extreme large values compared to the mean/median) or perhaps due to a potentially log-normal distribution (then the variable is better logarithm transformed), such as ferritin-1.

As a consequence of potential violation of normal distribution, for the test result between the levels for the released and deceased groups, we present both t-test p-values and Wilcoxon test p-values in Table 1. A striking observation of the result in Table 1 is that all factors are significantly different between the released and the deceased group. The deceased patients exhibit distinct level of many factors and form a cluster by itself. We have noticed this already in Cetin et al. and Ulgen et al. [15, 18, 19].

The mean/median time to either dead or release events is 10.3/8, the mean/median of time to release is 9.5/8, and that of time to death is 13.2/11. In other words, the time-to-death tends to be 2–3 days longer than the time-to-release. For both time-to-release and time-to-death, the logarithm-transformed time follows a normal distribution better than the non-log-transformed time. Though not included in this data, we also observed that ICU surviving patients have a hospitalization stay time longer than those of non-ICU surviving patients, and shorter than those of ICU-deceased patients [20].

Overview of the survival analysis and logistic regression analysis results: Large number of analysis run results are included in Tables 2 –4. Table 3 contains factors in a metabolic panel, Table 4 are factors in a complete blood count panel, and Table 2 lists the rest. The first column is the result from single-factor Cox regression survival analysis using death as event of interest and release as right-censored event. The second column is the Cox regression results by switching the two events. The listed results include the cause-specific hazard ratio (csHR) and its 95% confident interval (CI), and the p-value for testing csHR equal to 1. The third column is the logistic regression result comparing the dead and release samples, where the results include odds-ratio (OR) and its 95% CI, and p-value for testing OR = 1.

We also run the same group of analysis on log-transformed factor values, if that factor better follows a normal distribution after log transformation. Then, we run the same set of analysis, when the factor value is continuous, for binarized factors with optimally chosen threshold values. Both these runs will be discussed in detail later. The reason for running large number of analyses is not for “fishing expedition” in order to have a better chance to find statistically significant results, but to test the robustness of the results. Therefore, we do not do a multiple testing correction.

The proportional hazard ratio assumption is tested by the Schoenfeld residual method in implemented in the R function “cox.zph” in the “survival” package [17]. The test result sometimes depends on if the independent variable is log-transformed or not. We found that the proportional hazard ratio assumption is not rejected for time-to-death events, and mostly not rejected (at 0.01 level) for time-to-release events, except for age (pv = 0.0014), (log)d-dimer (pv = 3E-6), and (log)ALT (pv = 5E-5).

We mark those p-values that are smaller the 0.005 in Tables 2 –4. The reason to use 0.005 instead of 0.01 or 0.05 is explained in Ioannidis and Colquhoun [21, 22], and the practice of always adding a level when using the word “significant” is proposed in Wasserstein et al. [23]; see also Li et al. [24]. Strikingly, almost all factors significantly (at 0.005 level) influence the rate of event of COVID-19 patients, and are significantly different between the deceased and survived group.

Relationship between HR _d and HR _r for continuous factors: Another striking observation is the relationship between the two cause-specific hazard ratios. Denote csHR _d for csHR of the event of death from COVID-19 and csHR _r for csHR of the event of COVID-19 patient releasing from hospital. it can be easily seen from Tables 2 –4 that if csHR _d > 1, the corresponding csHR _r < 1, and vice versa. A simple mathematical expression of this fact is:

The only exception is the factor of gender. But the 95% CI is so large to have both < 1 and > 1 values, it should not really be considered as an exception. The opposite direction of csHR _d and csHR _r is understandable: a risk factor for faster death in a deceased patient would also make a surviving patient recover longer.

Tables 2 –4 seem to also show that the larger csHR _d , and smaller csHR _r , and vice versa. In order to check if there is a numerical relationship between csHR _d and csHR _r , we plot 1/HR _r as function of csHR _d in Figure 1. The line csHR _d × csHR _r = 1 is marked by the slope = 1 line in Figure 1. There are many factors clustered near the csHR _d = csHR _r = 1 point and a close-up plot is shown separately. A factor is labeled in red if p-values for both csHR is significantly (at level 0.001) different from 1, and in blue if one or both csHR is not significant. Our working hypothesis can be written as a reciprocal relation between csHR _d and csHR _r :

Figure 1:

The x-axis is the cause-specific hazard ratio csHR _d for death event, and y-axis is reciprocal of the hazard ratio for the release event (1/HR _r ), for the 18 blood test measurements as well as age and gender. The diagonal line indicates the exact relationship csHR _d = 1/csHR _r . A factor is in red if its p-values (for testing csHR _d and csHR _r = 1) are both smaller than 0.001; in blue if both p-values are larger than 0.001; and light-blue if one of the two p-values is smaller than 0.001. Because there are many factors having HR close to 1, the right subplot presents a close-up near csHR _d = csHR _r = 1.

Extension of the relationship between HR _d and HR _r in a multivariate context: One may wonder whether the approximate reciprocal relationship between HR _d and HR _r for a variable when other variables are also used in a multiple cause-specific Cox regression. Because adding more variables in a multiple regression often leads to unpredictable outcome due to collinearity between variables, and due to overfitting when the variable-per-sample ratio is too high, we investigate a simple situation: we run a three-variable Cox regression with the variable of interest, plus gender and age co-variates.

Table 5 shows the cause-specific hazard ratio of 18 continuous factors (log-transformed if necessary) conditional on gender and age, for either time-to-death or time-to-release. Without exception, if one of the csHR _d or csHR _r is larger than 1, the other is smaller than 1. The product of the two hazard ratios are close to 1 except for log(urea), log(creatine), and log(LYM). The urea and creatine factors are highly correlated, and both are significantly correlated with lymphocyte.

Relationship between hazard (rate) ratio and odds (risk) ratio: As emphasized in Sutradhar and Austin [1], survival analysis estimates the relative rates (of risk) whereas case-control type of analysis such as logistic regression estimates the relative (static or cumulative) risks. Tables 2 –4 show that OR seems to have a larger magnitude than csHR _d . Figure 2 shows y = OR as a function of x = csHR _d . Unlike Figure 1 where dots are scattered near the slope = 1 line, in Figure 2, dots systematically deviates from the diagonal line. In fact, if a factor is a risk (OR > 1 and csHR _d > 1), we have OR > csHR _d , and if a factor is a protection (OR < 1 and csHR _d < 1), then OR < csHR _d . We can summarize these into a working hypothesis:

Figure 2:

The x-axis is the cause-specific hazard ratio csHR _d for death event, and the y-axis the odds-ratio (OR) from logistic regression, for 20 factors. The right subplot presents a close-up near csHR _d = OR = 1.

Equation (8) might be proved in a simple approximation as follows: the risk function F(t) is known to be related to hazard rate function h(t):

Therefore, the odds-ratio is (the subscript 1, 2 refers to two states in a binary variable or a two numerical level of a continuous factor with one unit difference, and not refers to the two competing-risks):

In the proportional hazard assumption, h ₁(t) = α ₁ h ₀(t), h ₂(t) = α ₂ h ₀(t), where h ₀(t) is a baseline hazard function, and denote ∫ 0 t h 0 ( t ′ ) d t as H ₀(t), Eq. (10) becomes

One approach in getting the approximation in Eq. (11) is the Taylor expansion assuming small H ₀(t) [25].

Relationship between csHR _d and csHR _r for binarized continuous factors: The HR for continuous factors measures the ratio of two hazard rates when the unit of the factor increases by one. Therefore, when one unit change is negligible compared to the possible range, HR can be very close to 1. In order to see the true impact of a continuous variable, we discretize continuous factors into levels. Although one can choose three levels for below-normal, normal, and above-normal, the normal range of a factor may not be universally accepted.

We use binary levels (higher and lower than a threshold) where the threshold value is a compromise between two selections: the first threshold is chosen to maximize the Youden index [26], which is simply the sum of sensitivity and specificity (minus 1). The second threshold is chosen by providing the population prevalence of cases (which is set at 10%), which in turn gives weight to samples in the dataset (see, https://www.medcalc.org/manual/roc-curves.php and [27]). Both thresholds are obtained from the medcalc.org program. The final threshold is a geometric mean of the Youden-based threshold and that after considering the 10% population case prevalence. The resulting threshold values for all factor (except for age, neutrophil/lymphocyte ratio, platelet/lymphocyte ratio, where the thresholds are more intuitively selected) are given in Tables 2 –4. The resulting csHR _d and csHR _r have larger magnitude than the corresponding continuous value, because the “unit change” is much larger.

In order to study the numerical relationship between csHR _d and csHR _r for binarized factors, we examine other possible threshold values beyond that determined by the medcalc program and plot csHR _d (red) and 1/csHR _d (blue) as a function of threshold value, for 19 test measurements, in Figure 3. The 95% confidence interval (CI) of csHR _d or 1/csHR _r is marked with dash vertical lines. If the discretized factor is not significant at a corresponding threshold, red dots turn pink and blue dots turn light-blue. We also mark normal ranges of blood tests from two different sources (as grey horizontal lines) and the threshold used in Table 1 (as downward arrow in grey). We consistently found the 95% CI of csHR _d and 1/csHR _r overlap with each other at the chosen (optimal) threshold value. In other words, when a reasonable threshold value is used to convert a continuous factor to a binary factor, and running two survival analyses results in a roughly reciprocal relationship between the two cause-specific hazard ratios.

Figure 3:

The csHR _d (red) and 1/csHR _r (blue) for binarized 18 blood test measurements as a function of the threshold used to discretize these factors. The 95% CI of csHR _d or 1/csHR _r are shown in dashed vertical lines. If the discretized factor’s csHR is not significant (at 0.01 level) by the survival analysis, its color turns from red (blue) to pink (light-blue). The threshold used in Table 1 is shown as a downward arrow. The two horizontal lines represent the normal range of these blood test results (from two different sources), and horizontal dash line is csHR = 1.

Relationship between csHR _d and OR for binarized continuous factors: Similar to Figure 2 where we show the scatter plot for cause-specific hazard ratio for time-to-death (x-axis) and logistic regression odds-ratio (y-axis), Figure 4 shows the similar scatter plot for the discretized factors. The 95% CI for both are shown by horizontal and vertical segments. All points in the first quadrant are above the diagonal line, and those in the third quadrant below the diagonal line. Therefore, | log(OR)| ≥ | log(csHR _d )| is true for all binarized factors.

Figure 4:

Similar to Figure 2, but for discretized/binarized factors: the x-axis is the cause-specific hazard ratio csHR _d for death event, and the y-axis the odds-ratio (OR) from logistic regression.

Effect of log-transformation of factor values: Remember the meaning of HR for continuous variable is the ratio of two hazards evaluated at two factor values differing by one unit: x ₁ = c, x ₂ = c + 1. The dependence on c is supposed to be averaged out. If the factor is log-transformed, the two evaluation points are log(x ₁) = c′, log(x ₂) = c′ + 1, or x ₂ is x ₁ multiplied by a constant 2.718. Not only this one-unit change is much larger, but also, the change in the original scale x ₂ − x ₁ = e ^c ′⁺¹ − e ^c ′ = 1.718 ⋅ e ^c ′ depends on c′. Tables 2 –4 show that csHR _d or csHR _r are dramatically larger (in magnitude) than the un-logged factors. The effect of log-transformation on p-value is unclear, though in our examples, the test becomes more significant after the factor being log-transformed.

The concept of hazard ratio has been cautioned in Hernán [28]. In particular, the instantaneous incident rate for an event to occur may depend on the time, and HR is only an average over the potential time dependence. Our attempt to binarize a continuous factor or log-transform a factor illustrates a similar problem. If HR not only depends on the one-unit-change step, but also depends on which level this step is made, and depends on whether the step is in additive scale or multiplicative scale, then the average may not capture the whole spectrum of the behavior of hazard in its full range.

5 Discussion

The COVID-19 time-to-event data is unique not only the two events, mortality and discharge, are mutually exclusive, but also a susceptible factor might be behind faster death and slower release at the same time. Even though csHR and OR, i.e. instantaneous risk and cumulative risk, are not theoretically proven to always point to the same direction, in our data, they do. If we consider heart transplant and other open-heart operations as events that may have saved a patient’s life, thus are mutually exclusive with mortality, we can not say that the factors causing a longer waiting time until operation are the same ones causing a faster death without these operation. For this reason, we have the basis for the reciprocal cause-specific hazard (csHR) ratios hypothesis which can only be examined in data similar to COVID-19, but not in other survival data just because two events are mutually exclusive.

Besides cause-specific competing risk survival analysis proposed in Kalbfleisch and Prentice [12], there is another subdistribution hazard proposed in Fine and Gray [29]:

where ∪ (logical OR) means two groups of samples are considered in the calculation of hazard (h ^sd ): one is those who has not yet experienced the type-1 event at time t (e.g. still alive), and another is those who has already experienced the type-2 event (e.g. released) before time t. Why already released samples are still considered in the calculation for hazard (rate of risk) for death is not explained, and h ^sd does not have a good interpretation [11]. The reason that subdistribution hazard (Fine-Gray model) is still used is due to the fact that its hazard ratio (sdHR) always preserves the direction (larger or smaller than 1) in odds-ratio. The csHR does not have such a general proof for direction preserving. However, as seen in this paper, sdHR and OR are in the same direction without exception in our data. On the other hand, sdHR has its own problems [10, 30–33].

The possible links between the two csHR’s have at least two consequences. The first is on hospitalization stay time. If a patient has certain condition (e.g. high glucose or hyperglycemia), the larger-than-1 csHR _d implies that the patient has a higher death-rate than those with normal glucose level (remember again that although csHR and OR can not be proven to be in the same direction in theory, in practice such as our data, the two are pointing in the same direction); on the other hand, the larger-than-1 csHR _r value indicates the patients have a higher chance to be released. Therefore, whether hyperglycemia increases the hospital stay time or not depends on whether the patient survives or not. The hospital stay time is of interest because of the number of hospital beds is limited, and there is a need for bed management [34].

If a factor/condition causes the severity of a COVID-19 patient, intuitively we would conclude that patients with the condition will stay in hospital longer. In reality, if the disease is too severe, the patient will stay in hospital shorter, because the patient succumbed to death faster. Among the deceased patients, we would expect the co-existence of short and long stays, while less diverging in stay time for discharged patients. Indeed, though mean of log(stay time) between the deceased and released groups is not significantly different (t-test p-value = 0.15, though Wilcoxon test p-value is 0.00034), the variances are very different (F-test p-value = 1.1E-8).

The second consequence that two csHRs might be related is that if we focus on time-to-release events, we could collect much more samples simply because more patients being recovered/released than deceased. In a sense, this strategy examines which factor delays the recovering time in surviving patients. Larger sample sizes would help to detect more subtle causing factors. This strategy will become more relevant if life-saving drugs for COVID-19 are developed and nobody or almost nobody die from the disease. Even in that future event, we still have in possession surviving patients.

The fact that OR > csHR _d if OR > 1 (and OR < csHR _d if OR < 1) for both continuous factors and their discretized version, seem to be a consequence of the definition of the two quantities. Although one may use this result to obtain an upper limit of csHR _d , the result in Tables 2 –4 seems to indicate that the bound is not tight. In that case, if OR ≫ csHR _d , OR will not be very useful in estimating the csHR _d value.

As discussed thoroughly in the literature that we can not always assume csHR (unlike subdistribution HR) is in the same direction as OR [10, 14]. In other words, csHR > 1 does not universally imply OR > 1 . Individual csHR also can not determine the cumulative incidence function caused by multiple risks [35]. However, our results in Tables 2 –4 show that OR and csHR _d are always in the same direction (both larger than 1, or, both smaller than 1), indicating the difference between a theoretical possibility and reality. Drawing cumulative incidence function is also not a goal in our analysis. Considering all these, we consider the use of csHR better fitted for COVID-19 death/release survival data, than the subdistribution HR. In fact, we doubt subdistribution HR can be applied to this situation at all, because of the exclusive nature of the two events.

In conclusion, we draw attention to the connection between the two types of mutually exclusive events, mortality and discharge, in COVID-19 survival data. We also made three observations from COVID-19 data: the opposite direction between the two csHRs log(HR _d ) ⋅ log(HR _r ) < 0, approximately reciprocal link between them HR _d ⋅ HR _r ≈ 1, and odds-ratio as an upper limit of HR _d : |log(csHR _d )| ≤ |log(OR)|.