Power and sensitivity analyses [Regulatives / Guidelines]

posted by Helmut Homepage – Vienna, Austria, 2020-09-25 15:42 (1279 d 19:10 ago) – Posting: # 21939
Views: 4,079

Hi Pharma88,

in order to avoid surprises I recommend to perform a sensitivity analysis before designing the study (see also the vignette Power Analysis of PowerTOST).
In order to assess the impact of deviations from assumptions on power try this:

library(PowerTOST)
CV     <- 0.25  # assumed CV
theta0 <- 0.95  # assumed T/R-ratio
target <- 0.80  # target (desired) power
design <- "2x2" # any one given in known.designs()
# default BE limits: theta1 = 0.80, theta2 = 1.25
x      <- pa.ABE(CV = CV, theta0 = theta0,
                 targetpower = target, design = design)
plot(x, pct = FALSE, ratiolabel = "theta0")


[image]

Power depends on the CV, the T/R-ratio, and the final sample size of eligible subjects (dosed – dropouts). In each of the panels two parameters are kept constant and the third varied. We see that power is most sensitive to deviations of the T/R-ratio. Much less so to the CV followed by dropouts.

However, this is not the end of the story since potential deviations occur simultaneously. That’s a four-dimensional problem (power depends on theta0, CV, and n). A quick & dirty [image]-script at end.

[image]

The lower right quadrants of each panel show “nice” combinations (T/R-ratio > assumed and CV < assumed). Higher power than desired, great. NB sometimes you will see a statistically significant treatment effect. No worries about that (not clinically relevant).
The other combinations are tricky. Since power is most sensitive to the T/R-ratio, it would need a substantially lower CV to compensate for a worse T/R-ratio. Have a look at the 0.80 contour lines in the lower left quadrant of the first panel (no dropouts). Say, the T/R-ratio is just 0.92. Then with any CV > 0.2069 power will be below our target.
On the other hand, “better” T/R-ratios allow for higher CVs. That’s shown in the upper right quadrants. However, if the CV gets too large, even a T/R-ratio of 1 gives not the target power.
In the upper left quadrants are the worst case combinations (T/R-ratio < assumed and CV > assumed). It might still be possible to show BE though with a lower chance (power < 0.80).
Like in the Power Analysis above we see that dropouts don’t hurt that much.
Note that – since power curves are symmetrical in log-scale – you get the same power for \(\small{\theta_0}\) and \(\small{1/\theta_0}\).
With

sensitivity(CV.tgt = CV.tgt, do.rate = do.rate, theta0.lo = 0.8)

you see that power equals the nominal \(\small{\alpha}\) at the lower limit of the acceptance range.

But again, this should be be done before the study.
If you demonstrated BE with a post hoc power < target, all is good (answering the 3rd question of your post). If post hoc power is substantially lower than desired, you should reconsider your assumptions in designing the next study.

There is simple intuition behind results like these: If my car made it to the top of the hill, then it is powerful enough to climb that hill; if it didn’t, then it obviously isn’t powerful enough. Retrospective power is an obvious answer to a rather uninteresting question. A more meaningful question is to ask whether the car is powerful enough to climb a particular hill never climbed before; or whether a different car can climb that new hill. Such questions are prospective, not retrospective.
— Russell V. Lenth, Two Sample-Size Practices that I Don’t Recommend.



