目的

Kaplan-Meier 法による生命表を計算する

注：R のライブラリーにもある

使用法

km.surv(time, event)

引数

time      生存期間ベクトル
event     死亡なら 1，生存なら 0 の値をとるベクトル

ソース

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

# Kaplan-Meier 法による生命表
km.surv <- function( time,                                           # 生存期間ベクトル
                        event)                                          # エンドポイント（故障なら 1，そうでなければ 0）
{
        OK <- complete.cases(time, event)                            # 欠損値を持つデータを除く
        time <- time[OK]
        event <- event[OK]
        fraction <- min(time[time > 0])/1000                         # 打ち切りと故障の生存期間が同じときでも，ソートしたときに打ち切りの方が後になるような微小値を作る
        n <- length(time)                                            # 当初データ数
        truncate <- 1-event                                          # 打ち切りデータのとき 1 になるような変数を作る
        time <- time+truncate*fraction                                       # 打切りデータの生存期間に微小値を加える
        o <- order(time)                                             # 生存期間の順序づけ
        time <- time[o]                                                      # 生存期間の並べ替え
        truncate <- truncate[o]                                              # 打ち切りかどうかのデータの並べ替え
        time <- time-truncate*fraction                                       # 生存期間データから微小値を取り除く
        i <- 1:n                                                     # 整数ベクトル
        p <- ifelse(truncate, 1, (n-i)/(n-i+1))                              # 生存確率は，打ち切りなら 1，そうでないなら (n-i)/(n-i+1)
        P <- cumprod(p)                                                      # 累積生存率は累積
        se <- (1-truncate)/(n-i+1)/(n-i)                             # 標準誤差の計算に使用する
        SE <- ifelse(truncate == 0, P*sqrt(cumsum(se)), NA)          # 標準誤差は故障時点でのみ計算される
        result <- data.frame(time, truncate, p, P, SE)                       # 結果をデータフレームとして返す
        time <- c(0, time)                                           # 生存率グラフを描くために準備
        P <- c(1, P)                                                 # 同じく
        plot(time, P, xlim=c(0, time[n+1]), ylim=c(0, 1), type="s")     # type="s" で階段状のグラフが描ける
        points(time, P)                                                 # 打ち切りか故障があったポイントをマークする
        return(result)
}


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

# 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)
 
# A 群についての解析
> a.group <- group == 1	
> km.surv(time[a.group], event[a.group])
      time truncate         p         P         SE
 [1,]    1        1 1.0000000 1.0000000         NA
 [2,]    2        0 0.9333333 0.9333333 0.06440612
 [3,]    2        0 0.9285714 0.8666667 0.08777075
 [4,]    2        1 1.0000000 0.8666667         NA
 [5,]    3        1 1.0000000 0.8666667         NA	# NA は計算しない箇所を示す
 [6,]    3        1 1.0000000 0.8666667         NA
 [7,]    3        1 1.0000000 0.8666667         NA
 [8,]    5        1 1.0000000 0.8666667         NA
 [9,]    6        1 1.0000000 0.8666667         NA
[10,]    7        0 0.8571429 0.7428571 0.13710884
[11,]   11        1 1.0000000 0.7428571         NA
[12,]   11        1 1.0000000 0.7428571         NA
[13,]   13        0 0.7500000 0.5571429 0.19089707
[14,]   15        1 1.0000000 0.5571429         NA
[15,]   17        1 1.0000000 0.5571429         NA
[16,]   20        1 1.0000000 0.5571429         NA



# B 群についての解析
> b.group <- group == 2
> km.surv(time[b.group], event[b.group])
      time truncate         p         P         SE
 [1,]    1        0 0.9285714 0.9285714 0.06883029
 [2,]    1        0 0.9230769 0.8571429 0.09352195
 [3,]    2        0 0.9166667 0.7857143 0.10966421
 [4,]    2        0 0.9090909 0.7142857 0.12073632
 [5,]    3        0 0.9000000 0.6428571 0.12806021
 [6,]    3        1 1.0000000 0.6428571         NA
 [7,]    3        1 1.0000000 0.6428571         NA
 [8,]    4        1 1.0000000 0.6428571         NA
 [9,]    5        0 0.8333333 0.5357143 0.14475776
[10,]    8        0 0.8000000 0.4285714 0.15031551
[11,]    8        1 1.0000000 0.4285714         NA
[12,]   10        1 1.0000000 0.4285714         NA
[13,]   12        0 0.5000000 0.2142857 0.16913862
[14,]   14        1 1.0000000 0.2142857         NA



 解説ページ