Chapter 11 Appendix B: Calculation of the true causal quantities
11.1 True causal quantities without mediator-ouctome confounder affected by the exposure
11.1.1 Average total effects (ATE)
The following function true.ATE1 can be used to run the calculation for the average total effects (ATE).
true.ATE1 <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1), c(0,1)), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum")
for (n in 1:16) {
S[n,"sum"] <- ( ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) ) -
( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
ATE.death <- sum(S[,"sum"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1), c(0,1)), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum")
for (n in 1:16) {
S[n,"sum"] <- ( ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) ) -
( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
ATE.qol <- sum(S[,"sum"])
return(list(ATE.death = ATE.death, ATE.qol = ATE.qol))
}The average total effects \(\text{ATE} = \mathbb{E}(Y_1) - \mathbb{E}(Y_0)\) are:
- 0.058 for death and -4.9 for quality of life without interaction;
- 0.06955 for death and -6.825 for quality of life with interaction.
11.1.2 Controlled direct effects (CDE)
The following function true.CDE1 can be used to run the calculation for controlled direct effects (CDE).
true.CDE1 <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# binary outcome (death)
# we estimate both CDE, fixing do(M) = 0 et do(M) = 1 and
# using the corresponding lines in the S matrix
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^3))
colnames(S) <- list("male","soc_env","L1","M","sum")
for (n in 1:16) {
S[n,"sum"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
CDE.M0.death <- sum(S[1:8,"sum"])
CDE.M1.death <- sum(S[9:16,"sum"])
# quantitative outcome (QoL)
# we estimate both CDE, fixing do(M) = 0 et do(M) = 1 and using
# the corresponding lines in the S matrix
for (n in 1:16) {
S[n,"sum"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
CDE.M0.qol <- sum(S[1:8,"sum"])
CDE.M1.qol <- sum(S[9:16,"sum"])
return(list(CDE.M0.death = CDE.M0.death, CDE.M1.death = CDE.M1.death,
CDE.M0.qol = CDE.M0.qol, CDE.M1.qol = CDE.M1.qol))
}Setting \(do(M=0)\), the controlled direct effects \(\text{CDE}_{M=0} = \mathbb{E}\left(Y_{1,0} \right) - \mathbb{E}\left(Y_{0,0} \right)\) are:
- 0.05 for death and -4 for quality of life without interaction,
- 0.05 for death and -4 for quality of life with interaction.
Setting \(do(M=1)\), the controlled direct effects \(\text{CDE}_{M=1} = \mathbb{E}\left(Y_{1,1} \right) - \mathbb{E}\left(Y_{0,1} \right)\) are:
- 0.05 for death and -4 for quality of life without interaction,
- 0.08 for death and -9 for quality of life with interaction.
11.1.3 Pure natural direct effect and Total natural indirect effect
The following function true.PNDE.TNIE1 can be used to run the calculation for pure natural direct effects (PNDE) and total natural indirect effects (TNIE).
true.PNDE.TNIE1 <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.pnde", "sum.tnie")
for (n in 1:16) {
# PNDE
S[n,"sum.pnde"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# TNIE
S[n,"sum.tnie"] <- ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
( (( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) +
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) -
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) -
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
PNDE.death <- sum(S[,"sum.pnde"])
TNIE.death <- sum(S[,"sum.tnie"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.pnde", "sum.tnie")
for (n in 1:16) {
# PNDE
S[n,"sum.pnde"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# TNIE
S[n,"sum.tnie"] <- ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
( (( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) +
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) -
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) -
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
PNDE.qol <- sum(S[,"sum.pnde"])
TNIE.qol <- sum(S[,"sum.tnie"])
return(list(PNDE.death = PNDE.death, TNIE.death = TNIE.death,
PNDE.qol = PNDE.qol, TNIE.qol = TNIE.qol))
}true.PNDE.TNIE.no.inter <- true.PNDE.TNIE1(interaction = 0)
true.PNDE.TNIE.with.inter <- true.PNDE.TNIE1(interaction = 1)The \(\text{PNDE}=\mathbb{E}\left( Y_{1,M_0}\right) - \mathbb{E}\left(Y_{0,M_0}\right)\) and \(\text{TNIE}=\mathbb{E}\left( Y_{1,M_1}\right) - \mathbb{E}\left(Y_{1,M_0}\right)\) are respectively:
- 0.05 and 0.008 for death without interaction,
- 0.05855 and 0.011 for death with interaction,
- -4 and -0.9 for quality of life without interaction,
- -5.425 and -1.4 for quality of life with interaction.
11.1.4 Total natural direct effect and Pure natural indirect effect
The following function true.TNDE.PNIE1 can be used to run the calculation for total natural direct effects (TNDE) and pure natural indirect effects (PNIE).
true.TNDE.PNIE1 <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.tnde", "sum.pnie")
for (n in 1:16) {
# TNDE
S[n,"sum.tnde"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# PNIE
S[n,"sum.pnie"] <- ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
( (( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) +
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) -
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) -
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
TNDE.death <- sum(S[,"sum.tnde"])
PNIE.death <- sum(S[,"sum.pnie"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.tnde", "sum.pnie")
for (n in 1:16) {
# TNDE
S[n,"sum.tnde"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# PNIE
S[n,"sum.pnie"] <- ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
( (( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) +
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) -
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) -
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
TNDE.qol <- sum(S[,"sum.tnde"])
PNIE.qol <- sum(S[,"sum.pnie"])
return(list(TNDE.death = TNDE.death, PNIE.death = PNIE.death,
TNDE.qol = TNDE.qol, PNIE.qol = PNIE.qol))
}true.TNDE.PNIE.no.inter <- true.TNDE.PNIE1(interaction = 0)
true.TNDE.PNIE.with.inter <- true.TNDE.PNIE1(interaction = 1)The \(\text{TNDE}=\mathbb{E}\left( Y_{1,M_1}\right) - \mathbb{E}\left(Y_{0,M_1}\right)\) and \(\text{PNIE}=\mathbb{E}\left( Y_{0,M_1}\right) - \mathbb{E}\left(Y_{0,M_0}\right)\) are respectively:
- 0.05 and 0.008 for death without interaction,
- 0.06155 and 0.008 for death with interaction,
- -4 and -0.9 for quality of life without interaction,
- -5.925 and -0.9 for quality of life with interaction.
11.1.5 Vanderweele’s 3-way decomposition
The following function true.3way.decomp can be used to run the calculation for the 3-way decomposition of the total effect into a “pure natural direct effect” (PNDE), a “pure natural indirect effect” (PNIE) and a “mediated interactive effect” (MIE).
true.3way.decomp <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.pde", "sum.pie")
for (n in 1:16) {
# PDE
S[n,"sum.pde"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# PIE
S[n,"sum.pie"] <- ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
( (( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) +
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) -
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) -
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
# MI
S.MI <- cbind(expand.grid(c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4))
colnames(S.MI) <- list("male","soc_env","L1", "sum.mi")
for (n in 1:8) {
S.MI[n,"sum.mi"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S.MI[n,"male"] +
b["b_soc_env_Y"] * S.MI[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S.MI[n,"L1"] +
b["b_M_Y"] * 1 +
b["b_AM_Y"] * 1 * 1 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S.MI[n,"male"] +
b["b_soc_env_Y"] * S.MI[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S.MI[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S.MI[n,"male"] +
b["b_soc_env_Y"] * S.MI[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S.MI[n,"L1"] +
b["b_M_Y"] * 1 +
b["b_AM_Y"] * 0 * 1 * b["A.M.inter"] ) +
( b["b_Y"] +
b["b_male_Y"] * S.MI[n,"male"] +
b["b_soc_env_Y"] * S.MI[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S.MI[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 0 * 0 * b["A.M.inter"] )) *
( ( b["b_M"] +
b["b_male_M"] * S.MI[n,"male"] +
b["b_soc_env_M"] * S.MI[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) -
( b["b_M"] +
b["b_male_M"] * S.MI[n,"male"] +
b["b_soc_env_M"] * S.MI[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S.MI[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S.MI[n,"male"])) *
((b["p_L0_soc_env"])^(S.MI[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S.MI[n,"soc_env"]))
}
PDE.death <- sum(S[,"sum.pde"])
PIE.death <- sum(S[,"sum.pie"])
MI.death <- sum(S.MI[,"sum.mi"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.pde", "sum.pie")
for (n in 1:16) {
# PDE
S[n,"sum.pde"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# PIE
S[n,"sum.pie"] <- ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
( (( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) +
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) -
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) -
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
# MI
S.MI <- cbind(expand.grid(c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4))
colnames(S.MI) <- list("male","soc_env","L1","sum.mi")
for (n in 1:8) {
S.MI[n,"sum.mi"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S.MI[n,"male"] +
b["c_soc_env_Y"] * S.MI[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S.MI[n,"L1"] +
b["c_M_Y"] * 1 +
b["c_AM_Y"] * 1 * 1 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S.MI[n,"male"] +
b["c_soc_env_Y"] * S.MI[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S.MI[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S.MI[n,"male"] +
b["c_soc_env_Y"] * S.MI[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S.MI[n,"L1"] +
b["c_M_Y"] * 1 +
b["c_AM_Y"] * 0 * 1 * b["A.M.inter"] ) +
( b["mu_Y"] +
b["c_male_Y"] * S.MI[n,"male"] +
b["c_soc_env_Y"] * S.MI[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S.MI[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 0 * 0 * b["A.M.inter"] )) *
( ( b["b_M"] +
b["b_male_M"] * S.MI[n,"male"] +
b["b_soc_env_M"] * S.MI[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) -
( b["b_M"] +
b["b_male_M"] * S.MI[n,"male"] +
b["b_soc_env_M"] * S.MI[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S.MI[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S.MI[n,"male"])) *
((b["p_L0_soc_env"])^(S.MI[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S.MI[n,"soc_env"]))
}
PDE.qol <- sum(S[,"sum.pde"])
PIE.qol <- sum(S[,"sum.pie"])
MI.qol <- sum(S.MI[,"sum.mi"])
return(list(PDE.death = PDE.death, PIE.death = PIE.death, MI.death = MI.death,
PDE.qol = PDE.qol, PIE.qol = PIE.qol, MI.qol = MI.qol))
}true.3way.no.inter <- true.3way.decomp(interaction = 0)
true.3way.with.inter <- true.3way.decomp(interaction = 1)The \(\text{PNDE}=\mathbb{E}\left( Y_{1,M_0}\right) - \mathbb{E}\left( Y_{0,M_0}\right)\), the \(\text{PNIE}=\mathbb{E}\left(Y_{0,M_1}\right) - \mathbb{E}\left(Y_{0,M_0}\right)\) and the \(\text{MIE}=\mathbb{E}\left( (Y_{1,1} - Y_{1,0} - Y_{0,1} + Y_{0,0}) \times (M_1 - M_0) \right)\) are respectively:
- 0.05, 0.008 and 0.000 for death without interaction,
- 0.05855, 0.008 and 0.003 for death with interaction,
- -4, -0.9 and 0 for quality of life without interaction,
- -5.425, -0.9 and -0.5 for quality of life with interaction.
11.1.6 Vanderweele’s 4-way decomposition
The following function true.4way.decomp can be used to run the calculation for the 4-way decomposition of the total effect into a “controlled direct effect” (CDE), a “reference interaction effect” (RIE), a “mediated interaction effect” (MIE) and a “pure natural indirect effect” (PNIE).
true.4way.decomp <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1)), rep(NA,n=2^3), rep(NA,n=2^3),
rep(NA,n=2^3), rep(NA,n=2^3))
colnames(S) <- list("male","soc_env", "L1", "sum.cde", "sum.intref",
"sum.intmed", "sum.pie")
for (n in 1:8) {
# CDE
S[n,"sum.cde"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 0 * 0 * b["A.M.inter"] ) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# INTref
S[n,"sum.intref"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 1 +
b["b_AM_Y"] * 1 * 1 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 1 +
b["b_AM_Y"] * 0 * 1 * b["A.M.inter"] ) +
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 0 * 0 * b["A.M.inter"] )) *
( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# INTmed
S[n,"sum.intmed"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 1 +
b["b_AM_Y"] * 1 * 1 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 1 +
b["b_AM_Y"] * 0 * 1 * b["A.M.inter"] ) +
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 0 * 0 * b["A.M.inter"] )) *
( ( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) -
( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# PIE
S[n,"sum.pie"] <- ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 1 +
b["b_AM_Y"] * 0 * 1 * b["A.M.inter"] ) -
( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * 0 +
b["b_AM_Y"] * 0 * 0 * b["A.M.inter"] ) ) *
( ( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) -
( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
CDE.death <- sum(S[,"sum.cde"])
INTref.death <- sum(S[,"sum.intref"])
INTmed.death <- sum(S[,"sum.intmed"])
PIE.death <- sum(S[,"sum.pie"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1)), rep(NA,n=2^3), rep(NA,n=2^3),
rep(NA,n=2^3), rep(NA,n=2^3))
colnames(S) <- list("male","soc_env", "L1", "sum.cde", "sum.intref",
"sum.intmed", "sum.pie")
for (n in 1:8) {
# CDE
S[n,"sum.cde"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 0 * 0 * b["A.M.inter"] ) ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# INTref
S[n,"sum.intref"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 1 +
b["c_AM_Y"] * 1 * 1 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 1 +
b["c_AM_Y"] * 0 * 1 * b["A.M.inter"] ) +
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 0 * 0 * b["A.M.inter"] )) *
( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# INTmed
S[n,"sum.intmed"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 1 +
b["c_AM_Y"] * 1 * 1 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 1 * 0 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 1 +
b["c_AM_Y"] * 0 * 1 * b["A.M.inter"] ) +
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 0 * 0 * b["A.M.inter"] )) *
( ( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) -
( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
# PIE
S[n,"sum.pie"] <- ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 1 +
b["c_AM_Y"] * 0 * 1 * b["A.M.inter"] ) -
( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * 0 +
b["c_AM_Y"] * 0 * 0 * b["A.M.inter"] ) ) *
( ( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) -
( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
CDE.qol <- sum(S[,"sum.cde"])
INTref.qol <- sum(S[,"sum.intref"])
INTmed.qol <- sum(S[,"sum.intmed"])
PIE.qol <- sum(S[,"sum.pie"])
return(list(CDE.death = CDE.death, INTref.death = INTref.death,
INTmed.death = INTmed.death, PIE.death = PIE.death,
CDE.qol = CDE.qol, INTref.qol = INTref.qol,
INTmed.qol = INTmed.qol, PIE.qol = PIE.qol))
}true.4way.no.inter <- true.4way.decomp(interaction = 0)
true.4way.with.inter <- true.4way.decomp(interaction = 1)The \(\text{CDE}=\mathbb{E}\left(Y_{1,0} \right) - \mathbb{E}\left(Y_{0,0}\right)\), \(\text{RIE}=\left((Y_{1,1} - Y_{1,0} - Y_{0,1} + Y_{0,0}) \times M_0\right)\), \(\text{MIE}= \mathbb{E}\left((Y_{1,1} - Y_{1,0} - Y_{0,1} + Y_{0,0}) \times (M_1 - M_0)\right)\) and \(\text{PNIE}=\mathbb{E}\left(Y_{0,M_1}\right) - \mathbb{E}\left(Y_{0,M_0}\right)\) are respectively:
- 0.05, 0.000, 0.000 and 0.008 for death without interaction,
- 0.05, 0.00855, 0.003 and 0.008 for death with interaction,
- -4, 0, 0 and -0.9 for quality of life without interaction,
- -4, -1.425, -0.5 and -0.9 for quality of life with interaction.
11.1.7 Marginal randomized direct and indirect effects
The following function true.marg.random can be used to run the calculation for the marginal randomized natural direct (marginal MRDE) and indirect effects (marginal MRIE).
true.marg.random <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# marginal distribution of M
M.S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^5))
colnames(M.S) <- list("male","soc_env","L1","M","A","sum")
for (n in 1:32) {
M.S[n,"sum"] <- (( b["b_M"] +
b["b_male_M"] * M.S[n,"male"] +
b["b_soc_env_M"] * M.S[n,"soc_env"] +
b["b_L1_M"] * M.S[n,"L1"] +
b["b_A_M"] * M.S[n,"A"])^( M.S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * M.S[n,"male"] +
b["b_soc_env_M"] * M.S[n,"soc_env"] +
b["b_L1_M"] * M.S[n,"L1"] +
b["b_A_M"] * M.S[n,"A"]) )^( 1 - M.S[n,"M"] ))
}
M0.A0.L0_00.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M0.A0.L0_01.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M0.A0.L0_10.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M0.A0.L0_11.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M0.A0.L0_00.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M0.A0.L0_01.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
M0.A0.L0_10.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M0.A0.L0_11.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
M1.A0.L0_00.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M1.A0.L0_01.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M1.A0.L0_10.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M1.A0.L0_11.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M1.A0.L0_00.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M1.A0.L0_01.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
M1.A0.L0_10.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M1.A0.L0_11.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
M0.A1.L0_00.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M0.A1.L0_01.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M0.A1.L0_10.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M0.A1.L0_11.L1_0 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M0.A1.L0_00.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M0.A1.L0_01.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
M0.A1.L0_10.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M0.A1.L0_11.L1_1 <- M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
M1.A1.L0_00.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M1.A1.L0_01.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M1.A1.L0_10.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==0,"sum"]
M1.A1.L0_11.L1_0 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==0,"sum"]
M1.A1.L0_00.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M1.A1.L0_01.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
M1.A1.L0_10.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0 & M.S[,"L1"]==1,"sum"]
M1.A1.L0_11.L1_1 <- M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1 & M.S[,"L1"]==1,"sum"]
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4), rep(NA,n=2^4),
rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11", "sum.psi10", "sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A1.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) + # A'=1
M1.A1.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A1.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M1.A1.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A1.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A1.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M1.A1.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A1.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( S[n,"M"] )) *
((M0.A1.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A1.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A1.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A1.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A1.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A1.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M0.A1.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A1.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(M.S[n,"L1"])) *
((1 - b["p_L1"])^(1 - M.S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) + # A'=0
M1.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M1.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M1.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( S[n,"M"] )) *
((M0.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M0.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(M.S[n,"L1"])) *
((1 - b["p_L1"])^(1 - M.S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["b_Y"] + # A=0
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) + # A'=0
M1.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M1.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M1.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( S[n,"M"] )) *
((M0.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M0.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(M.S[n,"L1"])) *
((1 - b["p_L1"])^(1 - M.S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
mrNDE.death <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
mrNIE.death <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4),
rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11", "sum.psi10",
"sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A1.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) + # A'=1
M1.A1.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A1.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M1.A1.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A1.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A1.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M1.A1.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A1.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( S[n,"M"] )) *
((M0.A1.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A1.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A1.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A1.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A1.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A1.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M0.A1.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A1.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(M.S[n,"L1"])) *
((1 - b["p_L1"])^(1 - M.S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) + # A'=0
M1.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M1.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M1.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( S[n,"M"] )) *
((M0.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M0.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(M.S[n,"L1"])) *
((1 - b["p_L1"])^(1 - M.S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["mu_Y"] + # A=0
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) + # A'=0
M1.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M1.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M1.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M1.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M1.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( S[n,"M"] )) *
((M0.A0.L0_00.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_01.L1_0*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_10.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==0) +
M0.A0.L0_11.L1_0*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*(S[n,"L1"]==0) +
M0.A0.L0_00.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_01.L1_1*(S[n,"male"]==0)*(S[n,"soc_env"]==1)*(S[n,"L1"]==1) +
M0.A0.L0_10.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==0)*(S[n,"L1"]==1) +
M0.A0.L0_11.L1_1*(S[n,"male"]==1)*(S[n,"soc_env"]==1)*
(S[n,"L1"]==1) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(M.S[n,"L1"])) *
((1 - b["p_L1"])^(1 - M.S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
mrNDE.qol <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
mrNIE.qol <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
return(list(mrNDE.death = mrNDE.death, mrNIE.death = mrNIE.death,
mrNDE.qol = mrNDE.qol, mrNIE.qol = mrNIE.qol))
}true.marg.random.no.inter <- true.marg.random(interaction = 0)
true.marg.random.with.inter <- true.marg.random(interaction = 1)The marginal randomized direct effect \(\text{MRDE}=\mathbb{E}\left(Y_{1,G_{0|L(0)}}\right) - \mathbb{E}\left(Y_{0,G_{0|L(0)}}\right)\) and the marginal randomized indirect effect \(\text{MRIE}=\mathbb{E}\left(Y_{1,G_{1|L(0)}}\right) - \mathbb{E}\left(Y_{1,G_{0|L(0)}}\right)\) are respectively:
- 0.05 and 0.008 for death without interaction,
- 0.05855 and 0.011 for death with interaction,
- -4 and -0.9 for quality of life without interaction,
- -5.425 and -1.4 for quality of life with interaction.
11.1.8 Conditional randomized direct and indirect effects
The following function true.cond.random can be used to run the calculation for the conditional randomized natural direct (CRDE) and indirect effects (CRIE).
true.cond.random <- function(interaction = NULL) {
b <- param.causal.model.1(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4),
rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11", "sum.psi10",
"sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=1
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["b_Y"] + # A=0
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
crNDE.death <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
crNIE.death <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4),
rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11", "sum.psi10",
"sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=1
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["mu_Y"] + # A=0
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
((b["p_L1"])^(S[n,"L1"])) *
((1 - b["p_L1"])^(1 - S[n,"L1"])) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
crNDE.qol <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
crNIE.qol <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
return(list(crNDE.death = crNDE.death, crNIE.death = crNIE.death,
crNDE.qol = crNDE.qol, crNIE.qol = crNIE.qol))
}true.cond.random.no.inter <- true.cond.random(interaction = 0)
true.cond.random.with.inter <- true.cond.random(interaction = 1)The conditional randomized direct effect \(\text{CRDE}=\mathbb{E}\left(Y_{1,\Gamma_{0\mid L(0),L(1)}}\right) - \mathbb{E}\left(Y_{0,\Gamma_{0\mid L(0),L(1)}}\right)\) and conditional randomized indirect effect \(\text{CRIE}=\mathbb{E}\left(Y_{1,\Gamma_{1\mid L(0),L(1)}}\right) - \mathbb{E}\left(Y_{1,\Gamma_{0\mid L(0),L(1)}}\right)\) are respectively:
- 0.05 and 0.008 for death without interaction,
- 0.05855 and 0.011 for death with interaction,
- -4 and -0.9 for quality of life without interaction,
- -5.425 and -1.4 for quality of life with interaction.
| Effects | Without \(A*M\) interaction | with \(A*M\) interaction |
|---|---|---|
| Binary outcome | ||
| Average total effect (ATE) | 0.058 | 0.06955 |
| Controlled direct effect (CDE) | ||
| - CDE, setting do(M=0) | 0.05 | 0.05 |
| - CDE, setting do(M=1) | 0.05 | 0.08 |
| Pure NDE and Total NIE | ||
| - PNDE | 0.05 | 0.05855 |
| - TNIE | 0.008 | 0.011 |
| Total NDE and Pure NIE | ||
| - TNDE | 0.05 | 0.06155 |
| - PNIE | 0.008 | 0.008 |
| 3-way decomposition | ||
| - PDE | 0.05 | 0.05855 |
| - PIE | 0.008 | 0.008 |
| - MI | 0.000 | 0.003 |
| 4-way decomposition | ||
| - CDE | 0.05 | 0.05 |
| - INTref | 0.000 | 0.00855 |
| - INTmed | 0.000 | 0.003 |
| - PIE | 0.008 | 0.008 |
| Marginal randomized | ||
| - marginal rNDE | 0.05 | 0.05855 |
| - marginal rNIE | 0.008 | 0.011 |
| Conditional randomized | ||
| - conditional rNDE | 0.05 | 0.05855 |
| - conditional rNIE | 0.008 | 0.011 |
| Quantitative outcome | ||
| Average total effect (ATE) | -4.9 | -6.825 |
| Controlled direct effect (CDE) | ||
| - CDE, setting do(M=0) | -4 | -4 |
| - CDE, setting do(M=1) | -4 | -9 |
| Pure NDE and Total NIE | ||
| - PNDE | -4 | -5.425 |
| - TNIE | -0.9 | -1.4 |
| Total NDE and Pure NIE | ||
| - TNDE | -4 | -5.925 |
| - PNIE | -0.9 | -0.9 |
| 3-way decomposition | ||
| - PDE | -4 | -5.425 |
| - PIE | -0.9 | -0.9 |
| - MI | 0 | -0.5 |
| 4-way decomposition | ||
| - CDE | -4 | -4 |
| - INTref | 0.000 | -1.425 |
| - INTmed | 0.000 | -0.5 |
| - PIE | -0.9 | -0.9 |
| Marginal randomized | ||
| - marginal rNDE | -4 | -5.425 |
| - marginal rNIE | -0.9 | -1.4 |
| Conditional randomized | ||
| - conditional rNDE | -4 | -5.425 |
| - conditional rNIE | -0.9 | -1.4 |
11.2 True causal quantities with mediator-outcome confounder affected by the exposure
11.2.1 Average total effects (ATE)
The following function true.ATE.tv.conf can be used to run the calculation for the average total effects (ATE).
true.ATE.time.var.conf <- function(interaction = NULL) {
b <- param.causal.model.2(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1), c(0,1)), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum")
for (n in 1:16) {
S[n,"sum"] <- ( ( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) ) -
( ( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) ) ) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
ATE.death <- sum(S[,"sum"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1), c(0,1)), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum")
for (n in 1:16) {
S[n,"sum"] <- ( ( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) ) -
( ( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0) )^( 1 - S[n,"M"] )) ) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) ) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
ATE.qol <- sum(S[,"sum"])
return(list(ATE.death = ATE.death, ATE.qol = ATE.qol))
}true.ATE2.no.inter <- true.ATE.time.var.conf(interaction = 0)
true.ATE2.with.inter <- true.ATE.time.var.conf(interaction = 1)The average total effects \(\text{ATE} = \mathbb{E}(Y_1) - \mathbb{E}(Y_0)\) are:
- 0.0752 for death and -6.26 for quality of life without interaction
- 0.089282 for death and -8.607 for quality of life with interaction
11.2.2 Controlled direct effects (CDE)
The following function true.CDE.time.var can be used to run the calculation for controlled direct effects (CDE).
true.CDE.time.var <- function(interaction = NULL) {
b <- param.causal.model.2(A.M.interaction = interaction)
# binary outcome (death)
# we estimate both CDE, fixing do(M) = 0 et do(M) = 1 and using the
# corresponding lines in the S matrix
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^3))
colnames(S) <- list("male","soc_env","L1","M","sum")
for (n in 1:16) {
S[n,"sum"] <- ( (( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) ) -
(( b["b_Y"] +
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) ) ) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
CDE.M0.death <- sum(S[1:8,"sum"])
CDE.M1.death <- sum(S[9:16,"sum"])
# quantitative outcome (QoL)
# we estimate both CDE, fixing do(M) = 0 et do(M) = 1 and using the
# corresponding lines in the S matrix
for (n in 1:16) {
S[n,"sum"] <- ( (( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] ))) -
(( b["mu_Y"] +
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) ) ) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
CDE.M0.qol <- sum(S[1:8,"sum"])
CDE.M1.qol <- sum(S[9:16,"sum"])
return(list(CDE.M0.death = CDE.M0.death, CDE.M1.death = CDE.M1.death,
CDE.M0.qol = CDE.M0.qol, CDE.M1.qol = CDE.M1.qol))
}true.CDE2.no.inter <- true.CDE.time.var(interaction = 0)
true.CDE2.with.inter <- true.CDE.time.var(interaction = 1)Setting \(do(M=0)\), the controlled direct effects \(\text{CDE}_{M=0} = \mathbb{E}\left(Y_{1,0} \right) - \mathbb{E}\left(Y_{0,0} \right)\) are:
- 0.064 for death and -5 for quality of life without interaction
- 0.064 for death and -5 for quality of life with interaction
Setting \(do(M=1)\), the controlled direct effects \(\text{CDE}_{M=1} = \mathbb{E}\left(Y_{1,1} \right) - \mathbb{E}\left(Y_{0,1} \right)\) are:
- 0.064 for death and -5 for quality of life without interaction
- 0.094 for death and -10 for quality of life with interaction
11.2.3 Marginal randomized direct and indirect effects
The following function true.marg.random.time.var can be used to run the calculation for the marginal randomized natural direct (marginal MRDE) and indirect effects (marginal MRIE).
true.marg.random.time.var <- function(interaction = NULL) {
b <- param.causal.model.2(A.M.interaction = interaction)
# marginal distribution of M
M.S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^5))
colnames(M.S) <- list("male","soc_env","L1","M","A","sum")
for (n in 1:32) {
M.S[n,"sum"] <- (( b["b_M"] +
b["b_male_M"] * M.S[n,"male"] +
b["b_soc_env_M"] * M.S[n,"soc_env"] +
b["b_L1_M"] * M.S[n,"L1"] +
b["b_A_M"] * M.S[n,"A"])^( M.S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * M.S[n,"male"] +
b["b_soc_env_M"] * M.S[n,"soc_env"] +
b["b_L1_M"] * M.S[n,"L1"] +
b["b_A_M"] * M.S[n,"A"]) )^( 1 - M.S[n,"M"] )) *
(( b["b_L1"] +
b["b_male_L1"] * M.S[n,"male"] +
b["b_soc_env_L1"] * M.S[n,"soc_env"] +
b["b_A_L1"] * M.S[n,"A"])^( M.S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * M.S[n,"male"] +
b["b_soc_env_L1"] * M.S[n,"soc_env"] +
b["b_A_L1"] * M.S[n,"A"]))^( 1 - M.S[n,"L1"] ))
}
M0.A0.L00 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0,"sum"])
M0.A0.L01 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1,"sum"])
M0.A0.L10 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0,"sum"])
M0.A0.L11 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1,"sum"])
M1.A0.L00 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0,"sum"])
M1.A0.L01 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1,"sum"])
M1.A0.L10 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0,"sum"])
M1.A0.L11 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==0 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1,"sum"])
M0.A1.L00 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0,"sum"])
M0.A1.L01 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1,"sum"])
M0.A1.L10 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0,"sum"])
M0.A1.L11 <- sum(M.S[M.S[,"M"]==0 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1,"sum"])
M1.A1.L00 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==0,"sum"])
M1.A1.L01 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==0 &
M.S[,"soc_env"]==1,"sum"])
M1.A1.L10 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==0,"sum"])
M1.A1.L11 <- sum(M.S[M.S[,"M"]==1 & M.S[,"A"]==1 & M.S[,"male"]==1 &
M.S[,"soc_env"]==1,"sum"])
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4),
rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11",
"sum.psi10", "sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A1.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) + # A'=1
M1.A1.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M1.A1.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M1.A1.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( S[n,"M"] )) *
((M0.A1.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) +
M0.A1.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M0.A1.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M0.A1.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) + # A'=0
M1.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M1.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M1.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( S[n,"M"] )) *
((M0.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) +
M0.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M0.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M0.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["b_Y"] + # A=0
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) + # A'=0
M1.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M1.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M1.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( S[n,"M"] )) *
((M0.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) +
M0.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M0.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M0.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=0
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
mrNDE.death <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
mrNIE.death <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4),
rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11", "sum.psi10",
"sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A1.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) + # A'=1
M1.A1.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M1.A1.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M1.A1.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( S[n,"M"] )) *
((M0.A1.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) +
M0.A1.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M0.A1.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M0.A1.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) + # A'=0
M1.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M1.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M1.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( S[n,"M"] )) *
((M0.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) +
M0.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M0.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M0.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["mu_Y"] + # A=0
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
((M1.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) + # A'=0
M1.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M1.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M1.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( S[n,"M"] )) *
((M0.A0.L00*(S[n,"male"]==0)*(S[n,"soc_env"]==0) +
M0.A0.L01*(S[n,"male"]==0)*(S[n,"soc_env"]==1) +
M0.A0.L10*(S[n,"male"]==1)*(S[n,"soc_env"]==0) +
M0.A0.L11*(S[n,"male"]==1)*(S[n,"soc_env"]==1) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=0
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
mrNDE.qol <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
mrNIE.qol <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
return(list(mrNDE.death = mrNDE.death, mrNIE.death = mrNIE.death,
mrNDE.qol = mrNDE.qol, mrNIE.qol = mrNIE.qol))
}true.marg.random2.no.inter <- true.marg.random.time.var(interaction = 0)
true.marg.random2.with.inter <- true.marg.random.time.var(interaction = 1)The marginal randomized direct effect \(\text{MRDE}=\mathbb{E}\left(Y_{1,G_{0|L(0)}}\right) - \mathbb{E}\left(Y_{0,G_{0|L(0)}}\right)\) and the marginal randomized indirect effect \(\text{MRIE}=\mathbb{E}\left(Y_{1,G_{1|L(0)}}\right) - \mathbb{E}\left(Y_{1,G_{0|L(0)}}\right)\) are respectively:
- 0.064 and 0.0112 for death without interaction
- 0.073882 and 0.0154 for death with interaction
- -5 and -1.26 for quality of life without interaction
- -6.647 and -1.96 for quality of life with interaction
11.2.4 Conditional randomized direct and indirect effects
The following function true.cond.random.time.var can be used to run the calculation for the conditional randomized natural direct (CRDE) and the conditional randomized indirect effects (CRIE).
true.cond.random.time.var <- function(interaction = NULL) {
b <- param.causal.model.2(A.M.interaction = interaction)
# binary outcome (death)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4),
rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11", "sum.psi10",
"sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=1
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["b_Y"] + # A=1
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 1 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["b_Y"] + # A=0
b["b_male_Y"] * S[n,"male"] +
b["b_soc_env_Y"] * S[n,"soc_env"] +
b["b_A_Y"] * 0 +
b["b_L1_Y"] * S[n,"L1"] +
b["b_M_Y"] * S[n,"M"] +
b["b_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=0
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
crNDE.death <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
crNIE.death <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
# quantitative outcome (QoL)
S <- cbind(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1)), rep(NA,n=2^4),
rep(NA,n=2^4), rep(NA,n=2^4))
colnames(S) <- list("male","soc_env","L1","M","sum.psi11", "sum.psi10",
"sum.psi00")
for (n in 1:16) {
S[n,"sum.psi11"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=1
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 1 ) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi10"] <- ( b["mu_Y"] + # A=1
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 1 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 1 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=1
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 1))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
S[n,"sum.psi00"] <- ( b["mu_Y"] + # A=0
b["c_male_Y"] * S[n,"male"] +
b["c_soc_env_Y"] * S[n,"soc_env"] +
b["c_A_Y"] * 0 +
b["c_L1_Y"] * S[n,"L1"] +
b["c_M_Y"] * S[n,"M"] +
b["c_AM_Y"] * 0 * S[n,"M"] * b["A.M.inter"] ) *
(( b["b_M"] + # A'=0
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 )^( S[n,"M"] )) *
(( 1 - (b["b_M"] +
b["b_male_M"] * S[n,"male"] +
b["b_soc_env_M"] * S[n,"soc_env"] +
b["b_L1_M"] * S[n,"L1"] +
b["b_A_M"] * 0 ) )^( 1 - S[n,"M"] )) *
(( b["b_L1"] + # A=0
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0)^( S[n,"L1"] )) *
(( 1 - ( b["b_L1"] +
b["b_male_L1"] * S[n,"male"] +
b["b_soc_env_L1"] * S[n,"soc_env"] +
b["b_A_L1"] * 0))^( 1 - S[n,"L1"] )) *
((b["p_L0_male"])^(S[n,"male"])) *
((1 - b["p_L0_male"])^(1 - S[n,"male"])) *
((b["p_L0_soc_env"])^(S[n,"soc_env"])) *
((1 - b["p_L0_soc_env"])^(1 - S[n,"soc_env"]))
}
crNDE.qol <- sum(S[,"sum.psi10"]) - sum(S[,"sum.psi00"])
crNIE.qol <- sum(S[,"sum.psi11"]) - sum(S[,"sum.psi10"])
return(list(crNDE.death = crNDE.death, crNIE.death = crNIE.death,
crNDE.qol = crNDE.qol, crNIE.qol = crNIE.qol))
}true.cond.random2.no.inter <- true.cond.random.time.var(interaction = 0)
true.cond.random2.with.inter <- true.cond.random.time.var(interaction = 1)The conditional randomized direct effect \(\text{CRDE}=\mathbb{E}\left(Y_{1,\Gamma_{0\mid L(0),L(1)}}\right) - \mathbb{E}\left(Y_{0,\Gamma_{0\mid L(0),L(1)}}\right)\) and conditional randomized indirect effect \(\text{CRIE}=\mathbb{E}\left(Y_{1,\Gamma_{1\mid L(0),L(1)}}\right) - \mathbb{E}\left(Y_{1,\Gamma_{0\mid L(0),L(1)}}\right)\) are respectively:
- 0.0672 and 0.008 for death without interaction,
- 0.078282 and 0.011 for death with interaction,
- -5.36 and -0.9 for quality of life without interaction,
- -7.207 and -1.4 for quality of life with interaction.
| Effects | Without \(A*M\) interaction | with \(A*M\) interaction |
|---|---|---|
| Binary outcome | ||
| Average total effect (ATE) | 0.0752 | 0.089282 |
| Controlled direct effect (CDE) | ||
| - CDE, setting do(M=0) | 0.064 | 0.064 |
| - CDE, setting do(M=1) | 0.064 | 0.094 |
| Marginal randomized | ||
| - marginal rNDE | 0.064 | 0.073882 |
| - marginal rNIE | 0.0112 | 0.0154 |
| Conditional randomized | ||
| - conditional rNDE | 0.0672 | 0.078282 |
| - conditional rNIE | 0.008 | 0.011 |
| Quantitative outcome | ||
| Average total effect (ATE) | -6.26 | -8.607 |
| Controlled direct effect (CDE) | ||
| - CDE, setting do(M=0) | -5 | -5 |
| - CDE, setting do(M=1) | -5 | -10 |
| Marginal randomized | ||
| - marginal rNDE | -5 | -6.647 |
| - marginal rNIE | -1.26 | -1.96 |
| Conditional randomized | ||
| - conditional rNDE | -5.36 | -7.207 |
| - conditional rNIE | -0.9 | -1.4 |