gikoha’s blog

個人的メモがわり

RによるKaplan-Meier法

Prismで普段解析していますがGUIなだけに自動化しにくい．．．
特にグループ分けした場合それぞれを別のカラムに入力しなければいけないのがネック。
例えばgroup0-3まであった場合↓のようにexcelのマクロを使って整形しておかなければならない

そこでRの出番です

メモ代わりに．．
エクセルでデータを書く
1行目をDays event groupなどとしておく

範囲をコピー

R64を起動し

d <- read.table(pipe("pbpaste"), header=TRUE)

これで d$Days, d$event, d$groupとして参照可能になる
追加： read.csv("file.csv", fileEncoding="cp932")を使ってもよいっぽい

cmprskパッケージをインストールしていなければメニューからインストールしておく

library(cmprsk)

と入力、

 要求されたパッケージ survival をロード中です 
 要求されたパッケージ splines をロード中です

と出たら

result<-survfit(Surv(Days,event)~group, d, type="kaplan-meier")

これでKaplan-Meier解析完了

summary(result)   # 結果表示

Call: survfit(formula = Surv(Days, event) ~ group, data = d, type = "kaplan-meier")

                group=0 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
    200.0000      13.0000       1.0000       0.9231       0.0739       0.7890       1.0000 

                group=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  167     14       1    0.929  0.0688        0.803            1
  181     13       1    0.857  0.0935        0.692            1
  278     11       1    0.779  0.1129        0.587            1

                group=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
    267.0000      12.0000       1.0000       0.9167       0.0798       0.7729       1.0000 

                group=3 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     395.000        6.000        1.000        0.833        0.152        0.583        1.000

グラフだって書いてくれます

plot(result, col=1:4, xlab="Days", ylab="(%)")   # グラフ

ログランク検定

survdiff(Surv(Days,event)~group, d)

Call:
survdiff(formula = Surv(Days, event) ~ group, data = d)

         N Observed Expected (O-E)^2/E (O-E)^2/V
group=0 15        1     1.67     0.272    0.3782
group=1 19        3     1.59     1.259    1.7138
group=2 18        1     1.42     0.125    0.1667
group=3 14        1     1.32     0.076    0.0975

 Chisq= 1.7  on 3 degrees of freedom, p= 0.629

PDFへの出力

pdf(file="oge.pdf")
plot(result, lty=1:4, lwd=2, lheight=2, xlab="Days", ylab="(%)")
dev.off()