Desperate reader [Regulatives / Guidelines]

posted by Helmut Homepage – Vienna, Austria, 2021-08-27 17:51 (1202 d 23:09 ago) – Posting: # 22550
Views: 5,971

Dear all (and esp. ElMaestro),

with great interest I read two recent papers.1,2

On [image] the author posted:
I absolutely agree with this statement.
However, what can we do? We need the bloody confidence interval. At the end of the day we are [image] left out in the rain if we are dealing with a partial replicate design and the study has to be evaluated for average bioequivalence (according to the [image]’s wishful thinking if \(\small{s_\textrm{wR}<0.294}\)).
To obtain the variance components which are estimable and the point estimate is just not enough. As long as we can’t get the standard error of the treatment difference and its associated degrees of freedom, we are at a loss. Hence, I repeat my mantra:

 Avoid the partial replicate design 
 if intended for RSABE and if 
 there is even the slightest chance that you have to assess it for ABE. 


Of note: I used the supplementary [image] script and got for

res <- Run.REML()
PE  <- 100 * exp(res$Fixed.eff[["TrtT", 1]] - res$Fixed.eff[["TrtR", 1]])
print(res); cat("PE", PE, "\n")

$Outcome
[1] "The model converged."

$LogLikelihood
[1] 15.33728

$Fixed.eff
              [,1]
TrtT    8.02569184
TrtR    8.00330042
Per1   -0.04811877
Per2   -0.04682231
SeqRTR  0.00385480
SeqRRT -0.05092437

$Random.eff
  Var.Component  Ini.guess REML.Estimate
1          varT 0.06841615    0.07049570
2         varBR 0.03576077    0.03621596
3         varWR 0.01240137    0.01324664
4        varBRT 0.03576077    0.04563785

PE 102.2644


Pretty close to the modified script we discussed last year (given at the end). Much faster than mine. :-D

Nelder-Mead optimization (31.6 seconds)
  LogLikelihood 75.98722
  Evaluations   239
  Successful completion.

   Component    Initial   Estimate
   sigma²(T) 0.06841615 0.07049451
  sigma²(bR) 0.03576077 0.03621602
  sigma²(wR) 0.01240137 0.01324648
 covar(bTbR) 0.05208846 0.04563728

 Parameter  Estimate
  sigma(T) 0.2655080
     CV(T)   27.026%
 sigma(bR) 0.1903051
    CV(bR)   19.204%
 sigma(wR) 0.1150934
    CV(wR)   11.548%
        PE 102.2644%


[image]




  1. Fuglsang A. Professor Endrenyi’s Legacy: An Evaluation of the Regulatory Requirement “Fixed Effects, Rather Than Random Effects, Should Be Used for All Terms”. J Pharm Pharm Sci. 2021; 24: 413–20. doi:10.18433/jpps31872. [image] Open access.
  2. Fuglsang A. Mitigation of the convergence issues associated with semi-replicated bioequivalence data. Pharm Stat. 2021. Early View. doi:10.1002/pst.2142.


###########################
# Section 1: Housekeeping #
###########################
library(PowerTOST)
library(TeachingDemos)
GetData <- function(path, sep, dec) {
  # read data from a CSV-file
  if (missing(sep)) sep <- ","
  if (!sep %in% c(",", ";", "t"))
    stop("Column separator must be any of ',', ';', 't'.")
  if (sep == "t") sep <- "\t"
  if (missing(dec)) dec <- "."
  if (!dec %in% c(".", ","))
    stop("Decimal separator must be '.', or ','.")
  # strip quotes and skip eventual commentary lines
  D        <- read.csv(path, sep = sep, dec = dec,
                       quote = "", comment.char = "#")
  names(D) <- c("Subj", "Per", "Seq", "Trt", "Y")
  cols     <- c("Subj", "Per", "Seq", "Trt")
  D[cols]  <- lapply(D[cols], factor)
  return(invisible(D))
}

