Computing the G-estimation based on the Cox or Fine-Gray model : \[ \hat S(t,A=a) = n^{-1} \sum_i S(t,A=a,X_i) \] and this estimator has influence function \[ S(t,A=a,X_i) - S(t,A=a) + E( D_{A_0(t), \beta} S(t,A=a,X_i) ) \epsilon_i(t) \] where \(\epsilon_i(t)\) is the iid decomposition of \((\hat A(t) - A(t), \hat \beta- \beta)\).
set.seed(100)
data(bmt); bmt$time <- bmt$time+runif(nrow(bmt))*0.001
dfactor(bmt) <- tcell~tcell
bmt$event <- (bmt$cause!=0)*1
fg1 <- cifreg(Event(time,cause)~tcell+platelet+age,bmt,cause=1,
cox.prep=TRUE,propodds=NULL)
summary(survivalG(fg1,bmt,50))
#> risk:
#> Estimate Std.Err 2.5% 97.5% P-value
#> risk0 0.4331 0.02749 0.3793 0.4870 6.321e-56
#> risk1 0.2727 0.05863 0.1577 0.3876 3.313e-06
#>
#> Average Treatment effects (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 -0.1605 0.06353 -0.285 -0.03597 0.01153
#>
#> Average Treatment effect ratio (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> [p1] 0.6295004 0.139248 0.3565794 0.9024214 0.00779742
fg2 <- cifreg(Event(time,cause)~tcell+platelet+age,bmt,cause=2,
cox.prep=TRUE,propodds=NULL)
summary(survivalG(fg2,bmt,50))
#> risk:
#> Estimate Std.Err 2.5% 97.5% P-value
#> risk0 0.2127 0.02314 0.1674 0.2581 3.757e-20
#> risk1 0.3336 0.06799 0.2003 0.4668 9.281e-07
#>
#> Average Treatment effects (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 0.1208 0.07189 -0.02009 0.2617 0.09285
#>
#> Average Treatment effect ratio (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> [p1] 1.567915 0.3627528 0.8569321 2.278897 0.1174496
ss <- phreg(Surv(time,event)~tcell+platelet+age,bmt)
summary(survivalG(ss,bmt,50))
#> risk:
#> Estimate Std.Err 2.5% 97.5% P-value
#> risk0 0.3461 0.02709 0.2930 0.3992 2.252e-37
#> risk1 0.4360 0.05971 0.3189 0.5530 2.850e-13
#>
#> Average Treatment effects (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 0.08992 0.0629 -0.03337 0.2132 0.1529
#>
#> Average Treatment effect ratio (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> [p1] 1.259836 0.1894627 0.8884963 1.631176 0.1702385Comparing with binomial-regression ATE
br1 <- binregATE(Event(time,cause)~tcell+platelet+age,bmt,cause=1,
time=40,treat.model=tcell~platelet+age)
summary(br1)
#>
#> n events
#> 408 157
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.188365 0.130922 -0.444967 0.068237 0.1502
#> tcell1 -0.715361 0.352473 -1.406195 -0.024527 0.0424
#> platelet -0.537310 0.244804 -1.017117 -0.057502 0.0282
#> age 0.417814 0.107282 0.207545 0.628084 0.0001
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.82831 0.64085 1.0706
#> tcell1 0.48902 0.24507 0.9758
#> platelet 0.58432 0.36164 0.9441
#> age 1.51864 1.23065 1.8740
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 0.417795 0.027029 0.364819 0.470771 0.0000
#> treat1 0.266393 0.062041 0.144795 0.387991 0.0000
#> treat:1-0 -0.151402 0.067763 -0.284214 -0.018589 0.0255
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 0.417337 0.027120 0.364184 0.470491 0.0000
#> treat1 0.231224 0.060718 0.112218 0.350229 0.0001
#> treat:1-0 -0.186114 0.066117 -0.315700 -0.056527 0.0049
sr1 <- binregATE(Event(time,event)~tcell+platelet+age,bmt,cause=1,
time=40, treat.model=tcell~platelet+age)
summary(sr1)
#>
#> n events
#> 408 241
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.679693 0.138551 0.408138 0.951248 0.0000
#> tcell1 -0.032018 0.353415 -0.724698 0.660662 0.9278
#> platelet -0.504940 0.248245 -0.991492 -0.018387 0.0419
#> age 0.315033 0.117786 0.084178 0.545889 0.0075
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 1.97327 1.50401 2.5889
#> tcell1 0.96849 0.48447 1.9361
#> platelet 0.60354 0.37102 0.9818
#> age 1.37030 1.08782 1.7261
#>
#> Average Treatment effects (G-formula) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 0.6233534 0.0274214 0.5696085 0.6770983 0.000
#> treat1 0.6161006 0.0748225 0.4694512 0.7627499 0.000
#> treat:1-0 -0.0072528 0.0802736 -0.1645862 0.1500805 0.928
#>
#> Average Treatment effects (double robust) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 0.623341 0.027505 0.569433 0.677249 0.000
#> treat1 0.645159 0.085872 0.476853 0.813465 0.000
#> treat:1-0 0.021818 0.090254 -0.155076 0.198711 0.809sessionInfo()
#> 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