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

posted by martin  – Austria, 2011-11-06 00:43 (4578 d 10:50 ago) – Posting: # 7615
Views: 7,105

(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:My summary: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:

UA Flag
Activity
 Admin contact
23,029 posts in 4,834 threads, 1,644 registered users;
33 visitors (0 registered, 33 guests [including 6 identified bots]).
Forum time: 12:34 CEST (Europe/Vienna)

We must be careful not to confuse data with the abstractions
we use to analyze them.    William James

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