CreateX <- function(D) {
  # first the two treatments
  TrtT <- as.numeric(as.character(D$Trt) == "T")
  TrtR <- as.numeric(as.character(D$Trt) == "R")
  X <- cbind(TrtT, TrtR)
  R <- qr(X)$rank
  # now the periods
  for (p in unique(D$Per)) {
    v  <- as.numeric(D$Per == p)
    XX <- data.frame(X, v)
    names(XX)[ncol(XX)] <- paste0("P", p)
    rnk <- qr(XX)$rank
    if (rnk > R) {
      X <- XX
      R <- rnk
    }
  }
  for (q in unique(D$Seq)) {
    v  <- as.numeric(D$Seq == q)
    XX <- data.frame(X, v)
    names(XX)[ncol(XX)] <- paste0("Q", q)
    rnk <- qr(XX)$rank
    if (rnk > R) {
      X <- XX
      R <- rnk
    }
  }
  return(as.matrix(X))
}
Create.CovM <- function(Params) {
  # block diagonal covariance matrix
  varT  <- Params[1]
  varBR <- Params[2]
  varWR <- Params[3]
  varRT <- Params[4]
  Nobs  <- length(D$Y)
  V     <- matrix(0, ncol = Nobs, nrow = Nobs)
  for (iRow in 1:Nobs) {
    for (iCol in 1:Nobs) {
      if (iCol == iRow) { # diagonal elements
        if (D$Trt[iRow] == "T") {
          V[iRow, iCol] <- V[iRow, iCol] + varT
        } else {
          V[iRow, iCol] <- V[iRow, iCol] + varWR + varBR
        }
      } else {            # off diagonal elements
        if (D$Subj[iRow] == D$Subj[iCol]) {
          if (D$Trt[iCol] == D$Trt[iRow]) {
            V[iRow, iCol] <- V[iRow, iCol] + varBR
          } else {
            V[iRow, iCol] <- V[iRow, iCol] + varRT
          }
        }
      }
    }
  }
  return(V)
}

############################
# Section 2: Matrix things #
############################
Obj.F12 <- function(Pars) {
  j <<- j + 1
  int[[j]] <<- Pars
  # this is line 3 of page 10 of:
  # http://people.csail.mit.edu/xiuming/docs/tutorials/reml.pdf
  CovM  <- Create.CovM(Pars)
  A     <- -0.5*log(det(CovM))
  # Note: NaNs can be proceded in the next line.
  B     <- suppressWarnings(-0.5*log(det(t(X) %*% solve(CovM) %*% X)))
  est.b <- solve(t(X) %*% solve(CovM) %*% X) %*% t(X) %*% solve(CovM) %*% y
  tmp   <- y - X %*% est.b
  C     <- -0.5 *(t(tmp) %*% solve(CovM) %*% tmp)
  C     <- as.numeric(C) # make a scalar for other solvers
  return(A + B + C)
}
Final <- function(Pars) {
  CovM <- Create.CovM(Pars)
  tmp  <- solve(t(X) %*% solve(CovM) %*% X) %*% t(X) %*% solve(CovM) %*% y
  PE   <- exp(as.numeric(tmp["TrtT", ] - tmp["TrtR", ]))
  return(PE)
}

