Fit
in the context of mediation analysis using mediation weights as in the medFlex package. We thus fit natural effects models, that for example on the binary scale might state that \[\begin{align*} \mbox{logit}(P(Y(x,M(x^*))=1| Z) = \beta_0+ \beta_1 x + \beta_2 x^* + \beta_3^T Z, \end{align*}\] in this case the the Natural Direct Effect (NDE) for fixed covariates \(Z\) is \[\begin{align*} \mbox{OR}_{1,0|Z}^{\mbox{NDE}} = \frac{\mbox{odds}(Y(1,M(x))|Z)}{\mbox{odds}(Y(0,M(x))|Z)} = \exp(\beta_1), \end{align*}\] and the Natural Inderect Effect (NIE) for fixed covariates \(Z\) is \[\begin{align*} \mbox{OR}_{1,0|Z}^{\mbox{NIE}} = \frac{\mbox{odds}(Y(x,M(1))|Z)}{\mbox{odds}(Y(x,M(0))|Z)} = \exp(\beta_2). \end{align*}\] See the medFlex package for additional discussion of the parametrization.
The mediator can be
Both mediator and exposure must be coded as factors.
In the below example these are
and the outcome model is concerned with the risk/hazard of cause=2.
The key is that the standard errors are computed using the i.i.d influence functions and a Taylor expansion to deal with the uncertainty from the mediation weights.
First we simulate some data that mimics that of Kumar et al 2012. This is data from multiple myeloma patients treated with allogeneic stem cell transplantation from the Center for International Blood and Marrow Transplant Research (CIBMTR) Kumar et al (2012), “Trends in allogeneic stem cell transplantation for multiple myeloma: a CIBMTR analysis”. The data used in this paper consist of patients transplanted from 1995 to 2005, and we compared the outcomes between transplant periods: 2001-2005 (N=488) versus 1995-2000 (N=375). The two competing events were relapse (cause 2) and treatment-related mortality (TRM, cause 1)) defined as death without relapse. considered the following risk covariates: transplant time period (gp (main interest of the study): 1 for transplanted in 2001-2005 versus 0 for transplanted in 1995-2000), donor type (dnr: 1 for Unrelated or other related donor (N=280) versus 0 for HLA-identical sibling (N=584)), prior autologous transplant (preauto: 1 for Auto+Allo transplant (N=399) versus 0 for allogeneic transplant alone (N=465)) and time to transplant (ttt24: 1 for more than 24 months (N=289) versus 0 for less than or equal to 24 months (N=575))).
The interest is then on the effect of the period (gp) and the possible mediation via the amount of unrealted or related donors (dnr). A somewhat artificial example ! All adjusted for other important counfounders.
library(mets)
runb <- 0
options(warn=-1)
set.seed(1000) # to control output in simulatins for p-values below.
n <- 200; k.boot <- 10;
dat <- kumarsimRCT(n,rho1=0.5,rho2=0.5,rct=2,censpar=c(0,0,0,0),
beta = c(-0.67, 0.59, 0.55, 0.25, 0.98, 0.18, 0.45, 0.31),
treatmodel = c(-0.18, 0.56, 0.56, 0.54),restrict=1)
dfactor(dat) <- dnr.f~dnr
dfactor(dat) <- gp.f~gp
drename(dat) <- ttt24~"ttt24*"
dat$id <- 1:n
dat$ftime <- 1Then compute the mediation weights based on a mediation model
A simple multvariate regression of the probaibility of relapse at 50 months with both exposure and mediator (given the other covariates)
aaMss2 <- binreg(Event(time,status)~gp+dnr+preauto+ttt24+cluster(id),data=dat,time=50,cause=2)
summary(aaMss2)
#>
#> n events
#> 200 97
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -1.09545 0.33261 -1.74736 -0.44355 0.0010
#> gp 1.26427 0.36725 0.54447 1.98407 0.0006
#> dnr 0.45202 0.39524 -0.32263 1.22667 0.2528
#> preauto 0.35124 0.38217 -0.39780 1.10028 0.3581
#> ttt24 0.55819 0.41228 -0.24987 1.36624 0.1758
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.33439 0.17423 0.6418
#> gp 3.54050 1.72370 7.2722
#> dnr 1.57148 0.72424 3.4099
#> preauto 1.42083 0.67180 3.0050
#> ttt24 1.74750 0.77890 3.9206We first look at the probability of relapse at 50 months
### binomial regression ###########################################################
aaMss <- binreg(Event(time,status)~dnr.f0+dnr.f1+preauto+ttt24+cluster(id),data=wdata,
time=50,weights=wdata$weights,cause=2)
summary(aaMss)
#>
#> n events
#> 400 194
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.520113 0.259635 -1.028987 -0.011238 0.0452
#> dnr.f01 0.339934 0.376813 -0.398605 1.078473 0.3670
#> dnr.f11 0.274655 0.071476 0.134564 0.414747 0.0001
#> preauto 0.552689 0.365370 -0.163423 1.268800 0.1304
#> ttt24 0.300601 0.380049 -0.444282 1.045483 0.4290
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.59445 0.35737 0.9888
#> dnr.f01 1.40486 0.67126 2.9402
#> dnr.f11 1.31608 1.14404 1.5140
#> preauto 1.73792 0.84923 3.5566
#> ttt24 1.35067 0.64128 2.8448
ll <- mediatorSurv(aaMss,fit,data=dat,wdata=wdata)
summary(ll)
#>
#> n events
#> 400 194
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.520113 0.259239 -1.028211 -0.012014 0.0448
#> dnr.f01 0.339934 0.344113 -0.334515 1.014383 0.3232
#> dnr.f11 0.274655 0.117283 0.044785 0.504525 0.0192
#> preauto 0.552689 0.361565 -0.155966 1.261343 0.1264
#> ttt24 0.300601 0.381561 -0.447245 1.048447 0.4308
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.59445 0.35765 0.9881
#> dnr.f01 1.40486 0.71569 2.7577
#> dnr.f11 1.31608 1.04580 1.6562
#> preauto 1.73792 0.85559 3.5302
#> ttt24 1.35067 0.63939 2.8532
if (runb>0) { bll <- BootmediatorSurv(aaMss,fit,data=dat,k.boot=k.boot); summary(bll)}So the NDE is \(1.40 (0.72,2.76)\) and the NIE is \(1.32 (1.05,1.66)\).
We here also illustrate how to use the other models mentioned above.
### lin-ying model ################################################################
aaMss <- aalenMets(Surv(time/100,status==2)~dnr.f0+dnr.f1+preauto+ttt24+cluster(id),data=wdata,
weights=wdata$weights)
ll <- mediatorSurv(aaMss,fit,data=dat,wdata=wdata)
summary(ll)
#>
#> n events
#> 400 196
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> dnr.f01 1.169592 0.739323 -0.279454 2.618637 0.1137
#> dnr.f11 0.206757 0.131289 -0.050565 0.464078 0.1153
#> preauto 0.617537 0.504302 -0.370877 1.605950 0.2207
#> ttt24 0.457736 0.517822 -0.557175 1.472648 0.3767
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> dnr.f01 3.22068 0.75620 13.7170
#> dnr.f11 1.22968 0.95069 1.5905
#> preauto 1.85435 0.69013 4.9826
#> ttt24 1.58049 0.57282 4.3608
if (runb>0) { bll <- BootmediatorSurv(aaMss,fit,data=dat,k.boot=k.boot); summary(bll)}
### cox model ###############################################################################
aaMss <- phreg(Surv(time,status==2)~dnr.f0+dnr.f1+preauto+ttt24+cluster(id),data=wdata,
weights=wdata$weights)
summary(aaMss)
#>
#> n events
#> 400 196
#> coeffients:
#> Estimate S.E. dU^-1/2 P-value
#> dnr.f01 0.414565 0.213724 0.157231 0.0524
#> dnr.f11 0.100656 0.039308 0.144971 0.0104
#> preauto 0.284460 0.232166 0.162375 0.2205
#> ttt24 0.185561 0.226044 0.160886 0.4117
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> dnr.f01 1.51371 0.99568 2.3013
#> dnr.f11 1.10590 1.02389 1.1945
#> preauto 1.32904 0.84318 2.0949
#> ttt24 1.20389 0.77300 1.8750
ll <- mediatorSurv(aaMss,fit,data=dat,wdata=wdata)
summary(ll)
#>
#> n events
#> 400 196
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> dnr.f01 0.41456472 0.20869639 0.00552731 0.82360212 0.0470
#> dnr.f11 0.10065575 0.05121458 0.00027702 0.20103448 0.0494
#> preauto 0.28445952 0.23037280 -0.16706288 0.73598192 0.2169
#> ttt24 0.18556110 0.22549763 -0.25640614 0.62752835 0.4106
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> dnr.f01 1.51371 1.00554 2.2787
#> dnr.f11 1.10590 1.00028 1.2227
#> preauto 1.32904 0.84615 2.0875
#> ttt24 1.20389 0.77383 1.8730
if (runb>0) { bll <- BootmediatorSurv(aaMss,fit,data=dat,k.boot=k.boot); summary(bll)}
### Fine-Gray #############################################################3
aaMss <- cifreg(Event(time,status)~dnr.f0+dnr.f1+preauto+ttt24+cluster(id),data=wdata,
weights=wdata$weights,propodds=NULL,cause=2)
summary(aaMss)
#>
#> n events
#> 400 196
#>
#> 200 clusters
#> coeffients:
#> Estimate S.E. dU^-1/2 P-value
#> dnr.f01 0.18943 0.21986 0.15855 0.3889
#> dnr.f11 0.18730 0.04083 0.14503 0.0000
#> preauto 0.41452 0.22783 0.16098 0.0688
#> ttt24 0.17304 0.22892 0.16308 0.4497
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> dnr.f01 1.20856 0.78545 1.8596
#> dnr.f11 1.20599 1.11324 1.3065
#> preauto 1.51364 0.96849 2.3656
#> ttt24 1.18892 0.75910 1.8621
ll <- mediatorSurv(aaMss,fit,data=dat,wdata=wdata)
summary(ll)
#>
#> n events
#> 400 196
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> dnr.f01 0.189426 0.200088 -0.202740 0.581592 0.3438
#> dnr.f11 0.187298 0.075416 0.039486 0.335111 0.0130
#> preauto 0.414517 0.224290 -0.025083 0.854118 0.0646
#> ttt24 0.173042 0.229932 -0.277617 0.623701 0.4517
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> dnr.f01 1.20856 0.81649 1.7889
#> dnr.f11 1.20599 1.04028 1.3981
#> preauto 1.51364 0.97523 2.3493
#> ttt24 1.18892 0.75759 1.8658
if (runb>0) { bll <- BootmediatorSurv(aaMss,fit,data=dat,k.boot=k.boot); summary(bll)}
### logit model #############################################################3
aaMss <- cifreg(Event(time,status)~dnr.f0+dnr.f1+preauto+ttt24+cluster(id),data=wdata,
weights=wdata$weights,cause=2)
summary(aaMss)
#>
#> n events
#> 400 196
#>
#> 200 clusters
#> coeffients:
#> Estimate S.E. dU^-1/2 P-value
#> dnr.f01 0.357168 0.339848 0.158937 0.2933
#> dnr.f11 0.272392 0.064166 0.145076 0.0000
#> preauto 0.657010 0.326082 0.160361 0.0439
#> ttt24 0.191333 0.353606 0.167443 0.5884
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> dnr.f01 1.42928 0.73424 2.7822
#> dnr.f11 1.31310 1.15792 1.4891
#> preauto 1.92902 1.01806 3.6551
#> ttt24 1.21086 0.60549 2.4215
ll <- mediatorSurv(aaMss,fit,data=dat,wdata=wdata)
summary(ll)
#>
#> n events
#> 400 196
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> dnr.f01 0.357168 0.317645 -0.265405 0.979740 0.2608
#> dnr.f11 0.272392 0.111348 0.054154 0.490630 0.0144
#> preauto 0.657010 0.322612 0.024703 1.289317 0.0417
#> ttt24 0.191333 0.356870 -0.508119 0.890786 0.5919
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> dnr.f01 1.42928 0.76690 2.6638
#> dnr.f11 1.31310 1.05565 1.6333
#> preauto 1.92902 1.02501 3.6303
#> ttt24 1.21086 0.60163 2.4370
if (runb>0) { bll <- BootmediatorSurv(aaMss,fit,data=dat,k.boot=k.boot); summary(bll)}
### binomial outcome ############################
aaMss <- binreg(Event(ftime,status)~dnr.f0+dnr.f1+preauto+ttt24+cluster(id),data=wdata,
time=50,weights=wdata$weights,cens.weights=1,cause=2)
summary(aaMss)
#>
#> n events
#> 400 196
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.674433 0.235285 -1.135583 -0.213284 0.0042
#> dnr.f01 0.221834 0.318264 -0.401952 0.845620 0.4858
#> dnr.f11 0.262722 0.060281 0.144572 0.380871 0.0000
#> preauto 0.578077 0.319091 -0.047331 1.203484 0.0700
#> ttt24 0.214442 0.328183 -0.428784 0.857669 0.5135
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.50944 0.32123 0.8079
#> dnr.f01 1.24836 0.66901 2.3294
#> dnr.f11 1.30046 1.15555 1.4636
#> preauto 1.78261 0.95377 3.3317
#> ttt24 1.23917 0.65130 2.3577
ll <- mediatorSurv(aaMss,fit,data=dat,wdata=wdata)
summary(ll)
#>
#> n events
#> 400 196
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.674433 0.235022 -1.135069 -0.213798 0.0041
#> dnr.f01 0.221834 0.286717 -0.340122 0.783789 0.4391
#> dnr.f11 0.262722 0.107508 0.052011 0.473432 0.0145
#> preauto 0.578077 0.315260 -0.039822 1.195975 0.0667
#> ttt24 0.214442 0.329107 -0.430596 0.859480 0.5147
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.50944 0.32140 0.8075
#> dnr.f01 1.24836 0.71168 2.1898
#> dnr.f11 1.30046 1.05339 1.6055
#> preauto 1.78261 0.96096 3.3068
#> ttt24 1.23917 0.65012 2.3619
if (runb>0) { bll <- BootmediatorSurv(aaMss,fit,data=dat,k.boot=k.boot); summary(bll)}Also works with mediator with more than two levels
data(tTRACE)
dcut(tTRACE) <- ~.
weightmodel <- fit <- mlogit(wmicat.4 ~agecat.4+vf+chf,data=tTRACE,family=binomial)
wdata <- medweight(fit,data=tTRACE)
aaMss <- binreg(Event(time,status)~agecat.40+ agecat.41+ vf+chf+cluster(id),data=wdata,
time=7,weights=wdata$weights,cause=9)
summary(aaMss)
MultMed <- mediatorSurv(aaMss,fit,data=tTRACE,wdata=wdata)summary(MultMed)
#>
#> n events
#> 4000 2016
#>
#> 1000 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -1.837886 0.174671 -2.180236 -1.495537 0.0000
#> agecat.40(60.1,68.6] 0.911956 0.223683 0.473545 1.350366 0.0000
#> agecat.40(68.6,75.6] 1.367813 0.222395 0.931926 1.803699 0.0000
#> agecat.40(75.6,96.3] 2.284648 0.250033 1.794592 2.774704 0.0000
#> agecat.41(60.1,68.6] 0.121138 0.053348 0.016578 0.225698 0.0232
#> agecat.41(68.6,75.6] 0.119408 0.053222 0.015095 0.223722 0.0249
#> agecat.41(75.6,96.3] 0.095385 0.053904 -0.010265 0.201035 0.0768
#> vf 0.704175 0.295938 0.124147 1.284202 0.0173
#> chf 1.162855 0.154843 0.859368 1.466342 0.0000
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.15915 0.11301 0.2241
#> agecat.40(60.1,68.6] 2.48919 1.60568 3.8588
#> agecat.40(68.6,75.6] 3.92675 2.53940 6.0721
#> agecat.40(75.6,96.3] 9.82223 6.01702 16.0339
#> agecat.41(60.1,68.6] 1.12878 1.01672 1.2532
#> agecat.41(68.6,75.6] 1.12683 1.01521 1.2507
#> agecat.41(75.6,96.3] 1.10008 0.98979 1.2227
#> vf 2.02218 1.13218 3.6118
#> chf 3.19905 2.36167 4.3334sessionInfo()
#> R version 4.3.2 (2023-10-31)
#> Platform: aarch64-apple-darwin22.6.0 (64-bit)
#> Running under: macOS Sonoma 14.3.1
#>
#> Matrix products: default
#> BLAS: /Users/kkzh/.asdf/installs/R/4.3.2/lib/R/lib/libRblas.dylib
#> LAPACK: /Users/kkzh/.asdf/installs/R/4.3.2/lib/R/lib/libRlapack.dylib; LAPACK version 3.11.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] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mets_1.3.4 timereg_2.0.5 survival_3.5-7
#>
#> loaded via a namespace (and not attached):
#> [1] cli_3.6.2 knitr_1.45 rlang_1.1.3
#> [4] xfun_0.41 highr_0.10 jsonlite_1.8.8
#> [7] listenv_0.9.1 future.apply_1.11.1 lava_1.7.4
#> [10] htmltools_0.5.6.1 sass_0.4.7 rmarkdown_2.25
#> [13] grid_4.3.2 evaluate_0.23 jquerylib_0.1.4
#> [16] fastmap_1.1.1 mvtnorm_1.2-4 yaml_2.3.7
#> [19] numDeriv_2016.8-1.1 compiler_4.3.2 codetools_0.2-19
#> [22] ucminf_1.2.0 Rcpp_1.0.12 future_1.33.1
#> [25] lattice_0.22-5 digest_0.6.34 R6_2.5.1
#> [28] parallelly_1.37.0 parallel_4.3.2 splines_4.3.2
#> [31] bslib_0.5.1 Matrix_1.6-5 tools_4.3.2
#> [34] globals_0.16.2 cachem_1.0.8