臨床試験における生存解析では、参加登録した患者をランダムに2群に割り付け、新薬Xと旧薬Aをそれぞれの群に投与する。
新薬のMST(中央生存期間)は適当に17ヶ月、旧薬のMSTは10ヶ月。
カプランマイヤーの生存解析では、ハザード比がどの時刻においても新薬群と旧薬群で一定、という前提がある。つまり、任意の時刻において生存曲線の傾きは一定になる。
参加登録を開始した直後から登録人数は増えるのだが、ある程度のところで減り始めて、後は月に3人とかプラトーに達する。
参加登録したはいいけど、いろんな事情(患者が行方不明とかデータ紛失とか)で、年5%(とする)くらいの患者は臨床試験から脱落する。
という事情があるなかで、実際に17ヶ月と10ヶ月の差を有意水準でそれなりに検出しようと思ったら(power = 0.8)、参加登録する患者は何人必要で何ヶ月臨床試験を続けなければならないか、という相談を受けた。
お金があればEastというソフトウェアでできるらしいが、お金があればな!!
というわけで無料のRでシミュレーションすることで計算しよう。
28ヶ月の時点でこんな感じの参加登録があったとする。これを適当な関数で内挿して今後の動向を予測しておく。
n <- c(1,1,3,3,5,5,7,7,7,9,9,9,8,7,7,6,6,6,5,4,3,3,3,3,3,3,3,3) sp <- smooth.spline(n, df=20) ext.month <- 20 # 必要なフォローアップ月数を十分にとる plot(n, xlim=c(0, length(n)+ext.month), xlab="Month", ylab="No. patient in trial") lines(predict(sp, seq(0, length(n)+ext.month, length=200)), col=2, lwd=3) title("spline")
いまある28ヶ月の参加登録者数に、何ヶ月増やしたら十分な検出力で結果が得られるかをシミュレーションする。
ここで、脱落率は年単位なので、月単位でシミュレーションをする都合上、1ヶ月ごとに脱落する確率に変換して、誰が脱落するかをシミュレーションしている。
また、生存期間は傾きが一定ということで一様分布から作成している。
power <- 0.8 alpha <- 0.05 cencor <- 0.05 # 一年間での脱落率 m.cencor <- 1-(1-cencor)^(1/12) # 月での脱落率 MST1 <- 17 # 新薬の中央生存期間 MST2 <- 10 # 旧薬の中央生存期間 # ここからシミュレーション niter <- 1000 # フォローアップひと月ごとのシミュレーション回数 ext.month <- 20 # 必要なフォローアップ月数を十分にとる pval <- matrix(0, niter, ext.month) # p値 for(nit in seq(niter)){ pb <- txtProgressBar(min=1, max=niter, style=3) setTxtProgressBar(pb, nit) for(m0 in seq(ext.month)){ patient <- floor(predict(sp, seq(length(n) + m0-1))$y) N <- sum(patient) # 追加フォローアップ期間で見込まれる参加人数 n1 <- N %/% 2 # 2群の分け方はランダムで均等 n2 <- N - n1 # ということにする group <- replace(rep(1, N), sample(seq(N), size=n2), 2) # ランダム割り付け dropout <- replicate(length(n) + m0, sample(1:0, size=N, prob=c(m.cencor, 1-m.cencor), replace=TRUE)) # 1 なら観察期間中に脱落してしまう。月単位で見ている。 time.in <- unlist(mapply(rep, seq(patient), patient)) # 臨床試験に加わった時刻 x1 <- 2*MST1*(1-runif(n1)) x2 <- 2*MST2*(1-runif(n2)) X <- list(x1, x2) # 打ち切りや脱落がないときの、真の生存期間 sanka <- tapply(time.in, group, sort) # 割りつけたあとの参加月 time.true <- mapply(function(x, y) x+y, sanka, X, SIMPLIFY=FALSE) # 打ち切りや脱落がないときの、真の死亡時刻 #is.cencor <- mapply(function(x) sample(1:0, size=x, prob=c(m.cencor, 1-m.cencor), replace=TRUE), lapply(X, floor)) is.cencor <- mapply(function(y) mapply(function(x) sample(1:0, size=x, prob=c(m.cencor, 1-m.cencor), replace=TRUE), y, SIMPLIFY=FALSE), lapply(X, floor), SIMPLIFY=FALSE) Y <- X # 死亡まで観測できたか、それとも打ち切られたか event <- mapply(rep, 1, c(n1,n2), SIMPLIFY=FALSE) # とりあえずは、死亡まで観測できたとする for(j in seq(is.cencor)){ for(i in seq(is.cencor[[j]])){ event[[j]][i] <- ifelse(sum(is.cencor[[j]][[i]]) == 0, 1, 0) Y[[j]][i] <- ifelse(sum(is.cencor[[j]][[i]]) > 0, head(which(is.cencor[[j]][[i]] == 1), 1), X[[j]][i]) # 打ち切りのときは打ち切られた時刻までの観察期間を代入 } } time.cencor <- mapply(function(x, y) x+y, sanka, Y, SIMPLIFY=FALSE) # 脱落があるときの死亡時刻 time2 <- mapply(function(x) replace(x, x > length(patient)+m0, length(patient)+m0), time.cencor, SIMPLIFY=FALSE) # フォローアップ時に生きているものはその時刻 event1 <- mapply(function(x, y) replace(x, y > length(patient), 0), event, time.cencor, SIMPLIFY=FALSE) # 打ち切り時に生きているものをカウント s0 <- Surv(time=unlist(time2)-unlist(sanka), event=unlist(event1)) s0fit <- survfit(s0 ~ sort(group)) lr <- survdiff(s0 ~ sort(group)) #pval[nit, m0] <- 1 - pchisq(lr$chisq,1) # 2つの生存曲線が同じかどうかの検定 pval[nit, m0] <- s0summary$table[1, "0.95LCL"]-s0summary$table[2, "0.95UCL"] > 0 # 2つのMSTが同じかどうかの検定。0以上ならCIが被らず統計学的に有意 } }
plot(s0fit, col=1:2, lwd=3, xlab="month", ylab="survival") legend("topright", legend=c("新薬", "旧薬"), bty="n", col=1:2, lwd=4, cex=2) for(i in seq(nrow(s0summary$table))){ for(j in 6:7){ arrows(s0summary$table[i,"median"], 0.5, s0summary$table[i, j], length=0.05, angle=90, col=i, lwd=4) } }
powers <- colMeans(pval) plot(powers, xlab="Additional follow-up months", ylab="power", ylim=c(0, 1))
13ヶ月くらい追加フォローしたら検出力0.8でなんとか統計学的有意差を示せそう。このとき、2群で164人必要。
そもそも7ヶ月の差を検出するのはものすごい難病で治療法がないとかそういうのでないと多分試験自体が承認されそうにないかも知れないけどシミュレーション上、こうなる。