parameterize, estimate & interprete interaction terms between two factors in a Cox Proportional Hazard model

Question

Overview

I was trying to fit a cox proportional hazard model to look at interactions between two time-constant covariates, both of which are factors. I parameterized the model in two different but equivalent ways, but parameter estimates from the models were somehow different. I wonder what went wrong. Please see below for an elaboration of the issue.

Data

For reproducibility, I will use a public available dataset, 'lung', which comes with the library 'survival', to illustrate the issue. Suppose the key interest is age by sex interaction. Age is continuous in this example data set; for illustration purposes, I binned this variable into 4 groups as follows.

library(survival)
dat0=lung
dat0$age_group=1
dat0$age_group[is.na(dat0$age)]=NA
dat0$age_group[dat0$age>56 & dat0$age<=63]=2
dat0$age_group[dat0$age>63 & dat0$age<=69]=3
dat0$age_group[dat0$age>69]=4
table(dat0$age_group)

Models

I fitted two equivalent models. Model #1 is straightforward. But to fit model #2, I created a new variable, age_by_sex, based on a combination of age_group (4 levels) and sex (2 levels). They are coded as 0, 1, 2... up to 7.

# create a new variable: age_by_sex
age_by_sex = 0 # default value
age_by_sex[is.na(dat0$age_group)|is.na(dat0$sex) ]=NA

# for sex == 1
for(i in 2:4){
  age_by_sex[dat0$sex == 1 & dat0$age_group==i ]= i-1 
}

# for sex == 2
for(i in 1:4){
  age_by_sex[dat0$sex == 2 & dat0$age_group==i ]= i+3 
}

dat0$age_by_sex=age_by_sex
dat1=dat0[,c("time", "status", "sex","age_group", "age_by_sex")]
dat2=dat1[complete.cases(dat1),] # final dataset for model fitting

# Model #1
cph1= coxph(Surv(time, status) ~ as.factor(age_group)*as.factor(sex), data=dat2)

# Model #2
cph2= coxph(Surv(time, status) ~ as.factor(age_by_sex), data=dat2)

Output

# Model 1 -------------------------------------------------------------

> cph1
Call:
coxph(formula = Surv(time, status) ~ as.factor(age_group) * as.factor(sex), 
    data = dat2)

                                         coef exp(coef) se(coef)      z     p
as.factor(age_group)2                 -0.1040    0.9012   0.4699 -0.221 0.825
as.factor(age_group)3                  0.1725    1.1883   0.4743  0.364 0.716
as.factor(age_group)4                  0.3543    1.4252   0.3938  0.900 0.368
as.factor(sex)2                       -0.5445    0.5801   0.3367 -1.617 0.106
as.factor(age_group)2:as.factor(sex)2  0.3399    1.4049   0.4784  0.711 0.477
as.factor(age_group)3:as.factor(sex)2 -0.1607    0.8516   0.4777 -0.336 0.737
as.factor(age_group)4:as.factor(sex)2      NA        NA   0.0000     NA    NA

Likelihood ratio test=14.47  on 6 df, p=0.02477
n= 197, number of events= 140 

# Model 2 -------------------------------------------------------------

> cph2
Call:
coxph(formula = Surv(time, status) ~ as.factor(age_by_sex), data = dat2)

                           coef exp(coef) se(coef)      z     p
as.factor(age_by_sex)2  0.27652   1.31854  0.27900  0.991 0.322
as.factor(age_by_sex)3  0.45836   1.58148  0.25834  1.774 0.076
as.factor(age_by_sex)4 -0.44045   0.64375  0.32898 -1.339 0.181
as.factor(age_by_sex)5 -0.20455   0.81502  0.34017 -0.601 0.548
as.factor(age_by_sex)6 -0.42864   0.65140  0.32788 -1.307 0.191
as.factor(age_by_sex)7 -0.08613   0.91748  0.34785 -0.248 0.804

Likelihood ratio test=14.47  on 6 df, p=0.02477
n= 197, number of events= 140 

Problems

Information about the likelihood test, df, and p (down the bottom of each output) confirms the two models are equivalent. But I noticed two problems:

Problem # 1. The two models overlap in some parameters, but these overlapping parameters do not agree. For example, the exp(coef) of as.factor(age_group)2 in model 1 is the hazard of the group (sex ==1 & age_group==2) relative to the group (sex ==1 & age_group==1). The estimate is 0.9012. I was expecting the same estimate for as.factor(age_by_sex)2, which was 1.31854. I wonder what it is happening?

Problem # 2. What happened to the model estimate for as.factor(age_group)4:as.factor(sex)2? Why is it NA?

Many thanks for any insights into the problems.

Answer 1

You have an issue in your code that's causing this. Let's start from the end.

The reason the first model returns NA for one of the coefficients is that there is no record in your dat2 dataset that has sex = 1 and age_group = 1. These records do exist in dat1 but you have thrown them out via the complete.cases call. Because this combination does not exist the full interaction design is aliased, and coxph presumably throws out the last collinear predictor. This happens deep inside C code so I'm not exactly sure how this is implemented.

That issue does not occur in your second model, because there also isn't any record with age_by_sex = 0. The resulting combined factor has only 7 levels, the model automatically takes the first as reference - note that this is sex = 1 & age_group = 2! - and only 6 non-aliased parameters are fit.

You have records in dat1 where age_by_sex = NA because you created this variable as a scalar, and you seem to be unaware of what happens if you index a vector out of bounds:

age_by_sex = 0                      ## Only age_by_sex[1] is defined
age_by_sex[seq_len(nrow(dat0))]     ## You index beyond its bounds
> [1]  0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ...

These NAs carry forward into all next steps. An easy fix is to start with the correct dimensions: age_by_sex = rep(0, nrow(dat0)). No records will get dropped from dat1, and you'll see that some of the parameters match between models.