Bioequivalence and Bioavailability Forum

IUT and Dunnett: code for comparison of power [General Statistics]

posted by martin – Austria, 2011-11-06 00:43 (4552 d 09:43 ago) – Posting: # 7615
Views: 7,052
(edited by martin on 2011-11-06 11:04)

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)

Complete thread:

Two-stage (any method) and multiplicity Helmut 2011-10-30 17:36 [General Statistics]
- 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

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

Complete thread:

IUT and Dunnett: code for comparison of power [General Statistics]