Kaplan-Meyer 法による生命表     Last modified: Apr 16, 2004

目的

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

注: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-Meyer 法による生命表
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

graph

# 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 graph ・ 解説ページ


・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI

Made with Macintosh