Average treatment effect (ATE) for Restricted mean survival and years lost of Competing risks
Klaus Holst & Thomas Scheike
2025-01-11
Source:vignettes/rmst-ate.Rmd
rmst-ate.RmdRMST
Regression for rmst outcome 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) .
For competing risks the years lost can be computed via cumulative
incidence functions (cif.yearslost) or as IPCW estimator since
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)
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.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.0747SessionInfo
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=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mets_1.3.5 timereg_2.0.6 survival_3.7-0
#>
#> loaded via a namespace (and not attached):
#> [1] cli_3.6.3 knitr_1.49 rlang_1.1.4
#> [4] xfun_0.50 textshaping_0.4.1 jsonlite_1.8.9
#> [7] listenv_0.9.1 future.apply_1.11.3 lava_1.8.0
#> [10] htmltools_0.5.8.1 ragg_1.3.3 sass_0.4.9
#> [13] rmarkdown_2.29 grid_4.4.2 evaluate_1.0.3
#> [16] jquerylib_0.1.4 fastmap_1.2.0 mvtnorm_1.3-3
#> [19] numDeriv_2016.8-1.1 yaml_2.3.10 lifecycle_1.0.4
#> [22] compiler_4.4.2 codetools_0.2-20 fs_1.6.5
#> [25] Rcpp_1.0.13-1 future_1.34.0 systemfonts_1.1.0
#> [28] lattice_0.22-6 digest_0.6.37 R6_2.5.1
#> [31] parallelly_1.41.0 parallel_4.4.2 splines_4.4.2
#> [34] bslib_0.8.0 Matrix_1.7-1 tools_4.4.2
#> [37] globals_0.16.3 pkgdown_2.1.1 cachem_1.1.0
#> [40] desc_1.4.3