Time-Dependent Cox Regression Model

When the proportional hazard (PH) assumption for age is not met in a Cox regression model, it means that the effect of age on the hazard rate is not constant over time. For example, the hazard rate might increase at a slower rate at younger ages and then accelerate at older ages. The PH assumption test can be judged by the following four methods:

There is no intersection between the Kaplan-Meier survival curves of categorical covariates.
Draw the log cumulative hazard curve. The parallelism of the curve is a sufficient and necessary condition to meet the PH assumption.
Use the Schoenfeld residual method to view the P value and residual graph corresponding to the continuous variable.
The interaction term between each covariate and the logarithmic survival time can be put into the model. If the interaction term is not statistically significant, the proportional hazard condition is met.

The purpose of this case study is to explore the relationship between gender, age, and Karnofsky score and survival outcomes in patients with lung cancer.

Read data

setwd("C:\\Users\\hed2\\Downloads\\dell\\code-storage\\code")
mydata <- read.csv("Coxprop.csv")   
View(mydata)   
table(mydata$status)

## 
##   0   1 
##  42 158

prop.table(table(mydata$status))

## 
##    0    1 
## 0.21 0.79

Collinearity diagnosis

library(car)

## Loading required package: carData

lm.reg<-lm(status~age+sex+ph.karno,data=mydata)  
vif(lm.reg)

##      age      sex ph.karno 
## 1.044646 1.008268 1.036241

Cox regression

library(survival) 
fit<-coxph(Surv(time,status)~factor(sex)+age+ph.karno,data=mydata) 
summary (fit)

## Call:
## coxph(formula = Surv(time, status) ~ factor(sex) + age + ph.karno, 
##     data = mydata)
## 
##   n= 200, number of events= 158 
## 
##                   coef exp(coef)  se(coef)      z Pr(>|z|)   
## factor(sex)2 -0.470053  0.624969  0.170960 -2.749  0.00597 **
## age           0.011799  1.011869  0.009628  1.225  0.22042   
## ph.karno     -0.011885  0.988186  0.005975 -1.989  0.04670 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## factor(sex)2    0.6250     1.6001    0.4470    0.8737
## age             1.0119     0.9883    0.9930    1.0311
## ph.karno        0.9882     1.0120    0.9767    0.9998
## 
## Concordance= 0.626  (se = 0.026 )
## Likelihood ratio test= 15.48  on 3 df,   p=0.001
## Wald test            = 15.4  on 3 df,   p=0.002
## Score (logrank) test = 15.62  on 3 df,   p=0.001

Proportional Risk test

ph.test<-cox.zph(fit)
ph.test

##             chisq df      p
## factor(sex) 2.296  1 0.1297
## age         0.268  1 0.6044
## ph.karno    7.012  1 0.0081
## GLOBAL      8.804  3 0.0320

# schoenfeld test
library(survminer )

## Warning: package 'survminer' was built under R version 4.4.3

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 4.4.3

## Loading required package: ggpubr

## 
## Attaching package: 'survminer'

## The following object is masked from 'package:survival':
## 
##     myeloma

ggcoxzph(ph.test)

Survival function curves

# Plot the baseline survival function
ggsurvplot(survfit(fit,data=mydata), color = "#2E9FDF",
           ggtheme = theme_minimal())

## Warning: Now, to change color palette, use the argument palette= '#2E9FDF'
## instead of color = '#2E9FDF'

# Create the new data when age and  are are average
sex_df <- with(mydata,
               data.frame(sex = c(1, 2), 
                          age = rep(mean(age, na.rm = TRUE), 2),
                          ph.karno = rep(mean(ph.karno, na.rm = TRUE), 2)
                          )
               )
sex_df

##   sex    age ph.karno
## 1   1 62.335     81.9
## 2   2 62.335     81.9

# Survival curves by sex group
fit2 <- survfit(fit, newdata = sex_df)
ggsurvplot(fit2, data=mydata, conf.int = TRUE, legend.labs=c("Sex=1", "Sex=2"),
           ggtheme = theme_minimal() )

Time-Dependent Cox Regression

Add a interaction term x*log(t+20)

fit.ph<-coxph(Surv(time,status)~sex+age+ph.karno+tt(ph.karno),
              tt=function(x,t,...) x*log(t+20 ),data=mydata)
summary(fit.ph)

## Call:
## coxph(formula = Surv(time, status) ~ sex + age + ph.karno + tt(ph.karno), 
##     data = mydata, tt = function(x, t, ...) x * log(t + 20))
## 
##   n= 200, number of events= 158 
## 
##                   coef exp(coef)  se(coef)      z Pr(>|z|)   
## sex          -0.484654  0.615910  0.170935 -2.835  0.00458 **
## age           0.012503  1.012582  0.009682  1.291  0.19656   
## ph.karno     -0.093643  0.910608  0.039707 -2.358  0.01836 * 
## tt(ph.karno)  0.014687  1.014795  0.007091  2.071  0.03833 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## sex             0.6159     1.6236    0.4406    0.8610
## age             1.0126     0.9876    0.9935    1.0320
## ph.karno        0.9106     1.0982    0.8424    0.9843
## tt(ph.karno)    1.0148     0.9854    1.0008    1.0290
## 
## Concordance= 0.626  (se = 0.025 )
## Likelihood ratio test= 19.78  on 4 df,   p=6e-04
## Wald test            = 19.75  on 4 df,   p=6e-04
## Score (logrank) test = 20.18  on 4 df,   p=5e-04

The time-dependent coefficient of Karnofsky score β(t)= -0.094+0.015×ln(t+20), and the effect value HR=exp(-0.094+0.015×ln(t+20)). For example, when the time is 200 days, the HR corresponding to the Karnofsky score is exp(-0.094+0.015×ln(200+20))= 0.987. The "+20" is added to prevent the natural logarithm from becoming undefined when t = 0.

Time-Dependent Cox Regression Model

2025-05-23