Implementation of extended formulas when there are effect measure modifiers

Yi Li, Kazuki Yoshida

2024-01-06

In this document, we demonstrate including effect measure modification (EMM) terms in the mediator or the outcome models. The dataset used in this document is still vv2015.

library(regmedint)
library(tidyverse)
## Prepare dataset
data(vv2015)

No EMM by covariates

In the first model fit, we do not include any EMM term.

regmedint_obj1 <- regmedint(data = vv2015,
                            ## Variables
                            yvar = "y",
                            avar = "x",
                            mvar = "m",
                            cvar = c("c"),
                            eventvar = "event",
                            ## Values at which effects are evaluated
                            a0 = 0,
                            a1 = 1,
                            m_cde = 1,
                            c_cond = 3,
                            ## Model types
                            mreg = "logistic",
                            yreg = "survAFT_weibull",
                            ## Additional specification
                            interaction = TRUE,
                            casecontrol = FALSE)
summary(regmedint_obj1)
## ### Mediator model
## 
## Call:
## glm(formula = m ~ x + c, family = binomial(link = "logit"), data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.3545     0.3252  -1.090    0.276
## x             0.3842     0.4165   0.922    0.356
## c             0.2694     0.2058   1.309    0.191
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 138.59  on 99  degrees of freedom
## Residual deviance: 136.08  on 97  degrees of freedom
## AIC: 142.08
## 
## Number of Fisher Scoring iterations: 4
## 
## ### Outcome model
## 
## Call:
## survival::survreg(formula = Surv(y, event) ~ x + m + x:m + c, 
##     data = data, dist = "weibull")
##               Value Std. Error     z       p
## (Intercept) -1.0424     0.1903 -5.48 4.3e-08
## x            0.4408     0.3008  1.47    0.14
## m            0.0905     0.2683  0.34    0.74
## c           -0.0669     0.0915 -0.73    0.46
## x:m          0.1003     0.4207  0.24    0.81
## Log(scale)  -0.0347     0.0810 -0.43    0.67
## 
## Scale= 0.966 
## 
## Weibull distribution
## Loglik(model)= -11.4   Loglik(intercept only)= -14.5
##  Chisq= 6.31 on 4 degrees of freedom, p= 0.18 
## Number of Newton-Raphson Iterations: 5 
## n= 100 
## 
## ### Mediation analysis 
##              est         se         Z          p       lower      upper
## cde  0.541070807 0.29422958 1.8389409 0.06592388 -0.03560858 1.11775019
## pnde 0.505391952 0.21797147 2.3186151 0.02041591  0.07817572 0.93260819
## tnie 0.015988820 0.03171597 0.5041252 0.61417338 -0.04617334 0.07815098
## tnde 0.513662425 0.22946248 2.2385465 0.02518544  0.06392423 0.96340062
## pnie 0.007718348 0.02398457 0.3218047 0.74760066 -0.03929055 0.05472725
## te   0.521380773 0.22427066 2.3247837 0.02008353  0.08181835 0.96094319
## pm   0.039039346 0.07444080 0.5244348 0.59997616 -0.10686194 0.18494063
## 
## Evaluated at:
## avar: x
##  a1 (intervened value of avar) = 1
##  a0 (reference value of avar)  = 0
## mvar: m
##  m_cde (intervend value of mvar for cde) = 1
## cvar: c
##  c_cond (covariate vector value) = 3
## 
## Note that effect estimates can vary over m_cde and c_cond values when interaction = TRUE.

EMM by covariates

There is \(A\times C\) term in mediator model

Now suppose the covariate \(C\) modifies the treatment effect on the mediator. We add emm_ac_mreg = c("c") in regmedint(). Although there is only one covariate in our dataset, emm_ac_mreg can take a vector of multiple covariates. Please note that the covariates in emm_ac_mreg should be a subset of the covariates specified in cvar, i.e. if a covariate is an effect measure modifier included in emm_ac_mreg, it must be included in cvar, otherwise an error message will be printed.

regmedint_obj2 <- regmedint(data = vv2015,
                            ## Variables
                            yvar = "y",
                            avar = "x",
                            mvar = "m",
                            cvar = c("c"),
                            emm_ac_mreg = c("c"),
                            emm_ac_yreg = NULL,
                            emm_mc_yreg = NULL,
                            eventvar = "event",
                            ## Values at which effects are evaluated
                            a0 = 0,
                            a1 = 1,
                            m_cde = 1,
                            c_cond = 3,
                            ## Model types
                            mreg = "logistic",
                            yreg = "survAFT_weibull",
                            ## Additional specification
                            interaction = TRUE,
                            casecontrol = FALSE)
