Package {cox.rvph}

cox.rvph (Remedy for Violations of the Proportional Hazards Assumption in Cox Proportional Hazards Models)

Description

Stepwise or time-varying remedies for proportional hazards assumption violations in Cox Proportional Hazards Models

Usage

cox.rvph(
  data,
  time,
  event,
  covariate,
  adjust_vars = NULL,
  method = c("step", "timev"),
  g_candidates = NULL,
  max_K = 4,
  p_threshold = 0.05,
  verbose = TRUE
)

Arguments

Details

The event indicator specified by event must be binary, with 0 indicating a right-censored observation and 1 indicating that the event of interest occurred. Other event codings must be recoded before using cox.rvph().

Additionally, users should verify that the survival time variable does not contain negative values. For the step method, observations with a survival time of 0 should be handled using an appropriate preprocessing strategy based on the scientific context. For the timev method, the specified time-transformation functions should be defined over the observed range of survival times.

Users should first assess the proportional hazards (PH) assumption using cox.zph() before applying cox.rvph(). Variables showing evidence of non-proportional hazards may then be modeled using stepwise or time-varying remedies.

Currently, only a single continuous variable can be specified in covariate. Categorical variables may still be included in adjust_vars as adjustment covariates.

The step method performs segmented modeling by searching for optimal split points in time by maximizing the partial likelihood. The selected split points are incorporated into the Cox model through time-dependent interval-specific covariates using the counting-process formulation. For each number of segments, the PH assumption is assessed using the global test from cox.zph(). The search stops at the first model for which the global p-value exceeds p_threshold. If the resulting model satisfies the PH assumption with relatively few split points, hazard ratios (HRs) may be interpreted within each estimated time interval.

The resulting hazard ratio within each time interval k is given by:

HR(t) = \exp(\beta_k), \quad t \in I_k

However, if many split points are required or PH violations persist, a smooth time-varying effect may be more appropriate. In such cases, the timev method compares the base Cox model with several time-varying alternatives based on candidate time-transformation functions using AIC, and selects the one with the smallest AIC. By default, the following candidate functions are evaluated: linear (t), logarithmic (log) (\log(t+1)), square-root (sqrt) (\sqrt{t}), quadratic (t^2) and inverse (1/(t+1)). Users may alternatively provide their own candidate functions through the g_candidates argument.

The resulting hazard ratio at time t is given by:

HR(t) = \exp(\beta + \gamma g(t))

where \beta is the coefficient of the covariate and \gamma represents the time-varying interaction effect. Users may evaluate this expression at clinically relevant time points to interpret how the hazard ratio changes over time.

Example output (method: step):

$K
[1] 2

$tau
[1] 24.8

$p_value

[1] 0.2529169

$zph
                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

$fit
Call:
survival::coxph(formula = as.formula(f_str_final), data = final_data)

              coef exp(coef) se(coef)     z        p
age_seg1  0.27052   1.31065  0.07944  3.405 0.000661
age_seg2  0.05550   1.05707  0.09677  0.574 0.566284
.
.

Interpretation (method: step):

For the step method, hazard ratios are interpreted separately within each estimated time interval.

HR(t) = \begin{cases} \exp(0.27052) = 1.31, & t \le 24.8, \\ \exp(0.05550) = 1.06, & t > 24.8. \end{cases}

If the estimated split point is \tau = 24.8, this indicates that the hazard ratio associated with a one-unit increase in age is approximately 1.31 before time 24.8 and decreases to approximately 1.06 after time 24.8.

Example output (method: timev):

$selected_g
[1] "quadratic"

$AIC
    base    linear       log      sqrt quadratic   inverse    
1151.107  1139.675  1147.434  1142.833  1138.358  1152.789  

$fit
Call:
survival::coxph(formula = as.formula(f_str), data = data, tt = tt_fun)

               coef exp(coef)  se(coef)     z        p
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
.
.

Interpretation (method: timev):

For the timev method, hazard ratios vary continuously over time.

HR(t) = \exp(0.08036 + 3.563 \times 10^{-8} t^2)

In this example, the quadratic time transformation indicates that the effect of bili changes continuously over time. Users may substitute clinically meaningful values of t to estimate hazard ratios at specific time points.

Value

a list containing the fitted model and related results

Examples

if (requireNamespace("KMsurv", quietly = TRUE)) {
  data(psych, package = "KMsurv")
  
  psych$sex <- factor(psych$sex)

  cox.rvph(
    data = psych,
    time = "time",
    event = "death",
    covariate = "age",
    adjust_vars = "sex",
    method = "step"
  )
}

if (requireNamespace("survival", quietly = TRUE)) {
  data(pbc, package = "survival")
  
  pbc$status2 <- ifelse(pbc$status == 2, 1, 0)
  pbc$ascites <- factor(pbc$ascites)
  
  cox.rvph(
    data = pbc,
    time = "time",
    event = "status2",
    covariate = "bili",
    adjust_vars = c("ascites", "edema", "protime"),
    method = "timev"
  )
}