R/test.R

In PAsso: Assessing the Partial Association Between Ordinal Variables

#' @title Hypothesis testing of the partial association coefficients
#'
#' @description This function use bootstrapping to conduct hypothesis testing
#' for the partial association coefficients. It directly applies onto the "PAsso"
#' class of object generated by "PAsso".
#'
#' @param object An object of "PAsso" class, which is generated by "PAsso" function.
#'
#' @param bootstrap_rep The number of bootstrap replications. It may be slow.
#' @param H0 null hypothesis of partial correlation coefficient.
#' @param parallel logical argument whether conduct parallel for bootstrapping partial association.
#'
#' @importFrom progress progress_bar
#' @importFrom foreach foreach %dopar%
#'
#' @importFrom utils combn setTxtProgressBar txtProgressBar
#'
#' @export
#' @examples
#' # Import ANES2016 data in "PAsso"
#' data(ANES2016)
#' # Parial association:
#' PAsso_2v <- PAsso(responses = c("PreVote.num", "PID"),
#'                 adjustments = c("income.num", "age", "edu.year"),
#'                 data = ANES2016)
#'
#' summary(PAsso_2v, digits=4)
#'
#' PAsso_2v_test <- test(object = PAsso_2v, bootstrap_rep=20, H0=0, parallel=FALSE)
#' PAsso_2v_test
#'
test <- function(object, bootstrap_rep=300, H0=0, parallel=FALSE) {
  if (!inherits(object, "PAsso")) stop("Input object must be 'PAsso' class.")

  responses <- attr(object, "responses")
  adjustments <- attr(object, "adjustments")
  arguments <- attr(object, "arguments")
  models <- attr(object, "models")

  MatCor <- object$corr

  n_responses <- length(responses)

  mods_n <- object$mods_n

  n_corr <- n_responses*(n_responses-1)/2
  boot_Cor_temp <- matrix(NA, ncol = bootstrap_rep, nrow = n_corr)
  # A matrix to save all bootstrap correlation coefficients in upper triangle matrix.
  temp_pair <- combn(responses, 2) # result is a matrix, use each column!

  pair_list <- split(x = temp_pair, f = rep(1:ncol(temp_pair), each = nrow(temp_pair)))

  names(pair_list) <- apply(temp_pair, 2, FUN = function(x) paste(x, collapse = "-"))


  # Add progress bar --------------------------------------------------------
  pb <- progress_bar$new(
    format = "Replication = :letter [:bar] :percent :elapsed | eta: :eta",
    total = bootstrap_rep,    # 300
    width = 80)
  progress_repNo <- c(1:bootstrap_rep)  # token reported in progress bar


  # [FEATURE]: standard error from bootstrap: Done!
  if (arguments[1] == "partial") { # association=="partial"
    # StillBoot <- utils::menu(choices = c("Yes", "No"),
                             # title="Using bootstrap to conduct inference of partial association phi could be slow. \nDo you really want to continue?")
    message("The bootstrapping procedure may be slow when bootstrap_rep is large!")

    if (parallel) { # Use parallel for bootstrapping and "progress" package
      # Below is for .options.snow = opts: allowing progress bar to be used in foreach -----------------------------
      # progress <- function(n){
      #   pb$tick(tokens = list(letter = progress_repNo[n]))
      # }
      # opts <- list(progress = progress)

      data_temp <- object$data
      boot_Cor_temp <-
        foreach(i=1:bootstrap_rep,
                .packages = c('MASS', 'stats', 'PAsso'),
                # .export=c("mods_n", "data_temp", "arguments"),
                # Remove this for Warning message:
                  # In e$fun(obj, substitute(ex), parent.frame(), e$data) :
                  # already exporting variable(s): mods_n, data_temp, arguments
                .combine=cbind) %dopar% {
                  # ProgressBar
                  pb$tick(tokens = list(letter = progress_repNo[i]))

                  # tryCatch({
                    index <- sample(mods_n[1], replace=T)
                    boot_data <- as.data.frame(data_temp[index, ])
                    Pcor_temp <-
                      PAsso(data = boot_data,
                            responses = attr(object, "responses"),
                            adjustments = attr(object, "adjustments"),
                            models= attr(object, "models"), # in "PAsso" recover models arguments.
                            association = arguments[1],
                            method = arguments[2],
                            resids.type = arguments[3])
                    Pcor_temp$corr[upper.tri(Pcor_temp$corr)]
                  # }, error=function(e){
                  #   message("Error in bootstrap", i, ":\n")
                  #   message(error_message)
                  # })
                }

    } else { # Without using parallel.

      for(i in 1:bootstrap_rep){
        tryCatch({
          index <- sample(mods_n[1], replace=T)
          boot_data <- object$data[index, ]
          Pcor_temp <-
            PAsso(data = boot_data,
                  responses = attr(object, "responses"),
                  adjustments = attr(object, "adjustments"),
                  models = attr(object, "models"),
                  association = arguments[1],
                  method = arguments[2],
                  resids.type = arguments[3])
          boot_Cor_temp[,i] <- Pcor_temp$corr[upper.tri(Pcor_temp$corr)]
        }, error=function(e){
          message("Error in bootstrap", i, ":\n")
          # message(error_message)
        })
        # ProgressBar
        pb$tick(tokens = list(letter = progress_repNo[i]))
      }

    }

      # Calculate standard error of partial association coefficients!
      sd_MatCor <- matrix(NA, nrow = n_responses, ncol = n_responses)
      # sd_MatCor[upper.tri(sd_MatCor)] <- sqrt(matrixStats::rowVars(boot_Cor_temp))
      sd_MatCor[upper.tri(sd_MatCor)] <- apply(X = boot_Cor_temp, MARGIN = 1, FUN = sd)

    # Return values:
    boot_left <- sum(complete.cases(t(boot_Cor_temp)))
    boot_Cor_left <- matrix(data = boot_Cor_temp[, complete.cases(t(boot_Cor_temp))],
                            nrow = length(pair_list), ncol = boot_left)
    # Make above one as matrix to avoid bug when use apply! Keep complete column

    corr_stat <- corr_p.value <- matrix(NA, nrow = n_responses, ncol = n_responses)

    corr_stat_up <- corr_stat[upper.tri(corr_stat)] <-
      (MatCor[upper.tri(MatCor)] - H0)/sd_MatCor[upper.tri(sd_MatCor)]

    if (H0==0) {
      corr_p.value[upper.tri(corr_p.value)] <-
        2 * apply(
          X = cbind(
            # matrixStats::rowMeans2(boot_Cor_left>0),
            # matrixStats::rowMeans2(boot_Cor_left<0)
            apply(X = (boot_Cor_left > 0), MARGIN = 1, FUN = mean),
            apply(X = (boot_Cor_left < 0), MARGIN = 1, FUN = mean)
          ),
          MARGIN = 1,
          FUN = min)
    } else {
      corr_p.value[upper.tri(corr_p.value)] <-
        2 * apply(
          X = cbind(
            # matrixStats::rowMeans2(boot_Cor_left > -1*H0),
            # matrixStats::rowMeans2(boot_Cor_left < H0),
            apply(X = (boot_Cor_left > -1*H0), MARGIN = 1, FUN = mean),
            apply(X = (boot_Cor_left < H0), MARGIN = 1, FUN = mean),
            rep(0.5, n_corr)),
          MARGIN = 1,
          FUN = min)
    }

    # Save components: S.E.; statistics; p-value; CI_95.
    object$sd_MatCor <- sd_MatCor
    object$corr_stat <- corr_stat
    object$corr_p.value <- corr_p.value

    # Partial correlation test may be asymmetric, so use quantile.
    CI_temp1 <- apply(X = boot_Cor_left, MARGIN = 1,
                      function(x) quantile(x, probs = c(0.025, 0.975)))
      # rbind(MatCor[upper.tri(MatCor)] - qnorm(0.975) * MatCor[upper.tri(MatCor)] / corr_stat_up,
            # MatCor[upper.tri(MatCor)] + qnorm(0.975) * MatCor[upper.tri(MatCor)] / corr_stat_up)

    CI_temp <- split(x = CI_temp1, f = rep(1:ncol(CI_temp1), each = nrow(CI_temp1)))
    names(CI_temp) <- names(pair_list)

    object$CI_95 <- CI_temp

    # Return a class of "PAsso.test" inheriting from "PAsso".
    attr(object, "bootstrap_rep") <- bootstrap_rep
    attr(object, "H0") <- H0
    class(object) <- c("PAsso.test", class(object)[-1])
    object

  } else {
    #  association=="marginal", this can be replaced simply by "cor.test()"

    deal_pair <- function(pair) {
      resp_x <- as.numeric(object$data[,as.character(pair[1])])
      resp_y <- as.numeric(object$data[,as.character(pair[2])])
      cor.test(resp_x,resp_y, method = arguments[2])
      }
    pair_cortest <- lapply(pair_list, FUN = deal_pair) # Correlation tests of all pairs!
    object$pair_cortest <- pair_cortest

    corr_stat <- corr_p.value <- matrix(NA, nrow = n_responses, ncol = n_responses)
    corr_stat[upper.tri(corr_stat)] <- sapply(pair_cortest, FUN = function(x) x$statistic)
    corr_p.value[upper.tri(corr_p.value)] <- sapply(pair_cortest, FUN = function(x) x$p.value)

    sd_MatCor <- sapply(pair_cortest, FUN = function(x) as.vector(x$estimate/x$statistic))

    # Save components: S.E.; statistics; p-value; CI_95.
    object$sd_MatCor <- sd_MatCor
    object$corr_stat <- corr_stat
    object$corr_p.value <- corr_p.value

    # Marginal correlation test is conducted as symmetric.
    CI_temp <- lapply(pair_cortest,
                      FUN = function(x) {
                        as.vector(c(x$estimate - qnorm(0.975)*x$estimate/x$statistic,
                                    x$estimate + qnorm(0.975)*x$estimate/x$statistic))
                      })
    object$CI_95 <- CI_temp

    # Return an object of class "MAsso.test" inheriting from "MAsso".
    attr(object, "bootstrap_rep") <- c(bootstrap_rep)
    attr(object, "H0") <- H0
    class(object) <- c("MAsso.test", class(object)[-1])
    object

  }


}

