# Confidence intervals of rank correlation coefficients # 14-15 March 2016. Minato Nakazawa # See, http://minato.sip21c.org/ebhc-text.pdf for Japanese explanation # References: # http://support.sas.com/documentation/cdl/en/procstat/63104/HTML/default/corr_toc.htm spearman.ci.sas <- function(x, y, adj.bias=TRUE, conf.level=0.95) { n <- length(x) rx <- rank(x) ry <- rank(y) mx <- mean(rx) my <- mean(ry) rho <- sum((rx-mx)*(ry-my))/sqrt(sum((rx-mx)^2)*sum((ry-my)^2)) adj <- ifelse(adj.bias, rho/(2*(n-1)), 0) z <- 1/2*log((1+rho)/(1-rho)) gg <- qnorm(1-(1-conf.level)/2)*sqrt(1/(n-3)) ge <- z - adj gl <- ge - gg gu <- ge + gg rl <- (exp(2*gl)-1)/(exp(2*gl)+1) re <- (exp(2*ge)-1)/(exp(2*ge)+1) ru <- (exp(2*gu)-1)/(exp(2*gu)+1) print(sprintf("rho = %5.3f, %2d%% conf.int = [ %5.3f, %5.3f ]", re, conf.level*100, rl, ru)) return(list(rho=re, rho.ll=rl, rho.ul=ru, adj.bias=adj.bias)) } data(trees) plot(trees$Height, trees$Volume) cor.test(trees$Height, trees$Volume, method="spearman") cor.test(rank(trees$Height), rank(trees$Volume)) spearman.ci.sas(trees$Height, trees$Volume) spearman.ci.sas(trees$Height, trees$Volume, adj.bias=TRUE) if (require(RVAideMemoire)==FALSE) { install.packages("RVAideMemoire", dep=TRUE) library(RVAideMemoire) } spearman.ci(trees$Height, trees$Volume, nrep=1000) cor.test(trees$Height, trees$Volume, method="kendall") if (require(NSM3)==FALSE) { install.packages("NSM3", dep=TRUE) library(NSM3) } kendall.ci(trees$Height, trees$Volume, bootstrap=FALSE) kendall.ci(trees$Height, trees$Volume, bootstrap=TRUE, B=1000)