// looking for a difference of 5, between 10 and 15
// larger shape has more power (is more normal)
// larger n has more power

postfile nb group min max mean using nbinompower.dta, replace
postfile nbtest diff p using nbinomsig.dta, replace
quietly forvalues i = 1/5000 {
	clear
	*set obs 20  // power is ~ 78% for 5, ~41% for 3, ~5% for 0
	set obs 25  // power is ~ 86.5% for 5, ~50% for 3, ~5% for 0
	*set obs 30  // power is ~ 92% for 5, ~58% for 3, ~5% for 0
	generate y1 = rnbinomial(8, 8/18)
	generate y0 = rnbinomial(8, 8/18)
	generate y3 = rnbinomial(8, 8/21)
	generate y5 = rnbinomial(8, 8/23)
	summarize y1
	post nb (1) (r(min)) (r(max)) (r(mean))
	summarize y5
	post nb (2) (r(min)) (r(max)) (r(mean))
	ttest y1 == y0, unpaired
	post nbtest (0) (r(p))
	ttest y1 == y3, unpaired
	post nbtest (3) (r(p))
	ttest y1 == y5, unpaired
	post nbtest (5) (r(p))
}
postclose nb
postclose nbtest
use nbinompower.dta, clear
bysort group:  summarize
use nbinomsig.dta, clear
generate sig = (p < 0.05)
tab diff sig, row

postfile nb group min max mean using nbinompower.dta, replace
postfile nbtest p using nbinomsig.dta, replace
quietly forvalues i = 1/2500 {
	clear
	set obs 20
	generate y1 = rnbinomial(10, 0.5)
	generate y2 = rnbinomial(10, 0.4)
	summarize y1
	post nb (1) (r(min)) (r(max)) (r(mean))
	summarize y2
	post nb (2) (r(min)) (r(max)) (r(mean))
	ttest y1 == y2, unpaired
	post nbtest (r(p))
}
postclose nb
postclose nbtest
use nbinompower.dta, clear
bysort group:  summarize
use nbinomsig.dta, clear
generate sig = (p < 0.05)
tab sig