# ----------------------------------------------------------------------------
# power (or alpha) of 2-stage studies with a 2-group parallel desig according
# to Potvin et. al. # methods "B" and "C", modified to include a futility
# criterion Nmax and modified to use PE of stage 1 in sample size estimation
#
# variant of power.2stage.p() wich calculates power always via pooled t-test
# formulas according to the fuglsang 2014 paper
#
# author D.L.
# ----------------------------------------------------------------------------
power.2stage.pAF <- function(method=c("B","C"), alpha0=0.05, alpha=c(0.0294,0.0294),
n1, GMR, CV, targetpower=0.8,
pmethod=c("shifted", "nct", "exact"),
usePE=FALSE, Nmax=Inf, test=c("welch", "t-test", "anova"),
theta0, theta1, theta2, npct=c(0.05, 0.5, 0.95),
nsims, setseed=TRUE, details=FALSE)
{
# Check if called with .2stage. version
check2stage(fname=as.character(sys.call())[1])
if (missing(CV)) stop("CV(s) must be given!")
if (any(CV<=0)) stop("CV(s) must be >0!")
# equal CV's
if (length(CV)==1) {
CVT <- CVR <-CV
} else {
CV <- CV[1:2] # truncate if longer then 2
CVT <- CV[1]; CVR <- CV[2]
}
varT <- CV2mse(CVT)
varR <- CV2mse(CVR)
if (missing(n1)) stop("Number of subjects in stage 1 must be given!")
if (length(n1)>1) {
warning("n1 has to be a scalar. Sum of vector will be used.")
n1 <- sum(n1)
}
if (n1<=0) stop("Number of subjects in stage 1 must be >0!")
if (n1%%2!=0) warning("Number of subjects in stage 1 should be even.\n",
" Will be truncated to even.", immediate. = TRUE)
if (length(alpha) != 2) stop("alpha must have two elements")
if (missing(GMR)) GMR <- 0.95
if (missing(theta1) & missing(theta2)) theta1 <- 0.8
if (!missing(theta1) & missing(theta2)) theta2 <- 1/theta1
if (missing(theta1) & !missing(theta2)) theta1 <- 1/theta2
if (GMR<=theta1 | GMR>=theta2) stop("GMR must be within acceptance range!")
if (missing(theta0)) theta0 <- GMR
if (n1>Nmax) stop("n1>Nmax doesn\'t make sense!")
if(missing(nsims)){
nsims <- 1E5
if(theta0<=theta1 | theta0>=theta2) nsims <- 1E6
}
# check if Potvin B or C
method <- match.arg(method)
# check test
test <- match.arg(test)
test <- tolower(test)
# check if power calculation method is nct or exact
pmethod <- match.arg(pmethod)
if(details){
cat(nsims,"sims. Stage 1")
}
# start timer
ptm <- proc.time()
if (setseed) set.seed(1234567)
ltheta1 <- log(theta1)
ltheta2 <- log(theta2)
lGMR <- log(GMR)
mlog <- log(theta0)
bk <- 4 # 2-group parallel design constant
# reserve memory for the BE result
BE <- rep.int(NA, times=nsims)
# ----- stage 1 ----------------------------------------------------------
n1T <- n1R <- trunc(n1/2) # if not even
nT <- n1T; nR <- n1R
n1 <- n1T + n1R
Cfact <- bk/n1
df <- n1T + n1R -2
tval <- qt(1-alpha[1], df) # ANOVA and t-test
# simulate means via normal distributions
m1T <- rnorm(n=nsims, mean=mlog, sd=sqrt(varT/n1T))
m1R <- rnorm(n=nsims, mean=0, sd=sqrt(varR/n1R))
# point est.
pes <- m1T-m1R
# simulate variances via chi-squared distribution
varsT <- varT*rchisq(n=nsims, df=n1T-1)/(n1T-1)
varsR <- varR*rchisq(n=nsims, df=n1R-1)/(n1R-1)
Vpooled <- ((n1T-1)*varsT + (n1R-1)*varsR)/(n1-2)
if(method=="C"){
# if method=C then calculate power for alpha0=0.05 and plan GMR
# Here we use A.Fuglsangs settings, namely power monitoring steps
# and sample size via pooled t-test
pwr <- .calc.power(alpha=alpha0, ltheta1=ltheta1, ltheta2=ltheta2,
diffm=lGMR, sem=sqrt(bk*Vpooled/n1), df=df,
method=pmethod)
# calculate CIs at alpha0 for all
if (test=="t-test" || test=="anova"){
tval0 <- qt(1-alpha0, df)
hw <- tval0*sqrt(bk*Vpooled/n1)
} else {
# Welch's t-test
se <- sqrt(varsT/nT + varsR/nR)
dfs <- (varsT/nT + varsR/nR)^2/(varsT^2/nT^2/(nT-1)+varsR^2/nR^2/(nR-1))
# dfs <- trunc(dfs)
hw <- qt(1-alpha0,dfs)*se
}
lower <- pes - hw
upper <- pes + hw
# fail or pass
BE <- lower>=ltheta1 & upper<=ltheta2
# if power>0.8 then calculate CI for alpha=0.05
# i.e. if power<0.8 then
BE[pwr<targetpower] <- NA # not yet decided
}
# method "B" or power<=0.8 in method "C"
# evaluate BE with alpha[1]
Vpooled_tmp <- Vpooled[is.na(BE)]
pes_tmp <- pes[is.na(BE)]
varsT_tmp <- varsT[is.na(BE)]
varsR_tmp <- varsR[is.na(BE)]
BE1 <- rep.int(NA, times=length(Vpooled_tmp))
# calculate CI for alpha=alpha[1]
if (test=="t-test" || test=="anova"){
hw <- tval*sqrt(bk*Vpooled_tmp/n1)
} else {
# Welch's t-test
se <- sqrt(varsT_tmp/nT + varsR_tmp/nR)
dfs <- (varsT_tmp/nT + varsR_tmp/nR)^2/(varsT_tmp^2/nT^2/(nT-1) +
varsR_tmp^2/nR^2/(nR-1))
# dfs <- trunc(dfs)
hw <- qt(1-alpha[1],dfs)*se
}
rm(varsT_tmp, varsR_tmp)
lower <- pes_tmp - hw
upper <- pes_tmp + hw
BE1 <- lower>=ltheta1 & upper<=ltheta2
if (method=="C"){
#if BE met -> PASS and stop
#if not BE -> goto sample size estimation i.e flag BE1 as NA
BE1[!BE1] <- NA
} else {
# method B/E
# evaluate power at alpha[2] and planGMR
pwr <- .calc.power(alpha=alpha[2], ltheta1=ltheta1, ltheta2=ltheta2,
diffm=lGMR, sem=sqrt(bk*Vpooled_tmp/n1), df=df,
method=pmethod)
# Next is MSDBE scheme
# if BE met then decide BE regardless of power
# if not BE and power<0.8 then goto stage 2
#BE1[ !BE1 & pwr<targetpower] <- NA
# Xu et al., Potvin method E:
# if not BE and if power >= 0.8 (targetpower) make a second BE evaluation
# with alpha[2]
# only if alpha[1] != alpha[2] necessary, but works also without if(...)
BE12 <- BE1 # reserve memory for second BE decision with alpha[2]
BE11 <- BE1
# calculate CI for alpha=alpha2
if (test=="t-test" || test=="anova"){
hw <- qt(1-alpha[2], df)*sqrt(bk*Vpooled_tmp/n1)
} else {
# Welch's t-test
hw <- qt(1-alpha[2], dfs)*se
# now we are done with them
rm(dfs, se)
}
lower <- pes_tmp - hw
upper <- pes_tmp + hw
BE12 <- lower>=ltheta1 & upper<=ltheta2
# if BE(a1) then BE1=TRUE, regardless of power
BE1[BE11==TRUE] <- TRUE
# if not BE(a1) but power >= 0.8 then make BE decision at alpha2
BE1[BE11==FALSE & pwr>=targetpower] <- BE12[BE11==FALSE & pwr>=targetpower]
# if not BE(a1) and power<0.8 then not decided yet (marker NA)
# will be further decided by futility criterion
BE1[BE11==FALSE & pwr<targetpower] <- NA
# keep care of memory
rm(BE11, BE12)
}
# combine 'stage 0' from method C and stage 1
BE[is.na(BE)] <- BE1
# take care of memory
# done with them
rm(BE1, hw, lower, upper, pes_tmp, Vpooled_tmp)
# time for stage 1
if(details){
cat(" - Time consumed (secs):\n")
print(round((proc.time()-ptm),1))
}
# ------sample size for stage 2 -----------------------------------------
ntot <- rep(n1, times=nsims)
stage <- rep(1, times=nsims)
# filter out those were stage 2 is necessary
pes_tmp <- pes[is.na(BE)]
# Maybe we are already done with stage 1
if(length(pes_tmp)>0){
if(details){
cat("Keep calm. Sample sizes for stage 2 (", length(pes_tmp),
" studies)\n", sep="")
cat("will be estimated. May need some time.\n")
}
# preliminary setting stage=2 for those not yet decided BE
# may be altered for those with nts>Nmax or nts=Inf
# from sample size est. if pe outside acceptance range
# see below
stage[is.na(BE)] <- 2
Vpooled_tmp <- Vpooled[is.na(BE)]
m1T <- m1T[is.na(BE)]
m1R <- m1R[is.na(BE)]
varsT <- varsT[is.na(BE)]
varsR <- varsR[is.na(BE)]
BE2 <- rep.int(NA, times=length(Vpooled_tmp))
s2 <- rep.int( 2, times=length(Vpooled_tmp))
#------ sample size for stage 2 ---------------------------------------
ptms <- proc.time()
# degrees of freedom
if (test=="anova") dfc <- "n-3" else dfc <- "n-2"
if (usePE){
# use mse1 & pe1 like in the paper of Karalis/Macheras
# sample size function returns Inf if pe1 is outside acceptance range
nts <- .sampleN2(alpha=alpha[2], targetpower=targetpower, ltheta0=pes_tmp,
mse=Vpooled_tmp, ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod, bk=4, dfc=dfc)
} else {
# use mse1 & plan GMR to calculate sample size (original Potvin)
nts <- .sampleN2(alpha=alpha[2], targetpower=targetpower, ltheta0=lGMR,
mse=Vpooled_tmp, ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod, bk=4, dfc=dfc)
}
n2 <- ifelse(nts>n1, nts - n1, 0)
if(details){
cat("Time consumed (secs):\n")
print(round((proc.time()-ptms),1))
}
# futility rule: if nts > Nmax -> stay with stage 1 result not BE
# ntotal = n1 reasonable?
if (is.finite(Nmax) | any(!is.finite(nts))){
# sample size may return Inf if PE is used in ss estimation
# in that case we stay with stage 1
BE2[!is.finite(n2) | (n1+n2)>Nmax] <- FALSE
# and we are counting these for stage 1
s2[BE2==FALSE] <- 1
# debug print
# cat(sum(!BE2, na.rm=T)," cases with nts>Nmax or nts=Inf\n")
# save
stage[is.na(BE)] <- s2
# save the FALSE and NA in BE
BE[is.na(BE)] <- BE2
# filter out those were BE was yet not decided
pes_tmp <- pes_tmp[is.na(BE2)]
m1T <- m1T[is.na(BE2)]
m1R <- m1R[is.na(BE2)]
Vpooled_tmp <- Vpooled_tmp[is.na(BE2)]
varsT <- varsT[is.na(BE2)]
varsR <- varsR[is.na(BE2)]
n2 <- n2[is.na(BE2)]
} # end of futility Nmax
# ----- simulate stage 2 data ------
nsim2 <- length(pes_tmp)
n2T <- n2R <- n2/2
# to avoid warnings for n2X==0 in rnorm() and n2X<=1 in rchisq()
ow <- options("warn")
options(warn=-1)
m2T <- ifelse(n2T>0, rnorm(n=nsim2, mean=mlog, sd=sqrt(varT/n2T)), 0)
m2R <- ifelse(n2R>0, rnorm(n=nsim2, mean=0, sd=sqrt(varR/n2R)), 0)
# means T/R
nT <- n1T+n2T
nR <- n1R+n2R
mT <- (n1T*m1T+n2T*m2T)/nT
mR <- (n1R*m1R+n2R*m2R)/nR
# point est.
pes <- mT-mR
# means for s1, s2 (stages)
m1 <- (n1T*m1T + n1R*m1R)/(n1T+n1R)
m2 <- ifelse((n2T+n2R)>0, (n2T*m2T + n2R*m2R)/(n2T+n2R),0)
# total mean
mt <- (nT*mT+nR*mR)/(nT+nR)
# simulate variances via chi-squared distribution
# attention! in case of n2X==1 rchisq gives NaN!
# TODO: work out the correct way for n2X==1
vars2T <- ifelse(n2T>1, varT*rchisq(n=nsim2, df=n2T-1)/(n2T-1), 0)
vars2R <- ifelse(n2R>1, varR*rchisq(n=nsim2, df=n2R-1)/(n2R-1), 0)
# reset options
options(ow)
# vars T/R over stage 1, stage 2
# sy2 = sum of y squared
# SQ = sum((y-mean)^2) = sum(y^2) - sum(y)^2/n = sum(y^2) - n*mean^2
# var = SQ/(n-1)
# bug in pre 0.4-5 (was mean^2/n instead of n*mean^2) which lead to
# negativ residual variance and therefore BE = NA
# example:
# power.2stage.p(CV=0.4, n1=10, theta0=0.8, test="a")
sy2T <- (n1T-1)*varsT + n1T*m1T^2
sy2T <- ifelse(n2T>0, sy2T + (n2T-1)*vars2T + n2T*m2T^2, sy2T)
sy2R <- (n1R-1)*varsR + n1R*m1R^2
sy2R <- ifelse(n2R>0, sy2R + (n2R-1)*vars2R + n2R*m2R^2, sy2R)
varsT <- (sy2T - nT*mT^2)/(nT-1)
varsR <- (sy2R - nR*mR^2)/(nR-1)
sy2 <- sy2T+sy2R
sqtot <- (sy2-(nT+nR)*mt^2)
rm(sy2T, sy2R, sy2)
# calculate CI for stage 2 with alpha[2]
if (test=="t-test" || test=="anova"){
Vpooled <- ((nT-1)*varsT + (nR-1)*varsR)/(nT+nR-2)
dfs <- (nT+nR-2)
if (test=="anova"){
# no subjects in stage 2 ((n2R+n2T)==0) may occure in case of
# high n1 and/or haybittle-peto alpha's
# according to Anders paper
Vpooled <- ifelse(n2R+n2T>0,
(sqtot - (n1*(m1-mt)^2 + n2*(m2-mt)^2) # stage
-(nT*(mT-mt)^2 + nR*(mR-mt)^2))/(dfs-1), # treatment
Vpooled)
dfs <- ifelse(n2R+n2T>0, dfs-1, dfs)
}
hw <- qt(1-alpha[2],dfs)*sqrt(bk*Vpooled/(nT+nR))
rm(Vpooled, dfs)
} else {
# Welch's t-test
se <- sqrt(varsT/nT + varsR/nR)
dfs <- (varsT/nT + varsR/nR)^2/(varsT^2/nT^2/(nT-1)+varsR^2/nR^2/(nR-1))
#dfs <- trunc(dfs)
hw <- qt(1-alpha[2],dfs)*se
rm(se, dfs)
}
lower <- pes - hw
upper <- pes + hw
BE2 <- lower>=ltheta1 & upper<=ltheta2
# combine stage 1 & stage 2
ntot[is.na(BE)] <- n1+n2
BE[is.na(BE)] <- BE2
# done with them
rm(BE2, nts, lower, upper, hw, n2T, n2R, m2T, m2R)
} # end stage 2 calculations
# take care of memory
rm(pes_tmp, pes)
# the return list
res <- list(design="2 parallel groups", method=method,
alpha0=ifelse(method=="C",alpha0,NA), alpha=alpha, CV=CV, n1=n1,
GMR=GMR, test=test, targetpower=targetpower, pmethod=pmethod,
theta0=exp(mlog), theta1=theta1, theta2=theta2,
usePE=usePE, Nmax=Nmax, nsims=nsims,
# results
pBE=sum(BE)/nsims, pBE_s1=sum(BE[stage==1])/nsims,
# Dec 2014: meaning of pct_s2 changed
pct_s2=100*sum(ntot>n1)/nsims,
nmean=mean(ntot), nrange=range(ntot), nperc=quantile(ntot, p=npct))
# table object summarizing the discrete distri of ntot
# only if usePE=FALSE or if usePE=TRUE then Nmax must be finite
# or return it always?
if (usePE==FALSE | (usePE==TRUE & is.finite(Nmax))){
res$ntable <- table(ntot)
}
if (details){
cat("Total time consumed (secs):\n")
print(round((proc.time()-ptm),1))
cat("\n")
}
class(res) <- c("pwrtsd", "list")
return(res)
} #end function
# alias
power.tsd.pAF <- power.2stage.pAF
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.