summary(regmedint_obj2)
## ### Mediator model
## 
## Call:
## glm(formula = m ~ x + c + x:c, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.32727    0.34979  -0.936    0.349
## x            0.30431    0.56789   0.536    0.592
## c            0.24085    0.24688   0.976    0.329
## x:c          0.09216    0.44624   0.207    0.836
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 138.59  on 99  degrees of freedom
## Residual deviance: 136.04  on 96  degrees of freedom
## AIC: 144.04
## 
## Number of Fisher Scoring iterations: 4
## 
## ### Outcome model
## 
## Call:
## survival::survreg(formula = Surv(y, event) ~ x + m + x:m + c, 
##     data = data, dist = "weibull")
##               Value Std. Error     z       p
## (Intercept) -1.0424     0.1903 -5.48 4.3e-08
## x            0.4408     0.3008  1.47    0.14
## m            0.0905     0.2683  0.34    0.74
## c           -0.0669     0.0915 -0.73    0.46
## x:m          0.1003     0.4207  0.24    0.81
## Log(scale)  -0.0347     0.0810 -0.43    0.67
## 
## Scale= 0.966 
## 
## Weibull distribution
## Loglik(model)= -11.4   Loglik(intercept only)= -14.5
##  Chisq= 6.31 on 4 degrees of freedom, p= 0.18 
## Number of Newton-Raphson Iterations: 5 
## n= 100 
## 
## ### Mediation analysis 
##             est         se         Z          p       lower      upper
## cde  0.54107081 0.29422958 1.8389409 0.06592388 -0.03560858 1.11775019
## pnde 0.50404000 0.21666437 2.3263632 0.01999919  0.07938564 0.92869435
## tnie 0.02377050 0.05679639 0.4185213 0.67556602 -0.08754838 0.13508937
## tnde 0.51632801 0.23444392 2.2023519 0.02764046  0.05682637 0.97582965
## pnie 0.01148248 0.03882957 0.2957149 0.76744784 -0.06462208 0.08758704
## te   0.52781049 0.22811645 2.3137765 0.02067998  0.08071046 0.97491053
## pm   0.05727853 0.13042523 0.4391676 0.66054011 -0.19835021 0.31290728
## 
## Evaluated at:
## avar: x
##  a1 (intervened value of avar) = 1
##  a0 (reference value of avar)  = 0
## mvar: m
##  m_cde (intervend value of mvar for cde) = 1
## cvar: c
##  c_cond (covariate vector value) = 3
## 
## Note that effect estimates can vary over m_cde and c_cond values when interaction = TRUE.

There is \(A\times C\) term in both mediator and outcome models

Now suppose in addition to the EMM on mediator, the covariate \(C\) also modifies the treatment effect on the outcome We add emm_ac_yreg = c("c") in regmedint(). Please note that the covariates in emm_ac_yreg should be a subset of the covariates specified in cvar, i.e. if a covariate is an effect measure modifier included in emm_ac_yreg, it must be included in cvar, otherwise an error message will be printed.

regmedint_obj3 <- regmedint(data = vv2015,
                            ## Variables
                            yvar = "y",
                            avar = "x",
                            mvar = "m",
                            cvar = c("c"),
                            emm_ac_mreg = c("c"),
                            emm_ac_yreg = c("c"),
                            emm_mc_yreg = NULL,
                            eventvar = "event",
                            ## Values at which effects are evaluated
                            a0 = 0,
                            a1 = 1,
                            m_cde = 1,
                            c_cond = 3,
                            ## Model types
                            mreg = "logistic",
                            yreg = "survAFT_weibull",
                            ## Additional specification
                            interaction = TRUE,
                            casecontrol = FALSE)
summary(regmedint_obj3)
## ### Mediator model
## 
## Call:
## glm(formula = m ~ x + c + x:c, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.32727    0.34979  -0.936    0.349
## x            0.30431    0.56789   0.536    0.592
## c            0.24085    0.24688   0.976    0.329
## x:c          0.09216    0.44624   0.207    0.836
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 138.59  on 99  degrees of freedom
## Residual deviance: 136.04  on 96  degrees of freedom
## AIC: 144.04
## 
## Number of Fisher Scoring iterations: 4
## 
## ### Outcome model
## 
## Call:
## survival::survreg(formula = Surv(y, event) ~ x + m + x:m + c + 
##     x:c, data = data, dist = "weibull")
##                Value Std. Error     z       p
## (Intercept) -1.04148    0.19261 -5.41 6.4e-08
## x            0.43626    0.33092  1.32    0.19
## m            0.09138    0.26954  0.34    0.73
## c           -0.06844    0.10315 -0.66    0.51
## x:m          0.09681    0.43437  0.22    0.82
## x:c          0.00725    0.22300  0.03    0.97
## Log(scale)  -0.03473    0.08104 -0.43    0.67
## 
## Scale= 0.966 
## 
## Weibull distribution
## Loglik(model)= -11.4   Loglik(intercept only)= -14.5
##  Chisq= 6.31 on 5 degrees of freedom, p= 0.28 
## Number of Newton-Raphson Iterations: 5 
## n= 100 
## 
## ### Mediation analysis 
##             est         se         Z         p       lower      upper
## cde  0.55481498 0.51519657 1.0768996 0.2815251 -0.45495174 1.56458171
## pnde 0.51905802 0.51048271 1.0167984 0.3092493 -0.48146970 1.51958574
## tnie 0.02345019 0.05730239 0.4092359 0.6823666 -0.08886042 0.13576081
## tnde 0.53092081 0.50659327 1.0480218 0.2946285 -0.46198375 1.52382537
## pnie 0.01158740 0.03904582 0.2967641 0.7666466 -0.06494101 0.08811581
## te   0.54250821 0.50656000 1.0709654 0.2841850 -0.45033114 1.53534756
## pm   0.05535403 0.13968838 0.3962680 0.6919074 -0.21843016 0.32913822
## 
## Evaluated at:
## avar: x
##  a1 (intervened value of avar) = 1
##  a0 (reference value of avar)  = 0
## mvar: m
##  m_cde (intervend value of mvar for cde) = 1
## cvar: c
##  c_cond (covariate vector value) = 3
## 
## Note that effect estimates can vary over m_cde and c_cond values when interaction = TRUE.

