二種類の生存曲線を混ぜる。分子標的薬と細胞障害性抗癌剤。IPASS試験をRで科学する

二種類の生存曲線が混ざった状況だとどうなるか。分子標的薬と細胞障害性抗癌剤。

生存曲線の比較にはCOX回帰などが行われる。期間を通じて比例ハザード性が保たれていないと予期できない結果になる。

有名なのがIPASS試験の論文である。

http://medbio.utoronto.ca/students/courses/mbp1018/MOK_iPASS_nejm_M-S_Tsao.pdf

この例を説明しているのを聞く機会がこれまで２回あったが、いずれも釈然としない説明だった。ならば、Rで分析しようではないか。

A 分子標的薬群/遺伝子変異あり

B 分子標的薬群/遺伝子変異なし

C 通常の抗癌剤群/遺伝子変異あり

D 通常の抗癌剤群/遺伝子変異なし

前提は

遺伝子変異があると、ない人に比べて治療にかかわらず比較的予後良好。

遺伝子変異があると、分子標的薬は良く効くので予後良好

遺伝子変異がないと、分子標的薬はほぼ効かず予後不良

通常の抗癌剤は遺伝子変異にかかわらずそこそこ効く

といった感じだ。

上記論文の生存曲線から、見た目で各パラメータを決めた。

なお、試行錯誤は全く行わず、一発目のパラメータ設定での解析であることを明記しておく。

＃日本語が混じってもR_markdownは機能した。

＃ただし error=FALSE,warning=FALSEとしないと延々警告が出力されて残念

## 生存曲線の重ねあわせ

library(survival)

## Loading required package: splines

set.seed(1)
N <- 20000  #全症例
NAK <- N%/%4  #count N
NBK <- N%/%4
NCK <- N%/%4
NDK <- N - NAK * 3

AK <- 0.55  #survival rate at AKT
AKT <- 8  #x Months survive
BK <- 0.2
BKT <- 4
CK <- 0.5
CKT <- 6
DK <- 0.4
DKT <- 6

tinib <- c(rep(1, NAK + NBK), rep(0, NCK + NDK))
mut <- c(rep(1, NAK), rep(0, NBK), rep(1, NCK), rep(0, NDK))
data <- data.frame(tinib = tinib, mut = mut)

data$event <- rep(1, N)

AKk <- -log(AK)/AKT
BKk <- -log(BK)/BKT
CKk <- -log(CK)/CKT
DKk <- -log(DK)/DKT

AK <- rexp(NAK, rate = AKk)  #指数分布で
BK <- rexp(NBK, rate = BKk)
CK <- rexp(NCK, rate = CKk)
DK <- rexp(NDK, rate = DKk)

data$sur <- c(AK, BK, CK, DK)

## 分子標的薬ありなし。遺伝子変異考慮なし
all <- survfit(Surv(sur, event) ~ 1, data = data, conf.type = "none")
tinib_1 <- survfit(Surv(sur, event) ~ 1, subset(data, data$tinib == 1), conf.type = "log-log")
tinib_0 <- survfit(Surv(sur, event) ~ 1, subset(data, data$tinib == 0), conf.type = "log-log")

plot(all, xlim = c(0, 50), ylim = c(0, 1), col = "black")
par(new = TRUE)
plot(tinib_0, xlim = c(0, 50), ylim = c(0, 1), conf.int = TRUE, col = "red")
par(new = TRUE)
plot(tinib_1, xlim = c(0, 50), ylim = c(0, 1), conf.int = TRUE, col = "purple")
par(new = TRUE)
title(main = "遺伝子変異考慮せず。黒線は全員の曲線")
legend(25, 1, legend = c("All", "分子標的薬なし", "分子標的薬あり"), col = c("black", 
    "red", "purple"), lty = 1)

凄まじく面白い結果だ。ボランティアで無償で研修医に講義をしているが、こいつは使えるぜ。とはいっても、研修医にはちょっとマニアすぎるか。

９５％信頼区間は実線の上下の点線だ。２万人で全員完全追跡という、現実世界では実現不可能な超理想データなので信頼区間はめちゃくちゃ狭い。

生存曲線は、はじめは完全に分離、途中でクロスし、またその後は完全に分離する。ここまで見事に分離していても、何とP=0.58！！！！有意差なし。

曲線が異なっていても、比例ハザード性が崩れていると、今回の場合はどんなに非現実的なレベルにNを増やしても有意差は出ないようだ。

次に遺伝子変異などを考慮して解析する。

res1 <- coxph(formula = Surv(sur, event) ~ tinib, data = data, method = "breslow")
summary(res1)

