目的

分割表形式で与えられたデータに基づいて，ケンドールのτb（ケンドールの順位相関係数）を計算する。

使用法

tau.b(x)

引数

x	分割表（合計欄を含めない）

ソース

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

# 分割表形式で与えられたデータに基づいて，ケンドールのτb（ケンドールの順位相関係数）を計算する
tau.b <- function(f)                                                 # 分割表（合計欄を含めない）
{
        R <- nrow(f)                                                 # 行数
        C <- ncol(f)                                                 # 列数
        n <- sum(f)                                                  # 全数
        rt <- rowSums(f)                                             # 行和
        ct <- colSums(f)                                             # 列和
        dr <- n^2-sum(rt^2)
        dc <- n^2-sum(ct^2)
        g <- matrix(0, nr=R+2, nc=C+2)
        cada <- f
        g[2:(R+1), 2:(C+1)] <- f
        PQ <- 0
        for (i in 1:R) {
                for (j in 1:C) {
                        cada[i, j] <-        sum(g[1:i, 1:j], g[(i+2):(R+2), (j+2):(C+2)]) -
                                        sum(g[1:i, (j+2):(C+2)], g[(i+2):(R+2), 1:j])
                        PQ <- PQ+g[i+1, j+1]*cada[i, j]
                }
        }
        taub <- PQ/sqrt(dr*dc)                                               # τb
        ase0 <- 2*sqrt((sum(f*cada^2)-PQ^2/n)/(dr*dc))                       # 標準誤差
        Z <- taub/ase0                                                       # 検定統計量
        P <- pnorm(abs(Z), lower.tail=FALSE)*2                               # P 値
        for (i in 1:R) {
                for (j in 1:C) {
                        f[i, j] <- f[i, j]*(2*sqrt(dr*dc)*cada[i, j]+taub*(rt[i]*dc+ct[j]*dr))^2
                }
        }
        ase1 <- sqrt(sum(f)-n^3*taub^2*(dr+dc)^2)/(dr*dc)            # 別法（よく使われる方法）
        Z2 <- taub/sqrt((4*n+10)/(9*n*(n-1)))                                # ケンドールの順位相関係数の検定のときに出てくる式
        P2 <- pnorm(abs(Z2), lower.tail=FALSE)*2                     # P 値
        c("tau b"=taub, "ase1"=ase1, "ase0"=ase0, "Z value"=Z, "P value"=P, "Z value 2"=Z2, "P value 2"=P2)
}


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

# 分割表データから，元のデータを構成し，R にある cor.test を使う。
tau.b2 <- function(f)
{
        x <- rep(row(f), f)
        y <- rep(col(f), f)
        cor.test(x, y, method="kendall")
}


使用例

f  <- matrix(c(
	0, 3, 1, 1,
	1, 4, 4, 0,
	1, 3, 0, 1,
	1, 0, 1, 0), byrow=TRUE, nc=4)	# 4 行 4 列の分割表

tau.b(f)

出力結果例

> tau.b(f)
     tau b       ase1       ase0    Z value    P value  Z value 2  P value 2 
-0.1806515  0.2043195  0.2067817 -0.8736341  0.3823175 -1.1455752  0.2519710 

 tau b       # ケンドールのτb（ケンドールの順位相関係数）
 ase1        # 標準誤差（H0: τb = 0 の検定のときに使う）
 ase0        # 標準誤差
 Z value     # H0: τb = 0 の検定統計量
 P value     # P 値
 Z value 2   # H0: τb = 0 の検定統計量（よく使われる検定方式によるもの）
 value P 2   # P 値

別法

分割表データから，元のデータを構成し，R にある cor.test を使う tau.b2 を定義しておく。

出力結果例

> tau.b2(f)

	Kendall's rank correlation tau

data:  x and y 
z = -1.1456, p-value = 0.2520	# R は「よく使われる検定方式」を採用していることが分かる
alternative hypothesis: true tau is not equal to 0 
sample estimates:
       tau 
-0.1806515 

Warning message: 	# この警告は気にしなくてよい
Cannot compute exact p-value with ties in: cor.test.default(x, y, method = "kendall") 

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