Log in
RegisterIUT and Dunnett: code for comparison of power [General Statistics]
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:
Martin
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
- 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.
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)Complete thread:
- Two-stage (any method) and multiplicity Helmut 2011-10-30 17:36 [General Statistics]
![[Mix]](https://static.bebac.at/img/mix.png "open in Mix view")
- Two-stage (any method) and multiplicity ElMaestro 2011-10-30 18:26
- Dose-dependent PK Helmut 2011-11-02 22:30
- Dunnett not for continuous scales? Really? d_labes 2011-11-01 11:11
- Yes, yes – but another construction site Helmut 2011-11-02 22:18
- PowerTOST help d_labes 2011-11-03 13:18
- PowerTOST help (solved) Helmut 2011-11-03 14:28
- intersection-union test martin 2011-11-03 22:28
- IUT and Dunnett: code for comparison of powermartin 2011-11-05 23:43
- PowerTOST help d_labes 2011-11-03 13:18
- Yes, yes – but another construction site Helmut 2011-11-02 22:18
- Two-stage (any method) and multiplicity ElMaestro 2011-10-30 18:26
Ing. Helmut Schütz