## Call:
## coxph(formula = Surv(sur, event) ~ tinib, data = data, method = "breslow")
## 
##   n= 20000, number of events= 20000 
## 
##          coef exp(coef) se(coef)    z Pr(>|z|)
## tinib 0.00796   1.00799  0.01437 0.55     0.58
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## tinib      1.01      0.992      0.98      1.04
## 
## Concordance= 0.53  (se = 0.002 )
## Rsquare= 0   (max possible= 1 )
## Likelihood ratio test= 0.31  on 1 df,   p=0.58
## Wald test            = 0.31  on 1 df,   p=0.58
## Score (logrank) test = 0.31  on 1 df,   p=0.58


# 遺伝子変異考慮
all <- survfit(Surv(sur, event) ~ 1, data = data, conf.type = "none")
t1m0 <- survfit(Surv(sur, event) ~ 1, subset(data, data$tinib == 1 & data$mut == 
    0), conf.type = "log-log")
t0m0 <- survfit(Surv(sur, event) ~ 1, subset(data, data$tinib == 0 & data$mut == 
    0), conf.type = "log-log")
t1m1 <- survfit(Surv(sur, event) ~ 1, subset(data, data$tinib == 1 & data$mut == 
    1), conf.type = "log-log")
t0m1 <- survfit(Surv(sur, event) ~ 1, subset(data, data$tinib == 0 & data$mut == 
    1), conf.type = "log-log")

# 遺伝子変異なし
par(new = FALSE)
plot(all, xlim = c(0, 50), ylim = c(0, 1), col = "black")
par(new = TRUE)
plot(t0m0, xlim = c(0, 50), ylim = c(0, 1), conf.int = TRUE, col = "red")
par(new = TRUE)
plot(t1m0, xlim = c(0, 50), ylim = c(0, 1), conf.int = TRUE, col = "purple")
par(new = TRUE)
title(main = "遺伝子変異なし")
legend(25, 1, legend = c("All", "分子標的薬なし", "分子標的薬あり"), col = c("black", 
    "red", "purple"), lty = 1)

遺伝子変異がないので、分子標的薬の曲線は悲惨である。曲線はクロスせずP<0.001であり、激しく有意差ありだ。分子標的薬群が悪い。

res2 <- coxph(formula = Surv(sur, event) ~ tinib, subset(data, data$mut == 0), 
    method = "breslow")
summary(res2)

## Call:
## coxph(formula = Surv(sur, event) ~ tinib, data = subset(data, 
##     data$mut == 0), method = "breslow")
## 
##   n= 10000, number of events= 10000 
## 
##        coef exp(coef) se(coef)    z Pr(>|z|)    
## tinib 0.980     2.665    0.022 44.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## tinib      2.67      0.375      2.55      2.78
## 
## Concordance= 0.614  (se = 0.003 )
## Rsquare= 0.18   (max possible= 1 )
## Likelihood ratio test= 1986  on 1 df,   p=0
## Wald test            = 1992  on 1 df,   p=0
## Score (logrank) test = 2116  on 1 df,   p=0


# 遺伝子変異あり
par(new = FALSE)
plot(all, xlim = c(0, 50), ylim = c(0, 1), col = "black")
par(new = TRUE)
plot(t0m1, xlim = c(0, 50), ylim = c(0, 1), conf.int = TRUE, col = "red")
par(new = TRUE)
plot(t1m1, xlim = c(0, 50), ylim = c(0, 1), conf.int = TRUE, col = "purple")
par(new = TRUE)
title(main = "遺伝子変異あり")
legend(25, 1, legend = c("All", "分子標的薬なし", "分子標的薬あり"), col = c("black", 
    "red", "purple"), lty = 1)

これは遺伝子変異があるので分子標的薬が良好である。変異があると基本予後良好なので、どちらも黒線（全員統合）を上回っている。これも曲線はクロスしない。

res3 <- coxph(formula = Surv(sur, event) ~ tinib, subset(data, data$mut == 1), 
    method = "breslow")
summary(res3)

## Call:
## coxph(formula = Surv(sur, event) ~ tinib, data = subset(data, 
##     data$mut == 1), method = "breslow")
## 
##   n= 10000, number of events= 10000 
## 
##          coef exp(coef) se(coef)     z Pr(>|z|)    
## tinib -0.4442    0.6414   0.0205 -21.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## tinib     0.641       1.56     0.616     0.668
## 
## Concordance= 0.552  (se = 0.003 )
## Rsquare= 0.046   (max possible= 1 )
## Likelihood ratio test= 466  on 1 df,   p=0
## Wald test            = 468  on 1 df,   p=0
## Score (logrank) test = 475  on 1 df,   p=0

p<0.001で激しく有意差ありだ。

今度は多変量解析だ。

