Survival analysis - Cox proportional-hazards model problem with output

Question

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

Answer 1

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.

Answer 2

@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.