Package {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 from print() (method: 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

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 from print() (method: 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

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

An object of class "cox.rvph" containing the fitted model and method-specific results.

Examples

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"
  )

  print(fit_step)
  summary(fit_step)
}

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"
)

print(fit_timev)
summary(fit_timev)

Print a cox.rvph Object

Description

Print a cox.rvph Object

Usage

## S3 method for class 'cox.rvph'
print(x, ...)

Arguments

Summarize a cox.rvph Object

Description

Summarize a cox.rvph Object

Usage

## S3 method for class 'cox.rvph'
summary(object, ...)

Arguments

Value

The input "cox.rvph" object, returned invisibly.