Computing the standardized estimate (G-estimation) based on the Cox or Fine-Gray model : Ŝ(t, A = a) = n−1∑iS(t, A = a, Xi) and this estimator has influence function S(t, A = a, Xi) − S(t, A = a) + E(DA0(t), βS(t, A = a, Xi))ϵi(t) where ϵi(t) is the iid decomposition of (Â(t) − A(t), β̂ − β).
These estimates have a causal interpration under the assumption of no-unmeasured confounders, and even without the causal assumptions this standardization can still be a useful summary measure.
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
#> ps0 -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
#> [ps0] 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
#> ps0 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
#> [ps0] 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.6539 0.02709 0.6008 0.7070 9.218e-129
#> risk1 0.5640 0.05971 0.4470 0.6811 3.531e-21
#>
#> Average Treatment effects (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> ps0 -0.08992 0.0629 -0.2132 0.03337 0.1529
#>
#> Average Treatment effect ratio (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> [ps0] 0.8624974 0.09446477 0.6773499 1.047645 0.1455042We compare with the similar estimates using the Doubly Robust estimating equations using binregATE. The standardization from the G-computation can also be computed using a specialized function that takes less memory and is quicker (for large data).
## survival situation
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.809
## relative risk effect
estimate(coef=sr1$riskDR,vcov=sr1$var.riskDR,f=function(p) p[2]/p[1],null=1)
#> Estimate Std.Err 2.5% 97.5% P-value
#> [treat1] 1.035 0.1453 0.7503 1.32 0.8096
#>
#> Null Hypothesis:
#> [treat1] = 1
## competing risks
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.0049and using the specialized function
br1 <- binreg(Event(time,cause)~tcell+platelet+age,bmt,cause=1,time=40)
Gbr1 <- binregG(br1,data=bmt)
summary(Gbr1)
#> risk:
#> Estimate Std.Err 2.5% 97.5% P-value
#> risk0 0.4178 0.02703 0.3648 0.4708 6.752e-54
#> risk1 0.2664 0.06204 0.1448 0.3880 1.756e-05
#>
#> Average Treatment effects (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 -0.1514 0.06776 -0.2842 -0.01859 0.02546
#>
#> Average Treatment effect ratio (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> [p1] 0.6376167 0.1542628 0.3352673 0.9399661 0.01881733
## contrasting average age to +2-sd age, Avalues
Gbr2 <- binregG(br1,data=bmt,varname="age",Avalues=c(0,2))
summary(Gbr2)
#> risk:
#> Estimate Std.Err 2.5% 97.5% P-value
#> risk0 0.3935 0.02529 0.3439 0.4431 1.432e-54
#> risk2 0.5929 0.05580 0.4836 0.7023 2.261e-26
#>
#> Average Treatment effects (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 0.1994 0.05019 0.1011 0.2978 7.069e-05
#>
#> Average Treatment effect ratio (G-estimator) :
#> Estimate Std.Err 2.5% 97.5% P-value
#> [p1] 1.506855 0.132196 1.247756 1.765954 0.000126016sessionInfo()
#> 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