# R : 受け入れ期間参考：[新版]実用SAS生物統計ハンドブック[SAS®9.4/R3.2.0対応]（臨床評価研究会）
# Tr : 追跡期間
# TIME : TIME年生存率
# A : 100A(%)生存率(<1)
# B : 100B(%)生存率(<1)
# ALPHA : 有意水準(<1)
# BETA : 検出力（1-BETA）
# S : 例数比(Na/Nb)
R <- 3
Tr <- 5
TIME <- 5
A <- 0.2
B <- 0.05
ALPHA <- 0.05
BETA <- 0.2
S <- 1
LA <- -log(A)/TIME
LB <- -log(B)/TIME
LBER <- (LA + LB)/2
QA <- 1/(1 + S)
QB <- S/(1 + S)
ZA <- qnorm(1 - ALPHA/2)
ZB <- qnorm(1 - BETA)
P00 <- LBER * LBER
P01 <- LA * LA
P02 <- LB * LB
P10 <- LBER^2/(1 -exp(-LBER*Tr))
P11 <- LA^2/(1-exp(-LA*Tr))
P12 <- LB^2/(1-exp(-LB*Tr))
P20 <- LBER^3*Tr/(LBER*Tr-1+exp(-LBER*Tr))
P21 <- LA^3*Tr/(LA*Tr-1+exp(-LA*Tr))
P22 <- LB^3*Tr/(LB*Tr-1+exp(-LB*Tr))
P30 <- LBER^2*(1-(exp(-LBER*(Tr-R))-exp(-LBER*Tr))/(LBER*R))^(-1)
P31 <- LA^2*(1-(exp(-LA*(Tr-R))-exp(-LA*Tr))/(LA*R))^-1
P32 <- LB^2*(1-(exp(-LB*(Tr-R))-exp(-LB*Tr))/(LB*R))^-1
N <- ((ZA*sqrt(P00*(1/QB+1/QA)) + ZB*sqrt(P01/QA+P02/QB))/abs(LA-LB))^2/(1+S)
M <- ((ZA*sqrt(P10*(1/QB+1/QA)) + ZB*sqrt(P11/QA+P12/QB))/abs(LA-LB))^2/(1+S)
O <- ((ZA*sqrt(P30*(1/QB+1/QA)) + ZB*sqrt(P31/QA+P32/QB))/abs(LA-LB))^2/(1+S)
L <- ((ZA*sqrt(P20*(1/QB+1/QA)) + ZB*sqrt(P21/QA+P22/QB))/abs(LA-LB))^2/(1+S)
c(S=S, 受け入れ期間=R, 生存率=TIME, ハザードA=LA, ハザードB=LB,
ALPHA=ALPHA, ZA=ZA, BETA=BETA, ZB=ZB,
"R=0,T=Inf"=N, "R=0, T>0"=M, "0<R<T"=O, "O<R=T"=L)
No.22524 Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布 【青木繁伸】 2018/05/02(Wed) 23:23
> 打ち切り症例が発生する確率と，打ち切り例の場合の打ち切りまでの時間に，どのような分布とパラメータを仮定すべきなのか
シミュレーションデータの各例は，以下の順で生成されるのでは？
(1) 受入期間内に一様分布で研究に参入する
(2) エンドポイントまでの期間は指数分布で決まる
(3) エンドポイント発生が追跡期間後になるなら「打切り例」，前なら「故障例」
(4) 「打切り例」なら，研究対象であったのは研究参入から追跡期間終了までの期間
(5) 「故障例」なら，研究対象であったのは研究参入からエンドポイントまでの期間
これで，「打ち切り症例が発生する確率と，打ち切り例の場合の打ち切りまでの時間に，どのような分布とパラメータを仮定すべきなのか」というのは無関係，意味がないことになる。
なお，あなたがいっている「打ち切り症例」が「脱落例（追跡不能例）」を含んでいるなら，そもそもそういうものは考える必要がない。例えば，2群の平均値の差のサンプルサイズの決定などでも，測定不能，測定不全，拒否などによる例数減は考えないでしょう？受入開始 受入終了 追跡終了
+..........................+.......................+
o--------------x 受入終了前に死亡
o-------------------------------x 追跡終了前に死亡
o-------------------------------------@=======x 打切り例: 追跡終了後に死亡
- と = の和が指数分布で決定されるエンドポイントまでの期間
- が研究対象であった期間
No.22528 Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布 【kamo】 2018/05/03(Thu) 21:10
青木 繁伸 先生
ご返信頂きましてありがとうございます。
打ち切り例は，脱落例を含むという前提で質問申し上げておりました。
二群の平均値の検定の場合は，（補完などの特別な処理をしなければ）脱落例は解析対象に含められませんが，生存時間解析においては，そういった対象も情報として利用するようですので，
サンプルサイズの計算においても実験途中の脱落を考慮していると理解しておりました。
ご教授頂いた内容を踏まえ，Rコードを書いてみました。
しかし，各群57例の場合の検出力は70~75%程度にとどまり，上の計算結果と整合しませんでした。
どこか誤りがありますでしょうか。
恐れ入りますが，改めてご指導賜りたく，何卒よろしくお願い申し上げます。# 初期化
sign <- NULL
pval <- NULL
# ハザードを計算
B <- 0.05
A <- 0.2
LB <- -log(B)/(5*365.25)
LA <- -log(A)/(5*365.25)
n <- 57*2
accrualtime <- 3*365.25 # 登録期間
followuptime <- 5*365.25 # 観察期間※登録期間を含む
for(N in 1:1000){
# 初期化
x1 <- NULL # 群
time <- NULL # イベントまでの時間
timec <- NULL # 打ち切りまでの時間
randomized.at <- NULL
dt <- NULL
status <- NULL
i <- 0
for(cal in 1:accrualtime){# 登録期間の初日を1として，
entry <- rbinom(size=1,n=1,prob=n/accrualtime) # 当該時点に登録があるか（登録期間を通じて一様に集積されることを仮定する
if(entry==1){ # 当該時点に登録ありの場合
randomized.at <- c(randomized.at,cal) # 登録日
i <- i+1 # 例数カウンター
x1[i] <- rbinom(size=1,n=1,prob=0.5) # ランダム化
if(x1[i]==0){time[i] <- rexp(n=1,rate=LB)} # 対照群の場合の生存時間
else {time[i] <- rexp(n=1,rate=LA)} # 治療群の場合の生存時間
timec <- followuptime-randomized.at[i] # 試験終了による打ち切り
if(time[i] < timec){ # 打ち切り前にイベント発生
status[i] <- 1 # イベント
}
if(timec < time[i]){ # イベント発生前に打ち切り
status[i] <- 0 # 打ち切り
time[i] <- timec # 打ち切り症例には打ち切りまでの時間を代入置換
}
}
else {} # 登録なしの場合はレコード発生させない
if (i==n) break # 予定症例数に達したら終了
}
# ログランク検定
result <- survdiff(Surv(time)~x1)
# p値を直接取り出せないので，検定統計量を格納し，そこからp値を導出
pval[N] <- pchisq(result$chisq, length(result$n)-1, lower.tail = FALSE)
if (pval[N] < 0.05){sign[N] <- 1} else{sign[N] <- 0}
}
sum(sign)/N
hist(pval)
> sum(sign)/N
[1] 0.724
No.22529 Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布 【青木繁伸】 2018/05/04(Fri) 14:35
一番の間違いは，survdiff(Surv(time) ~ x1)としたところ。これは，survdiff(Surv(time, status) ~ x1)でしょう。
もう一つは，x1[i] <- rbinom(size=1,n=1,prob=0.5) # ランダム化ですが，これでは群ごとのサンプルサイズが同じになりません。
その他も，無駄を省き，分かりやすく書いて，ベクトル演算による高速化を図ったので，参考にしてみてください。f = function(n.each, B, A, accrualtime, followuptime, REPETITION = 1000) {
library(survival)
LB = -log(B) / followuptime
LA = -log(A) / followuptime
n = n.each * 2
x1 = rep(0:1, each = n.each)
pval = numeric(REPETITION)
for (N in 1:REPETITION) {
# 研究に登録される時刻（研究開始時刻を 0 とする）
entry = runif(n, min = 0, max = accrualtime) # 一様分布に従う
# 追跡可能如何に関わらず，登録時刻からエンドポイントまでの時間
# 前半に B 群，後半に A 群をまとめて生成しても一般性を損なわない
surv.time = c(rexp(n.each, rate = LB), rexp(n.each, rate = LA))
# 研究開始時刻（0）からエンドポイントまでの時間は entry + surv.time
# その時間が観察期間より短い場合は故障例（status = 1)
status = (entry + surv.time) < followuptime
# 故障例ならエンドポイントまでの時間は surv.time そのもの
# 打切り例 status = 0 なら観察期間終了までの時間
time = ifelse(status, surv.time, followuptime - entry)
# ログランク検定
# result = survdiff(Surv(time) ~ x1) # ここが間違い
result = survdiff(Surv(time, status) ~ x1) # こちらが正しい
# p値を直接取り出せないので，検定統計量を格納し，そこからp値を導出
pval[N] = pchisq(result$chisq, 1, lower.tail = FALSE)
}
mean(pval < 0.05)
}
n.each = 57
B = 0.05
A = 0.2
# 以下の2項目の単位は何でもよい
accrualtime = 3 * 365.25 # 登録期間
followuptime = 5 * 365.25 # 観察期間 ※登録期間を含む
REPETITION = 1000 # シミュレーション試行回数
> f(n.each, B, A, accrualtime, followuptime, REPETITION) # 実行結果は毎回異なる
[1] 0.823
> f(n.each, B, A, 3, 5, 10000)
[1] 0.8253
No.22533 Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布 【kamo】 2018/05/07(Mon) 18:36
青木 繁伸 先生
ご返信ありがとうございます。
>一番の間違いは，>survdiff(Surv(time) ~ x1)>としたところ。これは，>survdiff(Surv(time, status) ~ x1)>でしょう。
基本的なところを間違っておりました。大変失礼いたしました。
この点を修正したところ，理論式と概ね整合のとれる結果が得られました。
>もう一つは，>x1[i] <- rbinom(size=1,n=1,prob=0.5) #ランダム化>ですが，これでは群ごとのサンプルサイズが同じになりません。
乱数なので同じサンプルサイズになるとは限らない，といういことですね。ご指摘ありがとうございます。
シミュレーション回数10000回で，頂いたプログラムを何度か回してみましたが，検出力は全て82%と少し，という感じでございました。
80~80.5%程度の範囲に収まるのかなと勝手に想像しておりましたが，これはやはり，シミュレーションと理論式の違い（誤差?）と考えてよいのでしょうか？
No.22534 Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布 【青木繁伸】 2018/05/07(Mon) 21:06
他のサンプルサイズ決定については，
http://www.st.nanzan-u.ac.jp/info/ma-thesis/2008/TANAKA/m07mm028.pdf
が，参考になるのではないかと思います。
私はとても追い切れませんが，必要とするあなたらならそれぞれを試してみるとよいのではないでしょうか？
その他にも，いくつかあるかも知れませんね（追求する情熱が私にはない）
何かわかったら，ここで教えてください。
● 「統計学関連なんでもあり」の過去ログ--- 048 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る