Fit
in the context of mediation analysis using mediation weights as in the medFlex package.
The mediator can
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.
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 <- 1Compute mediation weights based on fitted model
weightmodel <- fit <- glm(gp.f~dnr.f+preauto+ttt24,data=dat,family=binomial)
wdata <- medweight(fit,data=dat)
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.9206
aaMss22 <- binreg(Event(time,status)~dnr+preauto+ttt24+cluster(id),data=dat,time=50,cause=2)
summary(aaMss22)
#>
#> n events
#> 200 97
#>
#> 200 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.445508 0.257170 -0.949552 0.058536 0.0832
#> dnr 0.678218 0.381764 -0.070027 1.426462 0.0756
#> preauto 0.409926 0.359574 -0.294826 1.114677 0.2543
#> ttt24 0.479037 0.375065 -0.256077 1.214150 0.2015
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.64050 0.38691 1.0603
#> dnr 1.97036 0.93237 4.1639
#> preauto 1.50671 0.74466 3.0486
#> ttt24 1.61452 0.77408 3.3674### 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)}
### 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
#>
#> 200 clusters
#> 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.234269 -1.133592 -0.215275 0.0040
#> dnr.f01 0.221834 0.297615 -0.361481 0.805149 0.4560
#> dnr.f11 0.262722 0.154635 -0.040358 0.565801 0.0893
#> preauto 0.578077 0.312548 -0.034507 1.190660 0.0644
#> ttt24 0.214442 0.328885 -0.430161 0.859045 0.5144
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.50944 0.32187 0.8063
#> dnr.f01 1.24836 0.69664 2.2370
#> dnr.f11 1.30046 0.96045 1.7609
#> preauto 1.78261 0.96608 3.2893
#> ttt24 1.23917 0.65040 2.3609
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.2.1 (2022-06-23)
#> Platform: x86_64-apple-darwin21.5.0 (64-bit)
#> Running under: macOS Monterey 12.5.1
#>
#> Matrix products: default
#> BLAS: /usr/local/Cellar/openblas/0.3.21/lib/libopenblasp-r0.3.21.dylib
#> LAPACK: /usr/local/Cellar/r/4.2.1_2/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] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mets_1.3.0 lava_1.7.0 timereg_2.0.2 survival_3.3-1
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.9 highr_0.9 bslib_0.3.1
#> [4] compiler_4.2.1 jquerylib_0.1.4 tools_4.2.1
#> [7] digest_0.6.29 jsonlite_1.8.0 evaluate_0.15
#> [10] lattice_0.20-45 ucminf_1.1-4 rlang_1.0.2
#> [13] Matrix_1.4-1 cli_3.3.0 yaml_2.3.5
#> [16] parallel_4.2.1 mvtnorm_1.1-3 xfun_0.30
#> [19] fastmap_1.1.0 stringr_1.4.0 knitr_1.37
#> [22] sass_0.4.0 globals_0.16.1 grid_4.2.1
#> [25] listenv_0.8.0 R6_2.5.1 future.apply_1.9.0
#> [28] parallelly_1.32.1 rmarkdown_2.13 magrittr_2.0.2
#> [31] codetools_0.2-18 htmltools_0.5.2 splines_4.2.1
#> [34] future_1.27.0 numDeriv_2016.8-1.1 stringi_1.7.6