Nothing
# --------------------------------------------------------------------------
# power (or alpha) of 2-stage studies according to Potvin et al.
# methods "B" and "C", modified to include a futility criterion Nmax
# modified to use PE of stage 1 in sample size estimation
#
# Author D.L.
# --------------------------------------------------------------------------
power.2stage <- function(method=c("B","C","B0"), alpha0=0.05, alpha=c(0.0294,0.0294),
n1, GMR, CV, targetpower=0.8,
pmethod=c("nct","exact", "shifted"),
usePE=FALSE, Nmax=Inf, min.n2=0, theta0, theta1, theta2,
npct=c(0.05, 0.5, 0.95), nsims, setseed=TRUE, details=FALSE)
{
check2stage(fname=as.character(sys.call())[1])
if (missing(CV)) stop("CV must be given.")
if (CV<=0) stop("CV must be >0.")
if (missing(n1)) stop("Number of subjects in stage 1 must be given.")
if (n1<=0) stop("Number of subjects in stage 1 must be >0.")
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
}
if(min.n2!=0 & min.n2<2) stop("min.n2 has to be at least +2.")
# make even (round up)
if( min.n2%%2 != 0) {
min.n2 <- min.n2 + min.n2%%2
message("min.n2 rounded up to next even ", min.n2)
}
# check if Potvin B or C
method <- match.arg(method)
# 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)
mse <- CV2mse(CV)
bk <- 2 # 2x2x2 crossover design const
# reserve memory
BE <- rep.int(NA, times=nsims)
# ----- stage 1 ----------------------------------------------------------
Cfact <- bk/n1
df <- n1-2
tval <- qt(1-alpha[1], df)
sdm <- sqrt(mse*Cfact)
# simulate point est. via normal distribution
pes <- rnorm(n=nsims, mean=mlog, sd=sdm)
# simulate mse via chi-squared distribution
mses <- mse*rchisq(n=nsims, df=df)/df
if(method=="C"){
# if method=C then calculate power for alpha0=0.05 and plan GMR
pwr <- .calc.power(alpha=alpha0, ltheta1=ltheta1, ltheta2=ltheta2,
diffm=lGMR, sem=sqrt(bk*mses/n1), df=df, method=pmethod)
tval0 <- qt(1-alpha0, df)
hw <- tval0*sqrt(Cfact*mses)
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"
# calculate CI for alpha=alpha1
mses_tmp <- mses[is.na(BE)]
pes_tmp <- pes[is.na(BE)]
BE1 <- rep.int(NA, times=length(mses_tmp))
hw <- tval*sqrt(Cfact*mses_tmp)
lower <- pes_tmp - hw
upper <- pes_tmp + hw
BE1 <- lower>=ltheta1 & upper<=ltheta2
if (method=="C"){
#if BE met -> PASS stop
#if not BE -> goto sample size estimation i.e flag BE1 as NA
BE1[!BE1] <- NA
} else {
# method E == B
pwr_alpha <- alpha[2]
# MSDBE == "B0"
if(method=="B0") pwr_alpha <- alpha[1]
pwr <- .calc.power(alpha=pwr_alpha, ltheta1=ltheta1, ltheta2=ltheta2,
diffm=lGMR, sem=sqrt(bk*mses_tmp/n1), df=df,
method=pmethod)
if(method=="B0"){
# B0 == MSDBE: if power > targetpower then STOP: FAIL
# 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
} else {
# Potvin method E == B:
# if not BE and if power >= 0.8 (targetpower) make a second BE evaluation
# with alpha[2]
# but only if alpha[1] != alpha[2] ?
# in case of alpha[1] == alpha[2] both BE12 and BE11 are identical
BE12 <- BE1 # reserve memory
BE11 <- BE1
# CI with alpha[2]
tval <- qt(1-alpha[2], df)
hw <- tval*sqrt(Cfact*mses_tmp)
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 (marker NA), goto stage 2
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)
# 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 nt>Nmax or nt=Inf
# from sample size est. if pe outside acceptance range
# see below
stage[is.na(BE)] <- 2
mses_tmp <- mses[is.na(BE)]
BE2 <- rep.int(NA, times=length(mses_tmp))
s2 <- rep.int(2, times=length(mses_tmp))
#------ sample size for stage 2 ---------------------------------------
ptms <- proc.time()
# Aug. 2017: .sampleN2() now uses N-3 as df in SSR
if (usePE){
# use mse1 & pe1 like in the paper of Karalis/Macheras
# sample size function returns Inf if pe1 is outside acceptance range
nt <- .sampleN2(alpha=alpha[2], targetpower=targetpower, ltheta0=pes_tmp,
mse=mses_tmp, ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod)
} else {
# use mse1 & GMR to calculate sample size (original Potvin)
nt <- .sampleN2(alpha=alpha[2], targetpower=targetpower, ltheta0=lGMR,
mse=mses_tmp, ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod)
}
n2 <- ifelse(nt>n1, nt - n1, 0)
# assure a min.n2
n2 <- ifelse(n2<min.n2, min.n2, n2)
# some upper bound due to numerics?
#n2 <- ifelse(n2>100000, 100000, n2) # may not necessary
if(details){
cat("Time consumed (secs):\n")
print(round((proc.time()-ptms),1))
}
# futility rule: if nt > Nmax -> stay with stage 1 result: not BE
if (is.finite(Nmax) | any(!is.finite(nt))){
# 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 nt>Nmax or nt=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)]
mses_tmp <- mses_tmp[is.na(BE2)]
n2 <- n2[is.na(BE2)]
}
# ---------- stage 2 evaluation --------------------------------------
m1 <- pes_tmp
SS1 <- (n1-2)*mses_tmp
nsim2 <- length(pes_tmp)
# to avoid warnings for n2=0 in rnorm() and rchisq()
ow <- options("warn")
options(warn=-1)
m2 <- ifelse(n2>0, rnorm(n=nsim2, mean=mlog, sd=sqrt(mse*bk/n2)), 0)
# ??? (n2-2) cancels out!
SS2 <- ifelse(n2>2, (n2-2)*mse*rchisq(n=nsim2, df=n2-2)/(n2-2), 0)
# reset options
options(ow)
SSmean <- ifelse(n2>0, (m1-m2)^2/(2/n1+2/n2), 0)
nt <- n1+n2
df2 <- ifelse(n2>0, nt-3, n1-2)
pe2 <- ifelse(n2>0, (n1*m1+n2*m2)/nt, pes_tmp)
mse2 <- ifelse(n2>0, (SS1+SSmean+SS2)/df2, mses_tmp)
# take care of memory
rm(m1, m2, SS1, SS2, SSmean)
# calculate CI for stage 2 with alpha[2]
tval2 <- qt(1-alpha[2], df2)
hw <- tval2*sqrt(mse2*bk/nt)
lower <- pe2 - hw
upper <- pe2 + hw
BE2 <- lower>=ltheta1 & upper<=ltheta2
# combine stage 1 & stage 2
ntot[is.na(BE)] <- nt
BE[is.na(BE)] <- BE2
# done with them
rm(BE2, nt, lower, upper, hw)
} # end stage 2 calculations
# take care of memory
rm(pes_tmp, mses_tmp)
# the return list
res <- list(design="2x2 crossover",
method=method, alpha0=ifelse(method=="C",alpha0,NA), alpha=alpha,
CV=CV, n1=n1, GMR=GMR, targetpower=targetpower, pmethod=pmethod,
theta0=exp(mlog), theta1=theta1, theta2=theta2, usePE=usePE,
Nmax=Nmax, min.n2=min.n2, nsims=nsims,
# results
pBE=sum(BE)/nsims,
# Dec 2014 changed to
pBE_s1=sum(BE[ntot==n1])/nsims,
# was
#pBE_s1=sum(BE[stage==2])/nsims,
# Dec 2014 changed the meaning of pct_s2, was
#pct_s2=100*length(BE[stage==2])/nsims,
# whereby in case of unsymmetric alpha's stage 2 was also assigned
# if n2=0 but not BE using alpha[1] to fascilate BE test with alpha[2]
# now it is
pct_s2=100*sum(ntot>n1)/nsims,
# which simply means all those with n2>0
nmean=mean(ntot), nrange=range(ntot), nperc=quantile(ntot, p=npct)
#, ntot=ntot # experimental: return also all sample sizes
)
# table object summarizing the discrete distri of ntot
# only if usePE=FALSE or if usePE=TRUE then Nmax must be finite?
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 for the function
power.tsd <- power.2stage
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.