RMST生存曲線下面積２ - r-statistics-fanの日記

#RMST生存曲線下面積２

前回の記事で(surv2sampleComp)パッケージを使ってみた
http://r-statistics-fan.hatenablog.com/entry/2014/08/04/225135

今回survRM2パッケージが出たことを知った。

使ってみる
前回と同じく、イレッサのシミュレーションデータを使用する。

#RMST2

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)

#######RMST2##########

library(survRM2)


res.RMST2 <- rmst2(time = data$sur, status = data$event, arm = data$tinib, tau=12)
print(res.RMST2)


############
library(surv2sampleComp)
# 調整なしRMST
res.RMST <- rmstaug(y = data$sur, x = data$cov, arm = data$tinib, delta = data$event, 
                    tau = 12, type = "difference")
res.RMST

結果

> res.RMST2 <- rmst2(time = data$sur, status = data$event, arm = data$tinib, tau=12)
> print(res.RMST2)
##survRM2 result
The truncation time: tau = 12  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 5.346 0.200     4.955     5.738
RMST (arm=0) 6.450 0.189     6.079     6.821


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 6.654 0.200     6.262     7.045
RMTL (arm=0) 5.550 0.189     5.179     5.921


Between-group contrast 
                       Est. lower .95 upper .95 p
RMST (arm=1)-(arm=0) -1.103    -1.643    -0.564 0
RMST (arm=1)/(arm=0)  0.829     0.755     0.910 0
RMTL (arm=1)/(arm=0)  1.199     1.097     1.310 0


##surv2sampleComp result
> res.RMST
$result.ini
               coef  se(coef)         z            p lower .95  upper .95
intercept  6.449593 0.1892441 34.080817 0.000000e+00  6.078674  6.8205112
arm       -1.103421 0.2751565 -4.010159 6.067775e-05 -1.642728 -0.5641147

$result.aug
               coef  se(coef)         z            p lower .95  upper .95
intercept  6.449593 0.1893388 34.063772 0.000000e+00  6.078489  6.8206968
arm       -1.103421 0.2752942 -4.008154 6.119524e-05 -1.642998 -0.5638448

Est. lower .95 upper .95 p
RMST (arm=1)-(arm=0) -1.103 -1.643 -0.564 0

coef se(coef) z p lower .95 upper .95
arm -1.103421 0.2752942 -4.008154 6.119524e-05 -1.642998 -0.5638448

結果は同じ感じになった。