Spline basis

Let \(x = \log t\). A restricted (natural) cubic spline with boundary knots \(k_{\min}, k_{\max}\) and interior knots \(k_1, \dots, k_m\) is written, in the Durrleman & Simon / Royston & Parmar parameterisation, as

\[ s(x) = \gamma_0 + \gamma_1 x + \sum_{j=1}^{m} \gamma_{1+j}\, v_j(x), \]

where the linear term keeps the domain \(x \in \mathbb{R}\) and each nonlinear basis term is

\[ v_j(x) = (x - k_j)_+^3 - \lambda_j (x - k_{\min})_+^3 - (1-\lambda_j)(x - k_{\max})_+^3, \qquad \lambda_j = \frac{k_{\max} - k_j}{k_{\max} - k_{\min}}, \]

with \((u)_+ = \max(u, 0)\). This is exactly the basis in rcs_basis() (R/splines.R); by construction \(s(x)\) is linear beyond the boundary knots, which keeps extrapolation well-behaved. Its derivative, needed for the hazard, is

\[ v_j'(x) = 3(x-k_j)_+^2 - 3\lambda_j (x-k_{\min})_+^2 - 3(1-\lambda_j)(x-k_{\max})_+^2 . \]

Boundary knots default to the min/max of \(\log t\) among events, and interior knots to equally spaced centiles of \(\log t\) among events (default_knots()), matching rstpm2/flexsurv defaults.

Linear predictor and scales

Covariates enter proportionally on the modelled scale:

\[ \eta(t \mid x) = s(\log t) + \beta^\top x . \]

rpsurv supports three transformations \(g\) of the survival function, selected via scale:

The hazard follows from differentiating \(S(t)=S(\eta(t))\) with respect to \(t\) via the chain rule through \(x=\log t\):

\[ h(t) = \frac{f(t)}{S(t)} = -\frac{d\eta}{d\log t}\cdot\frac{1}{t}\cdot\frac{S'(\eta)}{S(\eta)} . \]

Writing \(g(\eta) = \log h(t) - \log\!\big(\tfrac{d\eta}{d\log t}\big) + \log t\) (a function of \(\eta\) alone, on each scale), the three cases reduce to closed forms:

\[ g(\eta) = \begin{cases} \eta & \text{hazard scale} \\ \eta - \log(1+e^{\eta}) & \text{odds scale} \\ \log\phi(\eta) - \log\{1-\Phi(\eta)\} & \text{probit scale} \end{cases} \]

so that \(\log h(t) = \log\!\big(\tfrac{d\eta}{d\log t}\big) - \log t + g(\eta)\). This is exactly scale_terms() in src/rpsurv_loglik.cpp.

Log-likelihood

For a right-censored observation with event time \(t_i\), status \(d_i\), and (for left-truncated / counting-process data) entry time \(e_i\), the contribution is

\[ \ell_i(\beta) = d_i \log h(t_i) + \log S(t_i) - \log S(e_i), \]

with \(\log S(e_i) \equiv 0\) when \(e_i = 0\) (no truncation). Summing gives the full log-likelihood \(\ell(\beta) = \sum_i \ell_i(\beta)\), maximised by stats::optim(method = "BFGS") using the analytic gradient below.

Analytic gradient

Because \(\eta\) and \(\tfrac{d\eta}{d\log t}\) are both linear in \(\beta\) (through the design matrices \(X\) and \(\dot X = dX/d\log t\)), the gradient of \(\ell_i\) with respect to \(\beta_k\) is

\[ \frac{\partial \ell_i}{\partial \beta_k} = d_i\left[\frac{\dot X_{ik}}{\dot\eta_i} + g'(\eta_i) X_{ik}\right] + \frac{d\log S}{d\eta}(\eta_i)\, X_{ik} - \mathbb{1}[e_i>0]\,\frac{d\log S}{d\eta}(\eta_{e_i})\, X^{(e)}_{ik}, \]

where \(X^{(e)}\) is the design row evaluated at \(\log e_i\). The scale-specific pieces \(g'(\eta)\) and \(d\log S/d\eta\) are, again, closed form (hazard: \(g'=1\), \(d\log S/d\eta = -e^\eta\); odds: \(g'=S\), \(d\log S/d\eta=-(1-S)\); probit: \(g'=-\eta+\phi(\eta)/S\), \(d\log S/d\eta=-\phi(\eta)/S\)). Evaluating \(\ell\) and \(\nabla\ell\) this way, rather than by finite differences, avoids the \(O(p)\) extra likelihood evaluations per optimiser step that a numerical gradient would cost, which is the main source of rpsurv’s speed advantage together with parallel evaluation across observations (RcppParallel::parallelReduce, src/rpsurv_loglik.cpp).

An observation with \(d_i=1\) and \(\dot\eta_i \le 0\) (a locally decreasing, i.e. invalid, hazard) is given a large penalty rather than NaN, likewise any row where \(S(e_i) < S(t_i)\) (non-monotone survival across a truncation interval); both keep the optimiser inside the region where the model is a valid survival distribution.

Time-varying effect (non-proportional hazards)

A covariate’s value is fixed for a subject, but its coefficient \(\beta(t)\) is allowed to change over follow-up: instead of \(\beta x\), the linear predictor gets \(x \cdot s_{\text{tve}}(\log t)\), i.e. its own spline in log time multiplying the covariate. Request this with tve:

fit_tve <- rpsurv(Surv(rectime, censrec) ~ hormon, data = brcancer,
                   df = 4, tve = "hormon", tve.df = 2)
coef(fit_tve)
#>     (Intercept)        s(logt)1        s(logt)2        s(logt)3        s(logt)4 
#>    -19.19777720      2.93098413     -0.42391914      1.32597578     -0.87854493 
#>          hormon hormon:s(logt)1 hormon:s(logt)2 
#>     -2.34822990      0.32187562      0.03601483

The hormon:s(logt)* coefficients trace out how the hormonal-therapy hazard ratio evolves with time; a likelihood-ratio test against fit (both fitted by maximum likelihood, nested models) tests proportionality:

lrt_stat <- 2 * (fit_tve$loglik - fit$loglik)
lrt_df <- fit_tve$df - fit$df
pchisq(lrt_stat, lrt_df, lower.tail = FALSE)
#> [1] 0.8122301

Time-varying covariate

Here the covariate’s value itself changes during follow-up (e.g. a biomarker updated at clinic visits). This has nothing to do with tve; it is a data-representation question. Split each subject’s follow-up into consecutive intervals over which the covariate is constant, and supply counting-process data via Surv(start, stop, status):

d <- brcancer
d$id <- seq_len(nrow(d))
half <- d$rectime / 2
first  <- data.frame(id = d$id, start = 0, stop = half, status = 0,
                      hormon = d$hormon, x1 = d$x1)
second <- data.frame(id = d$id, start = half, stop = d$rectime, status = d$censrec,
                      hormon = d$hormon, x1 = d$x1 + 1)  # covariate value changes
long <- rbind(first, second)
long <- long[long$start < long$stop, ]

fit_cp <- rpsurv(Surv(start, stop, status) ~ hormon + x1, data = long, df = 4)
fit_cp$counting
#> [1] TRUE
coef(fit_cp)
#>   (Intercept)      s(logt)1      s(logt)2      s(logt)3      s(logt)4 
#> -19.826159484   2.991760891  -0.438400286   1.348304243  -0.878671736 
#>        hormon            x1 
#>  -0.385157266   0.004545335

rpsurv detects the 3-column Surv object and fits the left-truncated likelihood of Section 2 automatically: each interval contributes \(\log S(\text{stop}) - \log S(\text{start})\), correctly accounting for the subject having already survived to start under the covariate value(s) of the previous interval. The same mechanism (Surv(start, stop, status)) also handles ordinary left truncation (delayed entry) when the covariate values do not change.

What is delegated to C++

Everything on the hot path of fitting is C++, not just the likelihood sum:

• the log-likelihood and its analytic gradient are computed in a single parallel pass (RcppParallel::parallelReduce) over observations, rather than two separate passes;
• optim(method = "BFGS") calls the objective (fn) and the gradient (gr) separately at every trial point, which would otherwise mean two full data passes per iteration for the same point; rpsurv caches the fused C++ result keyed on the parameter vector, so each distinct point is evaluated exactly once;
• the restricted cubic spline basis itself (rcs_basis(), both the basis and its derivative) is computed by rcs_basis_cpp() in C++ rather than with R’s vectorised-but-still-interpreted arithmetic;
• the closed-form starting values (rp_start_values()) are computed from a bounded subsample (default 20,000 rows) of the Kaplan-Meier curve; they only need to be in the right neighbourhood for BFGS, so this avoids an \(O(n\log n)\) sort of the full data purely for initialisation.

What remains in R (model.frame()/model.matrix() formula parsing, one cbind() to assemble the design matrix) runs once per fit, not once per optimiser iteration, and is a small fraction of total time at the sample sizes below.

Benchmark

The figure and table compare wall-clock fit time against both canonical flexible parametric survival packages, rstpm2::stpm2 and flexsurv::flexsurvspline, on simulated Weibull-hazard data with two covariates (reproducible via data-raw/benchmark.R, shipped with the package source). flexsurv::flexsurvspline is materially slower per fit, so it is only run up to \(n=2\times10^4\) to keep the benchmark tractable; rstpm2::stpm2 is run up to \(n=2\times10^5\).

Fit time vs. sample size (log-log axes).

At \(n=2\times10^5\), rpsurv is about 11.9x faster than rstpm2::stpm2, and continues to scale to \(n=5\times10^5\) (rstpm2 was not run at that size). Against flexsurv::flexsurvspline, at \(n=2\times10^4\) rpsurv is about 29.1x faster. The advantage widens with \(n\) in both comparisons, consistent with rpsurv doing \(O(n)\) work per optimiser step against a single parallel pass, versus repeated R-level likelihood evaluation (numerical gradients, or gradients built from R-level matrix algebra) in the other two packages. The benchmark code itself:

simulate_data <- function(n, seed = 1) {
  set.seed(seed)
  x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5)
  lp <- 0.5 * x1 - 0.3 * x2
  u <- runif(n)
  event_time <- (-log(u) / (0.01 * exp(lp)))^(1 / 1.2)
  cens_time <- rexp(n, 0.008)
  data.frame(time = pmin(event_time, cens_time),
             status = as.numeric(event_time <= cens_time), x1 = x1, x2 = x2)
}
d <- simulate_data(200000)
system.time(rpsurv(Surv(time, status) ~ x1 + x2, data = d, df = 4))
system.time(stpm2(Surv(time, status) ~ x1 + x2, data = d, df = 4))
system.time(flexsurvspline(Surv(time, status) ~ x1 + x2, data = d, k = 3))