library(PowerTOST)
power.TOST.vectorized <- function(CV, theta0, ...) {
  # Supportive function, since only theta0 in power.TOST()
  # can be vectorized (not CV).
  # Returns a matrix of power with rows CV and columns theta0.

  power <- matrix(ncol = length(CV), nrow = length(theta0),
                  dimnames = list(c(paste0("theta0.", 1:length(theta0))),
                                  c(paste0("CV.", 1:length(CV)))))
  for (j in seq_along(CV)) {
    power[, j] <- suppressMessages( # for unbalanced cases
                    power.TOST(CV = CV[j], theta0 = theta0, ...))
  }
  return(power)
}
sensitivity <- function(alpha = 0.05, CV.tgt, CV.lo, CV.hi, theta0.tgt = 0.95,
                        theta0.lo, do.rate, target = 0.80, design = "2x2",
                        theta1, theta2, mesh = 25) {
  # alpha = 0.5 for assessing only the PE (Health Canada: Cmax)
  if (alpha <= 0 | alpha > 0.5)
    stop("alpha ", alpha, " does not make sense.")
  if (missing(CV.tgt))
    stop("CV.tgt must be given.")
  if (missing(CV.lo))
    CV.lo <- CV.tgt * 0.8
  if (missing(CV.hi))
    CV.hi <- CV.tgt * 1.2
  if (theta0.tgt >= 1)
    stop("theta0.tgt >=1 not implemented.")
  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 (missing(theta0.lo))
    theta0.lo <- theta0.tgt * 0.95
  if (theta0.lo < theta1) {
    message("theta0.lo ", theta0.lo, "< theta1 does not make sense. ",
            "Changed to theta1.")
    theta0.lo <- theta1
  }
  theta0.hi <- 1
  if (missing(do.rate))
    stop("do.rate must be given.")
  if (do.rate < 0)
    stop("do.rate", do.rate, " does not make sense.")
  if (target <= 0.5)
    stop("Target ", target, " does not make sense. Toss a coin instead.")
  if (target >= 1)
    stop("Target ", target, " does not make sense.")
  d.no <- PowerTOST:::.design.no(design)
  if (is.na(d.no))
    stop("design '", design, "' unknown.")
  if (mesh <= 10) {
    message("Too wide mesh is imprecise. Increased to 25.")
    mesh <- 25
  }
  CV      <- seq(CV.lo, CV.hi, length.out = mesh)
  theta0  <- seq(theta0.lo, theta0.hi, length.out = mesh)
  n       <- sampleN.TOST(alpha = alpha, CV = CV.tgt, theta0 = theta0.tgt,
                          theta1 = theta1, theta2 = theta2,
                          targetpower = target, design = design,
                          details = FALSE)[["Sample size"]]
  ns      <- n:floor(n * (1 - do.rate))
  windows(width = 6.5, height = 6.5)
  fig.col <- ceiling(sqrt(length(ns)))
  fig.row <- ceiling(length(ns)/fig.col)
  figs    <- c(fig.col, fig.row)
  op      <- par(no.readonly = TRUE)
  par(mar = c(3.5, 4, 0.1, 0.3))
  split.screen(figs)
  for (j in seq_along(ns)) {
    power <- power.TOST.vectorized(alpha = alpha, CV = CV, theta0 = theta0,
                                   n = ns[j], design = design)
    pwr   <- suppressMessages(
               power.TOST(alpha = alpha, CV = CV.tgt, theta0 = theta0.tgt,
                                   n = ns[j], design = design))
    screen(j)
    plot(theta0, CV, type = "n", xlab = "", ylab = "", cex.axis = 0.9, las = 1)
    axis(1, at = theta0.tgt, labels = FALSE)
    axis(2, at = CV.tgt, labels = FALSE)
    if (j %%figs[2] == 1) {       # y-label first columns
      mtext("CV", side = 2, line = 3)
    }
    if (j > prod(figs)-figs[2]) { # x-label last row
      mtext(expression(theta[0]), side = 1, line = 2.25)
    }
    grid(); box()
    pwrs <- as.vector(unlist(as.data.frame(power)))
    nl   <- length(pretty(pwrs, 20))
    clr  <- sapply(hcl.pals(type="sequential"),
                   hcl.colors, n = nl, rev = TRUE)[, "ag_Sunset"]
    contour(theta0, CV, power, col = clr, nlevels = nl, labcex = 0.8,
            labels = sprintf("%.2f", pretty(pwrs, nl)), add = TRUE)
    points(theta0.tgt, CV.tgt, cex = 1, pch = 21, col = "blue", bg = "#87CEFA")
    TeachingDemos::shadowtext(theta0.tgt, CV.tgt, col = "blue", bg = "white",
                              labels = paste0(signif(pwr, 3),
                                              " (n ", ns[j], ")"),
                              r = 0.25, adj = c(0.5, 1.6), cex = 0.8)
  }
  close.screen(all = TRUE)
  par(op)
}
#####################################################
# Specification of the study (mandatory values)     #
#####################################################

CV.tgt  <- 0.25 # assumed CV                        #
do.rate <- 0.10 # anticipated dropout-rate (10%)    #
#####################################################
# defaults (if not provided in named arguments)     #
#   alpha      = 0.05            common             #
#   CV.lo      = CV*0.80         best case          #
#   CV.hi      = CV*1.20         worst case         #
#   theta0.tgt = 0.95            assumed T/R-ratio  #
#   theta0.lo  = theta0.tgt*0.95 worst case         #
#   target     = 0.80            target power       #
#   theta1     = 0.80            lower BE limit     #
#   theta2     = 1.25            upper BE limit     #
#   design     = "2x2"           in known.designs() #
#   mesh       = 25              resolution         #
#####################################################

sensitivity(CV.tgt = CV.tgt, do.rate = do.rate)


Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

Complete thread:

UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,640 registered users;
84 visitors (0 registered, 84 guests [including 5 identified bots]).
Forum time: 09:53 CET (Europe/Vienna)

Nothing shows a lack of mathematical education more
than an overly precise calculation.    Carl Friedrich Gauß

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5