Power2Stage / power.2stage() [Power / Sample Size]
extending Detlew’s post. Try this code:
library(Power2Stage)
CV <- 0.4
n1 <- 12
method <- "C"
pmethod <- c("shifted", "nct", "exact")
st <- proc.time()[[3]]
res <- data.frame(method=method, CV=CV, n1=n1, pmethod=pmethod,
power=NA, TIE=NA, n.tot=NA, pct05=NA, pct50=NA,
pct95=NA, pct.stg2=NA)
for (j in seq_along(pmethod)) {
x <- power.2stage(alpha=rep(0.0294, 2), method=method, CV=CV, n1=n1,
GMR=0.95, theta0=0.95, pmethod=pmethod[j],
npct=c(0.05, 0.5, 0.95))
y <- power.2stage(alpha=rep(0.0294, 2), method=method, CV=CV, n1=n1,
GMR=0.95, theta0=1.25, pmethod=pmethod[j])
res[j, "power"] <- round(x$pBE, 4)
res[j, "TIE"] <- round(y$pBE, 4)
res[j, "n.tot"] <- round(x$nmean, 1)
res[j, 8:10] <- x$nperc
res[j, "pct.stg2"] <- round(x$pct_s2, 1)
}
et <- proc.time()[[3]] - st
cat("run time:", signif(et/60, 1), "minutes\n"); print(res, row.names=FALSE)
run time: 8.9 minutes
method CV n1 pmethod power TIE n.tot pct05 pct50 pct95 pct.stg2
C 0.4 12 shifted 0.7509 0.0342 79.0 34 74 142 99.0
C 0.4 12 nct 0.7506 0.0345 78.9 32 74 142 98.9
C 0.4 12 exact 0.7506 0.0345 78.9 32 74 142 98.9
power
overall power, TIE
empiric Type I Error (for true ratio 1.25), n.tot
expected average total sample size; pct05
, pct50
, pct95
5%, 50%, 95% percentiles of expected total sample size, pct.stg2
percent of studies expected to proceed to the 2nd stage.Reported by Potvin et al. in Tables I & II:
method CV n1 pmethod power TIE n.tot pct05 pct50 pct95 pct.stg2
C 0.4 12 shifted 0.7505 0.0346 79.0 32 74 142 99.0
Try this code with a different seed in each run:
runs <- 25
res <- data.frame(power=rep(NA, runs), TIE=rep(NA, runs))
st <- proc.time()[[3]]
for (j in 1:runs) {
res[j, "power"] <- power.2stage(alpha=rep(0.0294, 2), method=method,
CV=CV, n1=n1, GMR=0.95, theta0=0.95,
pmethod="shifted", setseed=FALSE)$pBE
res[j, "TIE"] <- power.2stage(alpha=rep(0.0294, 2), method=method,
CV=CV, n1=n1, GMR=0.95, theta0=1.25,
pmethod="shifted", setseed=FALSE)$pBE
}
et <- proc.time()[[3]] - st
cat("run time:", signif(et/60, 2), "minutes\n"); summary(res)
run time: 6.4 minutes
power TIE
Min. :0.7482 Min. :0.03403
1st Qu.:0.7503 1st Qu.:0.03443
Median :0.7511 Median :0.03455
Mean :0.7511 Mean :0.03451
3rd Qu.:0.7522 3rd Qu.:0.03460
Max. :0.7537 Max. :0.03485
If you want to estimate the sample size of the 2nd stage I recommend the function
sampleN2.TOST()
of Power2Stage
instead of sampleN.TOST()
of PowerTOST
and subtract n1. In the final analysis we have one degree of freedom less than in the conventional model (due to the additional stage-term). Comparison:x <- sampleN2.TOST(alpha=0.0294, CV=CV, n1=n1, theta0=0.95)
y <- sampleN.TOST(alpha=0.0294, CV=CV, theta0=0.95, details=FALSE, print=FALSE)
print(x, row.names=FALSE) # n2 directly with correct df
print(y, row.names=FALSE) # total n with incorrect df
Design alpha CV theta0 theta1 theta2 n1 Sample size Achieved power Target power
2x2 0.0294 0.4 0.95 0.8 1.25 12 68 0.8107469 0.8
Design alpha CV theta0 theta1 theta2 Sample size Achieved power Target power
2x2 0.0294 0.4 0.95 0.8 1.25 80 0.8105103 0.8
Let’s see whether we can reproduce Potvin’s examples for ‘Method C’:
potvin.examples <- function(alpha0=0.05, alpha=0.0294, n1,
log.GMR1, S.sq1, log.GMR2, S.sq2) {
require(PowerTOST)
require(Power2Stage)
method <- "C"
ltheta1 <- log(0.80)
ltheta2 <- log(1.25)
ltheta0 <- log(0.95)
pass.1 <- pass.2 <- FALSE
pwr.1 <- power.TOST(alpha=alpha0, CV=mse2CV(S.sq1), n=n1,
method="shifted")
if (pwr.1 >= 0.80) {
CI.stg1 <- round(100*CI.BE(alpha=alpha0, pe=exp(log.GMR1),
CV=mse2CV(S.sq1), n=n1), 2)
} else {
CI.stg1 <- round(100*CI.BE(alpha=alpha, pe=exp(log.GMR1),
CV=mse2CV(S.sq1), n=n1), 2)
}
if (CI.stg1[["lower"]] < 0.80 | CI.stg1[["upper"]] > 1.25 & pwr.1 < 0.80) {
n2 <- sampleN2.TOST(alpha=alpha, CV=mse2CV(S.sq1), n1=n1,
theta0=0.95, method="shifted")[["Sample size"]]
} else {
pass.1 <- TRUE
}
if (!pass.1) {
df2 <- n1 + n2 - 3
sem <- sqrt(2*S.sq2/(n1+n2))
CI.stg2 <- round(100*exp(log.GMR2 + c(-1, +1) *
qt(1-alpha, df2)*sem), 2) # balanced!
names(CI.stg2) <- c("lower", "upper")
if (CI.stg2[["lower"]] >= 0.80 | CI.stg2[["upper"]] <= 1.25) pass.2 <- TRUE
}
txt <- paste0("Power in the interim = ", round(100*pwr.1, 1), "%.",
"\nSince power is ")
if (pwr.1 >= 0.80) {
txt <- paste0(txt, "\u226580%, BE assessed with alpha ", alpha0, ".")
} else {
txt <- paste0(txt, "<80%, BE assessed with alpha ", alpha, ".")
}
if (pass.1) {
txt <- paste0(txt, "\nStudy passed in the 1st stage: ")
if (pwr.1 >= 0.80) {
txt <- paste0(txt, "90% CI: ")
} else {
txt <- paste0(txt, 100*(1-2*alpha), "% CI = ")
}
txt <- paste0(txt, CI.stg1[["lower"]], "% \u2013 ", CI.stg1[["upper"]], "%.")
} else {
txt <- paste0(txt, "\nStudy failed in the 1st stage: ")
if (pwr.1 >= 0.80) {
txt <- paste0(txt, "90% CI = ")
} else {
txt <- paste0(txt, 100*(1-2*alpha), "% CI = ")
}
txt <- paste0(txt, CI.stg1[["lower"]], "% \u2013 ", CI.stg1[["upper"]], "%.")
if (pwr.1 < 0.80)
txt <- paste0(txt, "\n2nd stage with ", n2, " subjects initiated.")
}
if (!pass.1 & pwr.1 < 0.80) {
if (pass.2) {
txt <- paste0(txt, "\nStudy passed in the 2nd stage: ")
} else {
txt <- paste0(txt, "\nStudy failed in the 2nd stage: ")
}
}
if (!pass.1 & pwr.1 < 0.80) {
pwr.21 <- Power2Stage:::.calc.power(alpha=alpha, ltheta1=ltheta1,
ltheta2=ltheta2, diffm=ltheta0,
sem=sem, df=df2, method="shifted")
pwr.22 <- Power2Stage:::.calc.power(alpha=alpha, ltheta1=ltheta1,
ltheta2=ltheta2, diffm=log.GMR2,
sem=sem, df=df2, method="shifted")
txt <- paste0(txt, 100*(1-2*alpha), "% CI = ",
CI.stg2[["lower"]], "% \u2013 ", CI.stg2[["upper"]], "%.")
txt <- paste0(txt, "\n(Irrelevant) post hoc powers:",
"\nFor GMR 0.95 = ", round(100*pwr.21, 1), "%, ",
"for GMR ", round(exp(log.GMR2), 4), " = ",
round(100*pwr.22, 1), "%.")
}
cat(txt, "\n")
}
potvin.examples(n1=12, log.GMR1=0.16785, S.sq1=0.020977,
log.GMR2=0.14401, S.sq2=0.0212240)
Power in the interim = 84.1%.
Since power is ≥80%, BE assessed with alpha 0.05.
Study failed in the 1st stage: 90% CI = 106.26% – 131.66%.
potvin.examples(n1=12, log.GMR1=0.08396, S.sq1=0.032634,
log.GMR2=0.014439, S.sq2=0.045896)
Power in the interim = 64.9%.
Since power is <80%, BE assessed with alpha 0.0294.
Study failed in the 1st stage: 94.12% CI = 92.93% – 127.28%.
2nd stage with 8 subjects initiated.
Study passed in the 2nd stage: 94.12% CI = 88.45% – 116.38%.
(Irrelevant) post hoc powers:
For GMR 0.95 = 66.3%, for GMR 1.0145 = 76.9%.
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
![[image]](https://static.bebac.at/pics/Blue_and_yellow_ribbon_UA.png)
Helmut Schütz
![[image]](https://static.bebac.at/img/CC by.png)
The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Complete thread:
- Power R WinNonLin Yura 2018-01-23 15:02 [Power / Sample Size]
- Power R WinNonLin d_labes 2018-01-23 15:25
- Power_TOST in PHX/WNL 6.4+ Helmut 2018-01-23 16:05
- Power_TOST in PHX/WNL 6.4+ d_labes 2018-01-23 16:54
- Power_TOST in PHX/WNL 6.4+ Helmut 2018-01-23 17:25
- Power_TOST in PHX/WNL 6.4+ Yura 2018-01-23 21:55
- Power_TOST in PHX/WNL 6.4+ d_labes 2018-01-23 16:54
- Power_TOST in PHX/WNL 6.4+ Helmut 2018-01-23 16:05
- Power R WinNonLin pjs 2018-03-23 10:41
- Power of two stage design d_labes 2018-03-23 17:02
- Power2Stage / power.2stage()Helmut 2018-03-24 15:49
- Power R WinNonLin d_labes 2018-01-23 15:25