目的

中央値の（差の）信頼区間を求める

使用法

CI4median(x, y = NULL, conf.level = 0.95, method = c("independent.sample", "paired.sample", "one.sample"))


引数

x, y		データベクトル
		一標本の場合には y は省略する
conf.level	信頼率。省略時は 0.95 を仮定する
method		データの種別（計算方法）の指定
		independent.sample, paired.sample, one.sample のいずれかを指定する
		最初の一文字だけでよい
		省略時はindependent.sample

ソース

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

# 中央値の（差の）信頼区間を求める
CI4median <- function(       x,                                              # データベクトル
                        y = NULL,                                       # 一標本の場合には y は省略する
                        conf.level = 0.95,                              # 信頼率
                        method = c(     "independent.sample",           # 独立二標本
                                        "paired.sample",                # 対応のある標本
                                        "one.sample"))                  # 一標本
{
        stopifnot(conf.level > 0, conf.level < 1)                    # 信頼率は割合で指定する
        method <- match.arg(method)                                  # 引数の補完
        if (is.null(y)) {                                               # y が NULL なら，
                method <- "one.sample"                                       # 一標本
        }

        if (method == "independent.sample") {                           # 独立二標本（中央値の差の信頼区間）
                x <- x[!is.na(x)]                                    # 欠損値を持つケースを除く
                y <- y[!is.na(y)]                                    # 欠損値を持つケースを除く
                n1 <- length(x)                                              # サンプルサイズ
                n2 <- length(y)                                              # サンプルサイズ
                vec <- sort(as.vector(outer(x, y, "-")))             # 全ての組み合わせで引き算し，小さい順に並べる
                k <- ifelse(n1 <= 100 && n2 <= 100,                      # 分位点を計算
                        qwilcox(0.5-conf.level/2, n1, n2),
                        n1*n2/2-qnorm(0.5+conf.level/2)*sqrt(n1*n2*(n1+n2+1)/12))
                n <- n1*n2                                           # 統計量の値の総数
        }
        else {
                if (method == "paired.sample") {                        # 対応のある標本（の中央値の差の信頼区間）
                        stopifnot(is.null(y) == FALSE)                  # y にもデータがあるはず
                        OK <- complete.cases(x, y)                   # 欠損値を持つケースを除く
                        x <- x[OK]
                        y <- y[OK]
                        x <- y-x                                     # 一標本の中央値の信頼区間に帰結する
                }
                n1 <- length(x)                                              # サンプルサイズ
                mean <- outer(x, x, "+")/2                           # あらゆる組み合わせで平均値を取る
                vec <- sort(mean[upper.tri(mean, diag=TRUE)])                # 小さい順に並べる
                k <- ifelse(n1 <= 300,                                    # 分位点を計算
                        qsignrank(0.5-conf.level/2, n1),
                        n1*(n1+1)/4-qnorm(0.5+conf.level/2)*sqrt(n1*(n1+1)*(2*n1+1)/24))
                n <- n1*(n1+1)/2                                     # 統計量の値の総数
        }
        result <- c("LCL"=vec[k], "UCL"=vec[n+1-k])                  # 信頼限界値
        return(list(name = method, result = result))
}


使用例

独立二標本（各 10 ケースずつ）のデータ（ファイルから読んでも良い）
x <- c(38, 26, 29, 41, 36, 31, 32, 30, 35, 33)
y <- c(45, 28, 27, 38, 40, 42, 39, 39, 34, 45)
CI4median(x, y, method="independent")

対応のある二標本（各 10 ケースずつ）のデータ（ファイルから読んでも良い）
a <- c(10.6, 5.2, 8.4, 9.0, 6.6, 4.6, 14.1, 5.2, 4.4, 17.4, 7.2)
b <- c(14.6, 15.6, 20.2, 20.9, 24.0, 25.0, 35.2, 30.2, 30.0, 46.2, 37.0)
CI4median(a, b, method="paired")

一標本のデータ（ファイルから読んでも良い）
CI4median(b-a, method="one")


出力結果例

> x <- c(38, 26, 29, 41, 36, 31, 32, 30, 35, 33)
> y <- c(45, 28, 27, 38, 40, 42, 39, 39, 34, 45)
> CI4median(x, y, method="independent")
$name
[1] "independent.sample"

$result
LCL UCL 
-10   1 

> a <- c(10.6, 5.2, 8.4, 9.0, 6.6, 4.6, 14.1, 5.2, 4.4, 17.4, 7.2)
> b <- c(14.6, 15.6, 20.2, 20.9, 24.0, 25.0, 35.2, 30.2, 30.0, 46.2, 37.0)
> CI4median(a, b, method="paired")
$name
[1] "paired.sample"

$result
 LCL  UCL 
11.9 25.1 

> CI4median(b-a, method="one")
$name
[1] "one.sample"

$result
 LCL  UCL 
11.9 25.1 

 解説ページ