Restricted mean survival time

This article covers ertte_rmst(), the more complicated of the two tools for scalar exposure-response (E-R) reduction of a time-to-event (TTE) endpoint. It’s the longer sibling of the landmark event probability article. Both of these tools are used to turn a full survival curve into a single number that can be plotted with the standard erplots grammar. If you haven’t already, start with the overview article for background on ertte_aft()/ ertte_coxph().

library(ertte)
library(survival)

Why RMST?

Most exposure-response summaries for a time-to-event (TTE) endpoint focus on either a hazard ratio (from a Cox model) or a survival probability at a fixed time (from ertte_predict()/ertte_landmark()). Both are useful, but neither directly answers a question that’s often more clinically intuitive: on average, how much event-free time does a subject with a given exposure profile experience, up to some horizon?

Restricted mean survival time (RMST) answers exactly that. For a survival function \(S(t \mid x)\) describing a subject with covariate profile \(x\), and a fixed horizon \(\tau\), RMST is defined as

\[ \text{RMST}(\tau \mid x) = \int_0^\tau S(t \mid x) \, dt. \]

Geometrically, this is the area under the survival curve between 0 and \(\tau\). It has units of time (e.g. “expected days event-free in the first 90 days”), doesn’t rely on a proportional-hazards assumption to interpret, and stays well-defined even when a study’s follow-up is too short to observe the survival curve reach zero – which is exactly why it’s “restricted” to \(\tau\) rather than defined over \([0, \infty)\).

ertte_rmst() computes RMST (with a confidence interval) for a fitted ertte_aft() or ertte_coxph() model, at one or more horizons tau, for one or more covariate profiles in newdata. This article works through how it does that for each model type, what assumptions the calculation relies on, and where the confidence interval is more approximate than you might expect.

This assumes you already have some familiarity with fitting TTE models via the survival package (survreg()/coxph()) and with ertte’s two model engines; if you’re new to ertte, start with ertte_aft()/ertte_coxph()’s documentation and ertte_predict() first.

Two engines, two survival curves

ertte fits two kinds of TTE models, and ertte_rmst() needs a genuinely different calculation for each, because they represent \(S(t \mid x)\) differently.

ertte_aft() fits a parametric accelerated failure time (AFT) model via survival::survreg(). All four supported distributions share the same log-location-scale structure,

\[ \log(T) = \mu(x) + \sigma W, \]

where \(\mu(x)\) is the linear predictor (intercept + covariate effects), \(\sigma\) is a scale parameter, and \(W\) follows a fixed “base” distribution (extreme-value for the Weibull/exponential case, standard normal for log-normal, standard logistic for log-logistic). This gives a fully parametric, closed-form survival function \(S(t \mid x) = 1 - F\!\left( \frac{\log t - \mu(x)}{\sigma}\right)\), where \(F\) is the base distribution’s CDF.

ertte_coxph() fits a semi-parametric Cox proportional-hazards model via survival::coxph(). There’s no parametric form for the baseline hazard \(h_0(t)\); it’s estimated nonparametrically (Breslow’s estimator, the same one behind survival::basehaz()), and the survival curve for a covariate profile \(x\) is

