R/apce.R

In aihuman: Experimental Evaluation of Algorithm-Assisted Human Decision-Making

#' Calculate APCE
#'
#' Calculate average principal causal effects (APCE) with ordinal decision. See Section 3.4 for more details.
#'
#' @param data A \code{data.frame} or \code{matrix} of which columns consists of pre-treatment covariates, a binary treatment (Z), an ordinal decision (D), and an outcome variable (Y). The column names of the latter three should be specified as "Z", "D", and "Y" respectively.
#' @param mcmc.re A \code{mcmc} object generated by \code{AiEvalmcmc()} function.
#' @param subgroup A list of numeric vectors for the index of each of the five subgroups.
#' @param name.group A list of character vectors for the label of five subgroups.
#' @param rho A sensitivity parameter. The default is  \code{0} which implies the unconfoundedness assumption (Assumption 4).
#' @param burnin A proportion of burnin for the Markov chain. The default is \code{0}.
#' @param out.length An integer to specify the progress on the screen. Every \code{out.length}-th iteration is printed on the screen. The default is \code{500}.
#' @param c0 The cost of an outcome. See Section 3.7 for more details. The default is \code{0}.
#' @param c1 The cost of an unnecessarily harsh decision. See Section 3.7 for more details. The default is \code{0}.
#' @param ZX The data matrix for interaction terms. The default is the interaction between Z and all of the pre-treatment covariates (X).
#' @param save.individual.optimal.decision A logical argument specified to save individual optimal decision rules. The default is \code{FALSE}.
#' @param parallel A logical argument specifying whether parallel computing is conducted. Do not change this argument manually.
#' @param optimal.decision.only A logical argument specified to compute only the optimal decision rule. The default is \code{FALSE}.
#' @param dmf A numeric vector of binary DMF recommendations. If \code{null}, use judge's decisions (0 if the decision is 0 and 1 o.w; e.g., signature or cash bond).
#' @param fair.dmf.only A logical argument specified to compute only the fairness of given DMF recommendations. The default is \code{FALSE}. Not used in the analysis for the JRSSA paper.
#'
#' @return An object of class \code{list} with the following elements:
#' \item{P.D1.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for quantity P(D(1)=d| R=r), dimension 1 is each posterior sample; dimension 2 is subgroup, dimension 3 is (k+1) values of D from 0 to k, dimension 4 is (k+2) values of R from 0 to k+1.}
#' \item{P.D0.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for quantity P(D(0)=d| R=r).}
#' \item{APCE.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for quantity P(D(1)=d| R=r)-P(D(0)=d| R=r).}
#' \item{P.R.mcmc}{An array with dimension n.mcmc by 5 by (k+2) for quantity P(R=r) for r from 0 to (k+1).}
#' \item{Optimal.Z.mcmc}{An array with dimension n.mcmc by 5 for the proportion of the cases where treatment (PSA provided) is optimal.}
#' \item{Optimal.D.mcmc}{An array with dimension n.mcmc by 5 by (k+1) for the proportion of optimal decision rule (average over observations). If \code{save.individual.optimal.decision = TRUE}, the dimension would be n by (k+1) (average over mcmc samples).}
#' \item{P.DMF.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for the proportion of binary DMF recommendations. Not used in the analysis for the JRSSA paper.}
#' \item{Utility.g_d.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n for the individual utility of judge's decisions.}
#' \item{Utility.g_dmf.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n for the individual utility of DMF recommendation.}
#' \item{Utility.diff.control.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n.mcmc for estimated difference in utility between judge's decisions and DMF recommendation among control group.}
#' \item{Utility.diff.treated.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n.mcmc for estimated difference in utility between judge's decisions and DMF recommendation among treated group.}
#'
#' @examples
#' data(synth)
#' sample_mcmc <- AiEvalmcmc(data = synth, n.mcmc = 2)
#' subgroup_synth <- list(
#'   1:nrow(synth), which(synth$Sex == 0), which(synth$Sex == 1),
#'   which(synth$Sex == 1 & synth$White == 0), which(synth$Sex == 1 & synth$White == 1)
#' )
#' sample_apce <- CalAPCE(data = synth, mcmc.re = sample_mcmc, subgroup = subgroup_synth)
#'
#' @references Imai, K., Jiang, Z., Greiner, D.J., Halen, R., and Shin, S. (2023).
#' "Experimental evaluation of algorithm-assisted human decision-making:
#' application to pretrial public safety assessment."
#' Journal of the Royal Statistical Society: Series A.
#' <DOI:10.1093/jrsssa/qnad010>.
#'
#' @useDynLib aihuman, .registration=TRUE
#' @export
#'
CalAPCE <- function(data,
                    mcmc.re,
                    subgroup,
                    name.group = c("overall", "Sex0", "Sex1", "Sex1 White0", "Sex1 White1"),
                    rho = 0,
                    burnin = 0,
                    out.length = 500,
                    c0 = 0,
                    c1 = 0,
                    ZX = NULL,
                    save.individual.optimal.decision = FALSE,
                    parallel = FALSE,
                    optimal.decision.only = FALSE,
                    dmf = NULL,
                    fair.dmf.only = FALSE) {
  ####### burn in
  mcmc.length <- dim(mcmc.re)[1]
  mcmc.matrix <- as.matrix(mcmc.re)

  ######## dimension of the data
  p <- dim(data)[2] - 3
  k <- length(unique(data$D)) - 1
  n <- dim(data)[1]
  s1 <- subgroup[[1]] - 1
  s2 <- subgroup[[2]] - 1
  s3 <- subgroup[[3]] - 1
  s4 <- subgroup[[4]] - 1
  s5 <- subgroup[[5]] - 1
  C <- matrix(1, k + 2, k + 1)
  C[upper.tri(C)] <- 1 - c1
  C[lower.tri(C)] <- -c0
  C[k + 2, ] <- rep(-c0, k + 1)
  C_binary <- matrix(1, k + 2, 2)
  C_binary[-1, 1] <- -c0
  C_binary[k + 2, 2] <- -c0
  C_binary[1, 2] <- 1 - c1

  ####### data
  Z <- data$Z
  D <- data$D
  D_binary <- D
  D_binary[D > 0] <- 1
  Y <- data$Y
  X <- as.matrix(subset(data, select = -c(Z, D, Y)))
  if (is.null(ZX)) {
    ZX <- as.matrix(subset(data, select = -c(Z, D, Y)))
  }
  lZX <- ncol(ZX)

  # If the data for DMF recommendation is not given, use optimal decision to compute g_dmf
  if (is.null(dmf)) {
    dmf <- rep(0, n)
    dmf.null <- 1
  } else {
    dmf.null <- 0
  }

  if (burnin != 0) {
    BETA <- mcmc.matrix[-(1:floor(burnin * mcmc.length)), 1:(p + 1 + lZX)]
    THETA <- mcmc.matrix[-(1:floor(burnin * mcmc.length)), (p + 2 + lZX):(p + 2 * k + 1 + lZX)]
    ALPHA <- mcmc.matrix[-(1:floor(burnin * mcmc.length)), (p + 2 * k + 2 + lZX):(2 * p + 2 * k + 1 + lZX)]
    DELTA <- mcmc.matrix[-(1:floor(burnin * mcmc.length)), (2 * p + 2 * k + 2 + lZX):(2 * p + 3 * k + 2 + lZX)]
  } else {
    BETA <- mcmc.matrix[, 1:(p + 1 + lZX)]
    THETA <- mcmc.matrix[, (p + 2 + lZX):(p + 2 * k + 1 + lZX)]
    ALPHA <- mcmc.matrix[, (p + 2 * k + 2 + lZX):(2 * p + 2 * k + 1 + lZX)]
    DELTA <- mcmc.matrix[, (2 * p + 2 * k + 2 + lZX):(2 * p + 3 * k + 2 + lZX)]
  }
  length <- dim(BETA)[1]
  if (!optimal.decision.only & !fair.dmf.only) {
    P.D1.mcmc <- array(0, dim = c(length, 5, k + 1, k + 2))
    P.D0.mcmc <- array(0, dim = c(length, 5, k + 1, k + 2))
    APCE.mcmc <- array(0, dim = c(length, 5, k + 1, k + 2))
    Optimal.Z.mcmc <- array(0, dim = c(length, 5))
    P.DMF.mcmc <- array(0, dim = c(length, 5, k + 1, k + 2))
  }
  P.R.mcmc <- array(0, dim = c(length, 5, k + 2))
  if (fair.dmf.only) {
    P.DMF.mcmc <- array(0, dim = c(length, 5, k + 1, k + 2))
    for (j in 1:length) {
      if (j %% out.length == 0) {
        print(paste(j, "/", length, sep = ""))
      }
      beta <- BETA[j, ]
      theta <- THETA[j, ]
      alpha <- ALPHA[j, ]
      delta <- DELTA[j, ]
      re <- .CalAPCEj_ordinal_rcpp(X, ZX, rho, beta, alpha, theta, delta, p, k, n, s1, s2, s3, s4, s5, lZX, C, C_binary, 0, D, D_binary, dmf, dmf.null, 1)
      P.R.mcmc[j, , ] <- re$Pj_R
      P.DMF.mcmc[j, , , ] <- re$Pj_dmf
    }
  } else {
    if (optimal.decision.only & save.individual.optimal.decision) {
      Optimal.D.mcmc <- array(0, dim = c(n, k + 1))
      Utility.g_d.mcmc <- rep(0, n)
      Utility.g_dmf.mcmc <- rep(0, n)
      Utility.diff.control.mcmc <- rep(0, length)
      Utility.diff.treated.mcmc <- rep(0, length)
      for (j in 1:length) {
        if (j %% out.length == 0) {
          print(paste(j, "/", length, sep = ""))
        }
        beta <- BETA[j, ]
        theta <- THETA[j, ]
        alpha <- ALPHA[j, ]
        delta <- DELTA[j, ]
        re <- .CalAPCEj_ordinal_rcpp(X, ZX, rho, beta, alpha, theta, delta, p, k, n, s1, s2, s3, s4, s5, lZX, C, C_binary, 1, D, D_binary, dmf, dmf.null, 0)
        P.R.mcmc[j, , ] <- re$Pj_R
        Optimal.D.mcmc <- Optimal.D.mcmc + re$Optimal_D_ind
        Utility.g_d.mcmc <- Utility.g_d.mcmc + re$Utility_g_d_ind
        Utility.g_dmf.mcmc <- Utility.g_dmf.mcmc + re$Utility_g_dmf_ind
        Utility.diff.control.mcmc[j] <- mean(re$Utility_g_d_ind[which(Z == 0)] - re$Utility_g_dmf_ind[which(Z == 0)])
        Utility.diff.treated.mcmc[j] <- mean(re$Utility_g_d_ind[which(Z == 1)] - re$Utility_g_dmf_ind[which(Z == 1)])
      }
    } else if (optimal.decision.only & !save.individual.optimal.decision) {
      Optimal.D.mcmc <- array(0, dim = c(length, 5, k + 1))
      for (j in 1:length) {
        if (j %% out.length == 0) {
          print(paste(j, "/", length, sep = ""))
        }
        beta <- BETA[j, ]
        theta <- THETA[j, ]
        alpha <- ALPHA[j, ]
        delta <- DELTA[j, ]
        re <- .CalAPCEj_ordinal_rcpp(X, ZX, rho, beta, alpha, theta, delta, p, k, n, s1, s2, s3, s4, s5, lZX, C, C_binary, 1, D, D_binary, dmf, dmf.null, 0)
        P.R.mcmc[j, , ] <- re$Pj_R
        Optimal.D.mcmc[j, , ] <- re$Optimal_D
      }
    } else if (!optimal.decision.only & !save.individual.optimal.decision) {
      Optimal.D.mcmc <- array(0, dim = c(length, 5, k + 1))
      for (j in 1:length) {
        if (j %% out.length == 0) {
          print(paste(j, "/", length, sep = ""))
        }
        beta <- BETA[j, ]
        theta <- THETA[j, ]
        alpha <- ALPHA[j, ]
        delta <- DELTA[j, ]
        re <- .CalAPCEj_ordinal_rcpp(X, ZX, rho, beta, alpha, theta, delta, p, k, n, s1, s2, s3, s4, s5, lZX, C, C_binary, 0, D, D_binary, dmf, dmf.null, 0)
        P.D1.mcmc[j, , , ] <- re$Pj_D1
        P.D0.mcmc[j, , , ] <- re$Pj_D0
        APCE.mcmc[j, , , ] <- re$APCE
        P.R.mcmc[j, , ] <- re$Pj_R
        Optimal.Z.mcmc[j, ] <- re$Optimal_Z
        Optimal.D.mcmc[j, , ] <- re$Optimal_D
        P.DMF.mcmc[j, , , ] <- re$Pj_dmf
      }
    } else if (!optimal.decision.only & save.individual.optimal.decision) {
      Optimal.D.mcmc <- array(0, dim = c(n, k + 1))
      Utility.g_d.mcmc <- rep(0, n)
      Utility.g_dmf.mcmc <- rep(0, n)
      for (j in 1:length) {
        if (j %% out.length == 0) {
          print(paste(j, "/", length, sep = ""))
        }
        beta <- BETA[j, ]
        theta <- THETA[j, ]
        alpha <- ALPHA[j, ]
        delta <- DELTA[j, ]
        re <- .CalAPCEj_ordinal_rcpp(X, ZX, rho, beta, alpha, theta, delta, p, k, n, s1, s2, s3, s4, s5, lZX, C, C_binary, 0, D, D_binary, dmf, dmf.null, 0)
        P.D1.mcmc[j, , , ] <- re$Pj_D1
        P.D0.mcmc[j, , , ] <- re$Pj_D0
        APCE.mcmc[j, , , ] <- re$APCE
        P.R.mcmc[j, , ] <- re$Pj_R
        P.DMF.mcmc[j, , , ] <- re$Pj_dmf
        Optimal.Z.mcmc[j, ] <- re$Optimal_Z
        Optimal.D.mcmc <- Optimal.D.mcmc + re$Optimal_D_ind
        Utility.g_d.mcmc <- Utility.g_d.mcmc + re$Utility_g_d_ind
        Utility.g_dmf.mcmc <- Utility.g_dmf.mcmc + re$Utility_g_dmf_ind
      }
    }
  }

  # ### name
  name.D <- numeric(k + 1)
  for (d in 0:k) {
    name.D[d + 1] <- paste("D=", d, sep = "")
  }
  name.R <- numeric(k + 2)
  for (r in 0:(k + 1)) {
    name.R[r + 1] <- paste("R=", r, sep = "")
  }

  if (!optimal.decision.only & !fair.dmf.only) {
    dimnames(P.D1.mcmc) <- list(NULL, name.group, name.D, name.R)
    dimnames(P.D0.mcmc) <- list(NULL, name.group, name.D, name.R)
    dimnames(APCE.mcmc) <- list(NULL, name.group, name.D, name.R)
    dimnames(P.DMF.mcmc) <- list(NULL, name.group, name.D, name.R)
  } else if (fair.dmf.only) {
    dimnames(P.DMF.mcmc) <- list(NULL, name.group, name.D, name.R)
  }

  if (!parallel & save.individual.optimal.decision) {
    Optimal.D.mcmc <- Optimal.D.mcmc / length
    Utility.g_d.mcmc <- Utility.g_d.mcmc / length
    Utility.g_dmf.mcmc <- Utility.g_dmf.mcmc / length
  }

  if (fair.dmf.only) {
    res <- list(
      P.R.mcmc = P.R.mcmc,
      P.DMF.mcmc = P.DMF.mcmc
    )
  } else {
    if (!optimal.decision.only & save.individual.optimal.decision) {
      res <- list(
        P.D1.mcmc = P.D1.mcmc,
        P.D0.mcmc = P.D0.mcmc,
        APCE.mcmc = APCE.mcmc,
        P.R.mcmc = P.R.mcmc,
        Optimal.Z.mcmc = Optimal.Z.mcmc,
        Optimal.D.mcmc = Optimal.D.mcmc,
        Utility.g_d.mcmc = as.vector(Utility.g_d.mcmc),
        Utility.g_dmf.mcmc = as.vector(Utility.g_dmf.mcmc),
        P.DMF.mcmc = P.DMF.mcmc
      )
    } else if (optimal.decision.only & save.individual.optimal.decision) {
      res <- list(
        P.R.mcmc = P.R.mcmc,
        Optimal.D.mcmc = Optimal.D.mcmc,
        Utility.g_d.mcmc = as.vector(Utility.g_d.mcmc),
        Utility.g_dmf.mcmc = as.vector(Utility.g_dmf.mcmc),
        Utility.diff.control.mcmc = Utility.diff.control.mcmc,
        Utility.diff.treated.mcmc = Utility.diff.treated.mcmc
      )
    } else if (!optimal.decision.only & !save.individual.optimal.decision) {
      res <- list(
        P.D1.mcmc = P.D1.mcmc,
        P.D0.mcmc = P.D0.mcmc,
        APCE.mcmc = APCE.mcmc,
        P.R.mcmc = P.R.mcmc,
        Optimal.Z.mcmc = Optimal.Z.mcmc,
        Optimal.D.mcmc = Optimal.D.mcmc,
        P.DMF.mcmc = P.DMF.mcmc
      )
    } else if (optimal.decision.only & !save.individual.optimal.decision) {
      res <- list(
        P.R.mcmc = P.R.mcmc,
        Optimal.D.mcmc = Optimal.D.mcmc
      )
    }
  }

  ##### P.D1.mcmc is an array with dimension n.mcmc*5*(k+1)*(k+2) for quantity P(D(1)=d| R=r), dimension 1 is each posterior sample; dimension 2 is  subgroup, dimension 3 is (k+1) values of D from 0 to k, dimension 4 is (k+2) values of R from 0 to k+1
  ##### P.D0.mcmc is an array with dimension n.mcmc*5*(k+1)*(k+2) for quantity P(D(0)=d| R=r)
  ##### APCE.mcmc is an array with dimension n.mcmc*5*(k+1)*(k+2) for quantity P(D(1)=d| R=r)-P(D(0)=d| R=r)
  return(res)
}

