Package {cox.rvph}

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

Description

Stepwise or time-varying remedies for proportional hazards assumption violations of a cox regression model

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

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 continuous variables 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 cut points in time by maximizing the partial likelihood. The selected cut points are incorporated into the Cox model through time-dependent interval-specific covariates using the counting-process formulation. If the resulting model satisfies the PH assumption with relatively few cut points, hazard ratios (HRs) may be interpreted within each estimated time interval.

However, if many cut points are required or PH violations persist, a smooth time-varying effect may be more appropriate. In such cases, the timev method fits several candidate time-transformation functions g(t) and selects the model with the smallest AIC. By default, the following candidate functions are evaluated: linear (t), log (\log(t+1)), sqrt (\sqrt{t}), quadratic (t^2), inverse (1/(t+1)), and scaled (t/\max(t)). Users may alternatively provide their own candidate functions through the g_candidates argument.

For the selected time-varying model, the hazard ratio at time t is given by:

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

where \beta is the baseline coefficient 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
covariate_seg1 0.0813  1 0.78
covariate_seg2 3.8528  1 0.05
variable        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
covariate_seg1  0.27052   1.31065  0.07944  3.405 0.000661
covariate_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.

⁠covariate_seg1 HR = 1.31⁠

⁠covariate_seg2 HR = 1.06⁠

If the estimated cut 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] "sqrt"

$AIC
    base    linear       log      sqrt quadratic   inverse    scaled 
1519.300  1508.827  1508.895  1507.789  1512.190  1516.242  1508.827

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

                  coef exp(coef)  se(coef)      z        p
covariate      0.693216  2.000137  0.132824  5.219 1.80e-07
tt(covariate) -0.013456  0.986634  0.003959 -3.399 0.000677
variable       0.118716  1.126050  0.012860  9.231  < 2e-16
.
.

Interpretation (method: timev):

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

Example model:

HR(t) = \exp(0.693 - 0.013 \sqrt{t})

This indicates that the hazard ratio decreases gradually over time. Users may substitute clinically meaningful values of t to estimate hazard ratios at specific time points.

Value

list containing model 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$event <- ifelse(pbc$status == 2, 1, 0)
  pbc <- na.omit(pbc[, c("time","event","bili","albumin","protime","age")])
  
  cox.rvph(
    data = pbc,
    time = "time",
    event = "event",
    covariate = "protime",
    adjust_vars = c("bili","albumin","age"),
    method = "timev"
  )
}