R/select_group_pocr.R

In causal.decomp: Causal Decomposition Analysis

# A function used in 'mmi' and 'smi' for comparing a pair of sel.lev.R vs ref.lev.group
select_group_pocr <- function(fit.x, fit.m, fit.y, sel.lev.group, ref.lev.group = 1, func){
  
  #:::::::::::::::::::::::::::#
  # Predict mediator::::::::::#
  #:::::::::::::::::::::::::::#
  y.data.new <- y.data
  y.data.new[, group] <- levels(y.data[, group])[sel.lev.group]
  
  # version 0.1.0 -> change!!!!
  #y.data.updated <- y.data
  #y.data.updated[, group] <- relevel(y.data.updated[, group], ref = sel.lev.group)
  #fit.y.updated <- update(fit.y, data = y.data.updated) ### correct?
  #meds <- all.vars(formula(fit.m))[[1]]
  #if(isGlm.m){
  #  meds <- which(meds == substring(colnames(vcov(fit.y.updated)), 1, nchar(meds)))[1]
  #}
  #var_gamma <- vcov(fit.y.updated)[meds, meds]
  #se_gamma <- sqrt(var_gamma)
  
  # Compute predicted values of mediator
  
  ### Case 1: GLM Mediator
  if(isGlm.m){
    
    if(FamilyM == "poisson" | FamilyM == "Gamma" |
       FamilyM == "gaussian" | FamilyM == "inverse.gaussian"){
      
      mod.m.f <- as.formula(paste(M, "~ ", paste(covariates, collapse = "+")))
      data.ref <- subset(y.data, y.data[, group] == levels(y.data[, group])[ref.lev.group])
      mod.m.ref <- glm(mod.m.f, data = data.ref, family = M.fun)
      data.sel <- subset(y.data, y.data[, group] == levels(y.data[, group])[sel.lev.group])
      mod.m.sel <- glm(mod.m.f, data = data.sel, family = M.fun)
      
      # alpha1
      if(LinkM == "identity") {
        alpha1 <- coef(mod.m.ref)[1] - coef(mod.m.sel)[1]
      } else if(LinkM == "log") {
        alpha1 <- exp(coef(mod.m.ref)[1]) - exp(coef(mod.m.sel)[1])
      } else if(LinkM == "sqrt") {
        alpha1 <- (coef(mod.m.ref)[1])^2 - (coef(mod.m.sel)[1])^2
      } else if(LinkM == "inverse") {
        alpha1 <- (coef(mod.m.ref)[1])^(-1) - (coef(mod.m.sel)[1])^(-1)
      } else if(LinkM == "1/mu^2") {
        alpha1 <- (coef(mod.m.ref)[1])^(-1/2) - (coef(mod.m.sel)[1])^(-1/2)
      } 
      
      # beta3, beta4
      name.beta3 <- paste(M, levels(y.data[, M])[2], sep = "")
      name.beta41 <- paste(paste(group, levels(y.data[, group])[sel.lev.group], sep = ""),
                           paste(M, levels(y.data[, M])[2], sep = ""), sep = ":")
      name.beta42 <- paste(paste(M, levels(y.data[, M])[2], sep = ""),
                           paste(group, levels(y.data[, group])[sel.lev.group], sep = ""), sep = ":")
      
      ind.beta3 <- which(names(coef(fit.y)) == name.beta3)
      ind.beta4 <- which(names(coef(fit.y)) == name.beta41 | names(coef(fit.y)) == name.beta42)
      
      beta3 <- coef(fit.y)[ind.beta3]
      beta4 <- coef(fit.y)[ind.beta4]
      
      # phi1
      mod.y.f <- as.formula(paste(outcome, "~ ", paste(c(group, covariates), collapse = "+")))
      mod.y <- lm(mod.y.f, data = y.data)
      name.phi1 <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
      ind.phi1 <- which(names(coef(mod.y)) == name.phi1)
      coef.phi1 <- coef(mod.y)[ind.phi1]
      
      # disparity reduction
      disp.red <- alpha1 * (beta3 + beta4)
      # disparity remaining for a continuous mediator
      disp.rem <- phi1 - alpha1 * (beta3 + beta4)
      
    } else if (FamilyM == "binomial"){
      
      mod.m.f <- as.formula(paste(M, "~ ", paste(covariates, collapse = "+")))
      data.ref <- subset(y.data, y.data[, group] == levels(y.data[, group])[ref.lev.group])
      mod.m.ref <- glm(mod.m.f, data = data.ref, family = M.fun)
      data.sel <- subset(y.data, y.data[, group] == levels(y.data[, group])[sel.lev.group])
      mod.m.sel <- glm(mod.m.f, data = data.sel, family = M.fun)
      
      # alpha1, alpha0
      if(LinkM == "logit"){
        alpha1 <- -(plogis(coef(mod.m.ref)[1]) - plogis(coef(mod.m.sel)[1]))
        alpha0 <- plogis(coef(mod.m.sel)[1])
      } else if(LinkM == "probit") {
        alpha1 <- pnorm(coef(mod.m.ref)[1]) - pnorm(coef(mod.m.sel)[1])
        alpha0 <- pnorm(coef(mod.m.sel)[1])
      } else if(LinkM == "cloglog") {
        alpha1 <- exp(- exp(coef(mod.m.sel)[1])) - exp(- exp(coef(mod.m.ref)[1]))
        alpha0 <- 1 - exp(- exp(coef(mod.m.sel)[1]))
      }
      
      # beta1, beta2, beta3, beta4
      name.beta1 <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
      name.beta2 <- X
      #name.X1 <- X1
      #name.X2 <- X2
      name.beta3 <- paste(M, levels(y.data[, M])[2], sep = "")
      name.beta41 <- paste(paste(group, levels(y.data[, group])[sel.lev.group], sep = ""),
                           paste(M, levels(y.data[, M])[2], sep = ""), sep = ":")
      name.beta42 <- paste(paste(M, levels(y.data[, M])[2], sep = ""),
                           paste(group, levels(y.data[, group])[sel.lev.group], sep = ""), sep = ":")
      
      ind.beta1 <- which(names(coef(fit.y)) == name.beta1)
      ind.beta2 <- which(names(coef(fit.y)) == name.beta2)
      #ind.X1 <- which(names(coef(fit.y)) == name.X1)
      #ind.X2 <- which(names(coef(fit.y)) == name.X2)
      ind.beta3 <- which(names(coef(fit.y)) == name.beta3)
      ind.beta4 <- which(names(coef(fit.y)) == name.beta41 | names(coef(fit.y)) == name.beta42)
      
      beta1 <- coef(fit.y)[ind.beta1]
      beta2 <- coef(fit.y)[ind.beta2]
      #coef.X1 <- coef(fit.y)[ind.X1]
      #coef.X2 <- coef(fit.y)[ind.X2]
      beta3 <- coef(fit.y)[ind.beta3]
      beta4 <- coef(fit.y)[ind.beta4]
      
      # gamma1
      name.gamma1 <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
      ind.gamma1 <- which(names(coef(fit.x)) == name.gamma1)
      gamma1 <- coef(fit.x)[ind.gamma1]
      #ind.gamma1.X1 <- which(names(coef(fit.x1)) == name.gamma1)
      #gamma1.X1 <- coef(fit.x1)[ind.gamma1.X1]
      #ind.gamma1.X2 <- which(names(coef(fit.X2)) == name.gamma1)
      #gamma1.X2 <- coef(fit.x2)[ind.gamma1.X2]
      
      # disparity reduction
      disp.red <- alpha1 * (beta3 + beta4)
      # disparity remaining for a binary mediator
      disp.rem <- beta1 + beta2 * gamma1 + beta4 * alpha0
      #disp.rem <- beta1 + coef.X1 * gamma1.X1 + coef.X2 * gamma1.X2 + beta4 * alpha0
      
    } else {
      stop("unsupported glm family")
    }
    
    ### Case 2: LM Mediator
  } else if(isLm.m & !isGlm.m){
    
    # alpha1
    name.R.in.M <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
    ind.R.in.M <- which(names(coef(fit.m)) == name.R.in.M)
    coef.R.in.M <- coef(fit.m)[ind.R.in.M] 
    
    # beta3
    ind.M.in.Y <- which(names(coef(fit.y)) == M)
    coef.M.in.Y <- coef(fit.y)[ind.M.in.Y]
    
    # beta4
    name.RM.in.Y1 <- paste(paste(group, levels(y.data[, group])[sel.lev.group], sep = ""), M, sep = ":")
    name.RM.in.Y2 <- paste(M, paste(group, levels(y.data[, group])[sel.lev.group], sep = ""), sep = ":")
    ind.RM.in.Y <- which(names(coef(fit.y)) == name.RM.in.Y1 | names(coef(fit.y)) == name.RM.in.Y2)
    coef.RM.in.Y <- coef(fit.y)[ind.RM.in.Y]
    
    # phi1
    mod.y.f <- as.formula(paste(outcome, "~ ", paste(c(group, covariates), collapse = "+")))
    mod.y <- lm(mod.y.f, data = y.data)
    name.R.in.y <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
    ind.R.in.y <- which(names(coef(mod.y)) == name.R.in.y)
    coef.R.in.y <- coef(mod.y)[ind.R.in.y]
    
    # disparity reduction = alpha1 * (beta3 + beta4)
    disp.red <- coef.R.in.M * (coef.M.in.Y + coef.RM.in.Y)
    # disparity remaining = phi1 - alpha1 * (beta3 + beta4)
    disp.rem <- coef.R.in.y - coef.R.in.M * (coef.M.in.Y + coef.RM.in.Y)
    
    ### Case 3: Nominal or Ordinal Mediator
  } else if (isNominal.m | isOrdinal.m) {
    
    mod.m.f <- as.formula(paste(M, "~ ", paste(covariates, collapse = "+")))
    m <- length(levels(y.data[, M]))
    prob.m.ref <- prob.m.sel <- rep(NA, m)
    
    if(isNominal.m){
      
      data.ref <- subset(y.data, y.data[, group] == levels(y.data[, group])[ref.lev.group])
      mod.m.ref <- multinom(mod.m.f, data = data.ref, trace = FALSE)
      data.sel <- subset(y.data, y.data[, group] == levels(y.data[, group])[sel.lev.group])
      mod.m.sel <- multinom(mod.m.f, data = data.sel, trace = FALSE)
      
      prob.m.ref[1] <- 1
      prob.m.sel[1] <- 1
      for(mm in 2:m){
        prob.m.ref[mm] <- exp(coef(mod.m.ref)[mm - 1, 1])
        prob.m.sel[mm] <- exp(coef(mod.m.sel)[mm - 1, 1])
      }
      prob.m.ref <- prob.m.ref/(1 + sum(prob.m.ref[2:m]))
      prob.m.sel <- prob.m.sel/(1 + sum(prob.m.sel[2:m]))
      
    } else if (isOrdinal.m){
      
      data.ref <- subset(y.data, y.data[, group] == levels(y.data[, group])[ref.lev.group])
      mod.m.ref <- polr(mod.m.f, data = data.ref, method = fit.m$method)
      data.sel <- subset(y.data, y.data[, group] == levels(y.data[, group])[sel.lev.group])
      mod.m.sel <- polr(mod.m.f, data = data.sel, method = fit.m$method)
      
      prob.m.ref[1] <- pfun(mod.m.ref$zeta[1])
      prob.m.sel[1] <- pfun(mod.m.sel$zeta[1])
      for(mm in 2:(m - 1)){
        prob.m.ref[mm] <- pfun(mod.m.ref$zeta[mm]) - pfun(mod.m.ref$zeta[mm - 1])
        prob.m.sel[mm] <- pfun(mod.m.sel$zeta[mm]) - pfun(mod.m.sel$zeta[mm - 1])
      }
      prob.m.ref[m] <- 1 - pfun(mod.m.ref$zeta[m - 1])
      prob.m.sel[m] <- 1 - pfun(mod.m.sel$zeta[m - 1])
      
    }
    #prob.m.ref <- colMeans(predict(mod.m.ref, type = "probs"))
    #prob.m.sel <- colMeans(predict(mod.m.sel, type = "probs"))
    
    # gamma1
    name.R <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
    ind.R <- which(names(coef(fit.y)) == name.R)
    ind.X <- which(names(coef(fit.y)) == X)
    ind.R.in.X <- which(names(coef(fit.x)) == name.R)
    
    disp.red <- disp.rem <- rep(NA, m)
    for(mm in 2:m){
      name.M <- paste(M, levels(y.data[, M])[mm], sep = "")
      name.RM1 <- paste(name.R, name.M, sep = ":")
      name.RM2 <- paste(name.M, name.R, sep = ":")
      ind.M <- which(names(coef(fit.y)) == name.M)
      ind.RM <- which(names(coef(fit.y)) == name.RM1 | names(coef(fit.y)) == name.RM2)
      
      # beta1, beta2, beta4
      name.beta1 <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
      name.beta2 <- X
      name.beta41 <- paste(paste(group, levels(y.data[, group])[sel.lev.group], sep = ""),
                           paste(M, levels(y.data[, M])[2], sep = ""), sep = ":")
      name.beta42 <- paste(paste(M, levels(y.data[, M])[2], sep = ""),
                           paste(group, levels(y.data[, group])[sel.lev.group], sep = ""), sep = ":")
      
      ind.beta1 <- which(names(coef(fit.y)) == name.beta1)
      ind.beta2 <- which(names(coef(fit.y)) == name.beta2)
      ind.beta4 <- which(names(coef(fit.y)) == name.beta41 | names(coef(fit.y)) == name.beta42)
      
      beta1 <- coef(fit.y)[ind.beta1]
      beta2 <- coef(fit.y)[ind.beta2]
      beta4 <- coef(fit.y)[ind.beta4]
      
      disp.red[mm] <- (coef(fit.y)[ind.M] + coef(fit.y)[ind.RM]) * (prob.m.sel[mm] - prob.m.ref[mm])
      disp.rem[mm] <- coef(fit.y)[ind.RM] * prob.m.ref[mm]
    }
    
    disp.red <- sum(disp.red, na.rm = TRUE)
    disp.rem <- coef(fit.y)[ind.R] + coef(fit.y)[ind.X] * coef(fit.x)[ind.R.in.X] + sum(disp.rem, na.rm = TRUE)
    
  } else {
    stop("mediator model is not yet implemented")
  }
  
  #:::::::::::::::#
  #:Results:::::::#
  #:::::::::::::::#
  # Initial disparity, Disparity remaining, and Disparity reduction in order
  out <- rep(NA, 5)
  out[2] <- disp.rem
  out[3] <- disp.red
  out[1] <- out[2] + out[3]
  
  # version 0.1.0
  #name.R.in.M <- paste(group, levels(y.data[, group])[sel.lev.group], sep = "")
  #ind.R.in.M <- which(names(coef(fit.m)) == name.R.in.M)
  out[4] <- 0 #as.numeric(coef(fit.m)[ind.R.in.M])
  out[5] <- 0 #se_gamma
  
  return(out)
}