R/combine.R

In causal.decomp: Causal Decomposition Analysis

# A function used in 'mmi', 'smi', and 'pocr' for comparing all pairs of R
combine <- function(fit.r, fit.x, fit.m, fit.y, func, weights){
  
  # Model frame Y model
  y.data <- model.frame(fit.y)
  # Number of observations
  n.y <- nrow(y.data)
  
  if(func == "mmi" | func == "smi"){
    ## Calculate weighted outcome from R model
    if(!is.null(fit.r) && isCBPS.r){
      y.data$w <- fit.r$weights
      prop.group <- table(y.data[, group])/nrow(y.data)
    } else if(!is.null(fit.r) && isSumStat.r){
      y.data$w <- fit.r$ps.weights[, 3]
      prop.group <- table(y.data[, group])/nrow(y.data)
    } else {
      y.data$w <- rep(1, n.y)
      prop.group <- rep(1, length(levels(y.data[, group])))
    }
    
    wy <- rep(NA, length(levels(y.data[, group])))
    if(isGlm.y){ ## fixed 8/27/2021 for conditional = T
      y.data$pred.prob <- predict(fit.y, newdata = y.data, type = "response")
      if(conditional){
        wy.f <- as.formula(paste("pred.prob", paste(covariates, collapse = "+"), sep = "~"))
      } else if (!conditional){
        wy.f <- as.formula(paste("pred.prob", "~ 1"))
      }
    } else if (!isGlm.y){
      if(conditional){
        wy.f <- as.formula(paste(outcome, paste(covariates, collapse = "+"), sep = "~"))
      } else if (!conditional){
        wy.f <- as.formula(paste(outcome, "~ 1"))
      }
    }
    for(i in 1:length(levels(y.data[, group]))){
      wy[i] <- lm(wy.f, weights = w * prop.group[i],
                  data = y.data[y.data[, group] == levels(y.data[, group])[i], ])$coef[1]
    }
  }
  
  if(func == "mmi" | func == "smi"){
    environment(select_group) <- environment()
    # Do all combinations, ref.lev.group vs. sel.lev.group
    res <- NULL
    for(l in 2:length(levels(y.data[, group]))){
      res <- c(res, select_group(fit.m = fit.m, fit.y = fit.y, sel.lev.group = l,
                                     ref.lev.group = 1,
                                     prop.group = prop.group, wy = wy, func = func, weights = weights))
    }
  } else if (func == "pocr"){
    environment(select_group_pocr) <- environment()
    # Do all combinations, ref.lev.group vs. sel.lev.group
    res <- NULL
    for(l in 2:length(levels(y.data[, group]))){
      res <- c(res, select_group_pocr(fit.x = fit.x, fit.m = fit.m, fit.y = fit.y, sel.lev.group = l,
                                      ref.lev.group = 1))
    }
  }
  
  return(res) 
}