r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

生存曲線下面積RMST(Restricted mean survival time)~比例ハザード性が仮定できない生存曲線での代替手法

生存曲線下面積RMST(Restricted mean survival time)というのを聞いた。 論文の多くは田舎病院では入手できなかったが、下記は読めた。

Royston, P. & Parmar, M.K., 2013. Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome. BMC medical research methodology, 13(1), p.152.

Uno, H.他, 2014. Moving Beyond the Hazard Ratio in Quantifying the Between-Group Difference in Survival Analysis. Journal of Clinical Oncology, 32(22), p.2380-2385.

Tian, L., Zhao, L. & Wei, L.J., 2013. On the restricted mean event time in survival analysis. Available at: http://biostats.bepress.com/cgi/viewcontent.cgi?article=1164&context=harvardbiostat

 

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

超有名なIPASS試験の論文について以前の記事で考察した。 

その他、早期には効いているが、後のほうでは曲線が重なる場合など 比例ハザード性が仮定できない実例が散見される。

今回の手法は、比例ハザード性は仮定せず、時間を区切って その時点までのカプランマイヤー生存曲線下面積で、2群の 差を検定しようというものである。 これは、比例ハザード性の仮定を要しない。

とりあえず、Rでやってみる。 以前の記事で作ったモデルが丁度はじめの方は悪く 後のほうで逆転する感じなのでそのまんま流用する。 本当はハザードが変わっていくモデルを作ればいい のだろうが、今回は流用で勘弁。

と思ったが、20000例だと計算が終わらなかった。 そのままではなく

1000例に減らす

library(survival)
set.seed(1)
N <- 1000  #全症例
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)
data$cov <- 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 = "none")
tinib_0 <- survfit(Surv(sur, event) ~ 1, subset(data, data$tinib == 0), conf.type = "none")

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)

 

f:id:r-statistics-fan:20140804224935p:plain

とりあえず、12ヶ月までのRMSTで検定してみよう。 regression 又は augmentationで、ベースラインのcovariateの調整ができるようだ。 ここではaugmentationを使用することにする。 covariateの調整なしでまずは解析したかったが、入力しないとエラーになる。 とりあえず全部1で投入する。

library(surv2sampleComp)
# 調整なしRMST
res.RMST <- rmstaug(y = data$sur, x = data$cov, arm = data$tinib, delta = data$event, 
    tau = 12, type = "difference")
res.RMST
## $result.ini
##             coef se(coef)     z         p lower .95 upper .95
## intercept  6.450   0.1892 34.08 0.000e+00     6.079    6.8205
## arm       -1.103   0.2752 -4.01 6.068e-05    -1.643   -0.5641
## 
## $result.aug
##             coef se(coef)      z        p lower .95 upper .95
## intercept  6.450   0.1893 34.064 0.00e+00     6.078    6.8207
## arm       -1.103   0.2753 -4.008 6.12e-05    -1.643   -0.5638

##Cox回帰
## Call: ## coxph(formula = Surv(sur, event) ~ tinib, data = data, method = "breslow") ## ## n= 1000, number of events= 1000 ## ## coef exp(coef) se(coef) z Pr(>|z|) ## tinib 0.1003 1.1055 0.0636 1.58 0.11 ## ## exp(coef) exp(-coef) lower .95 upper .95 ## tinib 1.11 0.905 0.976 1.25 ## ## Concordance= 0.537 (se = 0.009 ) ## Rsquare= 0.002 (max possible= 1 ) ## Likelihood ratio test= 2.49 on 1 df, p=0.115 ## Wald test = 2.49 on 1 df, p=0.115 ## Score (logrank) test = 2.49 on 1 df, p=0.115

###RMSTの各種解析
samp1<-surv2sample(time=data$sur, arm=data$tinib, status=data$event, tau=12)
samp1