########################
# Section 3: Execution #
########################
Some.Initial.Guesses <- function(check = FALSE) {
  # design (all possible ones, error handling, full or partial)
  subjs <- unique(as.numeric(D$Subj))
  trts  <- sort(unique(as.character(D$Trt)), decreasing = TRUE)
  if (!length(trts) == 2) stop("Only two treatments supported.")
  if (sum(!trts %in% c("T", "R")) !=0)
    stop("treatments must be coded as 'T' and 'R'.")
  seqs  <- sort(unique(as.character(D$Seq)), decreasing = TRUE) # T first
  type  <- paste(seqs, collapse = "|") # identifier
  nseq  <- length(unique(as.character(D$Seq)))
  nper  <- length(unique(as.numeric(D$Per)))
  nsubj <- length(subjs)
  if (nchar(type) == 19) {  # 4-period 4-sequence full replicate designs
    if (nper != 4) stop("4 periods required in this full replicate design.")
    if (nseq != 4) stop("4 sequences required in this full replicate design.")
    des <- "full"
  }
  if (nchar(type) == 9) {  # 4-period full replicate designs
    if (nper != 4) stop("4 periods required in this full replicate design.")
    if (nseq != 2) stop("2 sequences required in this full replicate design.")
    des <- "full"
  }
  if (nchar(type) == 7) {  # 3-period replicates
    if (type %in% c("TRT|RTR", "TRR|RTT")) {
      if (nper != 3) stop("3 periods required in this full replicate design.")
      if (nseq != 2) stop("2 sequences required in this full replicate design.")
      des <- "full"
    }
    if (type == "TRR|RTR") {
      if (nper != 3) stop("3 periods required in the extra-reference design.")
      if (nseq != 2) stop("2 sequences required in the extra-reference design.")
      des <- "partial"
    }
  }
  if (nchar(type) == 11) { # Balaam's design or 3-sequence partial
    if (!type == "TR|RT|TT|RR") { # Balamm's
      if (nper != 3) stop("3 periods required in this partial replicate design.")
      if (nseq != 3) stop("3 sequences required in this partial replicate design.")
      des <- "partial"
    } else {                      # TRT|RTR|RRT
      if (nper != 2) stop("2 periods required in Balaam's design.")
      if (nseq != 4) stop("4 sequences required in Balaam's design.")
      des <- "full"
    }
  }
  if (type == "TR|RT") stop("TR|RT design not supported yet.")
  gm.T <- gm.R <- numeric(nsubj) # vector of geometric means
  for (subj in seq_along(subjs)) {
    gm.T[subj] <- mean(log(D$Y[which(D$Subj == subj & D$Trt == "T")]))
    gm.R[subj] <- mean(log(D$Y[which(D$Subj == subj & D$Trt == "R")]))
  }
  # remove NAs
  gm.T <- gm.T[!is.na(gm.T)]  
  gm.R <- gm.R[!is.na(gm.R)]  
  if (des == "partial") { # guess varT (between + within)
    m    <- lm(log(Y) ~ Seq + Per, data = D[D$Trt == "T", ])
    varT <- anova(m)["Residuals", "Mean Sq"]
  } else {                # guess VarWT
    m     <- lm(log(Y) ~ Seq + Per + Subj, data = D[D$Trt == "T", ])
    varWT <- anova(m)["Residuals", "Mean Sq"]
  }
  # guess VarWR
  m     <- lm(log(Y) ~ Seq + Per + Subj, data = D[D$Trt == "R", ])
  varWR <- anova(m)["Residuals", "Mean Sq"]   # guess varBR
  m     <- lm(log(Y) ~ Seq + Per, data = D[D$Trt == "R", ])
  varBR <- anova(m)["Residuals", "Mean Sq"] - varWR
  # TODO: full replicates
  varRT <- 0.5 * (varT + varBR)
  rslt  <- c(varT, varBR, varWR, varRT)
  return(rslt)
}
MyREML <- function(maxit = 1000, method = "Nelder-Mead") {
  Ini <- Some.Initial.Guesses()
  ptm <- proc.time()[[3]]
  F   <- optim(par = Ini, fn = Obj.F12, method = method,
               control = list(reltol = 1E-12, maxit = maxit,
               fnscale = -1))
  et  <- proc.time()[[3]] - ptm
  Nms <- c("var_T", "var_bR", "var_wR", "covar_bTbR")
  X   <- data.frame(Nms, Ini, F$par)
  names(X) <- c("Component", "Initial", "Estimate")
  if (F$convergence == 0)  msg <- "Successful completion."
  if (F$convergence == 1)  msg <- paste("Iteration limit", maxit, "reached.")
  if (F$convergence == 10) msg <- "Degeneracy of the Nelder–Mead simplex."
  cat(paste0("\nNelder-Mead optimization (", signif(et, 3), " seconds)"),
             "\n  LogLikelihood", F$value,
             "\n  Evaluations  ", F$counts[["function"]], "\n ", msg, "\n")
  return(X)
}

