Implementation of formulas

Kazuki Yoshida

2021-05-11

The bodies of the following functions contain the R implementation of the formulas in V2015. Although many models are covered, there are only four patterns. See the package top page for further references. The code seen here is compacted and lacks the comments. For more redable code, please refer to the Github repo and search for the function names without the preceding regmedint:::. For the type-set mathematical expressions and LaTeX source, please see the supplement.

mreg linear yreg linear (V2015 p466 Proposition 2.3)

These functions are only used in the setting where both the mediator model and the outcome model are linear regression.

Point estimates

regmedint:::calc_myreg_mreg_linear_yreg_linear_est
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4)
##     fun_est <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         cde <- (theta1 + (theta3 * m_cde)) * (a1 - a0)
##         pnde <- (theta1 + (theta3 * beta0) + (theta3 * beta1 * 
##             a0) + (theta3 * beta2_c)) * (a1 - a0)
##         tnie <- ((theta2 * beta1) + (theta3 * beta1 * a1)) * 
##             (a1 - a0)
##         tnde <- (theta1 + (theta3 * beta0) + (theta3 * beta1 * 
##             a1) + (theta3 * beta2_c)) * (a1 - a0)
##         pnie <- ((theta2 * beta1) + (theta3 * beta1 * a0)) * 
##             (a1 - a0)
##         te <- pnde + tnie
##         pm <- tnie/te
##         c(cde = unname(cde), pnde = unname(pnde), tnie = unname(tnie), 
##             tnde = unname(tnde), pnie = unname(pnie), te = unname(te), 
##             pm = unname(pm))
##     }
##     return(fun_est)
## }
## <bytecode: 0x7f9f11272008>
## <environment: namespace:regmedint>

Standard error estimates

regmedint:::calc_myreg_mreg_linear_yreg_linear_se
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4, Sigma_beta, Sigma_theta) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4)
##     validate_myreg_vcovs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4, Sigma_beta = Sigma_beta, Sigma_theta = Sigma_theta)
##     Sigma <- Matrix::bdiag(Sigma_beta, Sigma_theta)
##     fun_se <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         Gamma_cde <- matrix(c(0, 0, rep(0, length(beta2)), 0, 
##             1, 0, m_cde, rep(0, length(theta4))))
##         Gamma_pnde <- matrix(c(theta3, (theta3 * a0), (theta3 * 
##             c_cond), 0, 1, 0, (beta0 + (beta1 * a0) + beta2_c), 
##             rep(0, length(theta4))))
##         Gamma_tnie <- matrix(c(0, (theta2 + (theta3 * a1)), rep(0, 
##             length(beta2)), 0, 0, beta1, (beta1 * a1), rep(0, 
##             length(theta4))))
##         Gamma_tnde <- matrix(c(theta3, (theta3 * a1), (theta3 * 
##             c_cond), 0, 1, 0, (beta0 + (beta1 * a1) + beta2_c), 
##             rep(0, length(theta4))))
##         Gamma_pnie <- matrix(c(0, (theta2 + (theta3 * a0)), rep(0, 
##             length(beta2)), 0, 0, beta1, (beta1 * a0), rep(0, 
##             length(theta4))))
##         Gamma_te <- Gamma_pnde + Gamma_tnie
##         pnde <- (theta1 + (theta3 * beta0) + (theta3 * beta1 * 
##             a0) + (theta3 * beta2_c)) * (a1 - a0)
##         tnie <- ((theta2 * beta1) + (theta3 * beta1 * a1)) * 
##             (a1 - a0)
##         d_pm <- grad_prop_med_yreg_linear(pnde = unname(pnde), 
##             tnie = unname(tnie))
##         Gamma_pm <- (d_pm[["pnde"]] * Gamma_pnde) + (d_pm[["tnie"]] * 
##             Gamma_tnie)
##         a1_sub_a0 <- abs(a1 - a0)
##         se_cde <- sqrt(as.numeric(t(Gamma_cde) %*% Sigma %*% 
##             Gamma_cde)) * a1_sub_a0
##         se_pnde <- sqrt(as.numeric(t(Gamma_pnde) %*% Sigma %*% 
##             Gamma_pnde)) * a1_sub_a0
##         se_tnie <- sqrt(as.numeric(t(Gamma_tnie) %*% Sigma %*% 
##             Gamma_tnie)) * a1_sub_a0
##         se_tnde <- sqrt(as.numeric(t(Gamma_tnde) %*% Sigma %*% 
##             Gamma_tnde)) * a1_sub_a0
##         se_pnie <- sqrt(as.numeric(t(Gamma_pnie) %*% Sigma %*% 
##             Gamma_pnie)) * a1_sub_a0
##         se_te <- sqrt(as.numeric(t(Gamma_te) %*% Sigma %*% Gamma_te)) * 
##             a1_sub_a0
##         se_pm <- sqrt(as.numeric(t(Gamma_pm) %*% Sigma %*% Gamma_pm)) * 
##             a1_sub_a0
##         c(se_cde = unname(se_cde), se_pnde = unname(se_pnde), 
##             se_tnie = unname(se_tnie), se_tnde = unname(se_tnde), 
##             se_pnie = unname(se_pnie), se_te = unname(se_te), 
##             se_pm = unname(se_pm))
##     }
##     return(fun_se)
## }
## <bytecode: 0x7f9f2623c550>
## <environment: namespace:regmedint>

