目的

κ統計量を求める。
vcd パッケージに，Kappa 関数が用意されている。

使用法

kappa.stat(x, conf.level=0.95)
kappa.stat(x, w, conf.level=0.95)

引数

x	    評定結果の二次元表を行列として与える。正方行列でなくてはならない。
w	    重みを表す行列。対角成分はどんな値でも良い（0 を入れておくとよい）
	    重み付きのκ統計量を求めないときには省略すればよい。
conf.level  信頼区間を求めるときの信頼率（デフォルトでは 95% 信頼区間を求める）

ソース

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

# カッパ統計量
kappa.stat <- function(      o,                                              # 評定結果の二次元表
                        w=FALSE,                                        # 重み行列（対角成分は何でも良い）
                        conf.level=0.95)                                # 信頼率
{
        data.name <- deparse(substitute(o))
        method <- "κ統計量"
        n <- sum(o)                                                  # サンプルサイズ
        e <- outer(rowSums(o), colSums(o))/n                         # 期待値
        if (is.matrix(w) == FALSE) {                                    # 重み付けしないときは，
                po <- sum(diag(o))/n
                qo <- 1-po
                pe <- sum(diag(e))/n
                qe <- 1-pe
                kappa <- 1-qo/qe
                sk <- sqrt(po*qo/(n*qe^2))                           # κの信頼区間を計算するときに使う標準誤差
                sk0 <- sqrt(pe/(n*qe))                                       # H0：κ = 0 の検定に使うκの標準誤差
        }
        else {                                                          # 重み付けするκを求めるときは，
                method <- paste(method, "（重み付け）", sep="")
                qo <- sum(w*o)/n
                qo2 <- sum(w*w*o)/n
                qe <- sum(w*e)/n
                qe2 <- sum(w*w*e)/n
                kappa <- 1-qo/qe
                sk <- sqrt((qo2-qo^2)/n/qe^2)                                # κの信頼区間を計算するときに使う標準誤差
                sk0 <- sqrt((qe2-qe^2)/n/qe^2)                               # H0：κ = 0 の検定に使うκの標準誤差
        }
        stopifnot(kappa >= 0)
        z <- kappa/sk0                                                       # H0：κ = 0 の検定
        P <- pnorm(z, lower.tail=FALSE)*2                            # P 値
        names(kappa) <- "kappa"
        names(z) <- "Z"
        cint <- kappa+c(-sk, sk)*qnorm(0.5+conf.level/2)
        attr(cint, "conf.level") <- conf.level
        return(structure(list(estimate=kappa, statistic=z, p.value=P,   # 結果は htest クラスのオブジェクト
                method=method, data.name=data.name, conf.int=cint, "sigma-kappa"=sk,
                "sigma-kappa-0"=sk0), class="htest"))
}


使用例

x <- matrix(c(	# ファイルから読んでも良い
	12,  6,  1,
	 3, 19,  4,
	 2,  5, 34
), byrow=TRUE, ncol=3)

kappa.stat(x)

w <- matrix(c(	# 重み行列の対角部分はどんな値でも良い
	0, 1, 3,
	1, 0, 1,
	3, 1, 0
), ncol=3, byrow=TRUE)
	
kappa.stat(x, w)

出力結果例

> kappa.stat(x)	# 重みを付けないκ

	κ統計量

data:  x 
Z = 7.5202, p-value = 5.467e-14	# 検定結果
95 percent confidence interval:	# 信頼区間
 0.4721927 0.7583143 
sample estimates:	# κ統計量
    kappa 
0.6152535 

> w <- matrix(c(
+ 	0, 1, 3,
+ 	1, 0, 1,
+ 	3, 1, 0
+ ), byrow=TRUE, ncol=3)

> kappa.stat(x, w)	# 重みを付けたκ

	κ統計量（重み付け）

data:  x 
Z = 6.1563, p-value = 7.447e-10	# 検定結果
95 percent confidence interval:	# 信頼区間
 0.5586969 0.8278289 
sample estimates:	# κ統計量
    kappa 
0.6932629 

# 作りつけの関数 Kappa による推定
> library(vcd)	# vcd ライブラリーを使うために必須
> Kappa(x)
               value        ASE
Unweighted 0.6152535 0.07299153
Weighted   0.6634051 0.14035169

 解説ページ