r-statistics-fanの日記

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

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

R

生存曲線下面積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試験の論文について以前の記事で考察した。 

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

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

とりあえず、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の調整なしでまずは解析したかったが、入力しないとエラーになる。 とりあえず全部１で投入する。

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の方が大きくなるようだ。