// Log-normal
twoway (function y=lnnormalden(x,0,1), range(0 4)) ///
	(function y=lnnormalden(x,0.5,1), range(0 4)) ///
	(function y=lnnormalden(x,1,1), range(0 4)), title("LogNormal, d = 0.5") ///
	legend(order(1 "mu=0" 2 "mu=0.5" 3 "mu=1"))

postfile kwmc rep d F dfm dfr X dfx using kwmclognorm.dta, replace
postfile kwperm rep d p3 p4 using kwpermlognorm.dta, replace

quietly forvalues d=0(0.1)1 {
	noisily display "difference " `d'
	forvalues i=1/1000 {
		noisily display "diff " `d' " rep " `i'
		clear
		set obs 60 // group size = 20
		generate group = mod(_n, 3)
		generate y = exp(rnormal(group*`d', 1))
		anova y group
		kwallis y, by(group)
		post kwmc (`i') (`d') (e(F)) (e(df_m)) (e(df_r)) (r(chi2_adj)) (r(df))
		
		permute group T=e(F), right reps(100) nodots: anova y group
		matrix A = r(p)
		scalar p3 = A[1,1]
		permute group T= r(chi2_adj), right reps(100) nodots nowarn: kwallis y, by(group)
		matrix A = r(p)
		scalar p4 = A[1,1]
		post kwperm (`i') (`d') (p3) (p4)
	}
}
postclose kwmc
postclose kwperm

use kwpermlognorm.dta, clear
sort rep d
save kwpermlognorm.dta, replace
use kwmclognorm.dta, clear
sort rep d
merge 1:1 rep d using kwpermlognorm.dta

generate sig1 = Ftail(dfm, dfr, F)<0.05
generate sig2 = chi2tail(dfx, X)<0.05
generate sig3 = p3<0.05
generate sig4 = p4<0.05
tabulate d sig1, row
tabulate d sig2, row
tabulate sig1 sig2

collapse (mean) sig1-sig4, by(d)
label variable sig1 "Anova"
label variable sig2 "Kruskal-Wallis"
label variable sig3 "permuted F"
label variable sig4 "permuted KW"
save kwlognormtable.dta, replace
line sig1 sig2 sig3 sig4 d, title("Log-Normal") ytitle("Power") ///
	xtitle("group differences, d")

// Log-normal arithmetic parameterization? Not! what van Hecke used.
twoway (function y=lnnormalden(x,ln(0.1)-0.5*ln(1/.1^2+1),sqrt(ln(1/.1^2+1))), range(0 4)) ///
	(function y=lnnormalden(x,ln(0.5)-0.5*ln(1/.5^2+1),sqrt(ln(1/.5^2+1))), range(0 4)) ///
	(function y=lnnormalden(x,ln(1)-0.5*ln(1/1^2+1),sqrt(ln(1/1^2+1))), range(0 4)), title("LogNormal, d ~= 0.5")

postfile kwpower test d T dfm dfr using kwpowermcln2.dta, replace
quietly forvalues d=0.1(0.05)0.5 {
	noisily display "difference " `d'
	forvalues i=1/2500 {
		clear
		set obs 60 // group size = 20
		generate group = mod(_n, 3)+1
		generate sigma = sqrt(log(1+1/(group*`d')^2))
		generate mu = log(group*`d')-0.5*sigma^2
		generate y = exp(rnormal(mu, sigma))
		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 kwpowermcln2.dta, clear

generate sig = Ftail(dfm, dfr, T)<0.05 if test==1
replace sig = chi2tail(dfm, T)<0.05 if test==2
*tabulate test sig, row
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 kwpowermcln2table.dta, replace
line sig1 sig2 d, title("Log-Normal, by MC") ytitle("Power") ///
	xtitle("group differences")