mreg linear yreg non-linear (V2015 p468 Proposition 2.4)

These functions are used in all cases where the mediator model is linear regression and the outcome model is any one of the non-linear models.

Point estimates

regmedint:::calc_myreg_mreg_linear_yreg_logistic_est
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4, sigma_sq) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4, sigma_sq = sigma_sq)
##     fun_est <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         cde <- (theta1 + (theta3 * m_cde)) * (a1 - a0)
##         pnde <- (theta1 + (theta3 * beta0) + (theta3 * beta1 * 
##             a0) + (theta3 * beta2_c) + (theta3 * theta2 * sigma_sq)) * 
##             (a1 - a0) + ((1/2) * theta3^2 * sigma_sq) * (a1^2 - 
##             a0^2)
##         tnie <- ((theta2 * beta1) + theta3 * beta1 * a1) * (a1 - 
##             a0)
##         tnde <- (theta1 + (theta3 * beta0) + (theta3 * beta1 * 
##             a1) + (theta3 * beta2_c) + (theta3 * theta2 * sigma_sq)) * 
##             (a1 - a0) + ((1/2) * theta3^2 * sigma_sq) * (a1^2 - 
##             a0^2)
##         pnie <- ((theta2 * beta1) + theta3 * beta1 * a0) * (a1 - 
##             a0)
##         te <- pnde + tnie
##         pm <- (exp(pnde) * (exp(tnie) - 1))/(exp(te) - 1)
##         c(cde = unname(cde), pnde = unname(pnde), tnie = unname(tnie), 
##             tnde = unname(tnde), pnie = unname(pnie), te = unname(te), 
##             pm = unname(pm))
##     }
##     return(fun_est)
## }
## <bytecode: 0x7f9f12ff34e8>
## <environment: namespace:regmedint>

Standard error estimates