samp1 $procedure
[1] "KM" $method [1] "Perturbation resampling" $survfit Call: survfit(formula = Surv(time, status) ~ arm) records n.max n.start events median 0.95LCL 0.95UCL arm=0 500 500 500 500 6.01 5.12 6.81 arm=1 500 500 500 500 3.75 3.15 4.41 $tau [1] 12 $npert [1] 1000 $timepoints [1] 12 24 36 40 $quanprobs [1] 0.10 0.15 0.20 $contrast.all Est. Lower 95% Upper 95% p-val RMST Group0-Group1 1.1034 0.56648 1.64037 5.63e-05 RMST Group1-Group0 -1.1034 -1.64037 -0.56648 5.63e-05 RMST Group0/Group1 1.2064 1.09836 1.32506 8.86e-05 RMST Group1/Group0 0.8289 0.75468 0.91045 8.86e-05 Loss time Group0-Group1 -1.1034 -1.64037 -0.56648 5.63e-05 Loss time Group1-Group0 1.1034 0.56648 1.64037 5.63e-05 Loss time Group0/Group1 0.8342 0.76338 0.91152 6.13e-05 Loss time Group1/Group0 1.1988 1.09707 1.30996 6.13e-05 Prob at 12 Group0-Group1 0.0020 -0.04798 0.05198 9.37e-01 Prob at 12 Group1-Group0 -0.0020 -0.05198 0.04798 9.37e-01 Prob at 12 Group0/Group1 1.0088 0.80731 1.26070 9.38e-01 Prob at 12 Group1/Group0 0.9912 0.79321 1.23868 9.38e-01 Prob at 24 Group0-Group1 -0.0240 -0.05535 0.00735 1.33e-01 Prob at 24 Group1-Group0 0.0240 -0.00735 0.05535 1.33e-01 Prob at 24 Group0/Group1 0.6923 0.43012 1.11432 1.30e-01 Prob at 24 Group1/Group0 1.4444 0.89741 2.32495 1.30e-01 Prob at 36 Group0-Group1 -0.0060 -0.02429 0.01229 5.20e-01 Prob at 36 Group1-Group0 0.0060 -0.01229 0.02429 5.20e-01 Prob at 36 Group0/Group1 0.7692 0.34354 1.72242 5.24e-01 Prob at 36 Group1/Group0 1.3000 0.58058 2.91089 5.24e-01 Prob at 40 Group0-Group1 -0.0020 -0.01493 0.01093 7.62e-01 Prob at 40 Group1-Group0 0.0020 -0.01093 0.01493 7.62e-01 Prob at 40 Group0/Group1 0.8333 0.25026 2.77490 7.66e-01 Prob at 40 Group1/Group0 1.2000 0.36037 3.99586 7.66e-01 Quantile at 10 % Group0-Group1 0.3051 -0.04546 0.65572 8.80e-02 Quantile at 10 % Group1-Group0 -0.3051 -0.65572 0.04546 8.80e-02 Quantile at 10 % Group0/Group1 1.6020 0.98627 2.60228 5.69e-02 Quantile at 10 % Group1/Group0 0.6242 0.38428 1.01393 5.69e-02 Quantile at 15 % Group0-Group1 0.6318 0.27800 0.98557 4.65e-04 Quantile at 15 % Group1-Group0 -0.6318 -0.98557 -0.27800 4.65e-04 Quantile at 15 % Group0/Group1 1.8615 1.37654 2.51741 5.45e-05 Quantile at 15 % Group1/Group0 0.5372 0.39723 0.72646 5.45e-05 Quantile at 20 % Group0-Group1 0.9475 0.41825 1.47683 4.50e-04 Quantile at 20 % Group1-Group0 -0.9475 -1.47683 -0.41825 4.50e-04 Quantile at 20 % Group0/Group1 1.9357 1.40877 2.65967 4.62e-05 Quantile at 20 % Group1/Group0 0.5166 0.37599 0.70984 4.62e-05 Ave of t-year event rates Group0-Group1 -0.0075 -0.02933 0.01433 5.01e-01 Ave of t-year event rates Group1-Group0 0.0075 -0.01433 0.02933 5.01e-01 Ave of t-year event rates Group0/Group1 0.9123 0.69844 1.19159 5.01e-01 Ave of t-year event rates Group1/Group0 1.0962 0.83921 1.43176 5.01e-01 Ave percentiles Group0-Group1 0.6282 0.25163 1.00468 1.08e-03 Ave percentiles Group1-Group0 -0.6282 -1.00468 -0.25163 1.08e-03 Ave percentiles Group0/Group1 1.8365 1.33737 2.52188 1.72e-04 Ave percentiles Group1/Group0 0.5445 0.39653 0.74774 1.72e-04 $group0 Est. Lower 95% Upper 95% SE RMST 6.450 6.10123 6.8196 0.18308 Loss time 5.550 5.18040 5.8988 0.18308 Prob at 12 0.228 0.19465 0.2644 0.01781 Prob at 24 0.054 0.03642 0.0770 0.01008 Prob at 36 0.020 0.01015 0.0330 0.00610 Prob at 40 0.010 0.00328 0.0203 0.00456 Quantile at 10 % 0.812 0.53226 1.1527 0.16614 Quantile at 15 % 1.365 1.08877 1.7279 0.16666 Quantile at 20 % 1.960 1.51232 2.4583 0.24890 Ave of t-year event rates 0.078 0.06424 0.0941 0.00738 Ave percentiles 1.379 1.05040 1.7363 0.17771 $group1 Est. Lower 95% Upper 95% SE RMST 5.3462 4.96701 5.7515 0.20606 Loss time 6.6538 6.24851 7.0330 0.20606 Prob at 12 0.2260 0.19073 0.2633 0.01867 Prob at 24 0.0780 0.05584 0.1038 0.01212 Prob at 36 0.0260 0.01368 0.0417 0.00717 Prob at 40 0.0120 0.00409 0.0225 0.00493 Quantile at 10 % 0.5068 0.34358 0.6405 0.07058 Quantile at 15 % 0.7333 0.62255 0.8777 0.06883 Quantile at 20 % 1.0127 0.81396 1.1796 0.10206 Ave of t-year event rates 0.0855 0.06902 0.1029 0.00838 Ave percentiles 0.7509 0.61131 0.8847 0.07349 $RMST Est. Lower 95% Upper 95% p-val RMST Group0-Group1 1.103 0.566 1.640 5.63e-05 RMST Group1-Group0 -1.103 -1.640 -0.566 5.63e-05 RMST Group0/Group1 1.206 1.098 1.325 8.86e-05 RMST Group1/Group0 0.829 0.755 0.910 8.86e-05 $RMLT Est. Lower 95% Upper 95% p-val Loss time Group0-Group1 -1.103 -1.640 -0.566 5.63e-05 Loss time Group1-Group0 1.103 0.566 1.640 5.63e-05 Loss time Group0/Group1 0.834 0.763 0.912 6.13e-05 Loss time Group1/Group0 1.199 1.097 1.310 6.13e-05 $contrast.diff10 Est. Lower 95% Upper 95% p-val RMST Group1-Group0 -1.1034 -1.64037 -0.5665 5.63e-05 Loss time Group1-Group0 1.1034 0.56648 1.6404 5.63e-05 Prob at 12 Group1-Group0 -0.0020 -0.05198 0.0480 9.37e-01 Prob at 24 Group1-Group0 0.0240 -0.00735 0.0553 1.33e-01 Prob at 36 Group1-Group0 0.0060 -0.01229 0.0243 5.20e-01 Prob at 40 Group1-Group0 0.0020 -0.01093 0.0149 7.62e-01 Quantile at 10 % Group1-Group0 -0.3051 -0.65572 0.0455 8.80e-02 Quantile at 15 % Group1-Group0 -0.6318 -0.98557 -0.2780 4.65e-04 Quantile at 20 % Group1-Group0 -0.9475 -1.47683 -0.4183 4.50e-04 Ave of t-year event rates Group1-Group0 0.0075 -0.01433 0.0293 5.01e-01 Ave percentiles Group1-Group0 -0.6282 -1.00468 -0.2516 1.08e-03 $contrast.diff01 Est. Lower 95% Upper 95% p-val RMST Group0-Group1 1.1034 0.5665 1.64037 5.63e-05 Loss time Group0-Group1 -1.1034 -1.6404 -0.56648 5.63e-05 Prob at 12 Group0-Group1 0.0020 -0.0480 0.05198 9.37e-01 Prob at 24 Group0-Group1 -0.0240 -0.0553 0.00735 1.33e-01 Prob at 36 Group0-Group1 -0.0060 -0.0243 0.01229 5.20e-01 Prob at 40 Group0-Group1 -0.0020 -0.0149 0.01093 7.62e-01 Quantile at 10 % Group0-Group1 0.3051 -0.0455 0.65572 8.80e-02 Quantile at 15 % Group0-Group1 0.6318 0.2780 0.98557 4.65e-04 Quantile at 20 % Group0-Group1 0.9475 0.4183 1.47683 4.50e-04 Ave of t-year event rates Group0-Group1 -0.0075 -0.0293 0.01433 5.01e-01 Ave percentiles Group0-Group1 0.6282 0.2516 1.00468 1.08e-03 $contrast.ratio01 Est. Lower 95% Upper 95% p-val RMST Group0/Group1 1.206 1.098 1.325 8.86e-05 Loss time Group0/Group1 0.834 0.763 0.912 6.13e-05 Prob at 12 Group0/Group1 1.009 0.807 1.261 9.38e-01 Prob at 24 Group0/Group1 0.692 0.430 1.114 1.30e-01 Prob at 36 Group0/Group1 0.769 0.344 1.722 5.24e-01 Prob at 40 Group0/Group1 0.833 0.250 2.775 7.66e-01 Quantile at 10 % Group0/Group1 1.602 0.986 2.602 5.69e-02 Quantile at 15 % Group0/Group1 1.862 1.377 2.517 5.45e-05 Quantile at 20 % Group0/Group1 1.936 1.409 2.660 4.62e-05 Ave of t-year event rates Group0/Group1 0.912 0.698 1.192 5.01e-01 Ave percentiles Group0/Group1 1.836 1.337 2.522 1.72e-04 $contrast.ratio10 Est. Lower 95% Upper 95% p-val RMST Group1/Group0 0.829 0.755 0.910 8.86e-05 Loss time Group1/Group0 1.199 1.097 1.310 6.13e-05 Prob at 12 Group1/Group0 0.991 0.793 1.239 9.38e-01 Prob at 24 Group1/Group0 1.444 0.897 2.325 1.30e-01 Prob at 36 Group1/Group0 1.300 0.581 2.911 5.24e-01 Prob at 40 Group1/Group0 1.200 0.360 3.996 7.66e-01 Quantile at 10 % Group1/Group0 0.624 0.384 1.014 5.69e-02 Quantile at 15 % Group1/Group0 0.537 0.397 0.726 5.45e-05 Quantile at 20 % Group1/Group0 0.517 0.376 0.710 4.62e-05 Ave of t-year event rates Group1/Group0 1.096 0.839 1.432 5.01e-01 Ave percentiles Group1/Group0 0.545 0.397 0.748 1.72e-04 attr(,"class") [1] "surv2sample"

というわけで、通常のCox回帰では有意差がなかったものが、 RMSTでは差がでた。 というわけで、はじめの方は効いているけど、後のほうでは効かない感じの 生存曲線で、この手法は有用であった。

文献をちら見すると、比例ハザード性が保たれる状況では大体Powerはlog-rankと同じくらい。 保たれない状況ではPowerはRMSTの方が大きくなるようだ。