Chapter slides here. (To convert html to pdf, press E \(\to\) Print \(\to\) Destination: Save to pdf)
R-code
Show the code
################################################################### This code generates the numerical results in chapter 3 #################################################################### install.packages("rmt") # if not already installedlibrary(rmt)library(tidyverse)##### Read in HF-ACTION DATA ########data(hfaction)head(hfaction)#> Displays the first few rows of the HF-ACTION dataset (patid, time, status, trt_ab)# TFE: take the first event per patient id# TFE = Time to First Event (death or hospitalization)hfaction_TFE <- hfaction |>arrange(patid, time) |>group_by(patid) |>slice_head() |>ungroup()###### Standard RMST analysis #################library(survRM2)# -------------------------------------------------------------------------# Mortality analysis# -------------------------------------------------------------------------## get mortality datahfaction_D <- hfaction |>filter(status !=1) # remove hospitalization records (status=1)## RMST (Restricted Mean Survival Time) analysis for overall survival# Here we compare arm=1 (training) vs arm=0 (usual care)rmst_obj <-rmst2(time = hfaction_D$time,status = hfaction_D$status >0,arm = hfaction_D$trt_ab,tau =3.97)rmst_obj#> Prints estimates of restricted mean survival (in years), differences, ratios, and p-values# Between-group contrast # Est. lower .95 upper .95 p# RMST (arm=1)-(arm=0) 0.238 0.013 0.464 0.039# RMST (arm=1)/(arm=0) 1.074 1.003 1.150 0.040# RMTL (arm=1)/(arm=0) 0.680 0.468 0.988 0.043# Extract the difference in RMST (in years) and convert to monthsrmst <- rmst_obj$unadjusted.result[1, 1] *12rmst_p <- rmst_obj$unadjusted.result[1, 4]# -------------------------------------------------------------------------# TFE analysis# -------------------------------------------------------------------------## Check how many of the first events are death (status=1) or hospitalization (status=2)hfaction_TFE |>count(status)## RMST analysis for hospitalization-free survivalrmest_obj <-rmst2(time = hfaction_TFE$time,status = hfaction_TFE$status >0,arm = hfaction_TFE$trt_ab,tau =3.97)rmest_obj#> Similar output for restricted mean event-free survival# Between-group contrast # Est. lower .95 upper .95 p# RMST (arm=1)-(arm=0) 0.198 -0.064 0.459 0.139# RMST (arm=1)/(arm=0) 1.145 0.957 1.370 0.139# RMTL (arm=1)/(arm=0) 0.924 0.832 1.027 0.141# Extract the difference in event-free survival (years) and convert to monthsrmest <- rmest_obj$unadjusted.result[1, 1] *12rmest_p <- rmest_obj$unadjusted.result[1, 4]# -------------------------------------------------------------------------# Mortality vs TFE (Kaplan-Meier plots with annotations)# -------------------------------------------------------------------------library(ggsurvfit)library(patchwork)# Plot for overall survivalpD <-survfit2(Surv(time, status >0) ~ trt_ab, data = hfaction_D) |>ggsurvfit(linewidth =1) +scale_ggsurvfit() +annotate("text", x =4, y =1, hjust =1, vjust =1,label =str_c("4y-RMST = ", round(rmst, 2), " months"," (P = ", round(rmst_p, 3), ")") ) +scale_color_discrete(labels =c("Usual care", "Training")) +scale_x_continuous("Time (years)", limits =c(0, 4)) +labs(y ="Overall survival")# Plot for hospitalization-free survivalpTFE <-survfit2(Surv(time, status >0) ~ trt_ab, data = hfaction_TFE) |>ggsurvfit(linewidth =1) +annotate("text", x =4, y =1, hjust =1, vjust =1,label =str_c("4y-RMEST = ", round(rmest, 2), " months"," (P = ", round(rmest_p, 3), ")") ) +scale_ggsurvfit() +scale_color_discrete(labels =c("Usual care", "Training")) +scale_x_continuous("Time (years)", limits =c(0, 4)) +labs(y ="Hospitalization-free survival",caption ="RMEST: restricted mean event-free survival time" )# Combine side-by-sidepD + pTFE +plot_layout(guides ="collect") &theme(legend.position ="top",legend.text =element_text(size =12) ) # # ggsave("images/est_hfaction_unis.png", width = 8, height = 4.6)#> Uncomment above to save the combined figure# -------------------------------------------------------------------------# RMT-IF analysis# -------------------------------------------------------------------------# rmtfit() from the rmt package fits the restricted mean time lost / # time free approach for recurrent events & deathobj <-rmtfit(rec(patid, time, status) ~ trt_ab, data = hfaction)summary(obj, Kmax =1, tau =3.97)#> Summarizes the model up to a maximum follow-up of 3.97 years, #> focusing on the first event (Kmax=1) or aggregated events############################################################## Graphical analysis of the HF-ACTION trial to# evaluate the effect of exercise training.###########################################################par(mfrow =c(1, 2))# bouquet() displays the "bouquet plot" for each k up to Kmax=4bouquet( obj, Kmax =4, cex.group =1.0,xlab ="Restricted mean win/loss time (years)",ylab ="Follow-up time (years)",group.label =FALSE, ylim =c(0, 4.2))text(-0.8, 4.15, paste("Usual care"))text(0.8, 4.15, paste("Exercise training"))# plot() visualizes the RMT-IF difference as a function of timeplot( obj, conf =TRUE, lwd =2,xlab ="Follow-up time (years)",ylab ="RMT-IF of training (years)", main ="")par(mfrow =c(1, 1))#> Reset the graphic device to a single panel### LaTeX table ###pval_fmt3 =function(x) {if (x <0.001) {return("$<$0.001") } else {return(round(x, 3)) }}ltable =NULL# aggregate the results for k=1,...,Khosp_sum =summary(obj, Kmax =1, tau =3.97)$tab# aggregate the results for k=4,...,Kall_sum =summary(obj, Kmax =4, tau =3.97)$tabltable =c("&", "&", "&",round(12* hosp_sum[1, 1], 2), "&", round(12* hosp_sum[1, 2], 2),"&", pval_fmt3(hosp_sum[1, 4]), "\\")for (i in1:6) { tmp =c("&", i, "&&",round(12* all_sum[i, 1], 2), "&", round(12* all_sum[i, 2], 2),"&", pval_fmt3(all_sum[i, 4]), "\\" ) ltable =rbind(ltable, tmp)}ltable[5, 2] ="4+"ltable[6:7, 2] =""rownames(ltable) =c("Hopitalization", "", "", "", "", "Death", "Overall")noquote(ltable)#> Produces a LaTeX-formatted table summarizing the RMT-IF results######################################################################## WA analysis ######################################################################### install.packages("WA") # if neededlibrary(WA)# load the hf-action study data (already loaded, but repeated here)head(hfaction)# descriptive analysis## death & hosp rates by treatment armhfaction |>group_by(trt_ab, patid) |>summarize(ND =sum(status ==2),NH =sum(status ==1) ) |>summarize(death_rate =mean(ND),avgNH =mean(NH),sdNH =sd(NH) )# weighted while-alive event rate analysis # with death weighted as 2 vs 1 for hospitalizationobj <-LRfit(id = hfaction$patid,time = hfaction$time,status = hfaction$status,trt = hfaction$trt_ab,Dweight =2)## print some descriptive information obj## summarize the inference results at tau=4 yearssummary(obj, tau =3.97, joint.test =TRUE)plot(obj)#> Plots the Weighted-Alive event rates over time for each arm# -------------------------------------------------------------------------# unadjusted cumulative mean (Ch 1)# -------------------------------------------------------------------------## fit proportional means model with death = 2 x hosplibrary(Wcompo)## change status codingstatus <- hfaction$statusstatus[status !=0] <-3- status[status !=0]#> This transforms status=1->2, status=2->1, effectively swapping themobj_ML <-CompoML( hfaction$patid, hfaction$time, status, hfaction$trt_ab,w =c(2, 1))## summary resultst <- obj_ML$t # time pointsmu0 <- obj_ML$y # baseline cumulative mean eventsmu1 <- mu0 *exp(as.numeric(obj_ML$beta))#> The exponentiated beta shifts the baseline function for the second arm## plot survival-adjusted cumulative event ## vs unadjusted under PM modelplot( obj,ylab ="Cumulative loss",xlab ="Time (years)")lines( t[t <=3.97], mu0[t <=3.97],lty =3, col ="red", lwd =2)lines( t[t <=3.97], mu1[t <=3.97],lty =3, col ="blue", lwd =2)legend("bottomright",col =c("red", "red", "blue", "blue"),c("Usual care (WA)", "Usual care (unadj)","Training (WA)", "Training (unadj)"),lty =c(1, 3, 1, 3),lwd =2)#> Compares Weighted-Alive-based event accumulation with the "unadjusted" PM model
The win ratio estimand depends on the censoring distribution, as each pairwise comparison is limited to the observed follow-up \([0, C_i^\t \wedge C_j^\c]\). This leads to a censoring-weighted average of time-specific win probabilities, which is trial-dependent and lacks generalizability. One solution is to pre-specify a restriction time \(\tau\) so all comparisons are evaluated over the same window.
For a two-tiered composite of death and a nonfatal event, the restricted win/loss probability is defined by \[\begin{align}
w_{a, 1-a}(\tau) &= \pr\{D^\b < \min(D^\a, \tau)\} \\
&\quad + \pr\{\min(D^\t, D^\c) > \tau, T^\b < \min(T^\a, \tau)\},
\end{align}\] with the restricted win ratio defined as \(\text{WR}(\tau) = w_{1,0}(\tau) / w_{0,1}(\tau)\). Estimation methods include inverse probability of censoring weighting (IPCW) and multiple imputation (MI), as implemented in the WINS package.
3.2 Restricted Mean Time in Favor (RMT-IF)
The restricted mean time in favor (RMT-IF) redefines the win-loss comparison in terms of time spent in more favorable states. For a multistate outcome \(Y(t)\), RMT-IF measures the net average time one subject remains in a better state than another: \[
\mu(\tau) = E\left[ \int_0^\tau I\{Y^{(1)}(t) < Y^{(0)}(t)\} \, \dd t \right]
- E\left[ \int_0^\tau I\{Y^{(0)}(t) < Y^{(1)}(t)\} \, \dd t \right].
\]
The process \(Y(t)\) is assumed progressive and represents a hierarchy of worsening states, terminating in death. The estimand \(\mu(\tau)\) admits a decomposition into component-specific effects: \[
\mu(\tau) = \sum_{k=1}^K \mu_k(\tau) + \mu_\infty(\tau),
\] where each \(\mu_k(\tau)\) captures net time gained in a specific nonfatal state and \(\mu_\infty(\tau)\) captures the difference in restricted mean survival time. Estimation proceeds by plugging in Kaplan–Meier estimators for transition times. The rmt::rmtfit() function provides estimation, inference, and graphics.
3.3 While-Alive Weighted Loss
Standard analyses of weighted total events may misrepresent treatment effects if survival time differs between groups. The while-alive loss framework defines a normalized event rate over time alive: \[
\ell^{(a)}(\tau) = \frac{E\left\{N_{\rm R}^{*(a)}(\tau)\right\}}{E\left(D^{(a)} \wedge \tau\right)}.
\]
More generally, the loss can be specified by a function \(\mathcal{L}(\mathcal{H}^*(t))\) that accumulates only during time alive. This leads to a class of estimands: \[
\ell^{(a)}(\tau) = \frac{E\left\{\mathcal L\left(\mathcal H^{*(a)}\right)(\tau)\right\}}{E\left(D^{(a)} \wedge \tau\right)},
\] which includes total events, mortality, and flexible weighted combinations. Estimation uses survival-weighted Aalen-type integrals and is implemented in the WA::LRfit() function.
3.4 HF-ACTION Illustrations
All three approaches were applied to the HF-ACTION dataset:
Restricted WR showed better interpretability than the original WR by standardizing follow-up to 4 years.
RMT-IF detected a significant average gain of 5.1 months in the training group, driven by both survival and reduced hospitalizations.
While-alive loss adjusted for differential follow-up and showed a 20% reduction in the average event rate, though the effect was attenuated compared to RMT-IF.
Each method offers a distinct perspective:
RMT-IF summarizes net benefit in time.
WA analyses account for differential exposure.
Restricted WR retains the pairwise framework but standardizes follow-up.
3.5 Example R code
# RMT-IFlibrary(rmt)obj <-rmtfit(id, time, status, trt, type ="recurrent")summary(obj, tau =4)# While-alive loss (death = 2 × nonfatal)library(WA)obj <-LRfit(id, time, status, trt, Dweight =2)summary(obj, tau =4)
3.6 Conclusion
Restricting follow-up time resolves many of the interpretability issues in standard composite endpoint analysis. The restricted win ratio, restricted mean time in favor, and while-alive loss rate each target clinically interpretable, time-standardized estimands. These tools are particularly helpful when survival times vary across groups, and should be considered complementary strategies in the analysis of composite outcomes.