regmedint:::calc_myreg_mreg_linear_yreg_logistic_se
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4, sigma_sq, Sigma_beta, Sigma_theta, Sigma_sigma_sq) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4, sigma_sq = sigma_sq)
##     validate_myreg_vcovs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4, sigma_sq = sigma_sq, Sigma_beta = Sigma_beta, 
##         Sigma_theta = Sigma_theta, Sigma_sigma_sq = Sigma_sigma_sq)
##     Sigma <- Matrix::bdiag(Sigma_beta, Sigma_theta, Sigma_sigma_sq)
##     fun_se <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         Gamma_cde <- matrix(c(0, 0, rep(0, length(beta2)), 0, 
##             1, 0, m_cde, rep(0, length(theta4)), 0))
##         Gamma_pnde <- matrix(c(theta3, (theta3 * a0), (theta3 * 
##             c_cond), 0, 1, (theta3 * sigma_sq), (beta0 + beta1 * 
##             a0 + beta2_c + theta2 * sigma_sq + theta3 * sigma_sq * 
##             (a1 + a0)), rep(0, length(theta4)), (theta3 * theta2 + 
##             (1/2) * theta3^2 * (a1 + a0))))
##         Gamma_tnie <- matrix(c(0, (theta2 + theta3 * a1), rep(0, 
##             length(beta2)), 0, 0, beta1, (beta1 * a1), rep(0, 
##             length(theta4)), 0))
##         Gamma_tnde <- matrix(c(theta3, (theta3 * a1), (theta3 * 
##             c_cond), 0, 1, (theta3 * sigma_sq), (beta0 + beta1 * 
##             a1 + beta2_c + theta2 * sigma_sq + theta3 * sigma_sq * 
##             (a1 + a0)), rep(0, length(theta4)), (theta3 * theta2 + 
##             (1/2) * theta3^2 * (a1 + a0))))
##         Gamma_pnie <- matrix(c(0, (theta2 + theta3 * a0), rep(0, 
##             length(beta2)), 0, 0, beta1, (beta1 * a0), rep(0, 
##             length(theta4)), 0))
##         Gamma_te <- Gamma_pnde + Gamma_tnie
##         pnde <- (theta1 + (theta3 * beta0) + (theta3 * beta1 * 
##             a0) + (theta3 * beta2_c) + (theta3 * theta2 * sigma_sq)) * 
##             (a1 - a0) + ((1/2) * theta3^2 * sigma_sq) * (a1^2 - 
##             a0^2)
##         tnie <- ((theta2 * beta1) + theta3 * beta1 * a1) * (a1 - 
##             a0)
##         d_pm <- grad_prop_med_yreg_logistic(pnde = unname(pnde), 
##             tnie = unname(tnie))
##         Gamma_pm <- d_pm[["pnde"]] * Gamma_pnde + d_pm[["tnie"]] * 
##             Gamma_tnie
##         a1_sub_a0 <- abs(a1 - a0)
##         se_cde <- sqrt(as.numeric(t(Gamma_cde) %*% Sigma %*% 
##             Gamma_cde)) * a1_sub_a0
##         se_pnde <- sqrt(as.numeric(t(Gamma_pnde) %*% Sigma %*% 
##             Gamma_pnde)) * a1_sub_a0
##         se_tnie <- sqrt(as.numeric(t(Gamma_tnie) %*% Sigma %*% 
##             Gamma_tnie)) * a1_sub_a0
##         se_tnde <- sqrt(as.numeric(t(Gamma_tnde) %*% Sigma %*% 
##             Gamma_tnde)) * a1_sub_a0
##         se_pnie <- sqrt(as.numeric(t(Gamma_pnie) %*% Sigma %*% 
##             Gamma_pnie)) * a1_sub_a0
##         se_te <- sqrt(as.numeric(t(Gamma_te) %*% Sigma %*% Gamma_te)) * 
##             a1_sub_a0
##         se_pm <- sqrt(as.numeric(t(Gamma_pm) %*% Sigma %*% Gamma_pm)) * 
##             a1_sub_a0
##         c(se_cde = unname(se_cde), se_pnde = unname(se_pnde), 
##             se_tnie = unname(se_tnie), se_tnde = unname(se_tnde), 
##             se_pnie = unname(se_pnie), se_te = unname(se_te), 
##             se_pm = unname(se_pm))
##     }
##     return(fun_se)
## }
## <bytecode: 0x7f9f13d39c80>
## <environment: namespace:regmedint>

mreg logistic yreg linear (V2015 p471 Proposition 2.5)

These functions are only used in the setting where the mediator model is logistic regression and the outcome model is non-linear regression.

