Using cox.rvph

Introduction

The cox.rvph package provides two approaches for handling violations of the proportional hazards assumption in Cox proportional hazards models:

a stepwise approach, and
a time-varying coefficient approach.

This vignette illustrates the basic workflow for both methods.

library(cox.rvph)
library(survival)

Stepwise approach

The stepwise approach partitions the follow-up timeline into discrete intervals using one or more split points and allows the effect of the selected covariate to remain constant within each interval.

The following example uses the psych data from the KMsurv package.

if (requireNamespace("KMsurv", quietly = TRUE)) {

  data("psych", package = "KMsurv")
  psych$sex <- factor(psych$sex)

  fit_step <- cox.rvph(
    data = psych,
    time = "time",
    event = "death",
    covariate = "age",
    adjust_vars = "sex",
    method = "step",
    verbose = FALSE
  )

  print(fit_step)
  summary(fit_step)
}
#> cox.rvph fit
#> Method: step
#> Number of segments: 2 
#> Split point(s): 24.8 
#> Global PH test p-value: 0.2529169 
#> PH assumption satisfied (p > 0.05)
#> 
#> Coefficients:
#> age_seg1 age_seg2     sex2 
#>   0.2705   0.0555  -0.4791 
#> Summary of cox.rvph fit
#> Method: step
#> Number of segments: 2 
#> Split point(s): 24.8 
#> 
#> Proportional hazards test:
#> 
#>           chisq df    p
#> age_seg1 0.0813  1 0.78
#> age_seg2 3.8528  1 0.05
#> sex      0.3214  1 0.57
#> GLOBAL   4.0804  3 0.25
#> 
#> Cox proportional hazards model:
#> 
#> Call:
#> survival::coxph(formula = as.formula(f_str_final), data = final_data)
#> 
#>   n= 44, number of events= 14 
#> 
#>              coef exp(coef) se(coef)      z Pr(>|z|)    
#> age_seg1  0.27052   1.31065  0.07944  3.405 0.000661 ***
#> age_seg2  0.05550   1.05707  0.09677  0.574 0.566284    
#> sex2     -0.47909   0.61935  0.75926 -0.631 0.528044    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>          exp(coef) exp(-coef) lower .95 upper .95
#> age_seg1    1.3106      0.763    1.1217     1.531
#> age_seg2    1.0571      0.946    0.8744     1.278
#> sex2        0.6193      1.615    0.1398     2.743
#> 
#> Concordance= 0.829  (se = 0.081 )
#> Likelihood ratio test= 24.26  on 3 df,   p=2e-05
#> Wald test            = 12.27  on 3 df,   p=0.007
#> Score (logrank) test = 24.2  on 3 df,   p=2e-05

The returned object contains the selected number of time intervals, estimated split-point locations, the proportional hazards test result, and the fitted Cox model.

Time-varying coefficient approach

The time-varying coefficient approach compares candidate time functions and selects the specification with the smallest AIC.

The following example uses the pbc data from the survival package.

data("pbc", package = "survival")

pbc$status2 <- ifelse(pbc$status == 2, 1, 0)
pbc$ascites <- factor(pbc$ascites)

fit_timev <- cox.rvph(
  data = pbc,
  time = "time",
  event = "status2",
  covariate = "bili",
  adjust_vars = c("ascites", "edema", "protime"),
  method = "timev",
  verbose = FALSE
)

print(fit_timev)
#> cox.rvph fit
#> Method: time-varying coefficient
#> Selected time function: quadratic 
#> 
#> Coefficients:
#>      bili  tt(bili)  ascites1     edema   protime 
#> 8.036e-02 3.563e-08 1.153e+00 1.128e+00 2.670e-01
summary(fit_timev)
#> Summary of cox.rvph fit
#> Method: time-varying coefficient
#> Selected time function: quadratic 
#> 
#> AIC values:
#>       Model      AIC
#> 1      base 1151.107
#> 2    linear 1139.675
#> 3       log 1147.434
#> 4      sqrt 1142.833
#> 5 quadratic 1138.358
#> 6   inverse 1152.789
#> 
#> Cox proportional hazards model:
#> 
#> Call:
#> survival::coxph(formula = as.formula(f_str), data = data, tt = tt_fun)
#> 
#>   n= 312, number of events= 125 
#> 
#>               coef exp(coef)  se(coef)     z Pr(>|z|)    
#> bili     8.036e-02 1.084e+00 1.807e-02 4.448 8.66e-06 ***
#> tt(bili) 3.563e-08 1.000e+00 8.323e-09 4.281 1.86e-05 ***
#> ascites1 1.153e+00 3.168e+00 2.746e-01 4.199 2.68e-05 ***
#> edema    1.128e+00 3.089e+00 2.968e-01 3.800 0.000145 ***
#> protime  2.670e-01 1.306e+00 7.679e-02 3.477 0.000507 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>          exp(coef) exp(-coef) lower .95 upper .95
#> bili         1.084     0.9228     1.046     1.123
#> tt(bili)     1.000     1.0000     1.000     1.000
#> ascites1     3.168     0.3157     1.849     5.427
#> edema        3.089     0.3237     1.727     5.527
#> protime      1.306     0.7657     1.124     1.518
#> 
#> Concordance= 0.819  (se = 0.02 )
#> Likelihood ratio test= 151.6  on 5 df,   p=<2e-16
#> Wald test            = 132.9  on 5 df,   p=<2e-16
#> Score (logrank) test = 310.7  on 5 df,   p=<2e-16

The returned object contains the selected time function, AIC values for the candidate time-function specifications, and the fitted Cox model.

Using cox.rvph

Hamin Kim

Introduction

Stepwise approach

Time-varying coefficient approach

Summary