#' Calculate APCE using parallel computing
#'
#' Calculate average principal causal effects (APCE) with ordinal decision using parallel computing. See Section 3.4 for more details.
#'
#' @param data A \code{data.frame} or \code{matrix} of which columns consists of pre-treatment covariates, a binary treatment (Z), an ordinal decision (D), and an outcome variable (Y). The column names of the latter three should be specified as "Z", "D", and "Y" respectively.
#' @param mcmc.re A \code{mcmc} object generated by \code{AiEvalmcmc()} function.
#' @param subgroup A list of numeric vectors for the index of each of the five subgroups.
#' @param name.group A list of character vectors for the label of five subgroups.
#' @param rho A sensitivity parameter. The default is  \code{0} which implies the unconfoundedness assumption (Assumption 4).
#' @param burnin A proportion of burnin for the Markov chain. The default is \code{0}.
#' @param out.length An integer to specify the progress on the screen. Every \code{out.length}-th iteration is printed on the screen. The default is \code{500}.
#' @param c0 The cost of an outcome. See Section 3.7 for more details. The default is \code{0}.
#' @param c1 The cost of an unnecessarily harsh decision. See Section 3.7 for more details. The default is \code{0}.
#' @param ZX The data matrix for interaction terms. The default is the interaction between Z and all of the pre-treatment covariates (X).
#' @param save.individual.optimal.decision A logical argument specified to save individual optimal decision rules. The default is \code{FALSE}.
#' @param optimal.decision.only A logical argument specified to compute only the optimal decision rule. The default is \code{FALSE}.
#' @param dmf A numeric vector of binary DMF recommendations. If \code{null}, use judge's decisions (0 if the decision is 0 and 1 o.w; e.g., signature or cash bond).
#' @param fair.dmf.only A logical argument specified to compute only the fairness of given DMF recommendations. The default is \code{FALSE}. Not used in the analysis for the JRSSA paper.
#' @param size The number of parallel computing. The default is \code{5}.
#'
#' @return An object of class \code{list} with the following elements:
#' \item{P.D1.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for quantity P(D(1)=d| R=r), dimension 1 is each posterior sample; dimension 2 is subgroup, dimension 3 is (k+1) values of D from 0 to k, dimension 4 is (k+2) values of R from 0 to k+1.}
#' \item{P.D0.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for quantity P(D(0)=d| R=r).}
#' \item{APCE.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for quantity P(D(1)=d| R=r)-P(D(0)=d| R=r).}
#' \item{P.R.mcmc}{An array with dimension n.mcmc by 5 by (k+2) for quantity P(R=r) for r from 0 to (k+1).}
#' \item{Optimal.Z.mcmc}{An array with dimension n.mcmc by 5 for the proportion of the cases where treatment (PSA provided) is optimal.}
#' \item{Optimal.D.mcmc}{An array with dimension n.mcmc by 5 by (k+1) for the proportion of optimal decision rule.}
#' \item{P.DMF.mcmc}{An array with dimension n.mcmc by 5 by (k+1) by (k+2) for the proportion of binary DMF recommendations. Not used in the analysis for the JRSSA paper.}
#' \item{Utility.g_d.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n for the individual utility of judge's decisions.}
#' \item{Utility.g_dmf.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n for the individual utility of DMF recommendation.}
#' \item{Utility.diff.control.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n.mcmc for estimated difference in utility between judge's decisions and DMF recommendation among control group.}
#' \item{Utility.diff.treated.mcmc}{Included if \code{save.individual.optimal.decision = TRUE}. An array with dimension n.mcmc for estimated difference in utility between judge's decisions and DMF recommendation among treated group.}
#'
#' @importFrom magrittr %>%
#' @importFrom purrr map
#' @importFrom purrr pluck
#' @importFrom abind abind
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#' @importFrom parallel detectCores
#' @importFrom doParallel registerDoParallel
#'
#'
#' @examples
#' \donttest{
#' data(synth)
#' sample_mcmc <- AiEvalmcmc(data = synth, n.mcmc = 10)
#' subgroup_synth <- list(
#'   1:nrow(synth), which(synth$Sex == 0), which(synth$Sex == 1),
#'   which(synth$Sex == 1 & synth$White == 0), which(synth$Sex == 1 & synth$White == 1)
#' )
#' sample_apce <- CalAPCEparallel(
#'   data = synth, mcmc.re = sample_mcmc,
#'   subgroup = subgroup_synth,
#'   size = 1
#' ) # adjust the size
#' }
#'
#' @references Imai, K., Jiang, Z., Greiner, D.J., Halen, R., and Shin, S. (2023).
#' "Experimental evaluation of algorithm-assisted human decision-making:
#' application to pretrial public safety assessment."
#' Journal of the Royal Statistical Society: Series A.
#' <DOI:10.1093/jrsssa/qnad010>.
#'
#' @useDynLib aihuman, .registration=TRUE
#' @export
#'
CalAPCEparallel <- function(data,
                            mcmc.re,
                            subgroup,
                            name.group = c("overall", "Sex0", "Sex1", "Sex1 White0", "Sex1 White1"),
                            rho = 0,
                            burnin = 0,
                            out.length = 500,
                            c0 = 0,
                            c1 = 0,
                            ZX = NULL,
                            save.individual.optimal.decision = FALSE,
                            optimal.decision.only = FALSE,
                            dmf = NULL,
                            fair.dmf.only = FALSE,
                            size = 5) {
  if (size == 1) {
    message("Increase the size or use unprallelized version instead.")
  } else {
    j <- NULL

    numCores <- detectCores()
    registerDoParallel(numCores)
    lb <- c(1, floor(nrow(mcmc.re) / size) * 1:(size - 1) + 1)
    ub <- c(floor(nrow(mcmc.re) / size) * 1:(size - 1), nrow(mcmc.re))
    apce_ls <- foreach(i = lb, j = ub) %dopar% {
      CalAPCE(data,
        mcmc.re[i:j, ],
        subgroup,
        name.group,
        rho,
        burnin,
        out.length,
        c0,
        c1,
        ZX,
        save.individual.optimal.decision,
        parallel = TRUE,
        optimal.decision.only,
        dmf,
        fair.dmf.only
      )
    }
    apce <- list()
    l <- names(apce_ls[[1]])
    for (i in l) {
      apce[[i]] <- apce_ls %>%
        map(pluck(i)) %>%
        abind(along = 1)
    }
    if (save.individual.optimal.decision) {
      apce[["Optimal.D.mcmc"]] <- apce_ls %>%
        map(pluck("Optimal.D.mcmc")) %>%
        Reduce("+", .data)
      apce[["Optimal.D.mcmc"]] <- apce[["Optimal.D.mcmc"]] / nrow(mcmc.re)
      apce[["Utility.g_d.mcmc"]] <- apce_ls %>%
        map(pluck("Utility.g_d.mcmc")) %>%
        Reduce("+", .data)
      apce[["Utility.g_d.mcmc"]] <- apce[["Utility.g_d.mcmc"]] / nrow(mcmc.re)
      apce[["Utility.g_dmf.mcmc"]] <- apce_ls %>%
        map(pluck("Utility.g_dmf.mcmc")) %>%
        Reduce("+", .data)
      apce[["Utility.g_dmf.mcmc"]] <- apce[["Utility.g_dmf.mcmc"]] / nrow(mcmc.re)
    }
    return(apce)
  }
}

