# Survival analysis - Cox proportional-hazards model problem with output

I have a time to event dataset were I’m looking at the time until individuals perform a specific action. I’m monitoring these individuals for a certain time and if that action is not performed during the monitoring time the individual gets censored in the sense that the monitored trait is neither “yes or no”.

I was thinking that a cox proportional hazards model is perfect to analyse my data and I’ve been able to fit the model. However, my problem is that the output of the model skips some of the independent variables and there might be something that I don’t grasp when it comes to the output. Let me explain this a bit further and I’m sure someone out there have an answer to my problem.

I’m working with R. The first line of code looks the following:

coxph_interaction <-coxph(Surv(Mating, No_mating) ~ A * B, data = dat)


Where “Mating” is the dependent variable, “No_mating” is the censored data and A and B are the independent variables where I’m also interested in their interaction. Both A and B consists of high, low and metabolite which gives us the following combinations:

A           B
High        High
Low         Low
Metabolite  Metabolite


With 9 different combinations of the two (3x3=9) When I’m running the model and then check:

Summary(coxph_interaction)


I get the following output:

> summary(coxph_interaction)
Call:
coxph(formula = Surv(Mating, No_mating) ~
A * B, data = dat)

n= 187, number of events= 187
(92 observations deleted due to missingness)

coef exp(coef) se(coef)      z Pr(>|z|)
A      high                  0.1400    1.1502   0.2784  0.503  0.61518
A      low                  -0.4099    0.6637   0.2813 -1.457  0.14509
B      high                 -0.2511    0.7780   0.3169 -0.792  0.42824
B      low                   -0.9005    0.4063   0.3001 -3.001  0.00269 **
A   high : B high           -0.1475    0.8628   0.4509 -0.327  0.74349
A   low  : B high            0.5164    1.6760   0.4315  1.197  0.23134
A   high : B low             0.4075    1.5030   0.4269  0.954  0.33988
A   low  : B low             0.3094    1.3625   0.4413  0.701  0.48327
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

exp(coef) exp(-coef) lower .95 upper .95
A      high                  1.1502     0.8694    0.6665    1.9852
A      low                   0.6637     1.5066    0.3824    1.1519
B      high                  0.7780     1.2854    0.4180    1.4479
B      low                    0.4063     2.4610    0.2257    0.7317
A   high : B high            0.8628     1.1590    0.3566    2.0878
A   low  : B high            1.6760     0.5967    0.7195    3.9043
A   high : B low             1.5030     0.6653    0.6510    3.4704
A   low  : B low             1.3625     0.7339    0.5738    3.2357

Concordance= 0.619  (se = 0.027 )
Likelihood ratio test= 22.58  on 8 df,   p=0.004
Wald test            = 20.9  on 8 df,   p=0.007
Score (logrank) test = 21.88  on 8 df,   p=0.005


What I can’t wrap my head around is why “metabolite” is missing. Are the other ones compared to that one? I would expect there to be:

A      high
A      low
A      metabolite
B      high
B      low
B      metabolite
A   high : B high
A   low  : B high
A metabolite : B high
A high : B low
A low : B low
A metabolite : B low
A high : b metabolite
A low : B metabolite
A Metabolite : B metabolite


Could someone give me an explanation to this. I would love to have output for the metabolite A, B and interactions (9 of them) as well.

library(survival)
library(lubridate)
library(ggsurvfit)
library(gtsummary)
library(tidycmprsk)
library(condsurv)
library(survminer)


These are the libraries loaded in R.

Here's the head of my dataset head(dat)

    A      B    Mating           No_mating
<chr>   <chr>   <chr>               <dbl>
1 High     High   NA                      0
2 High     High   2                       1
3 High     High   4                       1
4 High     High   2                       1
5 High     High   6                       1
6 High     High   2                       1


If mating did not occur within the timeframe there's a NA, I then put "No_mating" as 0 if mating did not occur to censor the data and 1 if it did occur.

EDIT:

Made the changes in Mating from NA to last event time (9) as suggested and now I get the following output:

    > summary(coxph_interaction)
Call:
coxph(formula = Surv(Mating, No_mating) ~
A * B, data = dat)

n= 279, number of events= 187

coef exp(coef) se(coef)      z Pr(>|z|)
A High                       0.11011   1.11640  0.27724  0.397  0.69125
A Low                       -0.01703   0.98312  0.27718 -0.061  0.95102
B High                      -0.60628   0.54537  0.31657 -1.915  0.05547 .
B Low                       -0.63799   0.52835  0.29368 -2.172  0.02983 *
A High : B High              0.12763   1.13613  0.44984  0.284  0.77663
A Low  : B High              1.12667   3.08536  0.43139  2.612  0.00901 **
A High : B Low               0.24185   1.27361  0.42380  0.571  0.56822
A Low  : B Low               0.11071   1.11707  0.43813  0.253  0.80051
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

exp(coef) exp(-coef) lower .95 upper .95
A High                         1.1164     0.8957    0.6484    1.9222
A Low                          0.9831     1.0172    0.5710    1.6926
B High                         0.5454     1.8336    0.2932    1.0143
B Low                          0.5284     1.8927    0.2971    0.9395
A High: B High                 1.1361     0.8802    0.4705    2.7437
A Low : B High                 3.0854     0.3241    1.3247    7.1863
A High: B Low                  1.2736     0.7852    0.5550    2.9226
A Low : B Low                  1.1171     0.8952    0.4733    2.6365