#' @keywords internal
#' This function is use to transfer the bootstrap pvalue to the confidence level it should be. For bootstrap_rep=20,
#' the actual confidence level is only 1/20.
format.pval.corr <-
  function(pv, digits = max(1L, getOption("digits") - 2L),
           eps = .Machine$double.eps,
           na.form = "NA", ..., bootstrap_rep # change .Machine$double.eps to 1/boot
           ) {
    # pv = x$corr_p.value[1,];

    eps <- 1/bootstrap_rep

    if ((has.na <- any(ina <- is.na(pv))))
      pv <- pv[!ina]
    r <- character(length(is0 <- pv < eps))
    if (any(!is0)) {
      rr <- pv <- pv[!is0]
      expo <- floor(log10(ifelse(pv > 0, pv, 1e-50)))
      fixp <- expo >= -3 | (expo == -4 & digits > 1)
      if (any(fixp))
        rr[fixp] <- format(pv[fixp], digits = digits, ...)
      if (any(!fixp))
        rr[!fixp] <- format(pv[!fixp], digits = digits, ...)
      r[!is0] <- rr
    }
    if (any(is0)) {
      digits <- max(1L, digits - 2L)
      if (any(!is0)) {
        nc <- max(nchar(rr, type = "w"))
        if (digits > 1L && digits + 6L > nc)
          digits <- max(1L, nc - 7L)
        sep <- if (digits == 1L && nc <= 6L)
          ""
        else " "
      }
      else sep <- if (digits == 1)
        ""
      else " "
      r[is0] <- paste("<", format(eps, digits = digits, ...),
                      sep = sep)
    }
    if (has.na) {
      rok <- r
      r <- character(length(ina))
      r[!ina] <- rok
      r[ina] <- na.form
    }
    r
  }


