*! 1.0.0 Univariate categorical GOF tests with Bonferroni post-hoc
* Author: Brent D. Hill (North Dakota State University)
* This program was adapted from the following packages:
* - tab_chi/chitest (N. J. Cox, 1999)
* - csgof (Statistical Consulting Group - UCLA, 2015)

program def ucgof, rclass
	version 10.0
	syntax varlist [, model(numlist min=2) freq ]

	preserve

	if ("`freq'" == "freq") {
		display _newline(0)
		display as text "--- Data in aggregate (tabular) format ---"
		sort `varlist'
		quietly tabulate `varlist' [fw=_freq]
		local J = `r(r)'
		local n = `r(N)'
	}

	if ("`freq'" != "freq") {
		quietly tabulate `varlist'
		local J = `r(r)'
		local n = `r(N)'
	}

	if ("`model'" == "") {
		local percexp = round(100/`J', .001)
		forvalues i = 1 / `J' {
			local model `percexp' `model'
		}
	}

	local expcat : word count `model'
	if (`expcat' != `J') {
		display in red "*** The variable `varlist' has `J' categories, but the specified model has `expcat' values, `model'."
		display in red "*** Please supply `J' values in the model() option or omit that option altogether."
		exit
	}

	if ("`freq'" != "freq") {
		quietly table `varlist', replace
		quietly rename table1 Obs_Freq
		sort `varlist'
	}
	if ("`freq'" == "freq") {
		//quietly table `varlist' [fw=Obs_Freq], replace
		quietly rename _freq Obs_Freq
		sort `varlist'
	}

	local i = 0
	local totperc = 0
	quietly generate Exp_Freq = .
	quietly generate model = .
	quietly generate Zres = .

	foreach percent of local model {
		local i = `i' + 1
		local totperc = `totperc' + `percent'
		quietly replace model = `percent' in `i'
		quietly replace Exp_Freq = round((`percent'/100)*`n', .001) in `i'
	}

	if abs(`totperc' - 100) > 1 {
		display in red "*** The specified model (null) percentages do not total to 100%."
		display in red "*** The sum of `model' is `totperc'."
		display in red "*** Please supply percentages that sum to 100%."
		exit
	}

	tempvar pearson1 pearson2 lrchisq1 lrchisq2
	quietly gen `pearson1' = (Obs_Freq - Exp_Freq)^2 / Exp_Freq
	quietly gen `pearson2' = sum(`pearson1')
	quietly gen `lrchisq1' = Obs_Freq*ln(Obs_Freq/Exp_Freq)
	quietly gen `lrchisq2' = 2*sum(`lrchisq1')

	local df = `J' - 1
	local osl_pearson = chi2tail(`df', `pearson2'[_N])
	local osl_lrchisq = chi2tail(`df', `lrchisq2'[_N])
	quietly replace Zres = (Obs_Freq-Exp_Freq)/sqrt(`n'*.01*model*(1-.01*model))
	quietly gen p = `J'*2*(1-normal(abs(Zres)))
	quietly replace p = 1 if p>1
	quietly gen Model = string(model) + "%"
	quietly drop model
	quietly format %9.2f Exp_Freq
	quietly format %9.0f Obs_Freq
	quietly format %9.2f Zres
	quietly format %4.3f p

	display _newline(0)
	display as text "  Chi-squared GOF tests (df = " as result "`df'" as text "):" _newline(1)
	display as text "  Likelihood-ratio = " as result %4.3f round(`lrchisq2'[_N], .001) as text ", p = " as result %5.4f round(`osl_lrchisq', 0.0001)
	display as text "           Pearson = " as result %4.3f round(`pearson2'[_N], .001) as text ", p = " as result %5.4f round(`osl_pearson', 0.0001)
	list `varlist' Model Exp_Freq Obs_Freq Zres p, noobs sep(0)
	display as text "  Note: P-values have been adjusted using the Bonferroni method."

end