cox.zph to Test the Proportional Hazards Assumption of a Cox Regression

1. Using cox.zph to test and visualize

## Load surival package
library(survival)

## Load Ovarian Cancer Survival Data
data(ovarian)

## Fit Cox regression model
fit <- coxph(Surv(futime, fustat) ~ age + ecog.ps,
              data = ovarian)

## Test the Proportional Hazards Assumption of a Cox Regression
temp <- cox.zph(fit)

## Print test results
temp

           rho chisq     p
age     -0.243 0.856 0.355
ecog.ps  0.520 2.545 0.111
GLOBAL      NA 3.195 0.202

## Show plots
plot(temp)

2. Using extended Cox model to look for interaction with time

See if the hazard ratio changes before and after 200 days.

## split follow-up interval at 200 days
ovarian.split <- survSplit(data    = ovarian,
                           cut     = 200,
                           end     = "futime",
                           event   = "fustat",
                           start   = "start",
                           id      = "id",
                           zero    = 0,
                           episode = NULL)

## Create interval indicator and survival object 
ovarian.split <- within(ovarian.split, {

    ## Interval indicator
    interval.indicator <- factor(start, c(0,200), c("<=200",">200"))

    ## Surv(start, stop, event) format
    SurvObj <- Surv(start, futime, fustat)
 })


## Fit model with interaction. Use cluster(id) to account for clustering by subject
fit2 <- coxph(SurvObj ~ age + ecog.ps + (age + ecog.ps):interval.indicator + cluster(id),
             data = ovarian.split)

## Show results: Both interaction terms are not significant.
summary(fit2)

Call:
coxph(formula = SurvObj ~ age + ecog.ps + (age + ecog.ps):interval.indicator + 
    cluster(id), data = ovarian.split)

  n= 49, number of events= 12 

                                  coef exp(coef) se(coef) robust se     z  Pr(>|z|)    
age                             0.2573    1.2934   0.1154    0.0507  5.07 0.0000004 ***
ecog.ps                        -1.3896    0.2492   1.4418    1.0538 -1.32      0.19    
age:interval.indicator>200     -0.1364    0.8725   0.1301    0.0934 -1.46      0.14    
ecog.ps:interval.indicator>200  1.8198    6.1707   1.6108    1.2528  1.45      0.15    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

                               exp(coef) exp(-coef) lower .95 upper .95
age                                1.293      0.773    1.1710      1.43
ecog.ps                            0.249      4.013    0.0316      1.97
age:interval.indicator>200         0.872      1.146    0.7266      1.05
ecog.ps:interval.indicator>200     6.171      0.162    0.5296     71.89

Concordance= 0.803  (se = 0.091 )
Rsquare= 0.289   (max possible= 0.76 )
Likelihood ratio test= 16.7  on 4 df,   p=0.00218
Wald test            = 33.3  on 4 df,   p=0.00000103
Score (logrank) test = 15.1  on 4 df,   p=0.00455,   Robust = 6.97  p=0.138

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).