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 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 compares the base Cox model with time-varying coefficient models based on several candidate time-transformation functions using AIC, 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 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
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 
30935.70  30917.65  30923.49  30916.70  30922.30  30935.77  30917.65

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

                     coef exp(coef)  se(coef)      z        p
covariate       0.0737807 1.0765707 0.0057407 12.852  < 2e-16
tt(covariate)   0.0005973 1.0005974 0.0001298  4.602 4.19e-06
variable        0.1640114 1.1782278 0.0194412  8.436  < 2e-16
.
.

Interpretation (method: timev):

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

Example model:

HR(t) = \exp(0.0737807 + 0.0005973 \times \sqrt{t})

In this example, the square-root time transformation indicates that the hazard ratio 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"
  )
}

## Not run: 
  
  data(flchain, package = "survival")
  
  cox.rvph(
    data = flchain,
    time = "futime",
    event = "death",
    covariate = "age",
    adjust_vars = c("lambda", "flc.grp", "creatinine"),
    method = "timev"
  )

## End(Not run)