3
$\begingroup$

I am undertaking some medical research using R. My outcome of interest is mortality in the intensive care unit.

Data

My data looks like this (there are ~15,000 rows).

      age   illness_score drug_y_use   var_x time mortality
1      48          35       1            1    5.88     0
2      84          56       0            0    1.78     0
3      54          25       1            0    6.15     0
4      67          24       1            0   13.71     1
5      71          20       0            0    2.53     0
6      81          28       0            0    1.02     0

I am interested in the effect of variable x (var_x) on icu mortality, correcting for the other variables: age, illness score and use of drug y.

logistic regression

Using the glm function:
glm(mortality ~ var_x + illness_score + drug_y_use + age,
              data = data.df, family=binomial()))

obtains the following:

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -7.074665   0.179259 -39.466  < 2e-16 ***
var_x            0.373837   0.085712   4.362 1.29e-05 ***
illness_score    0.104640   0.004004  26.131  < 2e-16 ***
drug_y_use       1.404151   0.066076  21.251  < 2e-16 ***
age              0.007041   0.001924   3.660 0.000252 ***

Exponentiating the coefficients reveals a positive association with var_x on mortality

exp(glm$coefficients)
var_x illness_score drug_y_use age 
1.45       1.11        4.07    1.007

Cox regression

Now asking a similar question with Cox Proportional Hazards (using the "time" column to determine the time of censoring (mortality = 0) or event (mortality = 1)

coxph(Surv(time, mortality) ~var_x + illness_score + drug_y_use + age , data = data.df

Gives the following results:

n= 15189, number of events= 1331 
                 coef    exp(coef)  se(coef)  z    Pr(>|z|)    
var_x         -0.285591  0.751570  0.069976 -4.081 4.48e-05 ***
illness_score  0.061434  1.063360  0.003378 18.185  < 2e-16 ***
age            0.017213  1.017362  0.001790  9.617  < 2e-16 ***
drug_y_use     0.435761  1.546139  0.061466  7.089 1.35e-12 ***

Now we have a negative coefficient and the exponentiated coefficient is below 0 suggesting a negative association with var_x

Summary / key question

Using the same data-

I have a significant positive association (exp coeff = 1.45) of var_x on mortality using glm

I have a significant negative association (exp coeff = 0.75) of var_x on mortality using coxph

How is this possible please? The only difference is the addition of timing data with coxph :how can this reverse an association?

$\endgroup$
3
$\begingroup$

Why would you expect the coefficients to be similar? After all, the likelihood function of logistic regression and Cox's proportional hazard function focus on different quantities. In logistic regression, we want to maximise the likelihood of mortality, given the other variables $$ \max \prod_{i=1}^n P(\text{mortality} \mid x_1, \ldots, x_d) . $$ In Cox's proportional hazards model the focus is time to mortality and the partial likelihood function is $$ \max \prod_{i=1}^n P (\text{subject experiences event at $y_i$} \mid \text{one event at $y_i$}, x_1, \ldots, x_d ) . $$ Accordingly, the meaning of coefficients is quite different. In binary logistic regression, they describe the log-odds-ratio, which indicate by how much the odds change when the value of the corresponding variable is increased by 1 \begin{align} \exp(\beta_j) &= \frac{P(y_i = 1 \mid x_{i1}, \ldots, x_{ij}, \ldots x_{id})}{P(y_i = 0 \mid x_{i1}, \ldots, x_{ij}, \ldots x_{id})} / \frac{P(y_i = 1 \mid x_{i1}, \ldots, x_{ij} + 1, \ldots x_{id})}{P(y_i = 0 \mid x_{i1}, \ldots, x_{ij} + 1, \ldots x_{id})} , \\ P(y_i = 1 \mid x_{i1}, \ldots, x_{id}) &= \frac{\exp \left(\beta_0 + \sum_{j=1}^d x_{ij} \beta_j \right)}{1 + \exp \left(\beta_0 + \sum_{j=1}^d x_{ij} \beta_j \right)} , \\ P(y_i = 0 \mid x_{i1}, \ldots, x_{id}) &= 1 - P(y_i = 1 \mid x_{i1}, \ldots, x_{id}) . \end{align}

In the Cox model, coefficients describe the change in hazard ratio if all variables are fixed and only the value of the $j$-th variable is incremented by 1: \begin{align} \exp ( \beta_j ) &= \frac{h(t | x_{i1}, \ldots, x_{ij}, \ldots x_{id})}{h(t | x_{i1}, \ldots, x_{ij} + 1, \ldots x_{id})} , \\ h(t | x_{i1}, \ldots, x_{id}) &= h_0(t) \exp \left( \sum_{j=1}^d x_{ij} \beta_j \right) , \end{align} with (unspecified) baseline hazard function $h_0(t)$ (cancels out).

You are correct, that the expressions for $\beta_j$ boil down to the same expression: $$ \exp ( \beta_j ) = \frac{\exp \left(x_{ij} \beta_j \right)}{\exp \left( (x_{ij} +1) \beta_j \right)} . $$ However, it is important to remember that semantically they aren't the same, because coefficients are determined according to different criteria (likelihood functions). In logistic regression, coefficients describe mortality (a binary variable), whereas in the Cox model they describe time to mortality (a continuous, non-negative variable). Moreover, Cox's proportional hazards model has no intercept.

Finally, it is important to remember that each model has different assumptions that need to be satisfied for the coefficients to be interpreted in this way. Most importantly, the proportional hazards assumption for the Cox model, which states the hazard ratio is a constant independent of time.

$\endgroup$
  • $\begingroup$ Thank you. Does this suggest var_x increases mortality (the binary variable), yet delays time to mortality? This just seemed an odd concept. Interestingly, instead of censoring at discharge, considering ICU discharge as a competing risk (using the cmprsk package in R) resulted in a positive association of var_x on mortality $\endgroup$ – Jon Aug 15 '18 at 12:25
  • 1
    $\begingroup$ Logistic regression is probably not what you want, because your negative class (no mortality) is actually a "I don't know" class. I believe a survival model is the preferred choice. Since you mentioned competing risks, another assumption of the Cox model is that censoring is random and independent of variables. This might not be true in your case. $\endgroup$ – sebp Aug 15 '18 at 13:35
  • $\begingroup$ Although because it is "ICU mortality", the fact that they have been discharged means that they haven't had the event. I am seeing logistic regression with odds ratios, and Cox regression with hazard rations being used in published studies investigating risk factors for ICU mortality. Still not certain why one would be chosen over another to be honest $\endgroup$ – Jon Aug 21 '18 at 7:34
  • 1
    $\begingroup$ Schoenfeld (of residuals fame) argues against survival methods, stating that mortality should be analysed as a binary variable because patients who die in the ICU do not benefit if the duration of their survival is prolonged: ncbi.nlm.nih.gov/pmc/articles/PMC1550820 $\endgroup$ – Jon Aug 21 '18 at 7:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.