Point estimates

regmedint:::calc_myreg_mreg_logistic_yreg_linear_est
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4)
##     fun_est <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         cde <- (theta1 + (theta3 * m_cde)) * (a1 - a0)
##         pnde <- (theta1 * (a1 - a0)) + (theta3 * (a1 - a0)) * 
##             (exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##                 (beta1 * a0) + beta2_c)))
##         tnie <- (theta2 + (theta3 * a1)) * ((exp(beta0 + (beta1 * 
##             a1) + beta2_c)/(1 + exp(beta0 + (beta1 * a1) + beta2_c))) - 
##             (exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##                 (beta1 * a0) + beta2_c))))
##         tnde <- (theta1 * (a1 - a0)) + (theta3 * (a1 - a0)) * 
##             (exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##                 (beta1 * a1) + beta2_c)))
##         pnie <- (theta2 + (theta3 * a0)) * ((exp(beta0 + (beta1 * 
##             a1) + beta2_c)/(1 + exp(beta0 + (beta1 * a1) + beta2_c))) - 
##             (exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##                 (beta1 * a0) + beta2_c))))
##         te <- pnde + tnie
##         pm <- tnie/te
##         c(cde = unname(cde), pnde = unname(pnde), tnie = unname(tnie), 
##             tnde = unname(tnde), pnie = unname(pnie), te = unname(te), 
##             pm = unname(pm))
##     }
##     return(fun_est)
## }
## <bytecode: 0x7f9f141c7b80>
## <environment: namespace:regmedint>

Standard error estimates

