## IUT and Dunnett: code for comparison of power [General Sta­tis­tics]

Dear HS !

Please find enclosed a modified (rather slow) code for evaluation of power and type I error for Dunnett and IUT for 2-to-1 comparisons (i.e. B vs. A and C vs. A) providing:
• Probability for showing equivalence in both comparisons with IUT
• Probability for showing equivalence in both comparisons with Dunnett
• Probability for showing equivalence in at least one comparison with Dunnett
My summary:
• When you are interested in showing equivalence in both comparisons use IUT
• When you are interested in showing equivalence in at least on comparison use Dunnett with the drawback of a loss in power regarding showing equivalence in both comparisons.
best regards

Martin

library(multcomp) library(gmodels) set.seed(24354) simfun <- function(mue, n, cv, nsim){     # 2-to-1 comparisons: B vs. A and C vs. A     grp <- factor(c(rep('A', n[1]), rep('B', n[2]), rep('C', n[3])))     count1 <- 0; count2 <- 0; count3 <- 0     # define logical variables necessary for evaluating the probability     # for showing equivalence in at least one comparison     comp1 <- TRUE; comp2 <- TRUE; scenario <- 'NA'     if(mue[2]/mue[1] < 0.8 | mue[2]/mue[1] > 1.25){comp1 <- FALSE} # comparison 1 is true inequivalent     if(mue[3]/mue[1] < 0.8 | mue[3]/mue[1] > 1.25){comp2 <- FALSE} # comparison 2 is true inequivalent     if((comp1==TRUE & comp2==TRUE) | (comp1==FALSE & comp2==FALSE)){scenario <- '0'}     if(comp1==FALSE & comp2==TRUE){scenario <- '1'}     if(comp1==TRUE & comp2==FALSE){scenario <- '2'}     print(scenario)       for(i in 1:nsim){        x <- c(rlnorm(n=n[1], meanlog=log(mue[1])-0.5*log(cv^2+1), sdlog=sqrt(log(cv^2+1))),               rlnorm(n=n[2], meanlog=log(mue[2])-0.5*log(cv^2+1), sdlog=sqrt(log(cv^2+1))),               rlnorm(n=n[3], meanlog=log(mue[3])-0.5*log(cv^2+1), sdlog=sqrt(log(cv^2+1))))         # contrasts of interest: ANOVA model to get unadjusted many-to-one comparisons         mod.raw <- lm(log(x) ~ 0 + grp)         res1.raw <- exp(estimable(mod.raw, c(-1,1,0), conf.int=0.9)[c(6,7)])         res2.raw <- exp(estimable(mod.raw, c(-1,0,1), conf.int=0.9)[c(6,7)])         # contrasts of interest: ANOVA model to get adjusted many-to-one comparisons using Dunnett         mod.dun <- summary(glht(aov(log(x) ~ grp), linfct=mcp(grp='Dunnett')))         res1.dun <- exp(confint(mod.dun, level=0.9)$confint[1,c(2,3)]) res2.dun <- exp(confint(mod.dun, level=0.9)$confint[2,c(2,3)])         # IUT: probability for equivalence in both comparisons         if(res1.raw[1] >= 0.8 & res1.raw[2] <= 1.25 &            res2.raw[1] >= 0.8 & res2.raw[2] <= 1.25){count1 <- count1 + 1}         # Dunnett: probability for equivalence in both comparisons         if(res1.dun[1] >= 0.8 & res1.dun[2] <= 1.25 &            res2.dun[1] >= 0.8 & res2.dun[2] <= 1.25){count2 <- count2 + 1}                 # Dunnett: probability for equivalence in at least one comparison         switch(scenario,            '0' = {if((res1.dun[1] >= 0.8 & res1.dun[2] <= 1.25) |                      (res2.dun[1] >= 0.8 & res2.dun[2] <= 1.25)){count3 <- count3 + 1}         }, '1' = {if(res1.dun[1] >= 0.8 & res1.dun[2] <= 1.25){count3 <- count3 + 1}         }, '2' = {if(res2.dun[1] >= 0.8 & res2.dun[2] <= 1.25){count3 <- count3 + 1}         })      }      res <- c(count1/nsim, count2/nsim, count3/nsim)      names(res) <- c('IUT: overall', 'Dunnett: overall', 'Dunnett: at least one')      return(res)  }  # define simulation parameters  n <- c(20,20,20)  cv <- 0.2  nsim <- 1E3  ## 2 equivalent scenarios: power  simfun(mue=c(1, 1, 1), n=n, cv=cv, nsim=nsim)  simfun(mue=c(1, 1.05, 1.05), n=n, cv=cv, nsim=nsim)  ## 3 inequivalent scenarios: type I error  # A and B & A and C are inequivalent  simfun(mue=c(1, 0.8-1E-16, 0.8-1E-16), n=n, cv=cv, nsim=nsim)  # A and B are equivalent & A and C are inequivalent  simfun(mue=c(1, 0.8, 0.8-1E-16), n=n, cv=cv, nsim=nsim)  # A and B are inequivalent & A and C are equivalent  simfun(mue=c(1, 0.8-1E-16, 0.8), n=n, cv=cv, nsim=nsim)