今回はRを使って、論文用のきれいなKaplan-Meier曲線を作成するコマンドを紹介します。最終的にはこんなグラフが出来上がります。
quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))
plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
XLIM <- 0:floor(max(m1$time) / 365.25)
XLABEL <- ifelse(XLIM%%2==0, XLIM, "")
axis(1, at=XLIM*365.25, labels=XLABEL)
axis(2, las=1)
title(main=" a", line=-1, adj=0, outer=TRUE)
title(xlab="Observational period (year)", line=2)
title(ylab="Cumulative probability of the incidence of outcome", line=2.5)
quartz.save(file="kaplan15.png", type="png")
StataでのKaplan-Meir曲線の作成方法は、「Stata: Kaplan-Meier曲線を作成する」をご覧下さい。
サンプルデータの読み込み
今回はsurvivalパッケージのサンプルデータflchainを利用します。data()コマンドを使って、flchainデータを読みこみます。さらにヘルプコマンドを使って、flchainデータの詳細を確認します。
> library(survival)
> data("flchain", package="survival"
> ?flchain
flchain package:survival R Documentation
Assay of serum free light chain for 7874 subjects.
Description:
This is a stratified random sample containing 1/2 of the subjects
from a study of the relationship between serum free light chain
(FLC) and mortality. The original sample contains samples on
approximately 2/3 of the residents of Olmsted County aged 50 or
greater.
Usage:
data(flchain)
Format:
A data frame with 7874 persons containing the following variables.
'age ' age in years
'sex' F=female, M=male
'sample.yr' the calendar year in which a blood sample was obtained
Usage:
data(flchain)
Format:
A data frame with 7874 persons containing the following variables.
'age ' age in years
'sex' F=female, M=male
'sample.yr' the calendar year in which a blood sample was obtained
'kappa' serum free light chain, kappa portion
'lambda' serum free light chain, lambda portion
'flc.grp' the FLC group for the subject, as used in the original
analysis
'creatinine' serum creatinine
'mgus' 1 if the subject had been diagnosed with monoclonal
gammapothy (MGUS)
'futime' days from enrollment until death. Note that there are 3
subjects whose sample was obtained on their death date.
'death' 0=alive at last contact date, 1=dead
'chapter' for those who died, a grouping of their primary cause of
death by chapter headings of the International Code of
Diseases ICD-9
Details:
'sample.yr' the calendar year in which a blood sample was obtained
'kappa' serum free light chain, kappa portion
'lambda' serum free light chain, lambda portion
'flc.grp' the FLC group for the subject, as used in the original
analysis
'creatinine' serum creatinine
'mgus' 1 if the subject had been diagnosed with monoclonal
gammapothy (MGUS)
'futime' days from enrollment until death. Note that there are 3
only found by chance from a test (serum electrophoresis) which is
ordered for other causes. Later work suggested that one component
of immunoglobulin production, the serum free light chain, might be
a possible marker for immune disregulation. In 2010 Dr. Angela
Dispenzieri and colleagues assayed FLC levels on those samples
from the original study for which they had patient permission and
from which sufficient material remained for further testing. They
found that elevated FLC levels were indeed associated with higher
death rates.
Patients were recruited when they came to the clinic for other
appointments, with a final random sample of those who had not yet
had a visit since the study began. An interesting side question
is whether there are differences between early, mid, and late
recruits.
This data set contains an age and sex stratified random sample
that includes 7874 of the original 15759 subjects. The original
subject identifiers and dates have been removed to protect patient
identity. Subsampling was done to further protect this
Proceedings 87:512-523.
R Kyle, T Therneau, SV Rajkumar, D Larson, M Plevak, J Offord, A
Dispenzieri, J Katzmann, and LJ Melton, III, 2006, Prevalence of
monoclonal gammopathy of undetermined significance, New England J
Medicine 354:1362-1369.
Examples:
data(flchain)
age.grp <- cut(flchain$age, c(49,54, 59,64, 69,74,79, 89, 110),
labels= paste(c(50,55,60,65,70,75,80,90),
c(54,59,64,69,74,79,89,109), sep='-'))
table(flchain$sex, age.grp)
qキーをクリックすれば、いつもの入力画面に復帰できます。とりあえず、head()、summary()、str()を使って、flchainデータセットの中身を確認します。
> head(flchain)
age sex sample.yr kappa lambda flc.grp creatinine mgus futime death
1 97 F 1997 5.70 4.860 10 1.7 0 85 1
2 92 F 2000 0.87 0.683 1 0.9 0 1281 1
3 94 F 1997 4.36 3.850 10 1.4 0 69 1
4 92 F 1996 2.42 2.220 9 1.0 0 115 1
5 93 F 1996 1.32 1.690 6 1.1 0 1039 1
6 90 F 1997 2.01 1.860 9 1.0 0 1355 1
chapter
1 Circulatory
2 Neoplasms
3 Circulatory
4 Circulatory
5 Circulatory
6 Mental
> summary(flchain)
age sex sample.yr kappa lambda
Min. : 50.00 F:4350 Min. :1995 Min. : 0.010 Min. : 0.040
1st Qu.: 55.00 M:3524 1st Qu.:1996 1st Qu.: 0.960 1st Qu.: 1.200
Median : 63.00 Median :1996 Median : 1.270 Median : 1.510
Mean : 64.29 Mean :1997 Mean : 1.431 Mean : 1.703
3rd Qu.: 72.00 3rd Qu.:1997 3rd Qu.: 1.680 3rd Qu.: 1.920
Max. :101.00 Max. :2003 Max. :20.500 Max. :26.600
flc.grp creatinine mgus futime
Min. : 1.000 Min. : 0.400 Min. :0.0000 Min. : 0
1st Qu.: 3.000 1st Qu.: 0.900 1st Qu.:0.0000 1st Qu.:2852
Median : 5.000 Median : 1.000 Median :0.0000 Median :4302
Mean : 5.471 Mean : 1.093 Mean :0.0146 Mean :3661
3rd Qu.: 8.000 3rd Qu.: 1.200 3rd Qu.:0.0000 3rd Qu.:4773
Max. :10.000 Max. :10.800 Max. :1.0000 Max. :5215
NA's :1350
death chapter
Min. :0.0000 Circulatory: 745
1st Qu.:0.0000 Neoplasms : 567
Median :0.0000 Respiratory: 245
Mean :0.2755 Mental : 144
3rd Qu.:1.0000 Nervous : 130
Max. :1.0000 (Other) : 338
NA's :5705
> str(flchain)
'data.frame': 7874 obs. of 11 variables:
$ age : num 97 92 94 92 93 90 90 90 93 91 ...
$ sex : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
$ sample.yr : num 1997 2000 1997 1996 1996 ...
$ kappa : num 5.7 0.87 4.36 2.42 1.32 2.01 0.43 2.47 1.91 0.791 ...
$ lambda : num 4.86 0.683 3.85 2.22 1.69 1.86 0.88 2.7 2.18 2.22 ...
$ flc.grp : num 10 1 10 9 6 9 1 10 9 6 ...
$ creatinine: num 1.7 0.9 1.4 1 1.1 1 0.8 1.2 1.2 0.8 ...
$ mgus : num 0 0 0 0 0 0 0 0 0 0 ...
$ futime : int 85 1281 69 115 1039 1355 2851 372 3309 1326 ...
$ death : num 1 1 1 1 1 1 1 1 1 1 ...
$ chapter : Factor w/ 16 levels "Blood","Circulatory",..: 2 13 2 2 2 11 11 14 15 2 ...
生存時間とアウトカムの設定
まずは、追跡期間とアウトカムをSurv()で設定します。アウトカム変数は、0=打ち切り、1=アウトカム発症で定義しておくのが無難です。
> s1 <- Surv(flchain$futime, flchain$death)
暴露因子としてflc.grpを設定します。
> m1 <- survfit(s1 ~ flc.grp, data=flchain)
m1オブジェクトにはいろいろな情報が詰まっています。str()で確認できます。
List of 18 $ n : int [1:10] 769 811 820 786 791 791 806 730 803 767 $ time : num [1:6191] 1 3 6 9 15 27 30 37 63 88 ... $ n.risk : num [1:6191] 769 768 767 766 765 763 762 761 760 759 ... $ n.event : num [1:6191] 0 0 0 1 2 1 0 1 0 1 ... $ n.censor : num [1:6191] 1 1 1 0 0 0 1 0 1 0 ... $ surv : num [1:6191] 1 1 1 0.999 0.996 ... $ std.err : num [1:6191] 0 0 0 0.00131 0.00227 ... $ cumhaz : num [1:6191] 0 0 0 0.00131 0.00392 ... $ std.chaz : num [1:6191] 0 0 0 0.00131 0.00226 ... $ start.time: num 0 $ strata : Named int [1:10] 530 585 609 611 596 615 651 631 672 691 ..- attr(*, "names")= chr [1:10] "flc.grp=1" "flc.grp=2" "flc.grp=3" "flc.grp=4" ... $ type : chr "right" $ logse : logi TRUE $ conf.int : num 0.95 $ conf.type : chr "log" $ lower : num [1:6191] 1 1 1 0.996 0.992 ... $ upper : num [1:6191] 1 1 1 1 1 ... $ call : language survfit(formula = s1 ~ flc.grp, data = flchain) - attr(*, "class")= chr "survfit"
累積アウトカム発症率の算出
暴露因子の各群における累積生存率を算出しましょう。研究開始から1年ごとに3年時までの累積生存率を算出します。
> summary(m1, time=seq(0, 3, 1)*365.25)
Call: survfit(formula = s1 ~ flc.grp, data = flchain)
flc.grp=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 769 0 1.000 0.00000 1.000 1.000
365 752 11 0.986 0.00431 0.977 0.994
730 744 3 0.982 0.00485 0.972 0.991
1096 730 9 0.970 0.00621 0.958 0.982
flc.grp=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 811 0 1.000 0.00000 1.000 1.000
365 797 9 0.989 0.00369 0.982 0.996
730 789 7 0.980 0.00491 0.971 0.990
1096 778 7 0.971 0.00587 0.960 0.983
flc.grp=3
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 820 0 1.000 0.00000 1.000 1.000
365 805 9 0.989 0.00365 0.982 0.996
730 791 11 0.975 0.00543 0.965 0.986
1096 778 11 0.962 0.00672 0.949 0.975
flc.grp=4
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 786 0 1.000 0.00000 1.000 1.000
365 768 14 0.982 0.00473 0.973 0.991
730 756 7 0.973 0.00578 0.962 0.985
1096 745 8 0.963 0.00677 0.950 0.976
flc.grp=5
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 791 0 1.000 0.00000 1.000 1.000
365 768 19 0.976 0.00546 0.965 0.987
730 758 5 0.970 0.00613 0.958 0.982
1096 741 13 0.953 0.00757 0.938 0.968
flc.grp=6
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 791 0 1.000 0.00000 1.000 1.000
365 774 12 0.985 0.00437 0.976 0.993
730 752 14 0.967 0.00640 0.954 0.979
1096 735 14 0.949 0.00789 0.933 0.964
flc.grp=7
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 806 0 1.000 0.00000 1.000 1.000
365 778 18 0.977 0.00524 0.967 0.988
730 758 17 0.956 0.00726 0.942 0.970
1096 740 16 0.936 0.00869 0.919 0.953
flc.grp=8
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 730 0 1.000 0.00000 1.000 1.000
365 696 27 0.963 0.00703 0.949 0.977
730 675 15 0.942 0.00870 0.925 0.959
1096 658 13 0.924 0.00988 0.905 0.943
flc.grp=9
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 803 0 1.000 0.0000 1.000 1.000
365 760 40 0.950 0.0077 0.935 0.965
730 728 31 0.911 0.0101 0.892 0.931
1096 700 27 0.877 0.0116 0.855 0.900
flc.grp=10
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 767 3 0.996 0.00225 0.992 1.000
365 651 105 0.858 0.01262 0.834 0.884
730 583 62 0.776 0.01512 0.747 0.807
1096 532 49 0.711 0.01648 0.679 0.744
各暴露群の累積発症率のデータフレームを作成してみましょう。names()コマンドでオブジェクトcp1 (cumulative probability)内のどこに累積生存率、95%信頼区間等が保存されているかを確認します。
> cp1 <- summary(m1, time=seq(0, 3, 1)*365.25) > names(cp1) [1] "n" "time" "n.risk" "n.event" [5] "n.censor" "surv" "std.err" "cumhaz" [9] "std.chaz" "start.time" "strata" "type" [13] "logse" "conf.int" "conf.type" "lower" [17] "upper" "call" "table" "rmean.endtime"
いろいろ確認すると、必要な情報は下記に保管されています。
- 累積生存率：surv
- 95％信頼区間：lower, upper
- 暴露群のカテゴリ名：strata
- 時間：time
上記の変数をまとめて、データフレームdfを作成します。
> df <- data.frame(group=cp1$strata, day=cp1$time, cp=cp1$surv, cpl=cp1$lower, cpu=cp1$upper)
> head(df)
group day cp cpl cpu
1 flc.grp=1 0.00 1.0000000 1.0000000 1.0000000
2 flc.grp=1 365.25 0.9856156 0.9772125 0.9940910
3 flc.grp=1 730.50 0.9816731 0.9722073 0.9912312
4 flc.grp=1 1095.75 0.9697498 0.9576508 0.9820017
5 flc.grp=2 0.00 1.0000000 1.0000000 1.0000000
6 flc.grp=2 365.25 0.9888582 0.9816461 0.9961233
累積生存率を累積アウトカム発症率に変換したければ、1からcp、cpl、cpuを引き算すればOKです。ただし、cplとcpuが逆になってしまうので、ご注意ください。
> df$cp <- 1 - df$cp
> df$cpl <- 1 - df$cpl
> df$cpu <- 1 - df$cpu
> head(df)
group day cp cpl cpu
1 flc.grp=1 0.00 0.00000000 0.00000000 0.000000000
2 flc.grp=1 365.25 0.01438439 0.02278754 0.005908989
3 flc.grp=1 730.50 0.01832687 0.02779273 0.008768846
4 flc.grp=1 1095.75 0.03025016 0.04234917 0.017998293
5 flc.grp=2 0.00 0.00000000 0.00000000 0.000000000
6 flc.grp=2 365.25 0.01114179 0.01835391 0.003876691
最後にdayをyearに変換します。
> df$yr <- df$day/365.25
> head(df)
group day cp cpl cpu yr
1 flc.grp=1 0.00 0.00000000 0.00000000 0.000000000 0
2 flc.grp=1 365.25 0.01438439 0.02278754 0.005908989 1
3 flc.grp=1 730.50 0.01832687 0.02779273 0.008768846 2
4 flc.grp=1 1095.75 0.03025016 0.04234917 0.017998293 3
5 flc.grp=2 0.00 0.00000000 0.00000000 0.000000000 0
6 flc.grp=2 365.25 0.01114179 0.01835391 0.003876691 1
Kaplan-Meier曲線の作成
とりあえず、デフォルト設定のままでKaplan-Meier曲線を作成します。3.7インチの正方形のデバイスの上にKM曲線を作成し、png形式で出力します。
> quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
> plot(m1)
> quartz.save(file="kaplan1.png", type="png")
累積生存率ではなく、累積アウトカム発症率のKM曲線を作成します。
> plot(m1, fun="event")
累積アウトカム発症率の最大値に計算して、小数点1桁レベルで切り上げてしまい、Y軸の最大値を調整します。累積アウトカム発症率の最大値が0.69なので、Y軸の最大値 YLIMを0.70に設定します。
> YLIM <- ceiling( (1 - min(m1$surv))*10 )/10 > plot(m1, fun="event", ylim=c(0, YLIM))
曲線がすべて黒実践なので、区別がつきません。Rでは下記の6種類のタイプを選択できます。
- lty=1: solid
- lty=2: dashed
- lty=3: dotted
- lty=4: dotdash
- lty=5: longdash
- lty=6: twodash
flc.grpのカテゴリ数は10なので、線の種類が足りません。そこで、2種類の色 (black = #000000とgray = #999999) を使うことによって、10種類の曲線を区別することにします。
%%の戻り値は割り算の余りであり (5%%3 = 2)、%/%の戻り値は商です (5%/%3 = 1) 。
> exname <- names(m1$strata) > ex <- 1:length(exname) > N <- 6 > exlty <- ex %% N > exlty <- ifelse(exlty==0, N, exlty) > > excol <- (ex - 1) %/% N > excol <- ifelse(excol==0, "#000000", excol) > excol <- ifelse(excol==1, "#999999", excol) > > plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol)
凡例の追加
凡例を追加します。
> plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol)
> legend("topleft", exname, lty=exlty, col=excol)
凡例の枠線を消します。また、凡例をrev()コマンドを使って降順ソートします。
> plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol)
> legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
軸の調整
軸の表示方法を調整します。まずは軸無しのグラフを作成します。
> plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
> legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
X軸の単位を年にします。追跡期間はm1$timeに保存されており、最大値は5215日＝ 14.3年です。X軸のラベルは0〜14にします。
> XLIM <- 0:floor(max(m1$time) / 365.25) > axis(1, at=XLIM*365.25, labels=XLIM)
Y軸を追加します。軸の数値を水平に表示します。
> axis(2, las=1)
デフォルト設定のまま、メインタイトルととxy軸のタイトルを追加するとこうなります。
title(main="a", xlab="Observational period (year)", ylab="Cumulative probability of the incidence of outcome")
メインタイトルをグラフの左上に記載するためには、figure regionの外側であるouter region (outer=TRUE) に、左揃え（adj=0) で、1行だけ内側に記載する設定が必要である。
title(main="a", line=-1, adj=0, outer=TRUE)
X軸のタイトルは、X軸から少し離れすぎですので、もう少し軸に近づけましょう。X軸から2行分だけ外向けに配置します。
title(xlab="Observational period (year)", line=2)
Y軸のラベルは長すぎるので途中で改行します。改行した場所に”\n”を挿入します。
> title(ylab="Cumulative probability of\nthe incidence of outcome", line=2.5)
余白の設定
figure regionの上下右の余白が大きいので、調整します。上下左右の余白のデフォルト値はpar()$marで確認できます。
> quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
> par()$mar
[1] 5.1 4.1 4.1 2.1
下5.1、左4.1、上4.1、右2.1に設定されています。
いろいろ微調整してみましたが、mar=c(3, 3.5, 1.5, 0.2)くらいが良さそうです。
quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))
plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
XLIM <- 0:floor(max(m1$time) / 365.25)
axis(1, at=XLIM*365.25, labels=XLIM)
axis(2, las=1)
title(main=" a", line=-1, adj=0, outer=TRUE)
title(xlab="Observational period (year)", line=2)
title(ylab="Cumulative probability of the incidence of outcome", line=2.5)
quartz.save(file="kaplan14.png", type="png")
graphics.off()
左右の余白が小さくなった分、figure regionが大きくなり、X軸のメモリの表示がおかしなことになっていますので、微調整します。偶数年だけ表示するようにX軸のラベルを設定します。2で割って余りが0ならば表示するように設定すればいいのです。
XLABEL <- ifelse(XLIM%%2==0, XLIM, "")
このXLABELをX軸のラベルに設定します。
quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))
plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
XLIM <- 0:floor(max(m1$time) / 365.25)
XLABEL <- ifelse(XLIM%%2==0, XLIM, "")
axis(1, at=XLIM*365.25, labels=XLABEL)
axis(2, las=1)
title(main=" a", line=-1, adj=0, outer=TRUE)
title(xlab="Observational period (year)", line=2)
title(ylab="Cumulative probability of the incidence of outcome", line=2.5)
quartz.save(file="kaplan15.png", type="png")
XY軸の目盛り線がちょっと長いのが気になりますが、今回はこれで完成とします。
今回のコマンドをplotKaplanMeier()にまとめてみました。
s1 <- Surv(flchain$futime, flchain$death)
m1 <- survfit(s1 ~ flc.grp, data=flchain)
quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))
MAIN <- " a"
XTITLE <- "Observational period (year)"
YTITLE <- "Cumulative probability of the incidence of outcome"
plotKaplanMeier <- function(m1, MAIN, XTITLE, YTITLE){
YLIM <- ceiling( (1 - min(m1$surv))*10 )/10
exname <- names(m1$strata)
ex <- 1:length(exname)
N <- 6
exlty <- ex %% N
exlty <- ifelse(exlty==0, N, exlty)
excol <- (ex - 1) %/% N
excol <- ifelse(excol==0, "#000000", excol)
excol <- ifelse(excol==1, "#999999", excol)
plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
XLIM <- 0:floor(max(m1$time) / 365.25)
XLABEL <- ifelse(XLIM%%2==0, XLIM, "")
axis(1, at=XLIM*365.25, labels=XLABEL)
axis(2, las=1)
title(main=MAIN, line=-1, adj=0, outer=TRUE)
title(xlab=XTITLE, line=2)
title(ylab=YTITLE, line=2.5)
}
plotKaplanMeir(m1, MAIN, XTITLE, YTITLE)
quartz.save(file=FILENAME, type="png")