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.
These functions are only used in the setting where both the mediator model and the outcome model are linear regression.
## 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: 0x7fcd49b32438>
## <environment: namespace:regmedint>## 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: 0x7fcd49300ac8>
## <environment: namespace:regmedint>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.
## 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: 0x7fcd6a2f30d0>
## <environment: namespace:regmedint>## 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: 0x7fcd2cec6c10>
## <environment: namespace:regmedint>These functions are only used in the setting where the mediator model is logistic regression and the outcome model is non-linear regression.
## 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: 0x7fcd2f14d258>
## <environment: namespace:regmedint>## 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: 0x7fcd2f310180>
## <environment: namespace:regmedint>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.
## 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: 0x7fcd3ac76a20>
## <environment: namespace:regmedint>## 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: 0x7fcd3acd7040>
## <environment: namespace:regmedint>