目的

ログランク検定を行う

使用法

logrank(group, event, time, method=c("SAS", "Tominaga"))

引数

gropu     群を識別するベクトル（1, 2 のいずれか）
event     死亡なら 1，生存なら 0 の値をとるベクトル
time      生存期間ベクトル
method    富永は一般的な検定法とちょっと違う

ソース

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

# ログランク検定を行う
logrank <- function( group,                                  # 群を識別するベクトル（1, 2 のいずれか）
                        event,                                  # 死亡なら 1，生存なら 0 の値をとるベクトル
                        time,                                   # 生存期間ベクトル
                        method=c("SAS", "Tominaga"))            # 一般的な SAS などの方法か，富永の方法か
{
        method <- match.arg(method)
        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]
        len <- length(group)
        stopifnot(length(event) == len, length(time) == len)
        tg <- table(c(time, rep(NA, 4)),                     # 後ろの4項目はダミー
                    c(group, 1, 1, 2, 2)*10+c(event, 1, 0, 1, 0))
        k <- nrow(tg)
        nia <- table(group)[1]
        nib <- len-nia
        na <- c(nia, (rep(nia, k)-cumsum(tg[,1]+tg[,2]))[-k])
        nb <- c(nib, (rep(nib, k)-cumsum(tg[,3]+tg[,4]))[-k])
        da <- tg[,2]
        db <- tg[,4]
        dt <- da+db
        nt <- na+nb
        d <- dt/nt
        O <- c(sum(da), sum(db))
        ea <- na*d
        eb <- nb*d
        E <- c(sum(ea), sum(eb))

        result <- data.frame(da, db, dt, na, nb, nt, d, ea, eb)

        if (method == "Tominaga") {                             # 富永による検定方式
                method <- "ログランク検定（富永）"
                chi <- sum((O-E)^2/E)
        }
        else {                                                  # SAS 等と同じ検定方式
                method <- "ログランク検定（一般的）"
                v <- sum(dt*(nt-dt)/(nt-1)*na/nt*(1-na/nt), na.rm=TRUE)
                chi <- (sum(da)-sum(na*d))^2/v
        }
        P <- pchisq(chi, 1, lower.tail=FALSE)
        return(structure(list(statistic=c("X-squared"=chi), parameter=c(df=1), p.value=P,
        method=method, data.name=data.name, result=result), class="htest"))
}


使用例

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

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

> logrank(group, event, time) # 検定結果　SAS 等と同じ検定方式によるもの

	ログランク検定（一般的）

data:  time: time, event: event, group: group 
X-squared = 3.3805, df = 1, p-value = 0.06597

> logrank(group, event, time, method="Tominaga") # 検定結果　上記，富永による検定方式によるもの

	ログランク検定（富永）

data:  time: time, event: event, group: group 
X-squared = 3.1528, df = 1, p-value = 0.0758

> logrank(group, event, time)$result

$result	# 計算表
   da db dt na nb nt          d        ea        eb
1   0  2  2 16 14 30 0.06666667 1.0666667 0.9333333
2   2  2  4 15 12 27 0.14814815 2.2222222 1.7777778
3   0  1  1 12 10 22 0.04545455 0.5454545 0.4545455
4   0  0  0  9  7 16 0.00000000 0.0000000 0.0000000
5   0  1  1  9  6 15 0.06666667 0.6000000 0.4000000
6   0  0  0  8  5 13 0.00000000 0.0000000 0.0000000
7   1  0  1  7  5 12 0.08333333 0.5833333 0.4166667
8   0  1  1  6  5 11 0.09090909 0.5454545 0.4545455
10  0  0  0  6  3  9 0.00000000 0.0000000 0.0000000
11  0  0  0  6  2  8 0.00000000 0.0000000 0.0000000
12  0  1  1  4  2  6 0.16666667 0.6666667 0.3333333
13  1  0  1  4  1  5 0.20000000 0.8000000 0.2000000
14  0  0  0  3  1  4 0.00000000 0.0000000 0.0000000
15  0  0  0  3  0  3 0.00000000 0.0000000 0.0000000
17  0  0  0  2  0  2 0.00000000 0.0000000 0.0000000
20  0  0  0  1  0  1 0.00000000 0.0000000 0.0000000

 解説ページ