Average treatment effect (ATE) for Restricted mean survival and years lost of Competing risks

Klaus Holst & Thomas Scheike

2025-01-11

RMST

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)

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.0747

SessionInfo

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