There are \(A\times C\) term in both mediator and outcome models, and \(M\times C\) term in outcome model

Now suppose in addition to the EMM of treatment effect, the covariate \(C\) also modifies the mediator effect on the outcome. We add emm_mc_yreg = c("c") in regmedint(). Please note that the covariates in emm_mc_yreg should be a subset of the covariates specified in cvar, i.e. if a covariate is an effect measure modifier included in emm_mc_yreg, it must be included in cvar, otherwise an error message will be printed.

regmedint_obj4 <- regmedint(data = vv2015,
                            ## Variables
                            yvar = "y",
                            avar = "x",
                            mvar = "m",
                            cvar = c("c"),
                            emm_ac_mreg = c("c"),
                            emm_ac_yreg = c("c"),
                            emm_mc_yreg = c("c"),
                            eventvar = "event",
                            ## Values at which effects are evaluated
                            a0 = 0,
                            a1 = 1,
                            m_cde = 1,
                            c_cond = 3,
                            ## Model types
                            mreg = "logistic",
                            yreg = "survAFT_weibull",
                            ## Additional specification
                            interaction = TRUE,
                            casecontrol = FALSE)
summary(regmedint_obj4)
## ### Mediator model
## 
## Call:
## glm(formula = m ~ x + c + x:c, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.32727    0.34979  -0.936    0.349
## x            0.30431    0.56789   0.536    0.592
## c            0.24085    0.24688   0.976    0.329
## x:c          0.09216    0.44624   0.207    0.836
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 138.59  on 99  degrees of freedom
## Residual deviance: 136.04  on 96  degrees of freedom
## AIC: 144.04
## 
## Number of Fisher Scoring iterations: 4
## 
## ### Outcome model
## 
## Call:
## survival::survreg(formula = Surv(y, event) ~ x + m + x:m + c + 
##     x:c + m:c, data = data, dist = "weibull")
##               Value Std. Error     z       p
## (Intercept) -0.9959     0.2071 -4.81 1.5e-06
## x            0.4185     0.3354  1.25    0.21
## m           -0.0216     0.3112 -0.07    0.94
## c           -0.1339     0.1405 -0.95    0.34
## x:m          0.0905     0.4265  0.21    0.83
## x:c          0.0327     0.2242  0.15    0.88
## m:c          0.1275     0.1861  0.69    0.49
## Log(scale)  -0.0406     0.0814 -0.50    0.62
## 
## Scale= 0.96 
## 
## Weibull distribution
## Loglik(model)= -11.1   Loglik(intercept only)= -14.5
##  Chisq= 6.78 on 6 degrees of freedom, p= 0.34 
## Number of Newton-Raphson Iterations: 5 
## n= 100 
## 
## ### Mediation analysis 
##             est         se         Z         p      lower     upper
## cde  0.60705735 0.52594922 1.1542128 0.2484129 -0.4237842 1.6378989
## pnde 0.57902523 0.51447701 1.1254638 0.2603926 -0.4293312 1.5873816
## tnie 0.05333600 0.10591830 0.5035579 0.6145721 -0.1542601 0.2609321
## tnde 0.58889505 0.51488644 1.1437377 0.2527324 -0.4202638 1.5980539
## pnie 0.04346618 0.09107534 0.4772552 0.6331804 -0.1350382 0.2219706
## te   0.63236123 0.52776615 1.1981845 0.2308452 -0.4020414 1.6667639
## pm   0.11082259 0.20960355 0.5287248 0.5969964 -0.2999928 0.5216380
## 
## Evaluated at:
## avar: x
##  a1 (intervened value of avar) = 1
##  a0 (reference value of avar)  = 0
## mvar: m
##  m_cde (intervend value of mvar for cde) = 1
## cvar: c
##  c_cond (covariate vector value) = 3
## 
## Note that effect estimates can vary over m_cde and c_cond values when interaction = TRUE.