カプランマイヤー生存曲線 を描こう（Rで）

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-17-2011

Kaplan-Meier estimates をグラフにするプログラムを書いた。Rだ。

仕事上、何度となくグラフにしていて、その都度せくせくコードを書いていたのだけど、今回、昔書いたプログラムを大幅に改善した。生存曲線を書いて、下のマージンに numbers at risk を書く。

例。セットアップから。

library(survival)
kma <- survfit( Surv(time, status) ~ rx + adhere, data=colon )

この kmaの中にグラフにするべき数字が隠されている。

僕が書いた function はkmplot()という名前で ここ にあるので、

source('https://biostat.app.vumc.org/wiki/pub/Main/TatsukiRcode/RFunctions1.R')

で読み込める。

例の続き：

シンプルにkmplot(kma)で、それなりのグラフは得られるんだけど、細かにいろいろと設定することが可能。

kmplot(kma, mark='',
   xaxis.at=c(0,.5,1:9)*365, xaxis.lab=c(0,.5,1:9),
   lty.surv=c(1,2), lwd.surv=1, col.surv=c(1,1,2,2,4,4),
   col.ci=0,
   group.names=c('Obs ','Obs tumor adh','Lev','Lev tumor adh',
    'Lev+5FU ','Lev+5FU tumor adh'),
   group.order=c(5,3,1,6,4,2),
   extra.left.margin=6, label.n.at.risk=FALSE, draw.lines=TRUE,
   cex.axis=.8, xlab='Years', ylab='Survival Probability',
   grid=TRUE, lty.grid=1, lwd.grid=1, col.grid=grey(.9),
   legend=TRUE, loc.legend='bottomleft',
   cex.lab=.8, xaxs='r', bty='L', las=1, tcl=-.2
)

これでこの図が描ける。

唯一絶対必要なインプットはsurvfit()で作ったオブジェクト（この例ではkma）だけ。それ以降のインプットは・・・

まずmarkは、打ち切りデータをどう表示するか。特に気にしなくて良い。いらない場合はmark=''。

次にxaxis.at。これはn at riskを計算して表示する x軸の値（つまり時）を指定する。何もしていないと適当な所になる。xaxis.labはxaxis.atで表示するラベル。例えば、時間は日にちで指定されているけれど、年を表示したいという時は、xaxis.at=(0:5)*365, xaxis.lab=0:5とすれば、０日、３６５日、７３０日のところに、０，１，２と表示することが出来る。

lty.surv, lwd.surv, col.survは、それぞれ線の種類、太さ、色だ。グループが４つあれば、４つ指定することもできるけれど、４つ未満の場合は繰り返して使われる。Rではいつものことだ。lty.surv=1:2, lty.lwd=1, lty.col=1:4とか。

lty.ci, lwd.ci, col.ciは confidence interval を表示する時に使う。表示したくなければ、lty.ci=0あるいはcol.ci=0とする。

group.namesは別に指定しなくていいんだけど、その場合、kmaオブジェクトで使われているものをそのまま引用する。それだと、たいていの場合うざい。
group.names=c('trt1','trt2','trt3','trt4')とかに書き換えてしまえば良い。 group.orderを指定しないと、ま、そのままグラフになるのだけど、n at riskの表示と、上のグラフがだいたい同じ順序だと見やすい。そこで、表示の順番を変えられる。例えば、group.order=c(2,4,1,3)とすると２番目のグループを１番上に書いて、その次が４番目のグループ、その下に、１番目のグループ、３番目のグループとなる。ま、例を見てください。

group.namesを短く指定してもそれでも見切れる場合はextra.left.marginを指定する。label.n.at.risk=TRUEとすると"Numbers at risk"と表示する。必要ないかな。

grid=TRUEで座標を表示する。細かな設定はcol.gridなどで。

legend=TRUEにすると、グループの名前と線の色、種類なんかを表示できる。場所はloc.legendで変えられる。

その他のインプットはplot()に行くので、cex.label、btyなども指定できる。

図だけで、下の数字はいらないや、という場合はsimple=TRUEで消せる。

他にも ggplot を使って似たようなことをした人は何人かいるけれど、ggplot なんか使わないで済むならそれにこしたことはない、と思う。

ブラヴォー！素晴らしいです

サンキュー！ありがとうございます

No. at riskをggplot2を使わずに作成する方法を探しており、このページにたどり着きました。
コードを拝見したいと思ったのですが、リンク切れのようでした。
可能でしたらコードをご教示いただければ幸いです。

リンクを訂正しておきました。リンク先の RFunctions1.R の中ほどに隠れています。

９年前に書いたコードですが、最終アップデートは2018年です。ggplot のもいいんですが、自分でゴテゴテ書くと細かくカスタマイズできる利点はあります。

コメントの投稿