Regression for rmst outcome \(E(T \wedge t | X) = exp(X^T \beta)\) based on IPCW adjustment for censoring, thus solving the estimating equation \[\begin{align*} & X^T [ (T \wedge t) \frac{I(C > T \wedge t)}{G_C(T \wedge t,X)} - exp(X^T \beta) ] = 0 . \end{align*}\] This is done with the resmeanIPCW function. For fully saturated model with full censoring model this is equal to the integrals of the Kaplan-Meier estimators as illustrated below.
We can also compute the integral of the Kaplan-Meier or Cox based survival estimator to get the RMST (with the resmean.phreg function) \[ \int_0^t S(s|X) ds \].
For competing risks the years lost can be computed via cumulative
incidence functions (cif.yearslost) or as IPCW estimator since
\[ E( I(\epsilon=1) ( t - T \wedge t ) | X) =
\int_0^t F_1(s) ds. \] For fully saturated model with full
censoring model these estimators are equivalent as illustrated
below.
set.seed(101)
data(bmt); bmt$time <- bmt$time+runif(nrow(bmt))*0.001
# E( min(T;t) | X ) = exp( a+b X) with IPCW estimation
out <- resmeanIPCW(Event(time,cause!=0)~tcell+platelet+age,bmt,
time=50,cens.model=~strata(platelet),model="exp")
summary(out)
#>
#> n events
#> 408 245
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 3.014321 0.065114 2.886701 3.141941 0.0000
#> tcell 0.106524 0.138268 -0.164475 0.377524 0.4410
#> platelet 0.247108 0.097337 0.056331 0.437885 0.0111
#> age -0.185939 0.043566 -0.271327 -0.100552 0.0000
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 20.37525 17.93404 23.1488
#> tcell 1.11240 0.84834 1.4587
#> platelet 1.28032 1.05795 1.5494
#> age 0.83032 0.76237 0.9043
### same as Kaplan-Meier for full censoring model
bmt$int <- with(bmt,strata(tcell,platelet))
out <- resmeanIPCW(Event(time,cause!=0)~-1+int,bmt,time=30,
cens.model=~strata(platelet,tcell),model="lin")
estimate(out)
#> Estimate Std.Err 2.5% 97.5% P-value
#> inttcell=0, platelet=0 13.60 0.8316 11.97 15.23 3.826e-60
#> inttcell=0, platelet=1 18.90 1.2696 16.41 21.39 3.997e-50
#> inttcell=1, platelet=0 16.19 2.4061 11.48 20.91 1.705e-11
#> inttcell=1, platelet=1 17.77 2.4536 12.96 22.58 4.463e-13
out1 <- phreg(Surv(time,cause!=0)~strata(tcell,platelet),data=bmt)
rm1 <- resmean.phreg(out1,times=30)
summary(rm1)
#> strata times rmean se.rmean years.lost
#> tcell=0, platelet=0 0 30 13.60295 0.8315418 16.39705
#> tcell=0, platelet=1 1 30 18.90127 1.2693263 11.09873
#> tcell=1, platelet=0 2 30 16.19121 2.4006185 13.80879
#> tcell=1, platelet=1 3 30 17.76610 2.4421987 12.23390
## competing risks years-lost for cause 1
out <- resmeanIPCW(Event(time,cause)~-1+int,bmt,time=30,cause=1,
cens.model=~strata(platelet,tcell),model="lin")
estimate(out)
#> Estimate Std.Err 2.5% 97.5% P-value
#> inttcell=0, platelet=0 12.105 0.8508 10.438 13.773 6.162e-46
#> inttcell=0, platelet=1 6.884 1.1741 4.583 9.185 4.536e-09
#> inttcell=1, platelet=0 7.261 2.3533 2.648 11.873 2.033e-03
#> inttcell=1, platelet=1 5.780 2.0925 1.679 9.882 5.737e-03
## same as integrated cumulative incidence
rmc1 <- cif.yearslost(Event(time,cause)~strata(tcell,platelet),data=bmt,times=30,cause=1)
summary(rmc1)
#> strata times intF11 intF12 se.intF11 se.intF12
#> tcell=0, platelet=0 0 30 12.105125 4.291930 0.8508102 0.6161439
#> tcell=0, platelet=1 1 30 6.884171 4.214556 1.1740988 0.9057028
#> tcell=1, platelet=0 2 30 7.260755 6.548034 2.3532867 1.9703340
#> tcell=1, platelet=1 3 30 5.780369 6.453535 2.0924946 2.0815225
#> total.years.lost
#> tcell=0, platelet=0 16.39705
#> tcell=0, platelet=1 11.09873
#> tcell=1, platelet=0 13.80879
#> tcell=1, platelet=1 12.23390
## plotting the years lost for different horizon's and the two causes
par(mfrow=c(1,3))
plot(rm1,years.lost=TRUE,se=1)
## cause refers to column of cumhaz for the different causes
plot(rmc1,cause=1,se=1)
plot(rmc1,cause=2,se=1)
Computes average treatment effect for restricted mean survival and years lost in competing risks situation
dfactor(bmt) <- tcell~tcell
bmt$event <- (bmt$cause!=0)*1
out <- resmeanATE(Event(time,event)~tcell+platelet,data=bmt,time=40,treat.model=tcell~platelet)
summary(out)
#>
#> n events
#> 408 241
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 2.852670 0.062493 2.730187 2.975153 0.0000
#> tcell1 0.021403 0.122907 -0.219489 0.262296 0.8618
#> platelet 0.303496 0.090758 0.125614 0.481379 0.0008
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 17.33400 15.33575 19.5926
#> tcell1 1.02163 0.80293 1.2999
#> platelet 1.35459 1.13384 1.6183
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 19.26228 0.95926 17.38217 21.14239 0.0000
#> treat1 19.67900 2.22777 15.31265 24.04536 0.0000
#> treat:1-0 0.41672 2.41067 -4.30810 5.14154 0.8628
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 19.3262 1.0516 17.2650 21.3873 0.0000
#> treat1 21.5621 3.8018 14.1106 29.0135 0.0000
#> treat:1-0 2.2359 4.1993 -5.9946 10.4664 0.5944
out1 <- resmeanATE(Event(time,cause)~tcell+platelet,data=bmt,cause=1,time=40,
treat.model=tcell~platelet)
summary(out1)
#>
#> n events
#> 408 157
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 2.806305 0.069618 2.669856 2.942754 0.0000
#> tcell1 -0.374071 0.247669 -0.859493 0.111351 0.1310
#> platelet -0.491745 0.164946 -0.815033 -0.168456 0.0029
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 16.54866 14.43789 18.9680
#> tcell1 0.68793 0.42338 1.1178
#> platelet 0.61156 0.44262 0.8450
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 14.53197 0.95707 12.65616 16.40778 0.0000
#> treat1 9.99695 2.37814 5.33588 14.65802 0.0000
#> treat:1-0 -4.53502 2.57516 -9.58224 0.51220 0.0782
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 14.520356 0.957707 12.643285 16.397428 0.0000
#> treat1 9.456750 2.405120 4.742802 14.170697 0.0001
#> treat:1-0 -5.063606 2.585664 -10.131415 0.004202 0.0502
out2 <- resmeanATE(Event(time,cause)~tcell+platelet,data=bmt,cause=2,time=40,
treat.model=tcell~platelet)
summary(out2)
#>
#> n events
#> 408 84
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 1.8266093 0.1312179 1.5694269 2.0837917 0.0000
#> tcell1 0.4751761 0.2403802 0.0040395 0.9463127 0.0481
#> platelet -0.0090727 0.2168458 -0.4340827 0.4159374 0.9666
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 6.21279 4.80389 8.0349
#> tcell1 1.60830 1.00405 2.5762
#> platelet 0.99097 0.64786 1.5158
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 6.19518 0.71372 4.79631 7.59405 0.0000
#> treat1 9.96369 2.09256 5.86236 14.06503 0.0000
#> treat:1-0 3.76851 2.21654 -0.57582 8.11285 0.0891
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 6.21830 0.71424 4.81842 7.61818 0.0000
#> treat1 10.38142 2.22242 6.02556 14.73728 0.0000
#> treat:1-0 4.16312 2.33581 -0.41500 8.74123 0.0747
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: aarch64-apple-darwin24.2.0
#> Running under: macOS Sequoia 15.2
#>
#> Matrix products: default
#> BLAS: /Users/klaus/.asdf/installs/R/4.4.2/lib/R/lib/libRblas.dylib
#> LAPACK: /Users/klaus/.asdf/installs/R/4.4.2/lib/R/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: Europe/Copenhagen
#> tzcode source: internal
#>
#> attached base packages:
#> [1] splines stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggplot2_3.5.1 cowplot_1.1.3 mets_1.3.5 timereg_2.0.6 survival_3.8-3
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.9 future_1.34.0 generics_0.1.3
#> [4] lattice_0.22-6 listenv_0.9.1 digest_0.6.37
#> [7] magrittr_2.0.3 evaluate_1.0.1 grid_4.4.2
#> [10] mvtnorm_1.3-2 fastmap_1.2.0 jsonlite_1.8.9
#> [13] Matrix_1.7-1 scales_1.3.0 isoband_0.2.7
#> [16] codetools_0.2-20 numDeriv_2016.8-1.1 jquerylib_0.1.4
#> [19] lava_1.8.1 cli_3.6.3 rlang_1.1.4
#> [22] parallelly_1.41.0 future.apply_1.11.3 munsell_0.5.1
#> [25] withr_3.0.2 cachem_1.1.0 yaml_2.3.10
#> [28] tools_4.4.2 parallel_4.4.2 ucminf_1.2.2
#> [31] dplyr_1.1.4 colorspace_2.1-1 globals_0.16.3
#> [34] vctrs_0.6.5 R6_2.5.1 lifecycle_1.0.4
#> [37] MASS_7.3-64 pkgconfig_2.0.3 bslib_0.8.0
#> [40] pillar_1.10.1 gtable_0.3.6 glue_1.8.0
#> [43] Rcpp_1.0.13-1 xfun_0.50 tibble_3.2.1
#> [46] tidyselect_1.2.1 knitr_1.49 farver_2.1.2
#> [49] htmltools_0.5.8.1 labeling_0.4.3 rmarkdown_2.29
#> [52] compiler_4.4.2