Nothing
# ----------------------------------------------------------------------------
# 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
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.