regmedint:::calc_myreg_mreg_logistic_yreg_linear_se
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4, Sigma_beta, Sigma_theta) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4)
##     validate_myreg_vcovs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4, Sigma_beta = Sigma_beta, Sigma_theta = Sigma_theta)
##     Sigma <- Matrix::bdiag(Sigma_beta, Sigma_theta)
##     fun_se <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         Gamma_cde <- matrix(c(0, 0, rep(0, length(beta2)), 0, 
##             1, 0, m_cde, rep(0, length(theta4))))
##         pnde_d1 <- theta3 * (exp(beta0 + (beta1 * a0) + beta2_c)/(1 + 
##             exp(beta0 + (beta1 * a0) + beta2_c))^2)
##         pnde_d2 <- a0 * pnde_d1
##         pnde_d3 <- c_cond * pnde_d1
##         pnde_d4 <- 0
##         pnde_d5 <- 1
##         pnde_d6 <- 0
##         pnde_d7 <- exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a0) + beta2_c))
##         pnde_d8 <- rep(0, length(theta4))
##         Gamma_pnde <- matrix(c(pnde_d1, pnde_d2, pnde_d3, pnde_d4, 
##             pnde_d5, pnde_d6, pnde_d7, pnde_d8))
##         tnie_Q <- exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a1) + beta2_c))^2
##         tnie_B <- exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a0) + beta2_c))^2
##         tnie_K <- exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a1) + beta2_c))
##         tnie_D <- exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a0) + beta2_c))
##         tnie_d1 <- (theta2 + (theta3 * a1)) * (tnie_Q - tnie_B)
##         tnie_d2 <- (theta2 + (theta3 * a1)) * ((a1 * tnie_Q) - 
##             (a0 * tnie_B))
##         tnie_d3 <- (theta2 + (theta3 * a1)) * c_cond * (tnie_Q - 
##             tnie_B)
##         tnie_d4 <- 0
##         tnie_d5 <- 0
##         tnie_d6 <- tnie_K - tnie_D
##         tnie_d7 <- a1 * (tnie_K - tnie_D)
##         tnie_d8 <- rep(0, length(theta4))
##         Gamma_tnie <- matrix(c(tnie_d1, tnie_d2, tnie_d3, tnie_d4, 
##             tnie_d5, tnie_d6, tnie_d7, tnie_d8))
##         tnde_d1 <- theta3 * (exp(beta0 + (beta1 * a1) + beta2_c)/(1 + 
##             exp(beta0 + (beta1 * a1) + beta2_c))^2)
##         tnde_d2 <- a1 * tnde_d1
##         tnde_d3 <- c_cond * tnde_d1
##         tnde_d4 <- 0
##         tnde_d5 <- 1
##         tnde_d6 <- 0
##         tnde_d7 <- exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a1) + beta2_c))
##         tnde_d8 <- rep(0, length(theta4))
##         Gamma_tnde <- matrix(c(tnde_d1, tnde_d2, tnde_d3, tnde_d4, 
##             tnde_d5, tnde_d6, tnde_d7, tnde_d8))
##         pnie_Q <- exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a1) + beta2_c))^2
##         pnie_B <- exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a0) + beta2_c))^2
##         pnie_K <- exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a1) + beta2_c))
##         pnie_D <- exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a0) + beta2_c))
##         pnie_d1 <- (theta2 + (theta3 * a0)) * (pnie_Q - pnie_B)
##         pnie_d2 <- (theta2 + (theta3 * a0)) * ((a1 * pnie_Q) - 
##             (a0 * pnie_B))
##         pnie_d3 <- (theta2 + (theta3 * a0)) * c_cond * (pnie_Q - 
##             pnie_B)
##         pnie_d4 <- 0
##         pnie_d5 <- 0
##         pnie_d6 <- pnie_K - pnie_D
##         pnie_d7 <- a0 * (pnie_K - pnie_D)
##         pnie_d8 <- rep(0, length(theta4))
##         Gamma_pnie <- matrix(c(pnie_d1, pnie_d2, pnie_d3, pnie_d4, 
##             pnie_d5, pnie_d6, pnie_d7, pnie_d8))
##         Gamma_te <- ((a1 - a0) * Gamma_pnde) + Gamma_tnie
##         pnde <- (theta1 * (a1 - a0)) + (theta3 * (a1 - a0)) * 
##             (exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##                 (beta1 * a0) + beta2_c)))
##         tnie <- (theta2 + (theta3 * a1)) * ((exp(beta0 + (beta1 * 
##             a1) + beta2_c)/(1 + exp(beta0 + (beta1 * a1) + beta2_c))) - 
##             (exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##                 (beta1 * a0) + beta2_c))))
##         d_pm <- grad_prop_med_yreg_linear(pnde = unname(pnde), 
##             tnie = unname(tnie))
##         Gamma_pm <- (d_pm[["pnde"]] * (a1 - a0) * Gamma_pnde) + 
##             (d_pm[["tnie"]] * Gamma_tnie)
##         a1_sub_a0 <- abs(a1 - a0)
##         se_cde <- sqrt(as.numeric(t(Gamma_cde) %*% Sigma %*% 
##             Gamma_cde)) * a1_sub_a0
##         se_pnde <- sqrt(as.numeric(t(Gamma_pnde) %*% Sigma %*% 
##             Gamma_pnde)) * a1_sub_a0
##         se_tnie <- sqrt(as.numeric(t(Gamma_tnie) %*% Sigma %*% 
##             Gamma_tnie))
##         se_tnde <- sqrt(as.numeric(t(Gamma_tnde) %*% Sigma %*% 
##             Gamma_tnde)) * a1_sub_a0
##         se_pnie <- sqrt(as.numeric(t(Gamma_pnie) %*% Sigma %*% 
##             Gamma_pnie))
##         se_te <- sqrt(as.numeric(t(Gamma_te) %*% Sigma %*% Gamma_te))
##         se_pm <- sqrt(as.numeric(t(Gamma_pm) %*% Sigma %*% Gamma_pm))
##         c(se_cde = unname(se_cde), se_pnde = unname(se_pnde), 
##             se_tnie = unname(se_tnie), se_tnde = unname(se_tnde), 
##             se_pnie = unname(se_pnie), se_te = unname(se_te), 
##             se_pm = unname(se_pm))
##     }
##     return(fun_se)
## }
## <bytecode: 0x7f9f14616448>
## <environment: namespace:regmedint>

