Regression for rmst outcome \(E(T \wedge t | X) = exp(X^T \beta)\) based on IPCW adjustment for censoring. For full censoring model this is equal to integral of Kaplan-Meier. Also possible to compute Kaplan-Meier or Cox based estimator or rmst \[ \int_0^t S(s|X) ds \].
For competing risks the years lost can be computed via cumulative incidence functions or as IPCW estimator since \[ E( I(\epsilon=1) ( t - T \wedge t ) | X) = \int_0^t F_1(s) ds. \] For full censoring model these estimators are also equivalent.
set.seed(100)
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.029039 0.068026 2.895710 3.162368 0.0000
#> tcell -0.084278 0.265484 -0.604616 0.436060 0.7509
#> platelet 0.264828 0.104368 0.060271 0.469386 0.0112
#> age -0.200527 0.082898 -0.363004 -0.038051 0.0156
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 20.67736 18.09635 23.6265
#> tcell 0.91918 0.54628 1.5466
#> platelet 1.30321 1.06212 1.5990
#> age 0.81830 0.69558 0.9627
### 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.999e-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.461e-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
#> 1 0 30 13.60294 0.8315422 16.39706
#> 2 1 30 18.90129 1.2693293 11.09871
#> 3 2 30 16.19119 2.4006192 13.80881
#> 4 3 30 17.76617 2.4421866 12.23383
## 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 total.years.lost
#> 1 0 30 12.105127 4.291933 0.8508102 0.6161447 16.39706
#> 2 1 30 6.884185 4.214527 1.1741034 0.9056995 11.09871
#> 3 2 30 7.260772 6.548042 2.3532912 1.9703328 13.80881
#> 4 3 30 5.780360 6.453473 2.0924927 2.0815028 12.23383Computes average treatment effect for restricted mean survival
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.8637403 0.0756653 2.7154390 3.0120416 0.0000
#> tcell1 0.0185112 0.1981777 -0.3699100 0.4069324 0.9256
#> platelet 0.2753483 0.1452274 -0.0092922 0.5599888 0.0580
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 17.52696 15.11124 20.3289
#> tcell1 1.01868 0.69080 1.5022
#> platelet 1.31699 0.99075 1.7507
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 19.26997 1.05138 17.20931 21.33064 0.0000
#> treat1 19.63001 3.43091 12.90554 26.35447 0.0000
#> treat:1-0 0.36003 3.87963 -7.24391 7.96397 0.9261
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 19.3253 1.0516 17.2642 21.3864 0.0000
#> treat1 21.5622 3.8017 14.1111 29.0133 0.0000
#> treat:1-0 2.2369 4.1991 -5.9932 10.4670 0.5942
out1 <- resmeanATE(Event(time,cause)~tcell+platelet,data=bmt,cause=1,outcome="rmst-cause",
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.807099 0.069703 2.670483 2.943715 0.0000
#> tcell1 -0.376884 0.247936 -0.862829 0.109061 0.1285
#> platelet -0.494384 0.165213 -0.818194 -0.170573 0.0028
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 16.56180 14.44695 18.9862
#> tcell1 0.68600 0.42197 1.1152
#> platelet 0.60995 0.44123 0.8432
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 14.53514 0.95673 12.65999 16.41029 0.0000
#> treat1 9.97105 2.37568 5.31479 14.62730 0.0000
#> treat:1-0 -4.56409 2.57161 -9.60437 0.47618 0.0759
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 14.5202792 0.9577035 12.6432148 16.3973436 0.0000
#> treat1 9.4567098 2.4051143 4.7427723 14.1706473 0.0001
#> treat:1-0 -5.0635694 2.5856565 -10.1313629 0.0042241 0.0502
out2 <- resmeanATE(Event(time,cause)~tcell+platelet,data=bmt,cause=2,outcome="rmst-cause",
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.8287103 0.1312246 1.5715149 2.0859057 0.0000
#> tcell1 0.4681132 0.2432252 -0.0085996 0.9448259 0.0543
#> platelet -0.0109542 0.2178062 -0.4378464 0.4159380 0.9599
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 6.22585 4.81394 8.0519
#> tcell1 1.59698 0.99144 2.5724
#> platelet 0.98911 0.64542 1.5158
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 6.20457 0.71391 4.80534 7.60381 0.0000
#> treat1 9.90857 2.10922 5.77458 14.04256 0.0000
#> treat:1-0 3.70399 2.23512 -0.67676 8.08474 0.0975
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 6.21823 0.71423 4.81836 7.61809 0.0000
#> treat1 10.38475 2.22282 6.02810 14.74140 0.0000
#> treat:1-0 4.16652 2.33621 -0.41237 8.74542 0.0745sessionInfo()
#> R version 4.2.2 (2022-10-31)
#> Platform: aarch64-apple-darwin22.1.0 (64-bit)
#> Running under: macOS Ventura 13.1
#>
#> Matrix products: default
#> BLAS: /opt/homebrew/Cellar/openblas/0.3.21/lib/libopenblasp-r0.3.21.dylib
#> LAPACK: /opt/homebrew/Cellar/r/4.2.2_1/lib/R/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] splines stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggplot2_3.4.0 cowplot_1.1.1 mets_1.3.2 timereg_2.0.5 survival_3.4-0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.0 xfun_0.36 bslib_0.4.2
#> [4] listenv_0.9.0 lattice_0.20-45 generics_0.1.3
#> [7] colorspace_2.0-3 vctrs_0.5.1 htmltools_0.5.4
#> [10] yaml_2.3.6 utf8_1.2.2 rlang_1.0.6
#> [13] isoband_0.2.7 jquerylib_0.1.4 pillar_1.8.1
#> [16] withr_2.5.0 glue_1.6.2 lifecycle_1.0.3
#> [19] lava_1.7.2 stringr_1.5.0 munsell_0.5.0
#> [22] gtable_0.3.1 future_1.30.0 mvtnorm_1.1-3
#> [25] codetools_0.2-18 evaluate_0.19 labeling_0.4.2
#> [28] knitr_1.41 fastmap_1.1.0 parallel_4.2.2
#> [31] fansi_1.0.3 highr_0.10 Rcpp_1.0.9
#> [34] scales_1.2.1 cachem_1.0.6 jsonlite_1.8.4
#> [37] farver_2.1.1 parallelly_1.34.0 digest_0.6.31
#> [40] stringi_1.7.8 dplyr_1.0.10 numDeriv_2016.8-1.1
#> [43] grid_4.2.2 cli_3.5.0 tools_4.2.2
#> [46] magrittr_2.0.3 sass_0.4.4 tibble_3.1.8
#> [49] ucminf_1.1-4.1 future.apply_1.10.0 pkgconfig_2.0.3
#> [52] MASS_7.3-58.1 Matrix_1.5-1 rmarkdown_2.19
#> [55] R6_2.5.1 globals_0.16.2 compiler_4.2.2