// Negative binomial
* Here r=8, the number of "successes" (the stopping number),
*    x is the number of "failures", and p=8/18=r/(r+mu) is
*  the probability of a "success".
twoway (function y=nbinomialp(8,x,8/18), range(0 20)) ///
	(function y=nbinomialp(8,x,8/(18+2.5)), range(0 20)) ///
	(function y=nbinomialp(8,x,8/(18+5)), range(0 20)), ///
	title("Negative binomail, d = 2.5")

postfile kwpower test d T dfm dfr using kwpowermcnb.dta, replace
quietly forvalues d=0(0.5)5 {
	noisily display "difference " `d'
	forvalues i=1/2500 {
		clear
		set obs 60 // group size = 20
		generate group = mod(_n, 3)
		generate y = rnbinomial(8, 8/(18+group*`d'))
		anova y group
		post kwpower (1) (`d') (e(F)) (e(df_m)) (e(df_r))
		kwallis y, by(group)
		post kwpower (2) (`d') (r(chi2_adj)) (r(df)) (.)
	}
}
postclose kwpower

use kwpowermcnb.dta, clear

generate sig = Ftail(dfm, dfr, T)<0.05 if test==1
replace sig = chi2tail(dfm, T)<0.05 if test==2
bysort test: tabulate d sig, row

collapse (mean) sig, by(test d)
separate sig, by(test)
label variable sig1 "Anova"
label variable sig2 "Kruskal-Wallis"
save kwpowermcnbtable.dta, replace
line sig1 sig2 d, title("Negative Binomial, by MC") ytitle("Power") ///
	xtitle("group differences")

// Negative binomial, more skewed
* Here r=2, the number of "successes" (the stopping number),
*    x is the number of "failures", and p=2/7=r/(r+mu) is
*  the probability of a "success".
twoway (function y=nbinomialp(2,x,2/7), range(0 20)) ///
	(function y=nbinomialp(2,x,2/(7+2.5)), range(0 20)) ///
	(function y=nbinomialp(2,x,2/(7+5)), range(0 20)), ///
	title("Negative binomial, d = 2.5")

postfile kwpower test d T dfm dfr using kwpowermcnb2.dta, replace
quietly forvalues d=0(0.5)5 {
	noisily display "difference " `d'
	forvalues i=1/2500 {
		clear
		set obs 60 // group size = 20
		generate group = mod(_n, 3)
		generate y = rnbinomial(2, 2/(7+group*`d'))
		anova y group
		post kwpower (1) (`d') (e(F)) (e(df_m)) (e(df_r))
		kwallis y, by(group)
		post kwpower (2) (`d') (r(chi2_adj)) (r(df)) (.)
	}
}
postclose kwpower

use kwpowermcnb2.dta, clear

generate sig = Ftail(dfm, dfr, T)<0.05 if test==1
replace sig = chi2tail(dfm, T)<0.05 if test==2
bysort test: tabulate d sig, row

collapse (mean) sig, by(test d)
separate sig, by(test)
label variable sig1 "Anova"
label variable sig2 "Kruskal-Wallis"
save kwpowermcnbtable2.dta, replace
line sig1 sig2 d, title("Negative Binomial, by MC") ytitle("Power") ///
	xtitle("group differences")