Bioequivalence and Bioavailability Forum

Dear HS !

thank you for your input and the references!

below some code to play around; it is an optimal contrast problem – the power to detect a deviation from dose-proportionality depend on the unknown monotone dose-response shape. Theoretically, the hypothesis regarding dose-proportionality should be specified in detail to “maximise power”. For example, the hypothesis

H0: all 4 doses are within the proportional range
H1: 3 highest doses are outside the proportional range

would suggest application of the reverse helmert contrast. Additionally, one-sided step-down testing allows a formal estimation of the minimum dose where we can not reject the null hypothesis of dose-proportionality – similar as estimation of the minimum effective dose.

best regards

martin

rm(list=ls(all=TRUE))
# function to simulate power for an overall assessment of dose-proportionality based on the power-law model and 4 dose levels
sim.dose4prop <- function(nsim, cv, n, a, b, mindose, maxdose, range, alpha=0.05){

# input parameters:
# nsim ... number of simulation runs
# cv ... coefficient of variation of random variables (e.g. AUCs)
# n ... n per dose level
# a, b ... parameters for power law model: AUC = a*dose^b
# mindose ... lowest dose
# maxdose ... highest dose
# range ... doses (index) for which parameter b of power law model is applied (e.g. b=1.1 and range=c(2:4) specifies that doses 2 to 4 are outside proportional range)
# alpha ... type I error, two-sided (default = 0.05)

   # function to define AUC using power law model
   powermod <- function(dose, a, b, range){return(c(a*dose[-range]^1, a*dose[range]^b))}

   # function to generate log-normal distributed RVs (e.g. AUCs) with identical CVs per dose. See e.g. Jaki et al. (2009) for more detail.
   # Jaki T, Wolfsegger MJ, Ploner M (2009). Confidence intervals for ratios of AUCs in the case of serial sampling: A comparison of seven methods.
   # Pharmaceutical Statistics, 8(1):12-24. DOI:10.1002/pst.321.
   data.gen <- function(dose, mue, n, cv){
      k <- length(dose)   
      sd.log <- sqrt(log(1+cv^2))
      res <- as.data.frame(matrix(nrow=k*n, ncol=6))
      res[,1] <- sort(rep(dose,n))
      res[,2] <- sort(rep(mue,n))
      res[,3] <- log(res[,2]) - 1/2 * log(1+cv^2)
      res[,4] <- sd.log
      res[,5] <- rlnorm(n=k*n, meanlog=res[,3], sdlog=res[,4])
      res[,6] <- res[,5]/res[,1]       
      colnames(res) <- c('dose', 'mue', 'mue.log', 'sd.log', 'auc', 'auc.adj')
      return(res)
   }

   # definition of contrasts see e.g. Bretz and Hothorn (2003) for more detail
   # Bretz F and Hothorn LA (2003). Statistical Analysis of Monotone or Non-monotone Dose–Response Data from In Vitro Toxicological Assays.
   # ATLA-ALTERN LAB ANIM 31:81-96 Suppl. 1 JUN 2003. URL: http://ecvam.jrc.it/publication/Bretz-Hothorn.pdf last visited at 2009-12-27
   numdose <- 4
   contr.lowhigh <- c(-1,0,0,1) # t-test of highest versus lowest dose (pairwise contrast)
   contr.helmert <- c(-1,-1,-1,3) # helmert contrast   
   contr.revhelm <- c(-3,1,1,1) # reverse helmert contrast   
   contr.linearc <- c(-3,-1,1,3) # linear contrast
   contr.jonckht <- c(-1, -0.67, -0.17, 1.83) # Jonkheere-Terpstra analogon contrast
   contr.maximum <- c(-0.87, -0.13, 0.13, 0.87) # Maximum contrast
   contr.linr24t <- c(-12, -2, 2, 12) # Linear 2–4 contrast

   # define doses: to be equidistant on log-scale
   step <- (log(maxdose)-log(mindose))/(numdose-1)
   log.dose <- c(log(mindose), rep(NA, numdose-1))
   for(i in 2:numdose){log.dose[i] <- log.dose[i-1]+step}
   dose <- exp(log.dose)
   mue <- powermod(dose=dose, a=a, b=b, range=range)
   k <- length(dose)
   power <- rep(0,9)

   for(i in 1:nsim){

      # generate log-normal distributed random variables
      data <- data.gen(dose=dose, mue=mue, n=n, cv=cv)

      # log-log regression approach: two-sided CI
      res1 <- confint(lm(log(data$auc) ~ log(data$dose)), level=1-alpha)[2,]
      
      # ANOVA approach
      temp <- summary(lm(log(data$auc.adj) ~ as.factor(data$dose)))
      res2 <- 1-pf(q=temp$fstatistic[1], df1=temp$fstatistic[2], df2=temp$fstatistic[3])
   
      # contrast approaches: two-sided hypothesis
      xq <- tapply(log(data$auc.adj), data$dose, mean)
      sd.pooled <- sqrt(sum((n-1)*tapply(log(data$auc.adj), data$dose, var))/(n*k-k))
      df <- n*k - k
   
      nom3 <- sum(contr.lowhigh*xq)
      den3 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.lowhigh)^2))

      nom4 <- sum(contr.helmert*xq)
      den4 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.helmert)^2))

      nom5 <- sum(contr.revhelm*xq)
      den5 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.revhelm)^2))

      nom6 <- sum(contr.linearc*xq)
      den6 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.linearc)^2))

      nom7 <- sum(contr.jonckht*xq)
      den7 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.jonckht)^2))

      nom8 <- sum(contr.maximum*xq)
      den8 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.maximum)^2))

      nom9 <- sum(contr.linr24t*xq)
      den9 <- (sd.pooled/sqrt(n))*sqrt(sum((contr.linr24t)^2))

      res3 <- 2*(1-pt(q=abs(nom3/den3), df=df))
      res4 <- 2*(1-pt(q=abs(nom4/den4), df=df))
      res5 <- 2*(1-pt(q=abs(nom5/den5), df=df))
      res6 <- 2*(1-pt(q=abs(nom6/den6), df=df))
      res7 <- 2*(1-pt(q=abs(nom7/den7), df=df))
      res8 <- 2*(1-pt(q=abs(nom8/den8), df=df))
      res9 <- 2*(1-pt(q=abs(nom9/den9), df=df))

      if(res1[1] > 1 | res1[2] < 1){power[1] <- power[1] + 1}
      if(res2 < alpha){power[2] <- power[2] + 1}
      if(res3 < alpha){power[3] <- power[3] + 1}
      if(res4 < alpha){power[4] <- power[4] + 1}
      if(res5 < alpha){power[5] <- power[5] + 1}
      if(res6 < alpha){power[6] <- power[6] + 1}
      if(res7 < alpha){power[7] <- power[7] + 1}
      if(res8 < alpha){power[8] <- power[8] + 1}
      if(res9 < alpha){power[9] <- power[9] + 1}
   }

   res <- rep(NA, 11)
   res[1] <- mindose
   res[2] <- maxdose
   res[c(3:11)] <- power/nsim
   names(res) <- c('mindose', 'maxdose', 'log-log reg', 'anova', 't-test', 'helmert', 'revhelm', 'linearc', 'jonckht', 'maximum', 'linr24t')
   return(res)
}


# set parameters for simulations
nsim <- 1E4
cv <- 0.5
n <- 10
a <- 1

# dose-proportionality for complete dose range (simulate type I error)
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1, mindose=10, maxdose=20, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1, mindose=10, maxdose=100, range=c(1:4))

# all 4 doses outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(1:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(1:4))

# last 3 doses outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(2:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(2:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(2:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(2:4))

# last 2 doses outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(3:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(3:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(3:4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(3:4))

# only highest dose outside proportional range
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=20, range=c(4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.1, mindose=10, maxdose=100, range=c(4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=100, range=c(4))
sim.dose4prop(nsim=nsim, cv=cv, n=n, a=a, b=1.15, mindose=10, maxdose=1000, range=c(4))