I was trying to figure out the power for a Dunnett many to one comparisons design, not just the omnibus F-test There is not a noncentral distribution understanding for the Dunnett t distribution that I know of, which is a multivariate normal type. So I did something that matched my simulation below, using Stata's power oneway command, using the contrast option, with an alpha determined from a critical Dunnett t statistic using the ttail function. The example gets the adjusted two tailed probability for an ANOVA design with k=7 treatment groups and 1 control group with n=4 subjects per group. I'm not aware of others approaching power for Dunnett's comparisons in this way, but it seems to have worked. I was going to translate this to a full power command, but wanted to know others opinions on how I've approached this.
Here is the full simulation and matching power oneway calculations:
Code:
local k = 7 local n = 4 local v = (`k'+1)*(`n'-1) local dunnett_t = invdunnettprob(`k', `v', 0.95) scalar adj_p = 2*(1-t(`v', `dunnett_t')) display adj_p
Code:
local k = 7 local n = 4 local v = (`k'+1)*(`n'-1) local dunnett_t = invdunnettprob(`k', `v', 0.95) scalar adj_p = 2*(1-t(`v', `dunnett_t')) display adj_p local kd "40 35 70 20 30 50 55" local t1 = 5 // control group dct local i = 2 local dct = "5" foreach l of local kd { local t`i' = `t1' + ln((100-`l')/100)/ln(2)*-1 display `t`i'' local dct = "`dct'" + " `t`i''" local ++i } numlist "`dct'" local mymeans `r(numlist)' display "`mymeans'" power oneway `mymeans', /// varerror(0.25) npergroup(4) // power .97 power oneway `mymeans', /// varerror(0.25) contrast(-1 1 0 0 0 0 0 0) /// npergroup(4) alpha(`=adj_p') // power .26 power oneway `mymeans', /// varerror(0.25) contrast(-1 0 1 0 0 0 0 0) /// npergroup(4) alpha(`=adj_p') // power .17 power oneway `mymeans', /// varerror(0.25) contrast(-1 0 0 1 0 0 0 0) /// npergroup(4) alpha(`=adj_p') // power .98 power oneway `mymeans', /// varerror(0.25) contrast(-1 0 0 0 1 0 0 0) /// npergroup(4) alpha(`=adj_p') // power .04 power oneway `mymeans', /// varerror(0.25) contrast(-1 0 0 0 0 1 0 0) /// npergroup(4) alpha(`=adj_p') // power .11 power oneway `mymeans', /// varerror(0.25) contrast(-1 0 0 0 0 0 1 0) /// npergroup(4) alpha(`=adj_p') // power .52 power oneway `mymeans', /// varerror(0.25) contrast(-1 0 0 0 0 0 0 1) /// npergroup(4) alpha(`=adj_p') // power .67 // compare power of F-test and all Dunnett p-value frame reset frame create results double(ftest dunnett2vs1 dunnett3vs1 dunnett4vs1 /// dunnett5vs1 dunnett6vs1 dunnett7vs1 dunnett8vs1) quietly { forvalues i = 1(1)5000 { frame create simdata frame change simdata set obs 4 local kd "40 35 70 20 30 50 55" local t1 = 5 // control group dct local j = 2 foreach l of local kd { local t`j' = `t1' + ln((100-`l')/100)/ln(2)*-1 local ++j } local ve = 0.5 local ntx :word count `kd' local groups = `ntx' + 1 generate y1 = rnormal(`t1', `ve') forvalues j = 2/`groups' { generate y`j' = rnormal(`t`j'', `ve') } generate row = _n reshape long y, i(row) j(group) regress y i.group scalar ftest = Ftail(e(df_m), e(df_r), e(F)) pwcompare group, effects mcompare(dunnett) forvalues j = 1/`ntx' { local g = `j' + 1 scalar dunnett`g'vs1 = el(r(table_vs),4,`j') // *1/2 if one-sided } frame post results /// (scalar(ftest)) /// (scalar(dunnett2vs1)) /// (scalar(dunnett3vs1)) /// (scalar(dunnett4vs1)) /// (scalar(dunnett5vs1)) /// (scalar(dunnett6vs1)) /// (scalar(dunnett7vs1)) /// (scalar(dunnett8vs1)) frame change default frame drop simdata if mod(`i',100) == 0 noisily display `i' } } frame change results save results_power_2.dta, replace generate sig_ftest = ftest < 0.05 // power = .97, above = .97 generate sig_2vs1 = dunnett2vs1 < 0.05 // power = .26, above = .26 generate sig_3vs1 = dunnett3vs1 < 0.05 // power = .18, above = .17 generate sig_4vs1 = dunnett4vs1 < 0.05 // power = .97, above = .98 generate sig_5vs1 = dunnett5vs1 < 0.05 // power = .04, above = .04 generate sig_6vs1 = dunnett6vs1 < 0.05 // power = .12, above = .11 generate sig_7vs1 = dunnett7vs1 < 0.05 // power = .52, above = .52 generate sig_8vs1 = dunnett8vs1 < 0.05 // power = .67, above = .67 summarize sig* frame change default frame drop results
Comment