Average treatment effect (ATE) for Restricted mean survival and years lost of Competing risks
RMST
Regression for rmst outcome E(T ∧ t|X) = exp(XTβ) based on IPCW adjustment for censoring, thus solving the estimating equation 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) ∫0tS(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(ϵ = 1)(t − T ∧ t)|X) = ∫0tF1(s)ds.
For fully saturated model with full censoring model these estimators are
equivalent as illustrated below.
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.014348 0.063466 2.889956 3.138740 0.0000
#> tcell 0.106223 0.142443 -0.172960 0.385406 0.4558
#> platelet 0.246994 0.098833 0.053285 0.440704 0.0125
#> age -0.185946 0.043533 -0.271270 -0.100622 0.0000
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 20.37580 17.99252 23.0748
#> tcell 1.11207 0.84117 1.4702
#> platelet 1.28017 1.05473 1.5538
#> age 0.83032 0.76241 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.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
#> tcell=0, platelet=0 0 30 13.60294 0.8315422 16.39706
#> tcell=0, platelet=1 1 30 18.90129 1.2693293 11.09871
#> tcell=1, platelet=0 2 30 16.19119 2.4006192 13.80881
#> tcell=1, platelet=1 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
#> tcell=0, platelet=0 0 30 12.105127 4.291933 0.8508102 0.6161447
#> tcell=0, platelet=1 1 30 6.884185 4.214527 1.1741034 0.9056995
#> tcell=1, platelet=0 2 30 7.260772 6.548042 2.3532912 1.9703328
#> tcell=1, platelet=1 3 30 5.780360 6.453473 2.0924927 2.0815028
#> total.years.lost
#> tcell=0, platelet=0 16.39706
#> tcell=0, platelet=1 11.09871
#> tcell=1, platelet=0 13.80881
#> tcell=1, platelet=1 12.23383
## 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)
Average treatment effect
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.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
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> 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.4 timereg_2.0.6 survival_3.7-0
#> [6] rmarkdown_2.29
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.9 utf8_1.2.4 future_1.34.0
#> [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 fansi_1.0.6 scales_1.3.0
#> [16] isoband_0.2.7 codetools_0.2-20 numDeriv_2016.8-1.1
#> [19] jquerylib_0.1.4 lava_1.8.0 cli_3.6.3
#> [22] rlang_1.1.4 parallelly_1.39.0 future.apply_1.11.3
#> [25] munsell_0.5.1 withr_3.0.2 cachem_1.1.0
#> [28] yaml_2.3.10 tools_4.4.2 parallel_4.4.2
#> [31] ucminf_1.2.2 colorspace_2.1-1 globals_0.16.3
#> [34] buildtools_1.0.0 vctrs_0.6.5 R6_2.5.1
#> [37] lifecycle_1.0.4 MASS_7.3-61 pkgconfig_2.0.3
#> [40] bslib_0.8.0 pillar_1.9.0 gtable_0.3.6
#> [43] glue_1.8.0 Rcpp_1.0.13-1 xfun_0.49
#> [46] tibble_3.2.1 sys_3.4.3 knitr_1.49
#> [49] farver_2.1.2 htmltools_0.5.8.1 labeling_0.4.3
#> [52] maketools_1.3.1 compiler_4.4.2