生存解析のサンプルサイズ

こんな感じの生存曲線を見た。
元論文では対照のマウスにある処置をすると生存期間が伸びる、ということで、それぞれ20匹のマウスを使っていた。

library(survival)
# 図から再現したらこんな感じ
n1 <- n2 <- 20
clp1 <- rep(c(24, 30, 42, 48), c(7,9,1,3))
clp2 <- rep(c(30, 42, 48, 54, 66, 96), c(6,7,1,3,1,2))
d1 <- rep(1, n1)
d2 <- replace(d1, clp2>=96, 0)

group <- rep(1:2, each=n1)
library(survival)
s1 <- survfit(Surv(time=c(clp1, clp2), event=c(d1, d2)) ~ group)
p <- survdiff(Surv(time=c(clp1, clp2), event=c(d1, d2)) ~ group)
pv <- pchisq(p$chisq, 1, lower.tail=FALSE)
plot(s1, col=c("blue", "red"), lwd=3, xaxp=c(0, 96, 12), las=1,
     xlab="Time [hrs]", ylab="Survival")
legend("topright", legend=c("vehicle", "treatment"), col=c("blue", "red"), lwd=3, bty="n", cex=2)
text(80, 0.3, paste("p =", round(pv,5)), cex=1.5)

　
あんまり差がないように見えたが、p値を計算すると割ときつい値が出てる。
となると、20匹も実験に使う必要があったのかというわけで、サンプルサイズと検出力をシミュレーションしてみると、事後解析的には検出力80%でいいなら10匹、15匹も使えば90%にはなるようである。

niter <- 300
mat <- matrix(NA, nr=n1, nc=niter)
for(j in seq(niter)){
	for(i in 3:n1){
		idx1 <- sample(seq(n1), size=i, replace=TRUE)
		idx2 <- sample(seq(n2), size=i, replace=TRUE)
		tmp.surv <- Surv(time=c(clp1[idx1], clp2[idx2]), event=c(d1[idx1], d2[idx2])) ~ rep(1:2, each=i)
		s1 <- survfit(tmp.surv)
		p <- try(survdiff(tmp.surv), silent=TRUE)
		pv <- ifelse(class(p)=="try-error", 1, pchisq(p$chisq, 1, lower.tail=FALSE))
		mat[i, j] <- pv
	}
}

po <- rowMeans(mat < 0.05)
plot(po, type="o", pch=16, xlab="sample size", ylab="power")