Nothing
#source(helper_dp.R)
# -----------------------------------------------------------------------
# large sample estimate of sample size
# -----------------------------------------------------------------------
sampleN0.dp <- function(alpha=0.05, CV, doses, targetpower=0.8, beta0=1,
theta1=0.8, theta2=1/theta1, design="crossover",
dm=NULL, CVb)
{
s2 <- CV2mse(CV)
rd <- max(doses)/min(doses)
grps <- length(doses)
prds <- ifelse(design=="parallel", 1, grps)
prds <- ifelse(design=="IBD",ncol(dm), prds)
grps <- ifelse(design=="IBD",nrow(dm), grps)
# all ni=1
n <- rep.int(1, times=grps)
# corrected sum of squares of log doses of design
css <- .css3(doses, design, dm, n, s02=CV2mse(CV), omega2=CV2mse(CVb))
# lower / upper acceptance range of slope
bl <- 1+log(theta1)/log(rd)
bu <- 1+log(theta2)/log(rd)
if (beta0<=bl | beta0>=bu){
warning("beta0 must lie between ",signif(bl,5)," and ",
signif(bu,5),". NaN returned.")
return(NaN)
}
# quantil for type2 error: U(power) or U(1-beta/2)
if (beta0!=1) Ub <- qnorm(targetpower) else Ub <- qnorm(1-(1-targetpower)/2)
# see Julious, Tan & Machin formula for beta0=1
# see also Sethuraman et al.
n <- s2*(Ub+qnorm(1-alpha))^2/css
if (beta0>=1) n <- n/(beta0-bu)^2 else n<- n/(beta0-bl)^2
n <- grps*round(n,0)
return(n)
}
# -----------------------------------------------------------------------
# sample size for dose prop., power model
# -----------------------------------------------------------------------
sampleN.dp <- function(alpha=0.05, CV, doses, targetpower=0.8, beta0,
theta1=0.8, theta2=1/theta1,
design=c("crossover", "parallel", "IBD"), dm=NULL,
CVb, print=TRUE, details=FALSE, imax=100)
{
design <- match.arg(design)
grps <- length(doses) # dose groups and periods in case of crossover
if (grps<=1) stop("At least two doses have to be given.")
if (design=="IBD"){
# check
if(!is.matrix(dm)) stop("Design matrix must be given.")
grps <- nrow(dm) # sequence groups
if(missing(CVb)) CVb <- 2*CV
} else {
if(missing(CVb)) CVb <- 0
}
if (CV<=0) stop("CV must be greater then zero.")
s2 <- CV2mse(CV)
# acceptance range for slope
rd <- max(doses)/min(doses)
bl <- 1+log(theta1)/log(rd)
bu <- 1+log(theta2)/log(rd)
if (missing(beta0)) beta0 <- 1+log(0.95)/log(rd)
if (beta0<=0) stop("beta0 must be greater then zero.")
if (beta0<=bl | beta0>=bu) stop("beta0 must lie between ",signif(bl,5)," and ",
signif(bu,5),".")
if (theta1<=0 | theta2<=0) stop("theta1/theta2 must be greater then zero.")
if (print){
cat("\n++++ Dose proportionality study, power model ++++\n")
cat(" Sample size estimation\n")
cat("-------------------------------------------------\n")
cat("Study design: ",design, sep="")
descomm <- paste0(" (",grps,"x",grps," Latin square)")
if (design=="parallel"){
descomm <- paste0(" (",grps," groups)")
}
if (design=="IBD"){
descomm <- paste0(" (",length(doses),"x",grps,"x",ncol(dm),")")
}
cat(descomm,"\n")
cat("alpha = ",alpha,", target power = ", targetpower,"\n", sep="")
cat("Equivalence margins of R(dnm) =",theta1,"...", theta2,"\n")
cat("Doses = "); cat(doses,"\n")
cat("True slope = ",beta0,", CV = ", CV, sep="")
if(design=="IBD"){
cat(", CVb = ", CVb, "\n", sep="")
} else {
cat("\n")
}
cat("Slope acceptance range =",signif(bl,5),"...", signif(bu,5),"\n")
}
# start sample size search
n <- sampleN0.dp(alpha=alpha, CV=CV, doses=doses, targetpower=targetpower,
beta0=beta0, theta1=theta1, theta2=theta2, design=design,
dm=dm, CVb=CVb)
# here for small CV < 0.1 n=0 may result
# is this reasonable?
if (n<grps) n <- grps
if (grps==2 & n<4) n <- 4 # else df=0 may result for "crossover" or "parallel"
pwr <- power.dp(alpha=alpha, CV=CV, doses=doses, n=n, beta0=beta0,
theta1=theta1, theta2=theta2, design=design, dm=dm, CVb=CVb)
if (details) {
cat("\nSample size search (ntotal)\n")
cat(" n power\n")
# do not print first too high
# this is for cases with only one step-down and than step up
if (pwr<=targetpower) cat( n," ", formatC(pwr, digits=6, format="f"),"\n")
}
iter <- 0
nmin <- 3
while(pwr>targetpower){
if (n<=nmin) {
if (details & iter==0) cat( n," ", formatC(pwr, digits=6, format="f"),"\n")
break
}
# step down
n <- n-grps
pwr <- power.dp(alpha=alpha, CV=CV, doses=doses, n=n, beta0=beta0,
theta1=theta1, theta2=theta2, design=design, dm=dm, CVb=CVb)
iter <- iter+1
if (details) cat( n," ", formatC(pwr, digits=6),"\n")
if (iter>imax) break
}
while(pwr<targetpower){
# step up
n <- n+grps
pwr <- power.dp(alpha=alpha, CV=CV, doses=doses, n=n, beta0=beta0,
theta1=theta1, theta2=theta2, design=design, dm=dm, CVb=CVb)
if (details) cat( n," ", formatC(pwr, digits=6),"\n")
iter <- iter+1
if (iter>imax) break
}
if(print && !details){
cat("\nSample size (total)\n")
cat(" n power\n")
cat( n," ", formatC(pwr, digits=6, format="f"), "\n")
}
#return results as data.frame
res <- data.frame(design=design, alpha=alpha, CV=CV, CVb=CVb,
doses=paste(doses, collapse=", "), beta0=beta0,
theta1=theta1, theta2=theta2, n=n, power=pwr,
targetpower=targetpower)
names(res) <-c("Design", "alpha", "CV", "CVb", "doses", "beta0", "theta1", "theta2",
"Sample size", "Achieved power", "Target power")
if (design!="IBD") res$CVb <- NULL
if (print) return(invisible(res)) else 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.