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