Mediation Analysis for survival data

Klaus Holst & Thomas Scheike

2024-02-16

Overview

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.

Simulated Data

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 <- 1

Mediation Weights

Then compute the mediation weights based on a mediation model

weightmodel <- fit <- glm(gp.f~dnr.f+preauto+ttt24,data=dat,family=binomial)
wdata <- medweight(fit,data=dat)

Binomial Regression

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

Binomial regression IPCW Mediation Analysis

We 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)\).

Mediation Analysis

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)}

Multinomial regression

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

SessionInfo

sessionInfo()
#> 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