# すべて
res4 <- coxph(formula = Surv(sur, event) ~ tinib + mut, data = data, method = "breslow")
summary(res4)

## Call:
## coxph(formula = Surv(sur, event) ~ tinib + mut, data = data, 
##     method = "breslow")
## 
##   n= 20000, number of events= 20000 
## 
##          coef exp(coef) se(coef)     z Pr(>|z|)    
## tinib  0.2521    1.2867   0.0151  16.6   <2e-16 ***
## mut   -0.9233    0.3972   0.0159 -58.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## tinib     1.287      0.777     1.249      1.33
## mut       0.397      2.518     0.385      0.41
## 
## Concordance= 0.625  (se = 0.002 )
## Rsquare= 0.156   (max possible= 1 )
## Likelihood ratio test= 3399  on 2 df,   p=0
## Wald test            = 3373  on 2 df,   p=0
## Score (logrank) test = 3543  on 2 df,   p=0

分子標的薬のみだと有意差はなかったが、遺伝子変異と両方投入するとどちらも激しく有意P<0.001だ。

しかし、これだと分子標的薬は単純に予後不良ってことになる。これは実態を表していない。

# 交互作用あり
res5 <- coxph(formula = Surv(sur, event) ~ tinib + mut + tinib:mut, data = data, 
    method = "breslow")
summary(res5)

## Call:
## coxph(formula = Surv(sur, event) ~ tinib + mut + tinib:mut, data = data, 
##     method = "breslow")
## 
##   n= 20000, number of events= 20000 
## 
##              coef exp(coef) se(coef)     z Pr(>|z|)    
## tinib      0.9773    2.6572   0.0209  46.8   <2e-16 ***
## mut       -0.2903    0.7481   0.0201 -14.4   <2e-16 ***
## tinib:mut -1.4222    0.2412   0.0296 -48.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## tinib         2.657      0.376     2.551     2.768
## mut           0.748      1.337     0.719     0.778
## tinib:mut     0.241      4.146     0.228     0.256
## 
## Concordance= 0.651  (se = 0.002 )
## Rsquare= 0.247   (max possible= 1 )
## Likelihood ratio test= 5684  on 3 df,   p=0
## Wald test            = 6110  on 3 df,   p=0
## Score (logrank) test = 6899  on 3 df,   p=0

交互作用も激しく有意差ありだ。

交互作用を入れてはじめて、分子標的薬の有用性がわかる。

つまり、遺伝子変異あり、かつ、分子標的薬使用アリで予後良好。

つまり、遺伝子変異なし、かつ、分子標的薬使用アリで予後不良とわかる。

最後に比例ハザード性を確認する。

## 比例ハザード性確認
cox.zph(res1)

##          rho chisq p
## tinib -0.198   810 0

cox.zph(res2)

##            rho chisq     p
## tinib -0.00607 0.366 0.545

cox.zph(res3)

##           rho chisq      p
## tinib -0.0172  2.97 0.0847

cox.zph(res4)

##            rho chisq p
## tinib  -0.1433   498 0
## mut     0.0845   167 0
## GLOBAL      NA   537 0

cox.zph(res5)

##                rho  chisq     p
## tinib     -0.00477 0.4531 0.501
## mut        0.00197 0.0776 0.781
## tinib:mut -0.00488 0.4782 0.489
## GLOBAL          NA 3.4239 0.331

res1 遺伝子変異を考慮しないと、激しく比例ハザード性が成り立たない。P<0.001 res2,res3 遺伝子変異で層別化すると、比例ハザード性は成り立たないとはいえない。OK

res4 交互作用を入れない多変量解析だと、激しく比例ハザード性が成り立たない。P<0.001

res5 交互作用を入れた多変量解析だと、比例ハザード性が成り立たないとは言えない。OK

と、いうわけでres5および遺伝子変異に基づく層別解析が有用であろう。

もし、実地臨床で遺伝子変異を計測できれば、変異があるときに分子標的薬を使うのが良さそう。

もし、実地臨床で遺伝子変異が計測できないなら、クロスする生存曲線で考えるしか無い。その場合、どちらが良いとも言えない。初期の予後を優先するか、長期の予後を優先するか、副作用なども考慮しつつの価値観の問題となるであろう。

今回のトリビア：生存曲線がクロスする場合、どんなに症例集めても有意差出ないかもよ。

r-statistics-fanの日記

統計好き人間の覚書のようなもの

二種類の生存曲線が混ざったら。分子標的薬と細胞障害性抗癌剤。IPASS試験のような状況をRで統計解析する

二種類の生存曲線を混ぜる。分子標的薬と細胞障害性抗癌剤。IPASS試験をRで科学する