mreg logistic yreg non-linear (V2015 p473 Proposition 2.6)

These functions are used in all cases where the mediator model is logistic regression and the outcome model is any one of the non-linear models.

Point estimates

regmedint:::calc_myreg_mreg_logistic_yreg_logistic_est
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4)
##     fun_est <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         cde <- (theta1 + (theta3 * m_cde)) * (a1 - a0)
##         pnde <- (theta1 * (a1 - a0)) + log(1 + exp(theta2 + (theta3 * 
##             a1) + beta0 + (beta1 * a0) + beta2_c)) - log(1 + 
##             exp(theta2 + (theta3 * a0) + beta0 + (beta1 * a0) + 
##                 beta2_c))
##         tnie <- log(1 + exp(beta0 + (beta1 * a0) + beta2_c)) - 
##             log(1 + exp(beta0 + (beta1 * a1) + beta2_c)) + log(1 + 
##             exp(theta2 + (theta3 * a1) + beta0 + (beta1 * a1) + 
##                 beta2_c)) - log(1 + exp(theta2 + (theta3 * a1) + 
##             beta0 + (beta1 * a0) + beta2_c))
##         tnde <- (theta1 * (a1 - a0)) + log(1 + exp(theta2 + (theta3 * 
##             a1) + beta0 + (beta1 * a1) + beta2_c)) - log(1 + 
##             exp(theta2 + (theta3 * a0) + beta0 + (beta1 * a1) + 
##                 beta2_c))
##         pnie <- log(1 + exp(beta0 + (beta1 * a0) + beta2_c)) - 
##             log(1 + exp(beta0 + (beta1 * a1) + beta2_c)) + log(1 + 
##             exp(theta2 + (theta3 * a0) + beta0 + (beta1 * a1) + 
##                 beta2_c)) - log(1 + exp(theta2 + (theta3 * a0) + 
##             beta0 + (beta1 * a0) + beta2_c))
##         te <- pnde + tnie
##         pm <- (exp(pnde) * (exp(tnie) - 1))/(exp(te) - 1)
##         c(cde = unname(cde), pnde = unname(pnde), tnie = unname(tnie), 
##             tnde = unname(tnde), pnie = unname(pnie), te = unname(te), 
##             pm = unname(pm))
##     }
##     return(fun_est)
## }
## <bytecode: 0x7f9f15043700>
## <environment: namespace:regmedint>

Standard error estimates