Concordance= 0.621  (se = 0.022 )
Likelihood ratio test= 22.92  on 8 df,   p=0.003
Wald test            = 23.83  on 8 df,   p=0.002
Score (logrank) test = 25.14  on 8 df,   p=0.001

• For the right-censored cases you need to use the last observation time in the Mating column. With the NAs the software is just ignoring those cases. You then should have a complete data set.
– EdM
Commented Jul 2, 2023 at 12:30
• So with last observation time, lets say we do 9 observations, you mean I'll add 9 in the ones that are supposed to be censored (if right censored)? Commented Jul 2, 2023 at 12:32
• Yes, that's correct, if you followed each of those individuals through time = 9. Note that you could also include individuals who dropped out before the end of the study this way, if you simply use each individual's own last observation time in that column. That might not be an issue for your study, but it's critical in many other applications.
– EdM
Commented Jul 2, 2023 at 13:02

This is the same problem as in other regression models with categorical predictors. One level of each predictor is chosen as the reference. It looks like “metabolite” was chosen as the reference for both A and B, In default coding with R, the regression coefficients represent differences in outcome from the reference condition that are associated with the other combinations of levels of those predictors.

In a Cox model, the baseline cumulative hazard is estimated to get the survival over time for the reference conditions. You can thus get survival estimates for each of your combinations of predictor values even though you don’t have as many coefficients as you might have originally expected,

Be sure that you are coding the survival data correctly; it’s hard to tell from your code. The first argument to Surv() should be the time to the event or censoring. The second argument should typically be 0 for right censoring at that time and 1 for an event.

You also seem to have a lot of missing data. Look carefully at the way your data are coded to make sure that isn’t due to some error. If it isn’t, you might consider multiple imputation to avoid bias in your results.

Three more thoughts, based on later information

The apparently missing data in the model output came from mis-coding the survival times. The survival times (called Mating in the code) were originally listed as NA for individuals who never experienced the event. Those data rows were thus completely ignored in the analysis, and analysis was restricted to those who experienced the event (total n analyzed was equal to the number of events). That's likely to bias the results. That can be fixed by properly coding the time for each individual as the last observation time, whether an event or no-event. There would be no need for imputation if that's the only source of NA values.

With only a few observation times this might better be handled with a discrete-time survival model. See this page and its links, for example. Cox models use some approximations if there are multiple individuals with the same event times, as you clearly have here.

I agree with Peter Flom's suggestion in another answer (+1): the best way to evaluate the model's results is to examine the predictions for each combination of predictor values. Those predictions will be the same regardless of the way that the categorical predictors are coded. For the Cox model you can display the entire predicted survival curves (ideally with confidence intervals). I tend to get into trouble when I try to interpret the results of other coding systems, so I typically rely on the default treatment coding for categorical predictors and do the comparisons post modeling.

• Thanks for the input. This is somewhat how I interpreted it. The problem is that its hard to grasp the output when you don't see all combinations. I've now made some edits with the head of my dataset. Is this the correct way as explained in the edit? The 92 missing observations are the NA's, but doesn't the formula assign those datapoints as censored since No_mating is 0? Thats how I have interpreted it. Otherwise, please feel free to suggest the correct way to do code. As for the multiple imputation, is this still necessary to do if i have N=20 per treatment (a*b combination) ? Commented Jul 2, 2023 at 11:54
• I don't quite get the last paragraph. What do you mean with examining predictors for each combination of predictor values? And doing comparisons post modelling, what do you mean with this? I think I might understand where this is going but I don't quite know how to get there. Could you please elaborate a bit. Commented Jul 2, 2023 at 18:38
• Also, on the comment about the discrete-time survival model, I've looked at the procedure and it seems good actually. Whats daunting is to change my dataset from the way its now to the way suggested by ngoodan in the link you've sent. I have hundreds of rows with individual observations and every observation should then be transformed to nine rows with binary values. I don't think reshape can handle that.. If I proceed this way, I would also have problems fitting the interaction term right (A*B)? Commented Jul 2, 2023 at 19:04
• @Blanca you certainly will end up with a very large data set if you use a discrete-time model, but with today's computers that's not a problem. (I grew up with punch cards for data input and magnetic tape for storing intermediate data values.) The R discSurv package has tools for reshaping data from your "short" format to the required "long" or "person-period" format for discrete survival models. See that package's dataLong() function.
– EdM
Commented Jul 2, 2023 at 19:17
• @Blanca with respect to post-modeling comparisons, I mean using the coefficient estimates from the model to examine different scenarios, with additional tools. For example, you have 9 combinations of the A and B predictor variables. You can use survfit.coxph() with 9 corresponding rows of the newdata data frame in the function call to get survival curves (with standard errors) for all those possibilities. You can use the tools in the emmeans package to evaluate pairwise differences among various A- and B`-variable combinations.
– EdM
Commented Jul 2, 2023 at 19:28

@EdM is correct, I'm just adding some thoughts.

First, there are different ways to operationalize categorical effects: Dummy coding, deviation coding, Helmert coding and more. See UCLA page for a full list of what R offers, with explanations. You can even make up your own. There are some operatinalizations where you will get a parameter estimate for each level (e.g. deviation coding), but be careful. The interpretation can be tricky.

Second, for interpretation of models with interactions, I like to look at the predicted values at different levels. This is often clearer (at least to me) than the parameter estimates. It also will be the same regardless of which operationalization you choose.