Log in
RegisterPower and sample size with 'uncertain' CV - part II [Power / Sample Size]
Stephen Senn in Statistical issues in drug development: "The sample size calculation is an excuse for a sample size and not a reason"
Dear All!
Here part II of the story: Sample size
BTW: If you copy this code take extremely care of the backslash-eating devil lurking round the corner
.
Edit: Code corrected according to an Erratum of D Labes in a later post to allow easy copying. [Helmut]
Dear All!
Here part II of the story: Sample size
#------------------------------------------------------------------------------
# helper function to calculate std err from CV of lognormal data
.CV2se <- function(CV) return(sqrt(log(1.+CV^2)))
# reverse helper function
.se2CV <- function(se) return(sqrt(exp(se*se)-1))
#------------------------------------------------------------------------------
# sample size start for Joulious "expected" power
.exp.sampleN0 <- function(alpha=0.05, targetpower=0.8, ltheta1, ltheta2, diffm,
se, dfse, steps=2, bk=2)
{
Z1 <- qnorm(1-alpha)
if (abs(diffm)>0.0001) tinv <- qt(targetpower, dfse, Z1) else
tinv <- qt(1-(1-targetpower)/2, dfse, Z1)
# is factor 2 in julious = bk?
n01 <- bk*(se*tinv/(ltheta1-diffm))^2
n02 <- bk*(se*tinv/(ltheta2-diffm))^2
n0 <- ceiling(max(n01,n02))
#make an even multiple of step (=2 in case of 2x2 cross-over)
n0 <- steps*trunc(n0/steps)
if (n0<4) n0 <- 6 # minimum sample size
return(n0)
}
#-------------------------------------------------------------------------
# Sample size for a desired "expected" power according to Julious:
# see known.designs() for covered experimental designs
# Only for log-transformed data
# leave upper BE margin (ltheta2) empty and the function will use 1/lower
# CV and dfCV can be vectors, if then a pooled CV, df will be calculated
exp.sampleN.TOST <- function(alpha=0.05, targetpower=0.8,
theta1=0.8, theta2, diff=0.95, CV, dfCV,
design="2x2", print=TRUE, details=FALSE)
{
#number of the design and check if design is implemented
d.no <- .design.no(design)
if (is.na(d.no)) stop("Err: design not known", call.=FALSE)
# design charcteristics
ades <-.design.props(d.no)
d.name <- ades$name # nice name of design
# get the df for the design as an unevaluated expression
dfe <- parse(text=ades$df,srcfile=NULL)
steps <- ades$steps # stepsize for sample size search
bk <- ades$bk # get design constant
if (missing(CV) | missing(dfCV)) {
stop("Err: CV and df must be given!", call.=FALSE)
}
# calculate pooled data if CV and dfCV are vectors
if (length(CV)>1){
if (length(dfCV)!=length(CV)) {
stop("Err: CV and df must have equal number of entries!", call.=FALSE)
}
dfse <- sum(dfCV)
CVp <- .CV2se(CV)^2
CVp <- CVp * dfCV
CVp <- sum(CVp)/dfse
CVp <- .se2CV(CVp)
} else {
dfse <- dfCV
CVp <- CV
}
# print the configuration:
if (print) {
cat("\n++++++++ Equivalence test - TOST ++++++++\n")
cat(" Sample size est. with uncertain CV\n")
cat("-----------------------------------------\n")
cat("Study design: ",d.name,"\n")
if (details) {
cat("Design characteristics:\n")
cat("df = ",ades$df,", design const. = ",bk,
", step = ",steps,"\n\n",sep="")
}
}
if (missing(theta2)) theta2=1/theta1
if ( (diff<=theta1) | (diff>=theta2) ) {
stop("Err: ratio not between margins!", call.=FALSE)
}
ltheta1 <- log(theta1)
ltheta2 <- log(theta2)
diffm <- log(diff)
se <- sqrt(log(1.+CVp^2))
if (print) cat("log-transformed data (multiplicative model)\n\n")
if (print) {
cat("alpha = ",alpha,", target power = ", targetpower,"\n", sep="")
cat("BE margins =",theta1,"...", theta2,"\n")
cat("Null (true) ratio = ",diff,"\n", sep="")
if (length(CV)>1){
cat("Variability data\\n")
print(data.frame(CV=CV,df=dfCV), row.names = FALSE)
cat("CV(pooled) = ", CVp, " with ", dfse," df\n", sep="")
} else cat("CV = ", CVp, " with ", dfse," df\n", sep="")
}
#start value from large sample approx.
n <- .exp.sampleN0(alpha, targetpower, ltheta1, ltheta2, diffm,
se, dfse, steps, bk)
df <- eval(dfe)
pow <- .exp.power.TOST(alpha, ltheta1, ltheta2, diffm, se, dfse, n=n, df, bk)
if (details) {
cat("\\nSample size search\n")
if (d.no == 0) cat("(n is sample size per group)\n") #parallel group design
cat(" n exp. power\n")
cat( n," ", formatC(pow, digits=6, format="f"),"\n")
}
# --- loop until power >= target power
iter <- 0
# iter>50 is emergency brake
# this is eventually not necessary, depends on quality of sampleN0
# in experimentation I have seen max of six steps
if (pow<targetpower) {
while (pow<targetpower) {
n <- n+steps
iter <- iter+1
df <- eval(dfe)
pow <- .exp.power.TOST(alpha, ltheta1, ltheta2, diffm, se, dfse,
n=n, df, bk)
if (details) cat( n," ", formatC(pow, digits=6, format="f"),"\n")
if (iter>50) break
}
} else {
while (pow>targetpower) {
if (n<=6) break # min. sample size
n <- n-steps # step down if start power is to high
iter <- iter+1
df <- eval(dfe)
pow <- .exp.power.TOST(alpha, ltheta1, ltheta2, diffm, se, dfse,
n=n, df, bk)
if (details) cat( n," ", formatC(pow, digits=6),"\n")
if (iter>50) break
}
if ( pow<targetpower ) {
n <- n+steps #step up once to have n with pow>=target
df <- eval(dfe)
pow <- .exp.power.TOST(alpha, ltheta1, ltheta2, diffm, se, dfse,
n=n, df, bk)
}
}
if (print && !details) {
cat("\nSample size\n")
if (d.no == 0) cat("(n is sample size per group)\n") #parallel group
cat(" n exp. power\\n")
cat( n," ", formatC(pow, digits=6, format="f"),"\n")
}
if (print) cat("\n")
#return results as data.frame
res <- data.frame(design=design, alpha=alpha, CV=CV, theta0=diff,
theta1=theta1, theta2=theta2, n=n, power=pow,
targetpower=targetpower)
names(res) <-c("Design","alpha","CV","theta0","theta1","theta2",
"Sample size", "Achieved power", "Target power")
if (print) return(invisible(res))
else return(res)
}
BTW: If you copy this code take extremely care of the backslash-eating devil lurking round the corner
.Edit: Code corrected according to an Erratum of D Labes in a later post to allow easy copying. [Helmut]
Complete thread:
- Power and sample size with 'uncertain' CV d_labes 2010-04-07 10:04
![[Mix]](https://static.bebac.at/img/mix.png "open in Mix view")
- Power and sample size with 'uncertain' CV - part IId_labes 2010-04-07 10:41
- Helper for rescue d_labes 2010-04-07 14:43
- Power and sample size with 'uncertain' CV yjlee168 2010-04-08 09:10
- Validation necessary d_labes 2010-04-08 10:04
- Scarcely validation d_labes 2010-04-08 13:36
- Power and sample size with 'uncertain' CV - part IId_labes 2010-04-07 10:41
Ing. Helmut Schütz