R Lecture lesson5
Ryohei Yamamoto
2020/1/22
動画
https://web.microsoftstream.com/video/eed1eb7f-822a-4d5c-a0a5-d3b20d0ecc1e
資料
https://bit.ly/2NQ7xyz
生存時間分析
今日はカプランマイヤー曲線を描くことが目標
検定の理解は統計で学んでください
生存時間分析の特徴
- 生存時間からイベントまでの時間（Time to event)での群間比較を行う
- 生存とはいうもののイベントはBinary variableであれば何でも良い
- 打ち切り（Cencoring)があってもアウトカムの評価が可能である
生存時間分析には観察期間と転帰の情報が必要
データの読み込み
生存時間分析をするにはsurvival packageが必要
データはパッケージに内装されている卵巣がんデータのovarianを使用
R markdownではInstall.package(“package_name”, repos=“http://cran.rstudio.com/”)とreposを付けないと読み込まれない
#survivalをインストール
install.packages("survival", repos="http://cran.rstudio.com/")
#libraryに読み込む
library(survival)なおデータはHelpに入れるとOvarian Cancer Survival Dataと内容が確認できる
#dataはovarianを使用する.
data(ovarian)
#ovarianを見てみる．
head(ovarian)## futime fustat age resid.ds rx ecog.ps
## 1 59 1 72.3315 2 1 1
## 2 115 1 74.4932 2 1 1
## 3 156 1 66.4658 2 1 2
## 4 421 0 53.3644 2 2 1
## 5 431 1 50.3397 2 1 1
## 6 448 0 56.4301 1 1 2str(ovarian)## 'data.frame': 26 obs. of 6 variables:
## $ futime : num 59 115 156 421 431 448 464 475 477 563 ...
## $ fustat : num 1 1 1 0 1 0 1 1 0 1 ...
## $ age : num 72.3 74.5 66.5 53.4 50.3 ...
## $ resid.ds: num 2 2 2 2 2 1 2 2 2 1 ...
## $ rx : num 1 1 1 2 1 1 2 2 1 2 ...
## $ ecog.ps : num 1 1 2 1 1 2 2 2 1 2 ...futime: survival or censoring time
fustat: censoring status
age: in years
resid.ds: residual disease present (1=no,2=yes) rx: treatment group
ecog.ps: ECOG performance status (1 is better, see reference)
カプランマイヤー曲線の書き方
カプランマイヤー曲線は以下の式でかける
fit <- survfit(Surv(futime, fustat)~rx, data=ovarian)
plot(fit)
survfit: 生存時間データの生成 Surv( , ): survival objectを生成
~ rx: どの群を比較するか，なにもないなら１
の３つをplotで生存曲線が書ける
まずはSurv()が何をしているか見てみる
Surv(ovarian$futime, ovarian$fustat)## [1] 59 115 156 421+ 431 448+ 464 475 477+ 563 638
## [12] 744+ 769+ 770+ 803+ 855+ 1040+ 1106+ 1129+ 1206+ 1227+ 268
## [23] 329 353 365 377+SurvはSurvival objcetを生成する関数
futime(フォローアップtime)とfustat(1,0のイベント)を組み合わせてベクトルを作ってくれた
このベクトルがSurvival object
Survのhelpを見てみる
Surv(time, time2, event,
type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'),
origin=0)timeは正しい打ち切りデータの場合，フォローアップ時間で，間隔データの場合，最初の引数は間隔の開始時間
time2は打ち切りされた間隔またはカウントプロセスデータのみの間隔の終了時刻.間隔は，左側が開いており．右側が閉じている（開始、終了]と見なされる．プロセスデータをカウントする場合，イベントは，間隔の終了時にイベントが発生したかどうかを示す．
今回はfutimeがすでに期間を表しているのでtime=futime, event=fustatである．
Survfitは生存曲線データを生成する関数である．
survival objectを入れるとデータが生成される． どんなデータが生成されているかsummaryでみてみる．
fit <-survfit(Surv(futime, fustat)~rx, data=ovarian)
summary(fit)## Call: survfit(formula = Surv(futime, fustat) ~ rx, data = ovarian)
##
## rx=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 59 13 1 0.923 0.0739 0.789 1.000
## 115 12 1 0.846 0.1001 0.671 1.000
## 156 11 1 0.769 0.1169 0.571 1.000
## 268 10 1 0.692 0.1280 0.482 0.995
## 329 9 1 0.615 0.1349 0.400 0.946
## 431 8 1 0.538 0.1383 0.326 0.891
## 638 5 1 0.431 0.1467 0.221 0.840
##
## rx=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 353 13 1 0.923 0.0739 0.789 1.000
## 365 12 1 0.846 0.1001 0.671 1.000
## 464 9 1 0.752 0.1256 0.542 1.000
## 475 8 1 0.658 0.1407 0.433 1.000
## 563 7 1 0.564 0.1488 0.336 0.946このデータをプロットすると生存曲線が描かれる
plot(fit)
これだとしょぼいので付け足してく
plot(fit, conf.int = T,col = c("red","blue"),xlab = "time",ylab = "event")
rxを１にすると群分けせずに全部の群での曲線を描く
fit2 <-survfit(Surv(futime, fustat)~1, data=ovarian)
plot(fit2, conf.int = T,xlab = "time",ylab = "event")
まだまだ心踊らないグラフ
美しいグラフを書くにはどうやったら良いのだろうと思った時は，以下のサイトを見る
CRAN task view: https://cran.r-project.org/web/views/
これのSurvivalを選ぶと，生存曲線に関する様々なパッケージが出てくる: https://cran.r-project.org/web/views/Survival.html
今回はggplot2を用いたggsurvplotというのが使えるsurvminerというパッケージを使ってみる．
install.packages("survminer", repos="http://cran.rstudio.com/")
install.packages("ggplot2", repos="http://cran.rstudio.com/")
library(ggplot2)
library(survminer)## Loading required package: ggpubr## Loading required package: magrittrggsurvplot(fit)
基本的にはggsurvplotのhelpを見ながらオプションを付けていく
conf.int: 信頼区間の表示
palette: 色の組み合わせ．雑誌を意識した組み合わせがある
risk.table: Nomber At riskを表示させる
ggsurvplot(fit, conf.int = T, palette = "Lancet", risk.table = T)
演習
PBCのRCTデータを使って生存分析曲線を書いてみる
pbc.csvはネット上に落ちているNEJMに掲載されたPBCに対するプラセボコントロールのRCTのData
pbc<-read.csv("pbc.csv")
head(pbc)## case time status drug age sex acitis hepatomegaly spiders edema
## 1 2 4500 0 penicillamine 56 F 0 1 1 0.0
## 2 3 1012 2 penicillamine 70 M 0 0 0 0.5
## 3 4 1925 2 penicillamine 54 F 0 1 1 0.5
## 4 5 1504 1 placebo 38 F 0 1 1 0.0
## 5 7 1832 0 placebo 55 F 0 1 0 0.0
## 6 8 2466 2 placebo 53 F 0 0 0 0.0
## bil chol alb copper_urine ALP SGOT triglicerides platelets
## 1 1.1 302 4.14 54 7394.8 113.52 88 221
## 2 1.4 176 3.48 210 516.0 96.10 55 151
## 3 1.8 244 2.54 64 6121.8 60.63 92 183
## 4 3.4 279 3.53 143 671.0 113.15 72 136
## 5 1.0 322 4.09 52 824.0 60.45 213 204
## 6 0.3 280 4.00 52 4651.2 28.38 189 373
## prothrombin Histologic_stage
## 1 10.6 3
## 2 12.0 4
## 3 10.3 4
## 4 10.9 3
## 5 9.7 3
## 6 11.0 3pbc$status<-pbc$status >=2 #死亡が2以上なので2以上のものとそうでないものの名義変数に変換
head(pbc$status)## [1] FALSE TRUE TRUE FALSE FALSE TRUEtimeが生存時間
statusが打ち切り情報 2死亡打ち切り，1肝障害打ち切り，０アウトカムなし
Drugがpenicillaminとplacebo
library(survival)
fit3 <- survfit(Surv(time, status)~drug, data=pbc)
ggsurvplot(fit3)
ggsurvplot(fit3, conf.int = T, palette = "jco", risk.table = T, pval=T, pval.method = T)