#' Summary of APCE
#'
#' Summary of average principal causal effects (APCE) with ordinal decision.
#'
#' @param apce.mcmc APCE.mcmc array generated from \code{CalAPCE} or \code{CalAPCEparallel}.
#'
#' @return A \code{data.frame} that consists of mean and quantiles (2.5%, 97.5%, 5%, 95%) of APCE (P(D(1)=d| R=r)-P(D(0)=d| R=r)) for each subgroup given specific value of D (decision) and R (principal strata).
#'
#' @importFrom stats quantile
#'
#' @examples
#' \donttest{
#' data(synth)
#' sample_mcmc <- AiEvalmcmc(data = synth, n.mcmc = 10)
#' subgroup_synth <- list(
#'   1:nrow(synth), which(synth$Sex == 0), which(synth$Sex == 1),
#'   which(synth$Sex == 1 & synth$White == 0), which(synth$Sex == 1 & synth$White == 1)
#' )
#' sample_apce <- CalAPCE(data = synth, mcmc.re = sample_mcmc, subgroup = subgroup_synth)
#' sample_apce_summary <- APCEsummary(sample_apce[["APCE.mcmc"]])
#' }
#'
#' @references Imai, K., Jiang, Z., Greiner, D.J., Halen, R., and Shin, S. (2023).
#' "Experimental evaluation of algorithm-assisted human decision-making:
#' application to pretrial public safety assessment."
#' Journal of the Royal Statistical Society: Series A.
#' <DOI:10.1093/jrsssa/qnad010>.
#'
#' @useDynLib aihuman, .registration=TRUE
#' @export
#'
APCEsummary <- function(apce.mcmc) {
  s <- dim(apce.mcmc)[2] # subgroup
  kp1 <- dim(apce.mcmc)[3] # k+1
  qoi <- data.frame()
  n <- m <- 1
  for (i in 1:s) {
    for (k in 1:(kp1 + 1)) {
      for (j in 1:kp1) {
        mcmc <- apce.mcmc[, i, j, k]
        qoi[n, 1] <- dimnames(apce.mcmc)[[2]][i]
        qoi[n, 2] <- dimnames(apce.mcmc)[[3]][j]
        qoi[n, 3] <- dimnames(apce.mcmc)[[4]][k]
        qoi[n, 4] <- mean(mcmc, na.rm = TRUE)
        qoi[n, 5:6] <- quantile(mcmc, na.rm = TRUE, probs = c(0.025, 0.975))
        qoi[n, 7:8] <- quantile(mcmc, na.rm = TRUE, probs = c(0.05, 0.95))
        n <- n + 1
      }
    }
  }
  colnames(qoi) <- c("subgroup", "D", "R", "Mean", "2.5%", "97.5%", "5%", "95%")

  qoi$D <- as.numeric(gsub("D=", "", qoi$D))
  qoi$R <- as.numeric(gsub("R=", "", qoi$R))
  return(qoi)
}