Nothing
# --------------------------------------------------------------------------
# power (or alpha) of 2-stage studies according to Potvin et al. methods "B"
# with 2 PK metrics
#
# Author D.L.
# --------------------------------------------------------------------------
power.tsd.2m <- function(alpha=c(0.0294,0.0294), CV, n1, rho=0, GMR,
targetpower=0.8, pmethod=c("nct","exact", "shifted"),
powerstep=FALSE, theta0, theta1, theta2,
npct=c(0.05, 0.5, 0.95), nsims, setseed=TRUE, details=FALSE)
{
if(missing(CV)) stop("CVs must be given.")
else {
if(length(CV)==1) CV <- rep(CV,2)
if(length(CV)!=2) stop("GMR must have 2 elements.")
if(any(CV<=0)) stop("CVs must be >0.")
}
if(missing(n1)) stop("Number of subjects in stage 1 must be given.")
if(length(n1)!=1) {
n1 <- n1[1]
message("n1 must be scalar. First element of n1 used.")
if(n1<=0) stop("Number of subjects in stage 1 must be >0.")
}
if(length(alpha)==1) {
alpha <- rep(alpha, 2)
message("Scalar alpha used for both stages.")
}
if(length(alpha) != 2) stop("alpha must have two elements.")
if(missing(GMR)) GMR <- rep(0.95, 2)
else {
if(length(GMR)==1) GMR <- rep(GMR, 2)
if(length(GMR)!=2) stop("GMR must have 2 elements.")
}
if (missing(theta1) & missing(theta2)) theta1 <- 0.8
if (!missing(theta1) & missing(theta2)) theta2 <- 1/theta1
if (missing(theta1) & !missing(theta2)) theta1 <- 1/theta2
stopifnot(length(theta1)==1, length(theta2)==1)
if(any(GMR<=theta1) | any(GMR>=theta2)) stop("GMRs must be within acceptance range.")
if (missing(theta0)) theta0 <- GMR
else {
if(length(theta0)==1) theta0 <- rep(theta0, 2)
if(length(theta0)!=2) stop("theta0 must have two elements.")
}
#if(rho!=0) warning("rho != 0 is only experimental.", call. = FALSE)
stopifnot(length(rho)==1, rho >= -1, rho <= 1)
if (rho < -1 || rho > 1)
stop("Correlation must be within {-1, +1}.")
# prevent error in rWishart(n, df, Sigma)
if (rho == -1) rho <- -1 + .Machine $double.eps
if (rho == +1) rho <- +1 - .Machine $double.eps
if(missing(nsims)){
nsims <- 1E5
if(any(theta0<=theta1) | any(theta0>=theta2)) nsims <- 1E6
}
# 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 <- matrix(0, nrow=nsims, ncol=2)
mses <- matrix(0, nrow=nsims, ncol=2)
if (rho==0){
pes[ , 1] <- rnorm(n=nsims, mean=mlog[1], sd=sdm[1]) # metric 1, f.i. AUC
pes[ , 2] <- rnorm(n=nsims, mean=mlog[2], sd=sdm[2]) # metric 2, f.i. Cmax
# simulate mse via chi-squared distribution
mses[ , 1] <- mse[1]*rchisq(n=nsims, df=df)/df
mses[ , 2] <- mse[2]*rchisq(n=nsims, df=df)/df
} else {
# multivariate normal with rho
# TODO: check this
#browser()
# (population) variance-covariance matrix of means
s_m <- diag(sdm^2)
s_m[1,2] <- s_m[2,1] <- rho*s_m[1,1]*s_m[2,2]
pes <- rmvnorm(nsims, mean=mlog, sigma=s_m)
# simulate covariance matrices via Wishart distribution
s_mse <- diag(mse)
s_mse[1,2] <- s_mse[2,1] <- rho*sqrt(s_mse[1,1]*s_mse[2,2])
#or do we have here df=n-1?
covm <- rWish2(n=nsims, df=df, Sigma=s_mse)/df
mses[, 1] <- covm[1, 1, ]
mses[, 2] <- covm[2, 2, ]
# take care of memory
rm(covm)
}
BE <- function(nu)
{
hw <- tval*sqrt(Cfact*mses[, nu])
loCL <- pes[, nu] - hw
upCL <- pes[, nu] + hw
BE <- loCL >= ltheta1 & upCL <= ltheta2
BE
}
# make BE decision for stage 1
BE_m1 <- BE(nu=1)
BE_m2 <- BE(nu=2)
# overall BE in stage 1
BE <- BE_m1 & BE_m2
#browser()
# set BE to NA if BE not decided yet
BE[!BE] <- NA
if(powerstep) {
# metric 1 power
pwr_m1 <- .calc.power(alpha=alpha[2], ltheta1=ltheta1, ltheta2=ltheta2,
diffm=lGMR[1], sem=sqrt(bk*mses[, 1]/n1), df=df,
method=pmethod)
BE[BE_m2 & !BE_m1 & pwr_m1>=targetpower] <- FALSE
# metric 2 power
pwr_m2 <- .calc.power(alpha=alpha[2], ltheta1=ltheta1, ltheta2=ltheta2,
diffm=lGMR[2], sem=sqrt(bk*mses[, 2]/n1), df=df,
method=pmethod)
BE[BE_m1 & !BE_m2 & pwr_m2>=targetpower] <- FALSE
# take care of memory
rm(pwr_m1, pwr_m2)
}
# take care of memory
rm(BE_m1, BE_m2)
# time for stage 1
if(details){
cat(" - Time consumed (secs):\n")
print(round((proc.time()-ptm),1))
}
# BE == TRUE/FALSE is decided yet
# BE == NA not yet decided
ntot <- rep(n1, nsims)
if(sum(is.na(BE)) > 0) {
ptms <- proc.time()
ind <- is.na(BE)
pes <- pes[ind, ]
mses <- mses[ind, ]
rm(ind)
# ------sample size for stage 2 -----------------------------------------
if(GMR[1]==GMR[2]){
if(details){
cat("Keep calm. Sample sizes for stage 2 (", sum(is.na(BE)),
" studies)\n", sep="")
cat("will be estimated. May need some time.\n")
}
# use the max. mse for sample size
mses_max <- pmax(mses[, 1], mses[, 2])
nt_max <- .sampleN2(alpha=alpha[2], targetpower=targetpower, ltheta0=lGMR[1],
mse=mses_max, ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod)
n2 <- ifelse(nt_max>n1, nt_max - n1, 0)
# keep care of memory
rm(mses_max, nt_max)
} else {
if(details){
cat("Keep calm. Sample sizes for stage 2 (2x ", sum(is.na(BE)),
" studies)\n", sep="")
cat("will be estimated. May need some time.\n")
}
# metric 1
nt_m1 <- .sampleN2(alpha=alpha[2], targetpower=targetpower, ltheta0=lGMR[1],
mse=mses[, 1], ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod)
n2_m1 <- ifelse(nt_m1>n1, nt_m1 - n1, 0)
# same for metric 2
nt_m2 <- .sampleN2(alpha=alpha[2], targetpower=targetpower, ltheta0=lGMR[2],
mse=mses[, 2], ltheta1=ltheta1, ltheta2=ltheta2,
method=pmethod)
n2_m2 <- ifelse(nt_m2>n1, nt_m2 - n1, 0)
n2 <- pmax(n2_m1, n2_m2 )
rm(nt_m1, n2_m1, nt_m2, n2_m2)
}
if(details){
cat("Time consumed (secs):\n")
print(round((proc.time()-ptms),1))
}
#browser()
# ------stage 2 evaluation ----------------------------------------------
# simulate stage 2 and evalute the combined data from stage 1 + 2
nsim2 <- nrow(pes)
pes2 <- matrix(0, ncol=2, nrow=nsim2 )
SS2 <- matrix(0, ncol=2, nrow=nsim2 )
ow <- options("warn")
options(warn=-1)
#browser()
if (rho!=0){
# we are simulating for mean=0, sigma = matrix 1 rho, rho 1
sigma2 <- diag(1,2,2)
sigma2[1,2] <- sigma2[2,1] <- rho
pes2 <- rmvnorm(nsim2, mean=rep(0,2), sigma=sigma2)
# now we transform to the actual variance and mean
# TODO: is this correct? Or have we to include rho in some way
# for the transformation?
pes2[, 1] <- ifelse(n2>0, pes2[ ,1]*sqrt(mse[1]*bk/n2) + mlog[1], 0)
pes2[, 2] <- ifelse(n2>0, pes2[ ,2]*sqrt(mse[2]*bk/n2) + mlog[2], 0)
# SS2 via Wishart distribution
#browser()
covm <- rWish2(n=nsim2, df=n2-2, Sigma=s_mse)
SS2[, 1] <- covm[1, 1, ]
SS2[, 2] <- covm[2, 2, ]
} else {
# pe's of stage 2 data via independent normal distri's
pes2[,1] <- ifelse(n2>0, rnorm(n=nsim2, mean=mlog[1], sd=sqrt(mse[1]*bk/n2)), 0)
pes2[,2] <- ifelse(n2>0, rnorm(n=nsim2, mean=mlog[2], sd=sqrt(mse[2]*bk/n2)), 0)
# SS2 via independent chi-squared distribution
SS2[, 1] <- ifelse(n2>2, (n2-2)*mse[1]*rchisq(n=nsim2, df=n2-2)/(n2-2), 0)
SS2[, 2] <- ifelse(n2>2, (n2-2)*mse[2]*rchisq(n=nsim2, df=n2-2)/(n2-2), 0)
}
options(ow)
BE2 <- function(nu)
{
# this function sees pes[], mses[], pes2[], SS2[] and n2[]
# nu is number of metric
#browser()
m1 <- pes[, nu]
SS1 <- (n1-2)*mses[, nu]
nsim2 <- length(m1)
m2 <- pes2[, nu]
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, m1)
mse2 <- ifelse(n2>0, (SS1+SSmean+SS2[, nu])/df2, mses[, nu])
# take care of memory
rm(m1, m2, SS1, 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
BE2
}
# browser()
BE2_m1 <- BE2(nu=1)
BE2_m2 <- BE2(nu=2)
# combine stage 1 & stage 2
ntot[is.na(BE)] <- n2 + n1
BE[is.na(BE)] <- BE2_m1 & BE2_m2
}
# the return list
res <- list(design="2x2 crossover",
method="B2m", alpha=alpha, CV=CV, n1=n1, GMR=GMR, rho=rho,
targetpower=targetpower, pmethod=pmethod,
theta0=exp(mlog), theta1=theta1, theta2=theta2, usePE=FALSE,
nsims=nsims,
# results
pBE=sum(BE)/nsims,
pBE_s1=sum(BE[ntot==n1])/nsims,
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
)
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
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.