# R/test.R In PAsso: Assessing the Partial Association Between Ordinal Variables

#### Documented in print.PAsso.testtest

```#' @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")
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 = "-"))

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"),
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"),
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)))
}
```

## Try the PAsso package in your browser

Any scripts or data that you put into this service are public.

PAsso documentation built on June 18, 2021, 5:09 p.m.