ハザード比

カプランマイヤーによる生存曲線をやっていてCox proportional hazard モデルによるハザード比を出そうとして、シミュレーションデータで遊んでいただけの話。
prodlim パッケージにあるSimSurv 関数は、変数X1とX2があるが、とりあえずこれらは考えず、まったく同じシミュレーターで2群作る。
これを複数回、生存曲線の解析をしたとして、ハザード比（HR）が1 をまたがない、すなわち2群の生存に差がある、のはシミュレーションしたらどれくらい生じるか、というと、これは
生存に差がない、というデータを用意しているので、帰無仮説は真、であるが、このときに有意な結果を生じる確率なので、0.05である。厳密にいうならば有意水準 は、「差がない」のに「差がある」という判断をしてしまう。
シミュレーション結果はこんな感じになる。左の図は、有意な差がある、という結果になってしまった生存解析だが、もともとは同じデータから2群が生成されているので、本当は差がない、はずなのに、サンプリングによるばらつきのせいで有意になってしまった。
これを複数回行い、有意な結果になってしまった場合のHRの95%信頼区間を赤で示しているのが右の図である。これは全体の5% 生じている。

library(survival)
library(prodlim)

niter <- 3000
res <- matrix(0, niter, 3)
n <- 100
g <- rep(1:2, n)
p0 <- 1

for(i in 1:niter){
  dat <- do.call(rbind, replicate(2, SimSurv(n), simplify=FALSE))
  dat <- cbind(dat, group=g)
  #dat <- subset(dat, X1==1)
  r <- coxph(Surv(time, event) ~ group, data = dat)
  p <- pchisq(survdiff(Surv(time, event) ~ group, data = dat)$chisq, df=1, lower.tail=FALSE)
  res[i, ] <- summary(r)$conf.int[c(3, 1, 4)]
  if(p < p0){ # もしp値がより小さい場合には更新する
    p0 <- p
    fit <- prodlim(Hist(time, event) ~ group, data = dat)
  }
}

t.risk <- seq(0, 15, by=3)
par(mfrow=c(1, 2))
plot(fit, atrisk.at=t.risk, axis1.at=t.risk, atrisk.title="No. at Risk", marktime=T)
title("Significant result under true null hypothesis")
legend("bottomleft", legend=sprintf("p = %.6f", p0), bty="o", bg="white", box.col=NA)

matplot(res, type="n", xlab="No. simulation", ylab="Hazard Ratio")
for(i in 1:nrow(res)){
  segments(i, res[i, 1], y1=res[i, 3], col=(any(res[i,1]>1, res[i,3]<1)+1))
}
title(sprintf("Significant simulation: %.3f", mean(res[,1]>1 | res[,3]<1)))