Nothing
# helper function to get CV from 'required sample size'
reqN2CV <- function(alpha=0.05, targetpower=0.8, theta0=1., theta1=0.8, n)
{
z1 <- qnorm(1-alpha)
if (log(theta0)!=0) z2 <- qnorm(targetpower) else z2 <- qnorm(1-(1-targetpower)/2)
s2e <- n*(log(theta0)-log(theta1))^2/2/(z1+z2)^2
return(mse2CV(s2e))
}
# ----------------------------------------------------------------------------
# power (or alpha) of 2-stage studies with (blinded) sample size re-estimation
# see Golkowski et al.
#
# author D.L.
# ----------------------------------------------------------------------------
# require(PowerTOST)
# source("C:/Users/dlabes/workspace/Power2Stage/R/sampsiz2.R")
# source("C:/Users/dlabes/workspace/Power2Stage/R/sampsiz_n0.R")
# source("C:/Users/dlabes/workspace/Power2Stage/R/power.R")
power.2stage.ssr <- function(alpha=0.05, n1, GMR, CV, targetpower=0.8,
pmethod=c("nct","exact", "shifted","ls"),
blind=FALSE, usePE=FALSE, min.n=0, max.n=Inf,
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])
# seems to give errors if alpha is a vector
alpha <- alpha[1L]
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 (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>max.n) stop("max.n < n1 doestn't make sense!")
if (min.n!=0 & min.n<n1) stop("min.n < n1 doestn't make sense!")
if(missing(nsims)){
nsims <- 1E5
if(theta0<=theta1 | theta0>=theta2) nsims <- 1E6
}
# check if power calculation method is nct, exact, shifted or large sample
pmethod <- match.arg(pmethod)
if(blind & usePE) {
usePE=FALSE
warning("usePE=TRUE does not make sense for blinded SSR. Reset to usePE=FALSE.")
}
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)
# ----- interim ----------------------------------------------------------
# simulate blinded est. of variance
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
#s2os <- mses + pes^2/2 # is this correct? example shows: it's not correct.
# s2os is smaller than mse! I would expect
# the other way round!
# SS of tmt formula is only valid if sequence balanced
# res. SS + SS of tmt
s2os <- (df*mses + n1*pes^2/2)/(n1-1)
# unblind
if (!blind) s2os <- mses
if(details){
# time for stage 1 sims
cat(" - Time consumed (secs):\n")
print(round((proc.time()-ptm),2))
cat("Keep calm. ")
}
# ------ recalculate sample size -----------------------------------------
ptms <- proc.time()
# first calc. power to avoid unnecessary time consuming sample size estimation
# this is not really part of the method (?) but a run-time goody
# will give some boost for small CV's and/or high n1
if(pmethod!="ls"){
# use GMR or pe1 in sample size re-est.
lpes <- lGMR
if(usePE) lpes <- pes
pwr <- .calc.power(alpha=alpha, ltheta1=ltheta1, ltheta2=ltheta2,
diffm=lpes, sem=sqrt(bk*s2os/n1), df=df, method=pmethod)
if(details){
cat("Sample sizes (", sum(pwr<targetpower),
" studies) will be re-estimated.\n", sep="")
cat("May need some time.\n")
}
} else {
# don't have power func for large sample approx.
# moreover ssr is so fast that we don't need the power step
pwr <- rep(0, times=nsims)
if(details){
cat("Sample sizes will be re-estimated.\n")
cat("May need some time.\n")
}
}
# total sample size
ntot <- rep(n1, times=nsims)
mse_tmp <- s2os[pwr<targetpower]
pes_tmp <- pes[pwr<targetpower]
# use GMR or pe1 in sample size re-est.
lpes <- lGMR
if(usePE) lpes <- pes_tmp
if(pmethod=="ls"){
# large sample approx. via normal distribution
nt <- .sampleN00(alpha=alpha, targetpower=targetpower, se=sqrt(mse_tmp),
diffm=lpes, ltheta1=ltheta1, ltheta2=ltheta2, bk=2,
steps=2, diffmthreshold=0.0)
} else {
# exact or approx. via t-distri
# Aug. 2017: .sampleN2() now uses N-3 as df in SSR due to extra term 'stage'
# TODO: Check df for sample size re-estimation
# I think N-2 is here appropriate, but Golkowski only uses pmethod=="ls"
nt <- .sampleN2(alpha=alpha, targetpower=targetpower, ltheta0=lpes,
mse=mse_tmp, ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod, bk=2, dfc="n-2")
}
# maybe we have enough power in all cases and thus no re-estimated sample size
if(length(nt)>0) ntot[pwr<targetpower] <- nt
# take care of memory
rm(nt, pwr, mse_tmp, lpes)
# sample size returns Inf if pe outside acceptance range, then stay with n1
# but this should not occure here since pe1 is not used
ntot <- ifelse(is.finite(ntot), ntot, n1)
# use max.n if nt > max.n
ntot <- ifelse(ntot>max.n, max.n, ntot)
# use min.n if nt < min.n
ntot <- ifelse(ntot<min.n, min.n, ntot)
n2 <- ntot-n1
# do not fall below n1
n2 <- ifelse(n2<0, 0, n2)
if(details){
cat("Time consumed (secs):\n")
print(round((proc.time()-ptms),1))
}
# ---------- stage 2 evaluation --------------------------------------
m1 <- pes
SS1 <- (n1-2)*mses
# to avoid warnings for n2=0 in rnorm() and rchisq()
ow <- options("warn")
options(warn=-1)
m2 <- ifelse(n2>0, rnorm(n=nsims, mean=mlog, sd=sqrt(mse*bk/n2)), 0)
# ??? (n2-2) cancels out!
SS2 <- ifelse(n2>2, (n2-2)*mse*rchisq(n=nsims, df=n2-2)/(n2-2), 0)
# reset options
options(ow)
ntot <- n1+n2
# do we need this 'stage' contribution?
# can't discover such in Golkowski et al., but ...
withStage <- TRUE
if(withStage){
SSmean <- ifelse(n2>0, (m1-m2)^2/(2/n1+2/n2), 0)
df2 <- ifelse(n2>0, ntot-3, n1-2)
} else {
SSmean <- 0
df2 <- ntot-2
}
pe2 <- ifelse(n2>0, (n1*m1+n2*m2)/ntot, pes)
mse2 <- ifelse(n2>0, (SS1+SSmean+SS2)/df2, mses)
# take care of memory
rm(m1, m2, SS1, SS2, SSmean)
# calculate CI and decide BE
hw <- qt(1-alpha, df2)*sqrt(mse2*bk/ntot)
lower <- pe2 - hw
upper <- pe2 + hw
BE <- (lower>=ltheta1) & (upper<=ltheta2)
# take care of memory
rm(lower, upper, hw)
# the return list
res <- list(method="SSR",
alpha=alpha, CV=CV, n1=n1, GMR=GMR, targetpower=targetpower,
pmethod=pmethod, theta0=theta0, theta1=theta1, theta2=theta2,
usePE=usePE, max.n=max.n, min.n=min.n, blind=blind, nsims=nsims,
# results
pBE=sum(BE)/nsims,
#pBE_s1 ?
pct_s2=100*length(ntot[ntot>n1])/nsims,
nmean=mean(ntot), nrange=range(ntot), nperc=quantile(ntot, p=npct))
# return a table object as summary of ntot distribution
# only if usePE==FALSE ? or always?
if(!usePE) res$ntable <- table(ntot)
if (details){
cat("Total time consumed (secs):\n")
print(round((proc.time()-ptm),1))
cat("\n")
}
# output is now via S3 print method
class(res) <- c("pwrtsd", "list")
return(res)
} #end function
# alias of the function
power.tsd.ssr <- power.2stage.ssr
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.