目的

スピアマンの順位相関係数行列，あるいは，ケンドールの順位相関係数行列を計算する。
ピアソンの積率相関係数行列も計算できる。
この rank.cor 関数は，無相関検定の結果も返す。

使用法

rank.cor(x, y = NULL, use = c("complete.obs", "pairwise.complete.obs"), method = c("pearson", "spearman", "kendall"))

引数

x, y      データフレーム，行列，ベクトル
use       欠損値の取り扱い方を表す。"complete.obs" か "pairwise.complete.obs"
method    "pearson" か "spearman" か "kendall"

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/rank-cor.R", encoding="euc-jp")

# 順位相関係数行列の計算
rank.cor <- function(   x,                                                              # ベクトルまたは行列（データフレーム）
                        y = NULL,                                                       # ベクトルまたは行列（データフレーム）
                        use = c("complete.obs", "pairwise.complete.obs"),               # 欠損値の取り扱い方
                        method = c("pearson", "spearman", "kendall"))                   # 求める相関係数の種類
{
        p.values <- function(method, r, n, square.matrix)                               # P 値を計算する関数
        {
                if (method != "kendall") {                                              # ピアソンの積率相関係数，スピアマンの順位相関係数
                        p <- ifelse(r == 1, NA, pt(abs(r)*sqrt((n-2)/(1-r^2)), n-2, lower.tail=FALSE)*2)
                }
                else {                                                                  # ケンドールの順位相関係数
                        p <- pnorm(abs(r)/sqrt((4*n+10)/(9*n*(n-1))), lower.tail=FALSE)*2
                        if (square.matrix) {
                                diag(p) <- NA
                        }
                }
                return(p)
        }

        to.matrix <- function(x)
        {
                if (is.matrix(x)) {
                        return(x)
                }
                else if (is.data.frame(x) || is.vector(x)) {
                        return(as.matrix(x))
                }
                else {
                        stopifnot(is.matrix(x))
                }
        }

        apply.rank <- function(method, x)
        {
                return(if (method == "pearson") x else rank(x))
        }

        use <- match.arg(use)
        na.method <- pmatch(use, c("complete.obs", "pairwise.complete.obs"))+1
        method <- match.arg(method)
        square.matrix <- FALSE

        x <- to.matrix(x)
        ncx <- ncol(x)

        if (is.null(y)) { # x only
                square.matrix <- TRUE
                if (use == "complete.obs") {
                        n <- nrow(x <- x[complete.cases(x),])
                        if (method != "pearson") {
                                x <- apply(x, 2, rank)
                        }
                        r <- .Internal(cor(x, NULL, na.method, method == "kendall"))
                }
                else { # pairwise.complete.obs
                        n <- r <- matrix(0, ncx, ncx)
                        for (i in 1:ncx) {
                                for (j in 1:i) {
                                        ok <- complete.cases(x[,i], x[,j])
                                        n[i, j] <- n[j, i] <- sum(ok)
                                        r[i, j] <- r[j, i] <- if (i == j) 1 else .Internal(cor(apply.rank(method, x[ok, i]), apply.rank(method, x[ok, j]), na.method, method == "kendall"))
                                }
                        }
                }
        }
        else { # x and y 
                ncy <- ncol(y <- to.matrix(y))
                if (ncx == 1 && ncy == 1) { # vector x vector
                        n <- sum(ok <- complete.cases(x, y)) # complete.obs == pairwise.complete.obs
                        r <- as.numeric(.Internal(cor(apply.rank(method, x[ok]), apply.rank(method, y[ok]), na.method, method == "kendall")))
                }
                else { # matrix x matrix
                        if (use == "complete.obs") {
                                n <- sum(ok <- complete.cases(x, y))
                                x <- x[ok,]
                                y <- y[ok,]
                                if (method != "pearson") {
                                        x <- if (ncx == 1) rank(x) else apply(x, 2, rank)
                                        y <- if (ncy == 1) rank(y) else apply(y, 2, rank)
                                }
                                r <- .Internal(cor(x, y, na.method, method == "kendall"))
                        }
                        else { # pairwise.complete.obs
                                n <- r <- matrix(0, ncx, ncy)
                                for (i in 1:ncx) {
                                        for (j in 1:ncy) {
                                                n[i, j] <- sum(ok <- complete.cases(x[,i], y[,j]))
                                                r[i, j] <- .Internal(cor(apply.rank(method, x[ok, i]), apply.rank(method, y[ok, j]), na.method, method == "kendall"))
                                        }
                                }
                        }
                }
        }
        p <- p.values(method, r, n, square.matrix)
        return(list(method=sprintf("%s(%s)", method, use), correlation=r, n=n, p=p))
}


使用例

> x <- c(7, 9, 8,  0, NA, NA)
> y <- c(2, 3, 4, NA,  4,  3)

> rank.cor(x, y, use="pairwise.complete.obs", method="spearman")
$method
[1] "spearman(pairwise.complete.obs)"

$correlation
[1] 0.5	# スピアマンの順位相関係数

$n
[1] 3	# データのペア数

$p	# 無相関検定の P 値
[1] 0.6666667

> a <- cbind(x, y, 1:6)
> a
      x  y  
[1,]  7  2 1
[2,]  9  3 2
[3,]  8  4 3
[4,]  0 NA 4
[5,] NA  4 5
[6,] NA  3 6

> rank.cor(a, use="pairwise.complete.obs", method="spearman")
$method
[1] "spearman(pairwise.complete.obs)"

$correlation
     x         y           
x  1.0 0.5000000 -0.4000000	# スピアマンの順位相関係数
y  0.5 1.0000000  0.5270463
  -0.4 0.5270463  1.0000000

$n	# データのペア数
  x y  
x 4 3 4
y 3 5 5
  4 5 6

$p	# 無相関検定の P 値
          x         y          
x        NA 0.6666667 0.6000000
y 0.6666667        NA 0.3614549
  0.6000000 0.3614549        NA

> rank.cor(a, use="complete.obs", method="spearman")
$method
[1] "spearman(complete.obs deletion)"

$correlation
    x   y    
x 1.0 0.5 0.5	# スピアマンの順位相関係数
y 0.5 1.0 1.0
  0.5 1.0 1.0

$n	# データのペア数
[1] 3

$p	# 無相関検定の P 値
          x         y          
x        NA 0.6666667 0.6666667
y 0.6666667        NA 0.0000000
  0.6666667 0.0000000        NA


 解説ページ1，解説ページ2