library(powerbydesign) # ----------------------------------------------------------------------- # -------------------------- CRD example -------------------------------- # ----------------------------------------------------------------------- # the first dialog screen is for within-subject designs, it is not # needed here, just close it to get to the dialog screen we need # to match our class example, mu1=2.82, mu2=3.89, mu3=3.04, sd=.274 designcrd <- design.anova(between = list(A = c("a1","a2","a3"))) # n_from and n_to are PER GROUP, not total power_by_samplesize_crd <- boot.power.anova(designcrd, alpha=.01, n_from = 2, n_to = 10, num_iterations_bootstrap = 1000) # use at least 1000 'bootstrap' samples for final results, # but you can try a smaller number first to estimate running time # ~20 sec with 1000 # the plot is n per group, but the listing is total N = n*g # list results power_by_samplesize_crd > power_by_samplesize_crd $A n_total power 1 6 0.162 2 9 0.644 3 12 0.933 4 15 0.983 5 18 1.000 6 21 1.000 7 24 1.000 8 27 1.000 9 30 1.000 plot(power_by_samplesize_crd$A$n_total,power_by_samplesize_crd$A$power,type="b", main="CRD Power Example", xlab="Total sample size",ylab="Power to detect Ha") abline(h=.9, col="blue") # ----------------------------------------------------------------------- # ----------------------- factorial example ----------------------------- # ----------------------------------------------------------------------- # the first dialog screen is for within-subject designs, it is not # needed here, just close it to get to the dialog screen we need # mu11=10, mu12=15, mu21=20, mu22=27, std=10 designcrf22 <- design.anova(between = list(A = c("a1","a2"), B = c("b1","b2"))) # n_from and n_to are PER GROUP, not total power_by_samplesize_crf22 <- boot.power.anova(designcrf22, alpha=.05, n_from = 5, n_to = 25, num_iterations_bootstrap = 1000) # use at least 1000 'bootstrap' samples for final results, # but you can try a smaller number first to estimate running time # ~50 sec with 1000 # list results power_by_samplesize_crf22 > power_by_samplesize_crf22 $A n_total power 1 20 0.637 2 24 0.703 3 28 0.797 4 32 0.866 5 36 0.892 6 40 0.934 ... 21 100 1.000 $B n_total power 1 20 0.239 2 24 0.277 3 28 0.304 4 32 0.373 ... 17 84 0.778 18 88 0.780 19 92 0.805 20 96 0.838 21 100 0.836 $`A:B` n_total power 1 20 0.044 2 24 0.056 ... 19 92 0.077 20 96 0.080 21 100 0.090 # plot results x11() par(mfrow=c(2,2)) plot(power_by_samplesize_crf22$A$n_total,power_by_samplesize_crf22$A$power,type="b", main="CRF22 A Main Effect",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.8, col="blue") plot(power_by_samplesize_crf22$B$n_total,power_by_samplesize_crf22$B$power,type="b", main="CRF22 B Main Effect",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.8, col="blue") plot(power_by_samplesize_crf22$'A:B'$n_total,power_by_samplesize_crf22$'A:B'$power, type="b",main="CRF22 AB Interaction",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.8, col="blue") par(mfrow=c(1,1)) # -------------------------------------------------------------------------- # ----------------------- split plot factorial ----------------------------- # -------------------------------------------------------------------------- # using the settings for means, std, and correlation structure as # shown in the SAS code: # Control 10 18 15 # Drug1 12 20 17 # Drug2 15 10 15 # std = 2.2 # matrix ("SPFCorr") = ( 1 # .5 1 # .5 .5 1 ) des_spf <- design.anova(between = list(TrtA = c("Control","Drug1","Drug2")), within = list(TrtB = c("b1","b2","b3")), default_within_correlation=.5 ) # n_from and n_to are PER WHOLE PLOT GROUP, not total power_sampsize_spf33 <- boot.power.anova(des_spf, alpha=.01, n_from = 3, n_to = 15, num_iterations_bootstrap = 1000) # use at least 1000 'bootstrap' samples for final results, # but you can try a smaller number first to estimate running time # ~ 50 sec with 1000 # list results power_sampsize_spf33 $TrtA n_total power 1 9 0.083 2 12 0.168 3 15 0.274 4 18 0.339 5 21 0.474 6 24 0.543 7 27 0.655 8 30 0.735 9 33 0.782 10 36 0.836 11 39 0.898 12 42 0.908 13 45 0.935 $TrtB n_total power 1 9 0.945 2 12 0.998 3 15 1.000 4 18 1.000 ... 13 45 1.000 $`TrtB:TrtA` n_total power 1 9 1 2 12 1 ... 13 45 1 # plot results # using a lowess line can help see the power function more clearly x11() par(mfrow=c(2,2)) plot(power_sampsize_spf33$TrtA$n_total,power_sampsize_spf33$TrtA$power,type="b", main="SPF33 A Main Effect",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.8, col="blue") lines(lowess(power_sampsize_spf33$TrtA$n_total,power_sampsize_spf33$TrtA$power)) plot(power_sampsize_spf33$TrtB$n_total,power_sampsize_spf33$TrtB$power,type="b", main="SPF33 B Main Effect",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.8, col="blue") lines(lowess(power_sampsize_spf33$TrtB$n_total,power_sampsize_spf33$TrtB$power)) plot(power_sampsize_spf33$'TrtB:TrtA'$n_total,power_sampsize_spf33$'TrtB:TrtA'$power, type="b",main="SPF33 AB Interaction",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.8, col="blue") lines(lowess(power_sampsize_spf33$'TrtB:TrtA'$n_total,power_sampsize_spf33$'TrtB:TrtA'$power)) par(mfrow=c(1,1)) # -------------------------------------------------------------------------- # -------------------- randomized complete block --------------------------- # -------------------------------------------------------------------------- # as mentioned in the SAS notes, this is an approximate method where we give # factorial design input but focus our attention just on the treatment A results # the first dialog screen is for within-subject designs, it is not # needed here, just close it to get to the dialog screen we need # mu11=10, mu12=14, mu13=16, mu21=12, mu22=16, mu18 = std=2 desrcb2blk <- design.anova(between = list(A = c("a1","a2","a3"), Blk = c("bl1","bl2"))) # n_from and n_to are PER GROUP, not total pow_by_sampsz_rcb2blk <- boot.power.anova(desrcb2blk, alpha=.05, n_from = 2, n_to = 15, num_iterations_bootstrap = 1000) # use at least 1000 'bootstrap' samples for final results, # but you can try a smaller number first to estimate running time # ~50 sec with 1000 # ----------------------------------------------------------------------- # ----------------------- repeated measures ----------------------------- # ----------------------------------------------------------------------- # using the settings for means, std, and correlation structure as # shown in the SAS code: # when filling in the correlation structure, remember to fill it in # for each level of the whole plot factor # SensoryFocus 2.40 2.38 2.05 1.90 # StandardOfCare 2.40 2.39 2.36 2.30 # std = .92 # matrix ("FourTimesCorr") = ( 1 # .6 1 # .491 .495 1 # .399 .402 .491 1 ) des_rptms <- design.anova(between = list(Trt = c("SensFocus","StdCare")), within = list(Time = c("Time0","Time1Wk","Time6Mo","Time12Mo")) ) # ------------------------------------------------------------------- # try 10 boots first just to time it ~ 14 min 30 sec with 10 going 50 to 200, power_rptms0 <- boot.power.anova(des_rptms,n_from = 50, n_to = 200, num_iterations_bootstrap = 10) # list results power_rptms0 # plot results x11() par(mfrow=c(2,2)) plot(power_rptms0$Time$n_total,power_rptms0$Time$power,type="b", main="Rpt Ms Example Time",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.9, col="blue") lines(lowess(power_rptms0$Time$n_total,power_rptms0$Time$power)) plot(power_rptms0$Trt$n_total,power_rptms0$Trt$power,type="b", main="Rpt Ms Example Trt",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.9, col="blue") lines(lowess(power_rptms0$Trt$n_total,power_rptms0$Trt$power)) plot(power_rptms0$'Time:Trt'$n_total,power_rptms0$'Time:Trt'$power,type="b", main="Rpt Ms Example Time x Trt",ylim=c(0.,1.), xlab="Total sample size",ylab="Power to detect Ha") abline(h=.9, col="blue") lines(lowess(power_rptms0$'Time:Trt'$n_total,power_rptms0$'Time:Trt'$power)) par(mfrow=c(1,1)) # initial results seem a bit to liberal on power estimates, but the # number of simulations is far too low