Log in
RegisterExact power and sample size, Part IV [🇷 for BE/BA]
Dear All!
Ok, ok, I will quit with this "Eierlegende Wollmilchsau".
The function .approx.power.TOST() not shown here is the code using non-central t-approximation for the power.
The whole code is released under the copy left (aka GPL).
Eventually there is someone out there to tying together all this snippets into a package, for christmas
.
Ok, ok, I will quit with this "Eierlegende Wollmilchsau".
# --------------------------------------------------------------------------
# Sample size for a desired power: see known.designs()
# for covered experimental designs
sampleN.TOST <- function(alpha=0.05, targetpower=0.8, logscale=TRUE,
ltheta1=0.8, ltheta2, ldiff=0.95, CV,
design="2x2", exact=TRUE, print=TRUE, details=FALSE)
{
#number of the design
d.no <- .design.no(design)
if (is.na(d.no)) stop("Err: design not known")
# design charcteristics
d.name <- .design.name(d.no)
# get the df for the design as an unevaluated expression
dfe <- .design.df(d.no,"Ns")
# get stepsize for search
steps <- .design.step(d.no)
# get design constant
bk <- .design.bk(d.no)
if (missing(CV)) stop("Err: CV must be given!")
# print the configuration:
if (print) {
cat("\n++++++++ Equivalence test - TOST ++++++++\n")
cat(" Sample size estimatio\n\n")
cat("-----------------------------------------\n")
cat("Study design: ",d.name,"\n")
if (details) {
cat("Design characteristics:\n")
cat("df = ",gsub("Ns","n",gsub(" ","",as.character(dfe))),
", design const = ",bk,", step = ",steps,"\n\n", sep="")
}
}
# handle the log transformation
if (logscale) {
if (missing(ltheta2)) ltheta2=1/ltheta1
if ( (ldiff<=ltheta1) | (ldiff>=ltheta2) ) {
cat("Err: ratio not between margins!\n")
return(NaN)
}
theta1 <- log(ltheta1)
theta2 <- log(ltheta2)
diffm <- log(ldiff)
se <- sqrt(log(1.+CV^2))
if (print) cat("log-transformed data (multiplicative model)\n\n")
} else {
if (missing(ltheta2)) ltheta2=-ltheta1
if ( (ldiff<=ltheta1) | (ldiff>=ltheta2) ) {
cat("Err: diff not between margins!\n")
return(NaN)
}
theta1 <- ltheta1
theta2 <- ltheta2
diffm <- ldiff
se <- CV
if (print) cat("untransformed data (additive model)\n\n")
}
if (print) {
cat("alpha = ",alpha,", target power = ", targetpower,"\n", sep="")
cat("BE margins =",ltheta1,"...", ltheta2,"\n")
if (logscale)
cat("Null (true) ratio = ",ldiff,", CV = ",CV,"\n", sep="")
else
cat("Null (true) diff. = ",ldiff,", CV = ",CV,"\n", sep="")
}
#start value from large sample approx. (hidden func.)
Ns <- .sampleN0(alpha, targetpower, theta1, theta2, diffm, se, steps, bk)
df <- eval(dfe)
if (exact) pow <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
else pow <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
if (details) {
cat("\nSample size search\n")
if (d.no == 0) cat("(n is sample size per group)\n")
cat(" n power\n")
cat( Ns," ", 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) {
Ns <- Ns+steps
iter <- iter+1
df <- eval(dfe)
if (exact)
pow <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
else
pow <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
if (details) cat( Ns," ", formatC(pow, digits=6, format="f"),"\n")
if (iter>50) break
}
} else {
while (pow>targetpower) {
if (Ns<=4) break # min. sample size
Ns <- Ns-steps # step down if start power is to high
iter <- iter+1
df <- eval(dfe)
if (exact)
pow <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
else
pow <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
if (details) cat( Ns," ", formatC(pow, digits=6),"\n")
if (iter>50) break
}
if ( pow<targetpower ) {
Ns <- Ns+steps #step up once to have n with pow>=target
df <- eval(dfe)
if (exact)
pow <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
else
pow <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)
}
}
if (print && !details) {
cat("\nSample size\n")
if (d.no == 0) cat("(n is sample size per group)\n")
cat(" n power\n")
cat( Ns," ", formatC(pow, digits=6, format="f"),"\n")
}
if (details && print) {
if (exact) cat("\nExact power calculation with\nOwen's Q functions.\n")
else cat("\nApproximate power calculation with\nnon-central t-distribution.\n")
}
# always print if approx.
if (print && !exact)
cat("\nApproximate power calculation with\nnon-central t-distribution.\n")
if (print) cat("\n")
res <- data.frame(design=design, alpha=alpha, CV=CV, theta0=ldiff,
theta1=ltheta1, theta2=ltheta2, n=Ns, 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)
}
The function .approx.power.TOST() not shown here is the code using non-central t-approximation for the power.
The whole code is released under the copy left (aka GPL).
Eventually there is someone out there to tying together all this snippets into a package, for christmas
.—
Regards,
Detlew
Regards,
Detlew
Complete thread:
- Exact power and sample size, Part I d_labes 2009-11-23 14:32 [🇷 for BE/BA]
![[Mix]](https://static.bebac.at/img/mix.png "open in Mix view")
- Exact power and sample size, Part II d_labes 2009-11-23 14:38
- Exact power and sample size, Part III d_labes 2009-11-23 14:50
- Exact power and sample size, Part IVd_labes 2009-11-23 14:57
- Backslashes in R-code Helmut 2009-11-23 15:13
- Backslashes eater d_labes 2009-11-23 15:19
- Backslashes eater Helmut 2009-11-23 15:28
- Backslashes eater d_labes 2009-11-23 15:19
- Backslashes in R-code Helmut 2009-11-23 15:13
- Exact power and sample size, Part IVd_labes 2009-11-23 14:57
- Exact power and sample size, Part III d_labes 2009-11-23 14:50
- Exact power and sample size, Part II d_labes 2009-11-23 14:38
Ing. Helmut Schütz