目的

Cox-Mantel 検定を行う

使用法

Cox.Mantel(group, event, time)

引数

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

ソース

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

# Cox-Mantel 検定
Cox.Mantel <- function(      group,                                  # 群を識別するベクトル（1, 2 のいずれか）
                        event,                                  # 死亡なら 1，生存なら 0 の値をとるベクトル
                        time)                                   # 生存期間ベクトル
{
        method <- "Cox-Mantel 検定"
        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]
        stopifnot(length(group) == length(event),               # 要素数が同じであること
                  length(group) == length(time))
        len <- length(group)                                 # 要素数
        tg <- table(time, group*10+event)                    # 群とエンドポイントごとに生存時間をまとめる
        k <- nrow(tg)                                                # 行数
        nia <- table(group)[1]                                       # 第 1 群のケース数
        nib <- len-nia                                               # 第 2 群のケース数
        na <- c(nia, (rep(nia, k)-cumsum(tg[,1]+tg[,2]))[-k])        # 第 1 群の総症例数
        nb <- c(nib, (rep(nib, k)-cumsum(tg[,3]+tg[,4]))[-k])        # 第 2 群の総症例数
        m <- tg[,2]+tg[,4]                                   # 両群の死亡数の合計
        s <- m > 0                                           # 死亡があった死亡期間のベクトル
        m <- m[s]                                            # 取り出す
        r <- (na+nb)[s]                                              # 各死亡期間における死亡リスク保有者数
        A <- nb[s]/r
        U <- sum(tg[,4])-sum(m*A)                            # 検定統計量
        I <- sum(m*(r-m)/(r-1)*A*(1-A))                              # 誤差分散
        Z <- U/sqrt(I)                                               # 検定統計量を標準得点に換算
        P <- pnorm(abs(Z), lower.tail=FALSE)*2                       # P 値
        return(structure(list(statistic=c(U=U, "V(U)"=I, Z=Z), p.value=P,
                method=method, data.name=data.name), class="htest"))
}


使用例

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

# 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)

> Cox.Mantel(group, event, time)

	Cox-Mantel 検定

data:  time: time, event: event, group: group 
U = 3.0298, V(U) = 2.7155, Z = 1.8386, p-value = 0.06597

 解説ページ