R/reg1_no_int.R

In causal.decomp: Causal Decomposition Analysis

reg1_no_int <- function(b.m = NULL, b.y = NULL, B, data, covariates, group, group.ref, intermediates, risk.factor, outcome, moderators, mods_formula = NULL, cluster){

  if (is.null(mods_formula)) {
    if (is.null(moderators) || length(moderators) == 0)
      stop("Provide either `mods_formula` or a non-empty `moderators` character vector.")
    mods_formula <- as.formula(paste0(" ~ ", paste(moderators, collapse = " + ")))
  }
  
  est1 <- function(data){

    data_R0 <- subset(data, data[[group]] == group.ref)
    data_R1 <- subset(data, data[[group]] != group.ref)

    # weighting
    moPropen1 <- buildModelObj(model = as.formula(paste0(" ~ ", paste(c(covariates,intermediates, "U"), collapse = " + "))),
                               solver.method = "glm",
                               solver.args = list("family" = "binomial"),
                               predict.method = 'predict.glm',
                               predict.args = list(type = "response"))

    moClass1 <- buildModelObj(model = mods_formula,
                              solver.method = rpart,
                              solver.args = list(method = "class"),
                              predict.args = list(type = 'class'))

    fitFS1 <- optimalClass(moPropen = moPropen1,
                           moClass = moClass1,
                           data = data.frame(data_R1), response = data_R1[[outcome]],
                           txName = risk.factor, verbose = FALSE)
    fitFS0 <- optimalClass(moPropen = moPropen1,
                           moClass = moClass1,
                           data = data.frame(data_R0), response = data_R0[[outcome]],
                           txName = risk.factor, verbose = FALSE)


    data$opt.M <- ifelse(data[[group]] == group.ref, (optTx(fitFS0)$optimalTx),(optTx(fitFS1)$optimalTx)) 

    ta <- table(data$opt.M)
    p.opt <- ta[2] / (ta[2] + ta[1]) * 100
    # Construct weight for mediator
    DATA1 <- data
    DATA1 <- DATA1 %>%
      mutate(Ind = (DATA1[[risk.factor]] == opt.M))

    fit.m <- glm(as.formula(paste0(risk.factor, " ~ ", paste(c(group,covariates,intermediates), collapse = " + "))), offset = b.m * DATA1$U,
                 family = binomial(logit), data = DATA1)
    coef.m <- c(fit.m$coef, U = b.m)
    pre.m <- 1 / (1 + exp(-(cbind(model.matrix(fit.m), U = DATA1$U) %*% coef.m)))
    p.med <- ifelse(DATA1[[risk.factor]] == 0, 1 - pre.m, pre.m)

    # Weights
    DATA1$w <- 1 / p.med

    # Calculate zeta_ICDE
    fit <- lm(paste0(outcome, " ~ ", paste(c(covariates,group), collapse = " + ")), data=DATA1, weights = Ind * w)
    estimated_zeta_ICDE <- coef(fit)[3] # coefficient of race

    # IIE: Method 2
    fit.lm1 <- lm(paste0(outcome, " ~ ", paste(c(covariates), collapse = " + ")), data=data_R1)
    W_Y1 <- coef(fit.lm1)[1]

    fit.lm0 <- lm(paste0(outcome, " ~ ", paste(c(covariates), collapse = " + ")), data=data_R0)
    W_Y0 <- coef(fit.lm0)[1]

    fit.lm <- lm(paste0(outcome, " ~ ", paste(c(covariates,group,"Ind", paste0(group, ":Ind")), collapse = " + ")), data = DATA1, weights = w)
    fit.Ind <- glm(paste0("Ind", " ~ ", paste(c(group), collapse = " + ")), data = DATA1, family = binomial(logit))

    alpha1 <- plogis(sum(coef(fit.Ind)[1:2])) - plogis(coef(fit.Ind)[1])

    coef_names <- names(coef(fit.lm))
    ind_coef <- coef_names[grepl(paste0("^", "Ind"), coef_names)]
    interaction_coef <- coef_names[grepl(paste0(":", "Ind"), coef_names)]

    reg_delta_IIE <- alpha1 * sum(coef(fit.lm)[c(ind_coef, interaction_coef)])
    reg_zeta_IIE <- (W_Y1 - W_Y0) - reg_delta_IIE

    return(list(fitFS1 = fitFS1, fitFS0 = fitFS0, p.opt = p.opt,
                estimated_zeta_ICDE = estimated_zeta_ICDE,
                reg_zeta_IIE = reg_zeta_IIE,
                reg_delta_IIE = reg_delta_IIE))
  }

  est.ori <- est1(data)

  reg_delta_IIE <- est.ori$reg_delta_IIE
  reg_zeta_IIE <- est.ori$reg_zeta_IIE
  estimated_zeta_ICDE <- est.ori$estimated_zeta_ICDE
  p.opt <- est.ori$p.opt
  fitFS1 <- est.ori$fitFS1
  fitFS0 <- est.ori$fitFS0

  regzeta_IIE.boot <- regdelta_IIE.boot <- zeta_ICDE.boot <- rep(NA, B)

  for(i in 1:B){
    if (is.null(cluster)) {
      data.boot <- data[sample(nrow(data), size = nrow(data), replace = TRUE), ]
    } else {
      clusters <- unique(data[,cluster])
      m <- length(clusters)
      units <- sample(clusters, size = length(clusters), replace = TRUE)
      df.bs <- sapply(units, function(x) which(data[,cluster] == x))
      sb <- unlist(df.bs)
      data.boot <- data[sb, ]
    }

    est.boot <- est1(data.boot)

    zeta_ICDE.boot[i] <- est.boot$estimated_zeta_ICDE
    regzeta_IIE.boot[i] <- est.boot$reg_zeta_IIE
    regdelta_IIE.boot[i] <- est.boot$reg_delta_IIE
  }

  SE_zeta_ICDE <- sd(zeta_ICDE.boot)
  SE_regdelta_IIE <- sd(regdelta_IIE.boot)
  SE_regzeta_IIE <- sd(regzeta_IIE.boot)

  results <- list(
    p.opt = p.opt, fitFS0 = fitFS0, fitFS1 = fitFS1,
    estimated_zeta_ICDE = as.numeric(estimated_zeta_ICDE), SE_zeta_ICDE = SE_zeta_ICDE,
    reg_delta_IIE = as.numeric(reg_delta_IIE), SE_regdelta_IIE = SE_regdelta_IIE,
    reg_zeta_IIE = as.numeric(reg_zeta_IIE), SE_regzeta_IIE = SE_regzeta_IIE
  )

  #print("results from reg1_no_int.R"); print(results)

  return(results)
}