\[ S(t \mid x) = S_0(t)^{\exp\left((x - \bar x)'\beta\right)}, \]

where \(S_0(t) = \exp(-H_0(t))\) is the fitted baseline survival curve and \(\bar x\) is the mean covariate profile the partial likelihood was centered on. Because \(H_0(t)\) is estimated from a finite set of observed event times, \(S(t \mid x)\) is a right-continuous step function – it only changes value at the observed event times in the fitting data, however far \(x\) is from \(\bar x\).

That structural difference – closed-form vs. step-function – is what drives most of the implementation differences below.

Computing the point estimate

AFT: numerical integration of a known curve

Since \(S(t \mid x)\) has a closed form for the AFT engine, ertte_rmst.ertte_aft() computes

\[ \widehat{\text{RMST}}(\tau \mid x) = \int_0^\tau \hat S(t \mid x) \, dt \]

via stats::integrate(). This is straightforward numerical quadrature on a smooth, fully known function – there’s no approximation of consequence here beyond ordinary quadrature error, which is negligible for a function this well-behaved.

Cox: an exact sum, not an approximation

For the Cox engine, it’s tempting to assume the same kind of numerical integration is needed, since the “curve” comes from an empirical baseline hazard rather than a formula. In fact the opposite is true: because \(\hat S(t \mid x)\) is a step function with jumps only at the observed event times \(t_{(1)} < t_{(2)} < \dots < t_{(m)}\), the area under it between two jumps is exactly a rectangle, and the whole integral reduces to an exact finite sum:

\[ \widehat{\text{RMST}}(\tau \mid x) = \sum_{k} \hat S\!\left(t_{(k)}^{-} \mid x\right) \cdot \Big(\min(t_{(k+1)}, \tau) - t_{(k)}\Big). \]

ertte_rmst.ertte_coxph() computes this directly from survival::survfit(object, newdata)’s own time/surv fields – there’s no stats::integrate() call anywhere in this path, and no quadrature error to worry about. The only genuine approximation in the point estimate shows up at the boundary, discussed under “Extrapolating beyond the data” below.

Confidence intervals: the harder half

Getting a point estimate is the easy part for both engines. Getting a defensible confidence interval turned out to need real design decisions, and the two engines’ methods are not analogous to each other.

AFT: a delta method that ignores scale uncertainty

ertte_predict.ertte_aft() already reports a Wald confidence interval on the survival probability by propagating the uncertainty in \(\mu(x)\) (the linear predictor) and not the uncertainty in \(\sigma\) (the scale parameter) – a simplification made throughout ertte, on the grounds that \(\sigma\)’s sampling variability is typically small relative to \(\mu(x)\)’s and that jointly propagating both would need a full delta-method gradient rather than a simple predict(..., se.fit = TRUE) call. ertte_rmst.ertte_aft() makes the identical simplification.

Differentiating the RMST integral with respect to \(\mu\) (differentiating under the integral sign) gives

\[ \frac{\partial}{\partial \mu} \text{RMST}(\tau \mid x) = \int_0^\tau \frac{f\!\left(\frac{\log t - \mu(x)}{\sigma}\right)}{\sigma} \, dt, \]

where \(f\) is the base distribution’s density. ertte_rmst.ertte_aft() evaluates this gradient numerically (another stats::integrate() call), then applies the standard delta-method variance,

\[ \widehat{\text{Var}}\!\left[\widehat{\text{RMST}}(\tau \mid x)\right] \approx \left(\frac{\partial \text{RMST}}{\partial \mu}\right)^{\!2} \widehat{\text{Var}}(\hat\mu), \]

with \(\widehat{\text{Var}}(\hat\mu)\) coming from predict(object, newdata, type = "linear", se.fit = TRUE), exactly as in ertte_predict.ertte_aft(). The resulting interval is a symmetric Wald interval, fit_rmst \(\pm\ z \cdot\) se_rmst.

Cox: a corrected delta method

For the Cox engine, the natural first idea is to reuse machinery survival already ships: print(survfit_object, rmean = tau) prints an rmean/se(rmean) pair for a fitted survfit object, and does run on survival::survfit(coxph_object, newdata = x). It’s the wrong standard error to use here, though: it’s built from n.risk/n.event, which are properties of the shared baseline hazard estimate, identical across any covariate profile passed through the same fitted model – not profile-specific uncertainty at all. For a covariate profile far from the average (e.g. an unusually high exposure value), that naive standard error can be dramatically too small, understating uncertainty for exactly the profiles an exposure-response analysis often cares most about (e.g. a high-exposure arm).

ertte_rmst.ertte_coxph() instead builds on survfit()’s own std.err field, which does vary correctly by covariate profile: it’s the estimated standard error of \(\hat H(t \mid x)\), the profile-specific cumulative hazard, coming from the same delta-method calculation survfit.coxph() uses internally for its own per-time confidence intervals. The construction otherwise matches the classic Greenwood-type sum – squared “area remaining beyond \(t_k\)” terms, each weighted by a variance increment at \(t_k\) – but with the increment itself replaced:

\[ \hat h_k = \widehat{\text{Var}}\!\left[\hat H(t_{(k)} \mid x)\right] - \widehat{\text{Var}}\!\left[\hat H(t_{(k-1)} \mid x)\right] = \text{std.err}(t_{(k)})^2 - \text{std.err}(t_{(k-1)})^2. \]

This is a genuine improvement, but still an approximation: it treats \(\hat H(t \mid x)\) as accumulating via statistically independent increments over time, which holds for the nonparametric part of the hazard estimator but not exactly for the part of its uncertainty coming from the regression coefficients \(\hat\beta\) – that’s really one shared random quantity affecting every time point together, not something that accumulates independently as \(t\) increases. In development testing against a nonparametric bootstrap for a couple of contrasting covariate profiles, this corrected standard error tracked the bootstrap substantially more closely than the naive survmean()-based version – but that check wasn’t exhaustive across sample sizes, censoring patterns, or covariate profiles, so treat it as a meaningfully better approximation, not a proven-exact one.

Extrapolating beyond the data

Both engines can, in principle, be asked for RMST at a horizon \(\tau\) that exceeds the longest follow-up time observed in the fitting data.

For ertte_coxph(), the baseline hazard has no information past the last observed event/censoring time. ertte_rmst.ertte_coxph() follows the same convention as ertte_predict.ertte_coxph(): survival is held flat at its last estimated value beyond that point. Because RMST integrates the entire curve up to \(\tau\), this flat-tail assumption has a much larger effect on an RMST value than it does on a single survival-probability prediction at one time point – a long flat tail can inflate the estimated area substantially. ertte_rmst() issues a warning whenever any requested tau exceeds the observed follow-up range, specifically because this matters more here than elsewhere in the package.
For ertte_aft(), there’s no equivalent hard boundary – the parametric survival curve is defined for all \(t > 0\) – but that comes with its own caveat: far beyond the observed data, the estimate is governed entirely by the assumed parametric family, with no data left to check whether that assumption still holds.

A worked example

mod_aft <- ertte_aft(Surv(time, event) ~ aucss, ertte_data)
mod_cox <- ertte_coxph(Surv(time, event) ~ aucss, ertte_data)

nd <- ertte_data[c(1, 50, 100), ]
nd[, c("id", "aucss")]
#>      id    aucss
#> 1     1 1114.089
#> 50   50  148.240
#> 100 100  427.802

ertte_rmst(mod_aft, nd, tau = c(60, 90, 120))
#> # A tibble: 9 × 15
#>      id sex     age weight  dose treatment aucss cmaxss  time event admin_censor
#>   <int> <fct> <int>  <dbl> <dbl> <fct>     <dbl>  <dbl> <dbl> <dbl>        <dbl>
#> 1     1 Fema…    27     70   200 Drug      1114.  187.   77.4     1          180
#> 2     1 Fema…    27     70   200 Drug      1114.  187.   77.4     1          180
#> 3     1 Fema…    27     70   200 Drug      1114.  187.   77.4     1          180
#> 4    50 Male     28     58   100 Drug       148.   10.4   0.3     0          180
#> 5    50 Male     28     58   100 Drug       148.   10.4   0.3     0          180
#> 6    50 Male     28     58   100 Drug       148.   10.4   0.3     0          180
#> 7   100 Fema…    26     48   200 Drug       428.   35.1 166.      0          180
#> 8   100 Fema…    26     48   200 Drug       428.   35.1 166.      0          180
#> 9   100 Fema…    26     48   200 Drug       428.   35.1 166.      0          180
#> # ℹ 4 more variables: tau <dbl>, fit_rmst <dbl>, ci_lower <dbl>, ci_upper <dbl>

ertte_rmst(mod_cox, nd, tau = c(60, 90, 120))
#> # A tibble: 9 × 15
#>      id sex     age weight  dose treatment aucss cmaxss  time event admin_censor
#>   <int> <fct> <int>  <dbl> <dbl> <fct>     <dbl>  <dbl> <dbl> <dbl>        <dbl>
#> 1     1 Fema…    27     70   200 Drug      1114.  187.   77.4     1          180
#> 2     1 Fema…    27     70   200 Drug      1114.  187.   77.4     1          180
#> 3     1 Fema…    27     70   200 Drug      1114.  187.   77.4     1          180
#> 4    50 Male     28     58   100 Drug       148.   10.4   0.3     0          180
#> 5    50 Male     28     58   100 Drug       148.   10.4   0.3     0          180
#> 6    50 Male     28     58   100 Drug       148.   10.4   0.3     0          180
#> 7   100 Fema…    26     48   200 Drug       428.   35.1 166.      0          180
#> 8   100 Fema…    26     48   200 Drug       428.   35.1 166.      0          180
#> 9   100 Fema…    26     48   200 Drug       428.   35.1 166.      0          180
#> # ℹ 4 more variables: tau <dbl>, fit_rmst <dbl>, ci_lower <dbl>, ci_upper <dbl>

The two engines’ point estimates and intervals should generally be in the same ballpark when both models fit the data reasonably well (as here, since ertte_data was simulated from a Weibull AFT ground truth that a Cox model should also approximate well) – but they are not expected to match exactly, and a persistent, large disagreement between the two is a useful diagnostic that one of the models may not fit the data well.

mod_cox_extrap <- ertte_rmst(mod_cox, nd[1, ], tau = 250)
#> Warning: `tau` exceeds the last observed follow-up time (180) for at least one
#> value. RMST beyond that point assumes survival stays flat at its last estimated
#> value, which may understate/overstate the true area.

The warning above fires because 250 exceeds the longest observed follow-up time in ertte_data; the returned value still assumes survival stays flat beyond that point.

Using RMST in a visual predictive check

ertte’s erplots integration also supports RMST-based visual predictive checks (VPCs): passing tau through simulate_args to erplots::er_vpc_add_simulated() reduces each simulated replicate to a per-subject restricted survival time, min(sim_time, tau), then compares its distribution against the observed data binned by exposure. This reuses the already-censored simulated event times – the same administrative- censoring convention described above – rather than the model’s raw, uncensored draws, since a VPC’s purpose is checking whether the model reproduces what you’d actually observe under the study’s real censoring pattern.

A consequence worth knowing about: a simulated replicate that’s censored before tau has a genuinely ambiguous restricted survival time (we don’t know what would have happened between the censoring time and tau), and is dropped (returned as NA) rather than imputed or reweighted. This is a complete-case convention, applied identically to the observed and simulated sides of the comparison – and, unlike a naive complete-case point estimate, isn’t obviously the wrong choice here, since it’s the observed study’s actual censoring pattern driving both sides equally. A more sophisticated correction (inverse-probability-of-censoring weighting, or pseudo-observations) was considered but isn’t offered: erplots’ VPC aggregation is currently an unweighted mean, so there’s nowhere for an IPCW-style weight to be applied downstream, and a pseudo-observations approach would need an expensive leave-one-out recomputation for every simulated replicate. If this convention drops an uncomfortably large fraction of replicates in a given analysis, supplying censor_time (a per-subject administrative follow-up horizon) to simulate_args changes how much ambiguity arises in the first place, rather than trying to correct for it after the fact.

Summary of assumptions and limitations

AFT confidence intervals ignore uncertainty in the scale parameter, propagating only the linear predictor’s uncertainty. This matches ertte_predict.ertte_aft()’s existing convention, not a special weakening for RMST specifically.
Cox confidence intervals rely on an independent-increments approximation for how coefficient uncertainty accumulates over time, which is not an exact result. It has been checked against a bootstrap for a small number of covariate profiles and found to track it substantially better than the naive alternative, but this is not exhaustive validation.
Neither engine’s interval is clipped to [0, tau]. Since a symmetric Wald interval isn’t automatically bounded the way a back-transformed survival-probability interval is, ci_lower/ci_upper can in principle fall outside [0, tau] for small samples or profiles near the edge of the observed covariate range. ertte_rmst() does not currently clip this.
RMST beyond the observed follow-up range assumes a flat survival tail (Cox engine) or relies entirely on extrapolating the fitted parametric family (AFT engine). Treat tau values near or beyond the longest observed follow-up time with real caution, and note the warning ertte_rmst() issues for the Cox engine when this happens.
RMST depends on the choice of tau. Unlike a hazard ratio, RMST is not a single summary of the whole exposure-response relationship – it’s specific to the chosen horizon, and different horizons can tell different (both valid) stories about the same fitted model.
RMST-based VPC simulations drop replicates censored before tau as ambiguous (complete-case), with no IPCW/pseudo-observation alternative currently offered – see “Using RMST in a visual predictive check” above for why, and for the censor_time mitigation.