目的

一般化 Wilcoxon 検定を行う

使用法

Gen.Wil(group, event, time)

引数

group     群を識別するベクトル（1, 2 のいずれか）
event     死亡なら 1，生存なら 0 の値をとるベクトル
time      生存期間ベクトル

ソース

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

# 一般化 Wilcoxon 検定を行う
Gen.Wil <- function( group,                                                  # 群を識別するベクトル
                        event,                                                  # 死亡なら 1，生存なら 0
                        time)                                                   # 生存期間ベクトル
{
        getU <- function(tx, sx, ty, sy)                                     # U の計算
        {
                if ((tx < ty && sx == 1 && sy == 1) | (tx <= ty && sx == 1 && sy == 0)) {
                        return(-1)
                }
                else if ((tx > ty && sx == 1 && sy == 1) | (tx >= ty && sx == 0 && sy == 1)) {
                        return(1)
                }
                else {
                        return(0)
                }
        }
        method <- "一般化 Wilcoxon 検定"
        data.name <- sprintf("time: %s, event: %s, group: %s",
                deparse(substitute(time)), deparse(substitute(event)), deparse(substitute(group)))
        OK <- complete.cases(group, event, time)                             # 欠損値を持つケースを除く
        group <- group[OK]
        event <- event[OK]
        time <- time[OK]
        n <- length(group)
        o <- order(group)                                                    # グループによって並べ替え
        time <- time[o]
        group <- group[o]
        event <- event[o]
        na <- table(group)[1]                                                        # 各群のデータ個数
        nb <- n-na                                                           # 各群のデータ個数

        W <- 0                                                                       # 検定統計量
        for (i in 1:na) {
                for (j in (na+1):n) {
                        W <- W+getU(time[i], event[i], time[j], event[j])
                }
        }

        Var.W <- 0                                                           # 分散
        for (i in 1:n) {
                temp <- 0
                for (j in 1:n) {
                        temp <- temp+getU(time[i], event[i], time[j], event[j])
                }
                Var.W <- Var.W+temp^2
        }
        Var.W <- na*nb*Var.W/n/(n-1)
        Z <- W/sqrt(Var.W)                                                   # Z 値
        P <- pnorm(abs(Z), lower.tail=FALSE)*2                                       # P 値
        return(structure(list(statistic=c(W=W, "Var(W)"=Var.W, Z=Z), p.value=P,
                method=method, data.name=data.name), class="htest"))
}


使用例

# 富永祐民「治療効果判定のための実用統計学 − 生命表法の解説 −」蟹書房(第3回改訂版）74 ページ，表 4.1 のデータ
# 解析結果は 107-109 ページ

# 1 は A 群，2 は B 群を表す
> group <- c(1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1)

# 1 は死亡，0 は 生存（打ち切り）を表す
> event <- c(1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0)

# 生存期間
> time <- c(2, 20, 5, 1, 3, 17, 2, 3, 15, 14, 12, 13, 11, 11, 10, 8, 8, 3, 7, 3, 6, 2, 5, 4, 2, 3, 1, 3, 2, 1)

> Gen.Wil(group, event, time)

	一般化 Wilcoxon 検定

data:  time: time, event: event, group: group 
W = 63.0000, Var(W) = 1320.3126, Z = 1.7338, p-value = 0.08295

 解説ページ