# R/corr.test.R In psych: Procedures for Psychological, Psychometric, and Personality Research

```#modified 9/18/21 to find Kendall and Spearman probabilities
"corr.test" <- "corTest" <-
cl <- match.call()

if(normal) {     #the normal case
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*(1 - pt(abs(t),(n-2)))
p <- -2 *  expm1(pt(abs(t),(n-2),log.p=TRUE))  #suggested by Nicholas Clark
se <- sqrt((1-r*r)/(n-2))
} else {temp <- prob.Kendall(x=x,y=y,method=method)
if(is.null(y)) { sym <- TRUE
n <- t(!is.na(x)) %*% (!is.na(x))} else {sym <- FALSE
n <- t(!is.na(x)) %*% (!is.na(y))}
r <- temp\$r
p <- temp\$p

se <- sqrt((1-r*r)/(n-2))
}

nvar <- ncol(r)

p[p>1] <- 1
pa <- p  #in case we don't do adjustments

if (is.null(y)) {lp <- upper.tri(p)  #the case of a symmetric matrix  (remember upper.tri =! lower tri)
pa <- p[lp]
p[upper.tri(p,diag=FALSE)] <- pa
} else {
} }
#find confidence intervals
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(sym) {ncors <- nvar * (nvar-1)/2} else ncors <- prod(dim(r))
if(adjust!="holm") {dif.corrected <- qnorm(1-alpha/(2* ncors)) } else {   # 1- alpha/2  /nvar *(nvar-1) /2)
ord <- order(abs(z),decreasing=FALSE)  #to find the HOlm correction, we need to order the size of the correlations
dif.corrected <- qnorm(1-alpha/(2*order(ord))) } #holm

alpha <- 1-alpha/2  #the raw alpha level for confidence intervals
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)])

cnR <- abbreviate(rownames(r),minlength=minlength)
cnC <- abbreviate(colnames(r),minlength=minlength)

k <- 1
for(i in 1:(nvar-1)) {for (j in (i+1):nvar) {
rownames(ci)[k] <- paste(cnC[i],cnR[j],sep="-")
k<- k +1 }}

} else { #non symmetric case
n.x <- NCOL(x)
n.y <- NCOL(y)
z <- fisherz(r)
} else {ord <- order(abs(z),decreasing=FALSE)  #to find the HOlm correction, we need to order the size of the correlations
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))
ci2 <- data.frame(lower=lower,r=as.vector(r),upper=upper,p=as.vector(p),pa=as.vector(pa))  #adding pa breaks two other packages
cnR <- abbreviate(rownames(r),minlength=minlength)  #added minlength as a parameter than fixed to 5  5/28/18
cnC <- abbreviate(colnames(r),minlength=minlength)
k <- 1
for(i in 1:NCOL(y)) {for (j in 1:NCOL(x)) {
rownames(ci)[k] <- paste(cnR[j],cnC[i],sep="-")
rownames(ci2)[k] <- paste(cnR[j],cnC[i],sep="-")
k<- k +1 }}
}
} else {ci <-  ci2 <-  sef <- ci.adj <- NULL
}
class(result) <- c("psych", "corr.test")
return(result)
}
#modified 1/4/14 to report sample size once if they are all equal
#modified 3/12/14 to report confidence intervals (suggested by Alexander Weiss)
#modified 3/27/14 to correct bug detected by Clemens Fell
#modified 3/27/14 to correct bug reported by Louis-Charles Vannier
#modified 2/21/15 to make confidence intervals an option (incredible decrease in speed if doing cis)
#modified 8/24/17 to include Bonferoni adjusted confidence intervals
#modified 3/27/21 to include the adjusted p values in the long output (requested by Alexander Weiss)

"corr.p" <-
cl <- match.call()
if(missing(n)) stop("The number of subjects must be specified")
sym <- FALSE
t <- (r*sqrt(n-2))/sqrt(1-r^2)
#p <- 2*(1 - pt(abs(t),(n-2)))
p <- -2 *  expm1(pt(abs(t),(n-2),log.p=TRUE))
pa <- p
p[p>1] <- 1
if(isSymmetric(unclass(p))) {sym <- TRUE
lp <- upper.tri(p)  #the case of a symmetric matrix
pa <- p[lp]
p[upper.tri(p,diag=FALSE)] <- pa
} else {
sym <- FALSE}
}
if(ci) {
nvar <- ncol(r)
if(sym) {z <- fisherz(r[lower.tri(r)])} else {z <- fisherz(r)
n.x <- NCOL(r)
n.y <- NROW(r)
if(adjust != "holm") {dif.corrected <- qnorm((1-alpha/2)/(n.x * n.y)) # adjust alpha by 2
} else {ord <- order(abs(z),decreasing=FALSE)  #to find the Holm correction, we need to order the size of the correlations
dif.corrected <- qnorm(1-alpha/(2*(order(ord))))
}}
if (min(n) < 4) {
warning("Number of subjects must be greater than 3 to find confidence intervals.")
}
if(sym & is.matrix(n)) {
se <- 1/sqrt(n[lower.tri(n)] - 3) } else { se <- 1/sqrt(n - 3)}
if(sym) {    dif.corrected <- qnorm(1-alpha/(nvar*(nvar-1)))  } # 1- alpha/2  /nvar *(nvar-1) /2
alpha <- 1-alpha/2
dif <- qnorm(alpha)
lower <- fisherz2r(z - dif * se)
upper <- fisherz2r(z + dif * se)
lower.corrected <- fisherz2r(z - dif.corrected * se)
upper.corrected <- fisherz2r(z + dif.corrected * se)
if(sym) {ci <- data.frame(lower=lower,r=r[lower.tri(r)],upper=upper,p=p[lower.tri(p)])
ci <- data.frame(lower=as.vector(lower),r=as.vector(r),upper=as.vector(upper),p=as.vector(p))
cnR <- abbreviate(colnames(r),minlength=minlength)
rnR <- abbreviate(rownames(r),minlength=minlength)
if(sym) {k <- 1
for(i in 1:(nvar-1)) {for (j in (i+1):nvar) {
rownames(ci)[k] <- paste(cnR[i],rnR[j],sep="-")
k<- k +1 } }

} else {k <- 1
for(i in 1:ncol(r)) {for (j in 1:nrow(r)) {
rownames(ci)[k] <- paste(cnR[i],rnR[j],sep="-")
k<- k +1 }}
}
result <- list(r=r,n=n,p=p,Call=cl)}
class(result) <- c("psych", "corr.p")
return(result)
}
#revised March 28, 2014 to be compatible with corr.test
#revised August 28, 2017 to include holm and bonferroini adjusted confidence intervals
#revised February 17, 2019 to allow cis not to be found
#could be replaced with the following
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)

result\$Call<- cl
class(result) <- c("psych", "corr.test")
return(result)
}

#added April 6, 2021 in response to a request by Uthpala Pinto
symp <- symnum(p, corr = FALSE,cutpoints = c(0,  .001,.01,.05, 1),
symbols = c("***","**","*"," "),legend=FALSE)
stars <- paste0(round(r,digits),symp)
nc <- NCOL(r)
nr <- NROW(r)
stars <- matrix(stars,nrow=nr,ncol=nc)
colnames(stars) <- colnames(r)
rownames(stars) <- rownames(r)
# print(stars,quote=FALSE)
return(stars)
}

#added 9/18/21 to find Spearman and Kendall probabilities
prob.Kendall <- function(x=x,y=NULL,method=c("kendall","spearman"),continuity=FALSE) {

exact <- FALSE
alternative <- "two.sided"

nvar <- NCOL(x)
if(!is.null(y)) {nvary <- NCOL(y)
sym <- FALSE} else {nvary<- nvar
sym <- TRUE}
r <- diag(.5,nvar,nvary)
p <- diag(0,nvar,nvary)
if(sym) {
for(i in 2:(nvar) ) {
xi <- x[,i]
for (j in 1:(i-1)){
yj <- x[,j]
temp <- cor.test(xi,yj,method,alternative=alternative,exact=exact)
r[i,j] <- temp\$estimate
p[i,j]<- temp\$p.value
}
}
r <- r + t(r)
p <- p + t(p)
} else {
for (i in 1:nvar) {
if(nvar > 1) {xi <- x[,i]} else {xi <- x}  #the weird case of a single variable
for (j in 1: nvary) {
if(nvary > 1) {yj <- y[,j]} else { yj <- y}
temp <- cor.test(xi,yj,method,alternative=alternative,exact=exact)
r[i,j] <- temp\$estimate
p[i,j]<- temp\$p.value
}}
}
rownames(r) <- rownames(p) <- colnames(x)
if(is.null (y) ) {colnames(r) <- colnames(p) <- colnames(x)} else {colnames(r) <- colnames(p) <- colnames(y)}
result <- list(r=r,p=p)
return(result)
}

```

## Try the psych package in your browser

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

psych documentation built on June 27, 2024, 5:07 p.m.