R Lecture lesson6
Ryohei Yamamoto
2020/1/22
動画
https://web.microsoftstream.com/video/963e3e8c-6959-40aa-86d8-5345f2bc636b
資料
https://goo.gl/vRcFR8
生存時間分析 2
- はじめに（前回の復習）
- 比例ハザードモデルとは？
- ハザード？
- ハザード比？
- ハザード？
- Rstudio 上での実行
- Coxの比例ハザードモデル
- 比例ハザード性の確認
- 生存時間解析において注意が必要なポイント
- 時間依存性共変量
- Competing Risk
- Informative Censoring
はじめに
生存時間分析とは
- 打ち切りあるデータを用いてイベント発生までにかかる時間を分析すること
- Kaplan-Meier曲線
- ログランク検定
前回の復習
df.csv をRに読み込んむ
治療の有無でKaplan-Meier曲線を描く
治療の有無でLog-Rank検定を行う
#インストール
install.packages("survival", repos = "http://cran.rstudio.com/")
install.packages("survminer", repos = "http://cran.rstudio.com/")
install.packages("ggplot2", repos = "http://cran.rstudio.com/")
install.packages("tableone", repos = "http://cran.rstudio.com/")
#libraryに読み込む
library(survival)
library(ggplot2)
library(survminer)## Loading required package: ggpubr## Loading required package: magrittrdf<-read.csv("df.csv")
head(df)## Treat Duration Outcome
## 1 Active 545 Died
## 2 Active 1826 Survived
## 3 Control 1826 Survived
## 4 Control 738 Died
## 5 Active 1826 Survived
## 6 Control 670 DiedDuration, Outcome, Treatのデータセット
library(survival)
KM<-survfit(Surv(Duration,Outcome=="Died")~Treat,data=df)
plot(KM, mark.time = T) #plotで書く
ggsurvplot(KM) #survminerのggsurvplotで書く
Log-Rank test
#２つの生存曲線のログランク検定
SVdif<-survdiff(Surv(Duration, Outcome == "Died")~Treat, data=df)
SVdif## Call:
## survdiff(formula = Surv(Duration, Outcome == "Died") ~ Treat,
## data = df)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Treat=Active 484 230 293 13.6 28.2
## Treat=Control 516 336 273 14.6 28.2
##
## Chisq= 28.2 on 1 degrees of freedom, p= 1e-07では次のような場合はどうしたら良いか * 先の治療有無に加えて年齢の影響を調べたい
* 二群の間で男女比が均一でなかったので調整したい
生存時間分析の多変量解析
- 線形回帰モデル
- 連続量の被説明変数を説明変数で回帰
- ロジスティック回帰モデル
- 離散量の被説明変数を説明変数で回帰
- 回帰係数はオッズ比で解釈できる。
- Cox の比例ハザードモデル
- 生存時間の被説明変数を説明変数で回帰
- 回帰係数はハザード比で解釈できる
ハザードとは？
ハザードとは、イベントの生じる“速度”である
2群間で異なっているのは
・イベントの数？
・イベントまでの時間？
・イベントが起こる速度？
いうなれば
「瞬間イベント発生率」
ハザード比
群間のイベントの起こる速度の比
- ハザードを計算するのは難しい
- 真の生存時間関数の形がわからない
- 比べられないのでは？ - 実は、ハザードの“比”に関しては計算することができる - 説明変数が一単位増加する際の比をハザード比とする
これがCoxの比例ハザードモデル
Sir David Roxbee Cox(born 15 July 1924)が作った
RでCoxの比例ハザードモデル
library(survival)
SVcox <- coxph(Surv(Duration, Outcome=="Died")~ Treat, data=df)
summary(SVcox)## Call:
## coxph(formula = Surv(Duration, Outcome == "Died") ~ Treat, data = df)
##
## n= 1000, number of events= 566
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TreatControl 0.45156 1.57076 0.08567 5.271 1.36e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TreatControl 1.571 0.6366 1.328 1.858
##
## Concordance= 0.554 (se = 0.011 )
## Likelihood ratio test= 28.34 on 1 df, p=1e-07
## Wald test = 27.79 on 1 df, p=1e-07
## Score (logrank) test = 28.26 on 1 df, p=1e-07Treatment の Control 群は Treat に比べて、1.25倍のハザード比
1.25倍の速度でイベントが起こっている
P値 ＜ 0.01
治療は効果がありそうだ
Cox比例ハザードモデルの仮定
モデルである以上仮定がある．
Coxの比例ハザードモデルはハザード比が時間によらず一定であることが前提にある
これを、外すと妥当な結論が得られない
比例ハザード性は
__Log-log plot__や __cox.zph__で確認できる
#cox比例ハザードモデルの結果をcox.zph()
cox.zph(SVcox)## chisq df p
## Treat 0.314 1 0.58
## GLOBAL 0.314 1 0.58#カプランマイヤー曲線を描く際のオプションにfun = “cloglog”
plot(KM,fun = "cloglog")
比例ハザード性が担保できない場合有意になる
cloglogの図でみると比例ハザード性が担保できない場合ラインが平行にならない。
時間依存性共変量
比例ハザード性が満たされないのなぜか？ * 例えば、治療の効果が短期間に出てその後は効果が弱くなる
* 例えば、特定の患者属性が時間によって変化してしまう
* 時間依存性共変量
* 現状の仮定は、患者属性は観察開始からずっと一緒
* 時間とともに変化する属性はそれを明示的に処理する必要がある
Competing Risk
競合リスク
* 例えばPrimary Outcome が癌による死亡
* 他の疾患による死亡は競合リスク（心筋梗塞死亡）
* 対応するには、特殊な生存時間分析を行う必要がある
* Cumulative Incidence Curve(KM-Curve の代わり）
* Gray’s test (Log-rank test の代わり）
* Fine&Gray Method(Cox Ph Model の代わり）
Informative Censoring
打ち切りは無情報でなければならない
* 例えば、治療群で打ち切りになった人が何らかの副作用（情報）で打ち切りになっている
* 打ち切った後で、イベントが多い可能性がある
まとめ
生存時間解析は奥が深い
* まずは、基本を押さえましょう。
* 落とし穴に注意しましょう。
* 専門家に相談しましょう。
演習
・Colon2.csv を読み込んでください
・rx:治療群
・status:イベント有無
・time:観察期間
・age:年齢
・治療群ごとのカプランマイヤー曲線を描いてください
・二群のカプランマイヤー曲線をLog-Rank検定で比較してください
・tableoneパッケージを用いて群間の特性の違いを表にしてください
・Coxの比例ハザードモデルを用いてObs 群とLev+5FU群間のハザード比、年齢と性別のハザード比を求めてください
wirte KM
colon_df<-read.csv("colon_df.csv")
head(colon_df)## id study rx sex age obstruct perfor adhere nodes status differ extent
## 1 1 1 Lev+5FU 1 43 0 0 0 5 1 2 3
## 2 2 1 Lev+5FU 1 63 0 0 0 1 0 2 3
## 3 3 1 Obs 0 71 0 0 1 7 1 2 2
## 4 4 1 Lev+5FU 0 66 1 0 0 6 1 2 3
## 5 5 1 Obs 1 69 0 0 0 22 1 2 3
## 6 6 1 Lev+5FU 0 57 0 0 0 9 1 2 3
## surg node4 time etype
## 1 0 1 1521 2
## 2 0 0 3087 2
## 3 0 1 963 2
## 4 1 1 293 2
## 5 1 1 659 2
## 6 0 1 1767 2KM1<-survfit(Surv(time,status)~rx,data = colon_df)
plot(KM1, mark.time = T)
Log-Rank test
SVdif1<-survdiff(Surv(time, status)~rx, data=colon_df)
SVdif1## Call:
## survdiff(formula = Surv(time, status) ~ rx, data = colon_df)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=Lev+5FU 304 123 150 4.82 9.97
## rx=Obs 315 168 141 5.12 9.97
##
## Chisq= 10 on 1 degrees of freedom, p= 0.002tableone
library(tableone)## Warning: package 'tableone' was built under R version 3.5.3table1<-CreateTableOne(vars= c("age", "sex","obstruct","perfor","adhere","nodes","status","extent","surg","node4","time"),strata = "rx", factorVars = c("sex","obstruct","perfor","adhere","extent","surg","node4"),data =colon_df)
table1## Stratified by rx
## Lev+5FU Obs p test
## n 304 315
## age (mean (SD)) 59.70 (12.26) 59.45 (11.97) 0.800
## sex = 1 (%) 141 (46.4) 166 (52.7) 0.136
## obstruct = 1 (%) 54 (17.8) 63 (20.0) 0.543
## perfor = 1 (%) 8 ( 2.6) 9 ( 2.9) 1.000
## adhere = 1 (%) 39 (12.8) 47 (14.9) 0.525
## nodes (mean (SD)) 3.49 (3.42) 3.79 (3.73) 0.313
## status (mean (SD)) 0.40 (0.49) 0.53 (0.50) 0.001
## extent (%) 0.367
## 1 10 ( 3.3) 8 ( 2.5)
## 2 32 (10.5) 38 (12.1)
## 3 251 (82.6) 249 (79.0)
## 4 11 ( 3.6) 20 ( 6.3)
## surg = 1 (%) 76 (25.0) 91 (28.9) 0.318
## node4 = 1 (%) 79 (26.0) 87 (27.6) 0.713
## time (mean (SD)) 1798.85 (862.52) 1599.98 (854.57) 0.004Cox PH model
fit_cox<-coxph(Surv(time,status)~rx+age+sex,data=colon_df)
summary(fit_cox)## Call:
## coxph(formula = Surv(time, status) ~ rx + age + sex, data = colon_df)
##
## n= 619, number of events= 291
##
## coef exp(coef) se(coef) z Pr(>|z|)
## rxObs 0.3758713 1.4562597 0.1189343 3.160 0.00158 **
## age -0.0009998 0.9990007 0.0048801 -0.205 0.83768
## sex -0.1292110 0.8787885 0.1175153 -1.100 0.27154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## rxObs 1.4563 0.6867 1.1535 1.839
## age 0.9990 1.0010 0.9895 1.009
## sex 0.8788 1.1379 0.6980 1.106
##
## Concordance= 0.554 (se = 0.017 )
## Likelihood ratio test= 11.25 on 3 df, p=0.01
## Wald test = 11.12 on 3 df, p=0.01
## Score (logrank) test = 11.24 on 3 df, p=0.01比例ハザード性の確認
cox.zph(fit_cox)## chisq df p
## rx 1.129 1 0.29
## age 0.129 1 0.72
## sex 2.873 1 0.09
## GLOBAL 4.018 3 0.26plot(KM1,fun="cloglog")