-

 

解析応用編

7. カプランマイヤー法(Kaplan-Meier method)によるイベント発生率の点推定とNumber at riskの算出

このサイトは無料の統計ソフトである「R」を用いて
カプランマイヤー法(Kaplan-Meier method)によるイベント発生率の点推定とNumber at riskの算出
が実行できるように説明したサイトです。
臨床医のためのRコマンダーによる医学統計解析マニュアル」の読者を対象に作成しており、本書の内容を理解しているものとして解説していきます。

その他、臨床研究や英語論文執筆にご興味のある方は無料のメールマガジン
「臨床研究の立ち上げから英語論文発表までを最速最短で行うための極意」への登録も御考慮下さい。

イベント発生率の点推定

それでは早速やってみましょう。
まずは「臨床医のためのRコマンダーによる医学統計解析マニュアル」のp102以降を参考に
練習用データを用いて糖尿病の有無による生存率の違いをカプランマイヤー法(Kaplan-Meier method)を用いて描出してみて下さい。

手順に従い実行するとRコマンダーのスクリプトウィンドウには下記のように表示されます。

.Survfit <- survfit(Surv(death_time, death_censor) ~ diabetes, conf.type="log", conf.int=0.95,
type="kaplan-meier", error="greenwood", data=practice)
.Survfit
plot(.Survfit, col=1:2, lty=1:2, conf.int=FALSE, mark.time=TRUE)
legend("bottomleft", legend=c("No","Yes"), title="diabetes", col=1:2, lty=1:2, bty="n")
quantile(.Survfit, quantiles=c(.25,.5,.75))
remove(.Survfit)

やることは一つだけです。
最後の
remove(.Survfit)のremovesummaryに書き換えて下さい。
そして再度スクリプトを実行するとRコマンダーの出力ウィンドウ下記のような結果が得られます。



もうおわかりだと思いますが、この結果をみるとある時点における推定生存率とその95%信頼区間、number at risk(=その時点で生存している患者の数)がすべてわかります。

それでは詳しく見て行きましょう。
まずこのデータはタイトルにdiabetes=Noと書いてあることから糖尿病のない患者におけるデータであることがわかります。
次に赤四角で囲った部分を左から見て行きます。
timeは時点を、n.riskはnumber at risk(=その時点で生存している患者の数)を、survivalは生存率を、lower 95% CIはその95%信頼区間下限値を、upper 95% CIは上限値を意味しています。

死亡率が知りたい場合は1から生存率を引いた値を用いればよいだけです。

Number at riskの算出

次に、任意の時点におけるNumber at riskの点推定を行いたいと思います。
残念ながらこちらはRのスクリプトを使用しなければ計算できないですが、せっかくの機会ですのでやり方をご説明したいと思います。
準備としてsurvivalパッケージを読み込んでおきます。
library(survival)

次にデータセット名は「Dataset」とし、時間情報を「Time」という変数名で連続変数で、イベント情報を「Censor」という変数名で0 or 1で、グループ情報をGroupという変数名(2グループのみに対応)でデータシートを作成し、以下をR Consoleにコピーアンドペーストして下さい。

#####コピーアンドペースト_ここから################
names(Dataset)
NoatRisk = function(dat, time, event, group, point){

mat = survfit(Surv(as.numeric(time), event=='1') ~ group, data=dat)
a = mat$time
b = mat$n.risk
c = mat$n.event
d = mat$n.censor
e = c(rep(names(mat$strata)[1], mat$strata[[1]]), rep(names(mat$strata)[2], mat$strata[[2]]))
table = data.frame(time=a, no.at.risk=b, event=c, censor=d, group=e)

Group = unique(table$group)
sub1 = subset(table, group==Group[[1]] & time>=point)
sub2 = subset(table, group==Group[[2]] & time>=point)
nar1 = head(sub1,1)[,"no.at.risk"]
nar2 = head(sub2,1)[,"no.at.risk"]
res = data.frame(Group, NoatRisk=c(nar1,nar2), Day=point)
return(res)

}
#####コピーアンドペースト_ここまで################

最後に下記の一文の最後の数字の部分を推定したいTimeの数字に置き換えて実行して下さい。
180日時点であれば

NoatRisk(Dataset, Dataset$Time, Dataset$Censor, Dataset$Group, 180)

1群しかない場合も、適当にGroupという変数を名をデータシートに作成して、A or Bといった情報を付与した上で出力された結果の数を合計すればNo at riskの点推定が可能です。

以上です。





このページの先頭へ