#' @title Correlation matrix
#' @description Based on corr.test (package: psych)
#' @export
corr.test_v2 <- function (x,
y = NULL,
use = "pairwise",
method = "pearson",
adjust = "holm",
alpha = 0.05,
ci = TRUE)
{
cl <- match.call()
if (is.null(y)) {
r <- cor(x, use = use, method = method)
sym <- TRUE
n <- t(!is.na(x)) %*% (!is.na(x))
}
else {
r <- cor(x, y, use = use, method = method)
sym = FALSE
n <- t(!is.na(x)) %*% (!is.na(y))
}
if ((use == "complete") | (min(n) == max(n)))
n <- min(n)
t <- (r * sqrt(n - 2))/sqrt(1 - r^2)
p <- -2 * expm1(pt(abs(t), (n - 2), log.p = TRUE))
se <- sqrt((1 - r * r)/(n - 2))
nvar <- ncol(r)
p[p > 1] <- 1
if (adjust != "none") {
if (is.null(y)) {
lp <- upper.tri(p)
pa <- p[lp]
pa <- p.adjust(pa, adjust)
p[upper.tri(p, diag = FALSE)] <- pa
}
else {
p[] <- p.adjust(p, adjust)
}
}
z <- fisherz(r[lower.tri(r)])
if (ci) {
if (min(n) < 4) {
warning("Number of subjects must be greater than 3 to find confidence intervals.")
}
if (adjust != "holm") {
dif.corrected <- qnorm(1 - alpha/(nvar * (nvar -
1)))
}
else {
ord <- order(abs(z), decreasing = FALSE)
dif.corrected <- qnorm(1 - alpha/(2 * order(ord)))
}
alpha <- 1 - alpha/2
dif <- qnorm(alpha)
if (sym) {
if (is.matrix(n)) {
sef <- 1/sqrt(n[lower.tri(n)] - 3)
}
else {
sef <- 1/sqrt(n - 3)
}
lower <- fisherz2r(z - dif * sef)
upper <- fisherz2r(z + dif * sef)
lower.corrected <- fisherz2r(z - dif.corrected *
sef)
upper.corrected <- fisherz2r(z + dif.corrected *
sef)
ci <- data.frame(lower = lower, r = r[lower.tri(r)],
upper = upper, p = p[lower.tri(p)])
ci.adj <- data.frame(lower.adj = lower.corrected,
upper.adj = upper.corrected)
cnR <- colnames(r)
k <- 1
for (i in 1:(nvar - 1)) {
for (j in (i + 1):nvar) {
rownames(ci)[k] <- paste(cnR[i], cnR[j], sep = "-")
k <- k + 1
}
}
}
else {
n.x <- NCOL(x)
n.y <- NCOL(y)
z <- fisherz(r)
if (adjust != "holm") {
dif.corrected <- qnorm(1 - (1 - alpha)/(n.x *
n.y))
}
else {
ord <- order(abs(z), decreasing = FALSE)
dif.corrected <- qnorm(1 - (1 - alpha)/(order(ord)))
}
sef <- 1/sqrt(n - 3)
lower <- as.vector(fisherz2r(z - dif * sef))
upper <- as.vector(fisherz2r(z + dif * sef))
lower.corrected <- fisherz2r(z - dif.corrected *
sef)
upper.corrected <- fisherz2r(z + dif.corrected *
sef)
ci <- data.frame(lower = lower, r = as.vector(r),
upper = upper, p = as.vector(p))
ci.adj <- data.frame(lower.adj = as.vector(lower.corrected),
r = as.vector(r), upper.adj = as.vector(upper.corrected))
cnR <- rownames(r)
cnC <- colnames(r)
k <- 1
for (i in 1:ncol(y)) {
for (j in 1:ncol(x)) {
rownames(ci)[k] <- paste(cnR[j], cnC[i], sep = "-")
k <- k + 1
}
}
}
}
else {
ci <- sef <- ci.adj <- NULL
}
result <- list(r = r, n = n, t = t, p = p, se = se, sef = sef,
adjust = adjust, sym = sym, ci = ci, ci.adj = ci.adj,
Call = cl)
class(result) <- c("psych", "corr.test")
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.