# Data has to have five columns in this order:
#   Subject (numeric), Period (numeric), Sequence (character),
#   Treatment (character), Response (numeric)
# Column separator 'sep' comma (default), semicolon, or t (tab)
# Decimal separator 'dec' period (default) or comma
path <- "https://bebac.at/downloads/ds02.csv" # EMA partial replicate
D    <- res <- PE <- int <- NULL # nullify eventual previous data and results
D    <- GetData(path)
X    <- CreateX(D)
y    <- log(D$Y)
# capture intermediate results:
# https://stackoverflow.com/questions/23975101/getting-more-details-from-optim-function-from-r
j    <- 0
int  <- list()
# the workhorse
#  optional arguments:
#  maxit = X to reduce the number of iterations from the default 1000
#  plot = TRUE for a plot of the progress of optimization
#  print = TRUE to print the progress of optimization
res <- MyREML()
res[, 1]     <- c("sigma\u00B2(T)", "sigma\u00B2(bR)",
                  "sigma\u00B2(wR)", "covar(bTbR)")
fin          <- data.frame(Parameter = c("sigma(T)", "CV(T)",
                                         "sigma(bR)", "CV(bR)",
                                         "sigma(wR)", "CV(wR)", "PE"),
                           Estimate = NA_real_)
for (j in 1:3) {
  fin$Estimate[j * 2 - 1] <- sprintf("%.7f", sqrt(res[j, 3]))
  fin$Estimate[j * 2]     <- sprintf("%3.3f%%", 100 * mse2CV(res[j, 3]))
}
fin$Estimate[7] <- sprintf("%3.4f%%", 100 * Final(res$Estimate))
int          <- data.frame(matrix(unlist(int), nrow = length(int),
                                  byrow = TRUE))
if (is.null(attr(dev.list(), "names"))) windows(width = 6.5, height = 4)
# variance components
lbls <- c(expression(sigma[T]^2),
          expression(sigma[bR]^2),
          expression(sigma[wR]^2),
          expression(italic(cov)(bTbR)))
cols <- c("darkgreen", "blue", "magenta", "red")
op   <- par(no.readonly = TRUE)
par(mar = c(3.3, 3.1, 1.6, 4.1),  mgp = c(2.2, 0.5, 0), lend = 2,
    ljoin = 1, tcl = 0.25, cex.axis = 0.8)
plot(1:nrow(int), int$X1, type = "n", ylim = range(int), log = "x",
     xlab = "Iteration", ylab = "Estimate", main = "Progress of optimization",
     cex.lab = 0.95, cex.main = 1, font.main = 1, axes = FALSE)
grid()
axis(1)
axis(2, las = 1)
box()
ax <- axTicks(1, log = TRUE)
for (j in 1:4) {
  lines(1:nrow(int), int[, j], type="s", lwd = 2, col = cols[j])
  shadowtext(sqrt(2), int[1, j], cex = 0.65, pos = 3,
             labels = signif(int[1, j], 4), col = cols[j],
             bg = "white", r = 0.2)
  shadowtext(nrow(int)*0.9, tail(int[, j], 1), cex = 0.65, pos = 3,
             labels = signif(tail(int[, j], 1), 4), col = cols[j],
             bg = "white", r = 0.2)
  mtext(lbls[j], 4, at = tail(int[, j], 1), cex = 0.8, line = 0.5, las = 1)
}
par(op)
print(res, row.names = FALSE); cat("\n"); print(fin, row.names = FALSE)

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
23,336 posts in 4,902 threads, 1,669 registered users;
22 visitors (1 registered, 21 guests [including 10 identified bots]).
Forum time: 16:01 CET (Europe/Vienna)

Biostatistician. One who has neither the intellect for mathematics
nor the commitment for medicine but likes to dabble in both.    Stephen Senn

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