regmedint:::calc_myreg_mreg_logistic_yreg_logistic_se
## function (beta0, beta1, beta2, theta0, theta1, theta2, theta3, 
##     theta4, Sigma_beta, Sigma_theta) 
## {
##     validate_myreg_coefs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4)
##     validate_myreg_vcovs(beta0 = beta0, beta1 = beta1, beta2 = beta2, 
##         theta0 = theta0, theta1 = theta1, theta2 = theta2, theta3 = theta3, 
##         theta4 = theta4, Sigma_beta = Sigma_beta, Sigma_theta = Sigma_theta)
##     Sigma <- Matrix::bdiag(Sigma_beta, Sigma_theta)
##     fun_se <- function(a0, a1, m_cde, c_cond) {
##         if (is.null(beta2)) {
##             assertthat::assert_that(is.null(c_cond))
##             beta2_c <- 0
##         }
##         else {
##             assertthat::assert_that(!is.null(c_cond))
##             assertthat::assert_that(length(c_cond) == length(beta2))
##             beta2_c <- sum(t(matrix(beta2)) %*% matrix(c_cond))
##         }
##         Gamma_cde <- matrix(c(0, 0, rep(0, length(beta2)), 0, 
##             (a1 - a0), 0, (a1 - a0) * m_cde, rep(0, length(theta4))))
##         pnde_Q <- exp(theta2 + (theta3 * a1) + beta0 + (beta1 * 
##             a0) + beta2_c)/(1 + exp(theta2 + (theta3 * a1) + 
##             beta0 + (beta1 * a0) + beta2_c))
##         pnde_B <- exp(theta2 + (theta3 * a0) + beta0 + (beta1 * 
##             a0) + beta2_c)/(1 + exp(theta2 + (theta3 * a0) + 
##             beta0 + (beta1 * a0) + beta2_c))
##         pnde_d1 <- pnde_Q - pnde_B
##         pnde_d2 <- a0 * pnde_d1
##         pnde_d3 <- c_cond * pnde_d1
##         pnde_d4 <- 0
##         pnde_d5 <- (a1 - a0)
##         pnde_d6 <- pnde_Q - pnde_B
##         pnde_d7 <- (a1 * pnde_Q) - (a0 * pnde_B)
##         pnde_d8 <- rep(0, length(theta4))
##         Gamma_pnde <- matrix(c(pnde_d1, pnde_d2, pnde_d3, pnde_d4, 
##             pnde_d5, pnde_d6, pnde_d7, pnde_d8))
##         tnie_Q <- exp(theta2 + (theta3 * a1) + beta0 + (beta1 * 
##             a1) + beta2_c)/(1 + exp(theta2 + (theta3 * a1) + 
##             beta0 + (beta1 * a1) + beta2_c))
##         tnie_B <- exp(theta2 + (theta3 * a1) + beta0 + (beta1 * 
##             a0) + beta2_c)/(1 + exp(theta2 + (theta3 * a1) + 
##             beta0 + (beta1 * a0) + beta2_c))
##         tnie_K <- exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a1) + beta2_c))
##         tnie_D <- exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a0) + beta2_c))
##         tnie_d1 <- (tnie_D + tnie_Q) - (tnie_K + tnie_B)
##         tnie_d2 <- (a0 * (tnie_D - tnie_B)) + (a1 * (tnie_Q - 
##             tnie_K))
##         tnie_d3 <- c_cond * ((tnie_D + tnie_Q) - (tnie_K + tnie_B))
##         tnie_d4 <- 0
##         tnie_d5 <- 0
##         tnie_d6 <- tnie_Q - tnie_B
##         tnie_d7 <- a1 * (tnie_Q - tnie_B)
##         tnie_d8 <- rep(0, length(theta4))
##         Gamma_tnie <- matrix(c(tnie_d1, tnie_d2, tnie_d3, tnie_d4, 
##             tnie_d5, tnie_d6, tnie_d7, tnie_d8))
##         tnde_Q <- exp(theta2 + (theta3 * a1) + beta0 + (beta1 * 
##             a1) + beta2_c)/(1 + exp(theta2 + (theta3 * a1) + 
##             beta0 + (beta1 * a1) + beta2_c))
##         tnde_B <- exp(theta2 + (theta3 * a0) + beta0 + (beta1 * 
##             a1) + beta2_c)/(1 + exp(theta2 + (theta3 * a0) + 
##             beta0 + (beta1 * a1) + beta2_c))
##         tnde_d1 <- tnde_Q - tnde_B
##         tnde_d2 <- a1 * tnde_d1
##         tnde_d3 <- c_cond * tnde_d1
##         tnde_d4 <- 0
##         tnde_d5 <- (a1 - a0)
##         tnde_d6 <- tnde_Q - tnde_B
##         tnde_d7 <- (a1 * tnde_Q) - (a0 * tnde_B)
##         tnde_d8 <- rep(0, length(theta4))
##         Gamma_tnde <- matrix(c(tnde_d1, tnde_d2, tnde_d3, tnde_d4, 
##             tnde_d5, tnde_d6, tnde_d7, tnde_d8))
##         pnie_Q <- exp(theta2 + (theta3 * a0) + beta0 + (beta1 * 
##             a1) + beta2_c)/(1 + exp(theta2 + (theta3 * a0) + 
##             beta0 + (beta1 * a1) + beta2_c))
##         pnie_B <- exp(theta2 + (theta3 * a0) + beta0 + (beta1 * 
##             a0) + beta2_c)/(1 + exp(theta2 + (theta3 * a0) + 
##             beta0 + (beta1 * a0) + beta2_c))
##         pnie_K <- exp(beta0 + (beta1 * a1) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a1) + beta2_c))
##         pnie_D <- exp(beta0 + (beta1 * a0) + beta2_c)/(1 + exp(beta0 + 
##             (beta1 * a0) + beta2_c))
##         pnie_d1 <- (pnie_D + pnie_Q) - (pnie_K + pnie_B)
##         pnie_d2 <- (a0 * (pnie_D - pnie_B)) + (a1 * (pnie_Q - 
##             pnie_K))
##         pnie_d3 <- c_cond * ((pnie_D + pnie_Q) - (pnie_K + pnie_B))
##         pnie_d4 <- 0
##         pnie_d5 <- 0
##         pnie_d6 <- pnie_Q - pnie_B
##         pnie_d7 <- a0 * (pnie_Q - pnie_B)
##         pnie_d8 <- rep(0, length(theta4))
##         Gamma_pnie <- matrix(c(pnie_d1, pnie_d2, pnie_d3, pnie_d4, 
##             pnie_d5, pnie_d6, pnie_d7, pnie_d8))
##         Gamma_te <- Gamma_pnde + Gamma_tnie
##         pnde <- (exp(theta1 * a1) * (1 + exp(theta2 + (theta3 * 
##             a1) + beta0 + (beta1 * a0) + beta2_c)))/(exp(theta1 * 
##             a0) * (1 + exp(theta2 + (theta3 * a0) + beta0 + (beta1 * 
##             a0) + beta2_c)))
##         tnie <- ((1 + exp(beta0 + (beta1 * a0) + beta2_c)) * 
##             (1 + exp(theta2 + (theta3 * a1) + beta0 + (beta1 * 
##                 a1) + beta2_c)))/((1 + exp(beta0 + (beta1 * a1) + 
##             beta2_c)) * (1 + exp(theta2 + (theta3 * a1) + beta0 + 
##             (beta1 * a0) + beta2_c)))
##         d_pm <- grad_prop_med_yreg_logistic(pnde = unname(pnde), 
##             tnie = unname(tnie))
##         Gamma_pm <- (d_pm[["pnde"]] * Gamma_pnde) + (d_pm[["tnie"]] * 
##             Gamma_tnie)
##         se_cde <- sqrt(as.numeric(t(Gamma_cde) %*% Sigma %*% 
##             Gamma_cde))
##         se_pnde <- sqrt(as.numeric(t(Gamma_pnde) %*% Sigma %*% 
##             Gamma_pnde))
##         se_tnie <- sqrt(as.numeric(t(Gamma_tnie) %*% Sigma %*% 
##             Gamma_tnie))
##         se_tnde <- sqrt(as.numeric(t(Gamma_tnde) %*% Sigma %*% 
##             Gamma_tnde))
##         se_pnie <- sqrt(as.numeric(t(Gamma_pnie) %*% Sigma %*% 
##             Gamma_pnie))
##         se_te <- sqrt(as.numeric(t(Gamma_te) %*% Sigma %*% Gamma_te))
##         se_pm <- sqrt(as.numeric(t(Gamma_pm) %*% Sigma %*% Gamma_pm))
##         c(se_cde = unname(se_cde), se_pnde = unname(se_pnde), 
##             se_tnie = unname(se_tnie), se_tnde = unname(se_tnde), 
##             se_pnie = unname(se_pnie), se_te = unname(se_te), 
##             se_pm = unname(se_pm))
##     }
##     return(fun_se)
## }
## <bytecode: 0x7f9f150b4948>
## <environment: namespace:regmedint>

Bibliography