#' @rdname print
#' @method print PAsso.test
#'
#' @export
print.PAsso.test <- function(x, digits = max(3L, getOption("digits") - 3L), ...) {

  x$corr[lower.tri(x$corr)] <- NA

  cat("-------------------------------------------- \n")
  cat("The partial association analysis: \n\n")
  n_rep <- dim(x$corr)[1]
  mat_tem <- matrix(NA, nrow = n_rep*4, ncol = n_rep)

  bootstrap_rep <- attr(x, "bootstrap_rep")

  for (i in 1:(n_rep-1)) {
    mat_tem[(4*i-3),i:n_rep] <-
      format(round(x$corr[i, i:n_rep], digits = max(2, digits)),
             digits = max(2, digits), ...)

    mat_tem[(4*i-2),(i+1):n_rep] <-
      format(round(x$sd_MatCor[i,(i+1):n_rep], digits = max(2, digits)),
             digits = max(2, digits), ...)

    temp_p <- x$corr_p.value[i,(i+1):n_rep]
    # Need to compare temp_p with boot_p!
    temp_p[temp_p < (1/bootstrap_rep)] <- 1/bootstrap_rep

    Cf_p <- format.pval.corr(round(temp_p, digits=max(2, digits)),
                             bootstrap_rep = bootstrap_rep,
                             digits = max(2, digits), ...)

    Signif <- symnum(temp_p, corr = FALSE, na = FALSE,
                     cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                     symbols = c("***", "**", "*", ".", " "))
    mat_tem[(4*i-1),(i+1):n_rep] <- paste(Cf_p, format(Signif), sep = "")
    # S.E.
  }
  mat_tem[(4*n_rep-3),n_rep] <- format(x$corr[n_rep, n_rep],
                                       digits = max(2, digits),
                                       nsmall= max(2, digits), ...)

  colnames(mat_tem) <- colnames(x$corr)
  temp_rowname <- rep(rownames(x$corr), each=4)

  for (i in 1:(n_rep)) {
    temp_rowname[(4*i-2)] <- "S.E."
    temp_rowname[(4*i-1)] <- "Pr"
    temp_rowname[(4*i)] <- "---"
  }

  rownames(mat_tem) <- temp_rowname

  print.default(mat_tem, na.print = "",
                print.gap = 2, quote = FALSE, ...)

  # Add stars of significance level --------------------------------------------------
  if ((w <- getOption("width")) <
      nchar(sleg <- attr(Signif, "legend"))) {
    sleg <- strwrap(sleg, width = w - 2, prefix = "  ")
  }
  cat("Signif. codes:  ", sleg, sep = "",
      fill = w + 4 + max(nchar(sleg, "bytes") - nchar(sleg)))
}