*! MLB 1.1.7 07 Apr 2011
*above shift to version number of betafit
* MLB 1.3.0 06 Apr 2011
* MLB 1.2.0 01 Nov 2009
* MLB 1.1.0 20 Jan 2007
* MLB 1.0.0 08 Aug 2006
program dbetafit, rclass
	if c(stata_version) >= 11 {
		version 11
	}
	else {
		version 8.2
	}
	syntax [, at(string) ]
	if "`e(cmd)'" != "betafit" {
		di as err "results for betafit not found"
		exit 301
	}
	if "`e(title)'" != "ML fit of beta (mu, phi)" {
		di as err "dbetafit only possible with alternative parameterization"
		exit 301
	}

	/* get variable names */
	local vars: colnames(e(b_mu))

	if "`e(munocons)'" == "" {
		local cons "_cons"
		local vars : list vars - cons
	}
	else {
		foreach var of local vars {
			capture assert `var' == 1 if e(sample)
			if !_rc {
				local hascons "`var'"
			}
		}
	}
	if c(stata_version) >= 11 {
		foreach var of local vars {
			_ms_parse_parts `var'
			if inlist("`r(type)'", "interaction", "product") {
				di as err "dbetafit cannot meaningfully handle interactions, use margins instead"
				exit 198
			}
		}
	}

	/* make matrix with values at which effects are to be evaluated
	   and summary stats */
	tempname partat atx
	if "`at'" != "" {
		tokenize `at'
		local length : word count `at'
		local length = `length'/2
		local end = 2*`length'

		forvalues i = 1(2)`end' {
			local atvars "`atvars' ``i''"
		}

		local i = 1
		matrix `partat' = J(`length',5,0)

		foreach var of local atvars {
			local valid  : list var in vars
			if `valid' == 0 {
				di as err "`var' is not an explanatory variable in the last model"
				exit 198
			}
			local j = `i' *2
			if real("``j''") != . {
				quietly sum `var' if e(sample)
				matrix `partat'[`i',1] = ``j'', r(mean), r(sd), r(min), r(max)
			}
			else {
				quietly sum `var' if e(sample), detail
				local stat = cond("``j''"=="median", "p50", "``j''")
				local stats "mean min max p1 p5 p10 p25 p50 p75 p90 p95 p99"
				local valid : list stat in stats
				if `valid' {
					matrix `partat'[`i',1] = r(`stat'), /*
					   */ r(mean), r(sd), r(min), r(max)
				}
				else {
					di as err "You requested that the variable `var' is evaluated at `stat'."
					di as err "Variables can only be evaluated at a number or at"
					di as err "mean, median, min, max, p1, p5, p10, p25, p50, p75, p90, p95, and p99."
					exit 198
				}
			}
			local atvarlist "`atvarlist' `var'"
			local `++i'
		}
		matrix rownames `partat' = `atvarlist'
		matrix colnames `partat' = x mean sd min max
	}



	local k : word count `vars'
	matrix `atx' = J(`k',5,0)
	local i = 1
	foreach var of local vars {
		local part : list var in atvarlist
		if `part' == 1 matrix `atx'[`i',1] = `partat'[rownumb(`partat',"`var'"),1..5]
		else {
			quietly sum `var' if e(sample)
			matrix `atx'[`i',1] = r(mean), r(mean), r(sd), r(min), r(max)
		}
		local `++i'
	}

	matrix rownames `atx' = `vars'
	matrix colnames `atx' = x mean sd min max

	tempname summary b v xb xb1 xb2
	matrix `b' = e(b)
	matrix `b' = `b'[1, "mu:"]
	matrix `v' = e(V)

	matrix `summary' = `atx'
	if "`e(munocons)'" == ""  {
		matrix `atx' = `atx'[1...,1] \ 1
	}

	matrix `xb' = `b'*`atx'
	scalar `xb' = el(`xb', 1, 1)

	local vars : list vars - hascons
	local k : word count `vars'

	tempname minmax unit sd marg maxmarg b
	matrix `minmax'  = J(`k',2,.)
	matrix `unit'    = J(`k',2,.)
	matrix `sd'      = J(`k',2,.)
	matrix `marg'    = J(`k',2,.)
	matrix `maxmarg' = J(`k',2,.)


	local i = 0
	foreach var of local vars {
		quietly {
			local `i++'
			local varn : subinstr local var "." ""
			tempvar unique`varn'
			capture bys `var': gen byte `unique`varn'' =  _n==1
			if _rc {
				tempvar newvar
				gen `newvar' = `var'
				bys `newvar': gen byte `unique`varn'' =  _n==1
			}
			count if `unique`varn''
			if r(N) > 2 {
				scalar `xb1' = `xb' - `atx'[`i',1]*_b[`var'] + _b[`var']*`summary'[`i',5]
				scalar `xb2' = `xb' - `atx'[`i',1]*_b[`var'] + _b[`var']*`summary'[`i',4]
				scalar `b' = invlogit(`xb1') - invlogit(`xb2')
				local dg1db "invlogit(`xb1')*invlogit(-1*`xb1') *`summary'[`i',5]"
				local dg2db "invlogit(`xb2')*invlogit(-1*`xb2') *`summary'[`i',4]"
				matrix `minmax'[`i',1] = `b', sqrt(`v'[`i',`i']*(`dg1db' - `dg2db')^2)

				scalar `xb1' = `xb' + _b[`var']*`summary'[`i',3]/2
				scalar `xb2' = `xb' - _b[`var']*`summary'[`i',3]/2
				scalar `b' = invlogit(`xb1') - invlogit(`xb2')
				local dg1db "invlogit(`xb1')*invlogit(-1*`xb1') *(`summary'[`i',1]+`summary'[`i',3]/2)"
				local dg2db "invlogit(`xb2')*invlogit(-1*`xb2') *(`summary'[`i',1]-`summary'[`i',3]/2)"
				matrix `sd'[`i',1] = `b', sqrt(`v'[`i',`i']*(`dg1db' - `dg2db')^2)

				scalar `xb1' = `xb' + _b[`var']*.5
				scalar `xb2' = `xb' - _b[`var']*.5
				scalar `b' = invlogit(`xb1') - invlogit(`xb2')
				local dg1db "invlogit(`xb1')*invlogit(-1*`xb1') *(`summary'[`i',1]+.5)"
				local dg2db "invlogit(`xb2')*invlogit(-1*`xb2') *(`summary'[`i',1]-.5)"
				matrix `unit'[`i',1] = `b', sqrt(`v'[`i',`i']*(`dg1db' - `dg2db')^2)

				local p "invlogit(`xb')"
				local pp "invlogit(-1*`xb')"
				scalar `b' = `p'*`pp'*_b[`var']
				local dgdb "`p'*`pp'* (_b[`var']*`summary'[`i',1] - 2*_b[`var']*`summary'[`i',1]*`p'+1)"
				matrix `marg'[`i',1] = `b', sqrt(`v'[`i',`i']*(`dgdb')^2)

				matrix `maxmarg'[`i',1] = _b[`var']/4, _se[`var']/4
			}
			else {
				scalar `xb1' = `xb' - `atx'[`i',1]*_b[`var'] + _b[`var']*`summary'[`i',5]
				scalar `xb2' = `xb' - `atx'[`i',1]*_b[`var'] + _b[`var']*`summary'[`i',4]
				scalar `b' = invlogit(`xb1') - invlogit(`xb2')
				local dg1db "invlogit(`xb1')*invlogit(-1*`xb1') *`summary'[`i',5]"
				local dg2db "invlogit(`xb2')*invlogit(-1*`xb2') *`summary'[`i',4]"
				matrix `minmax'[`i',1] = `b', sqrt(`v'[`i',`i']*(`dg1db' - `dg2db')^2)
			}
		}
	}


	di in text "{hline 14}{c TT}{hline 64}
	di in text "discrete" _col(15) "{c |}  Min --> Max               +-SD/2                  +-1/2"
	di in text "change"   _col(15) "{c |}  coef.     se           coef.     se            coef.     se"
	di in text "{hline 14}{c +}{hline 64}
	local i = 1
	foreach var of local vars{
		local bminmax  = el(`minmax',`i',1)
		local seminmax = el(`minmax',`i',2)
		if el(`sd',`i',1) != . {
			local bsd      = el(`sd',`i',1)
			local sesd     = el(`sd',`i',2)
			local b1       = el(`unit',`i',1)
			local se1      = el(`unit',`i',2)
		}
		else {
			local bsd      = ""
			local sesd     = ""
			local b1       = ""
			local se1      = ""
		}
		output_line `var' `bminmax' `seminmax' `bsd' `sesd' `b1' `se1'
		local `i++'
	}
	di in text "{hline 14}{c BT}{hline 64}

	di _newline(1)

	di in text "{hline 14}{c TT}{hline 41}
	di in text "Marginal"  _col(15) "{c |}    MFX at x                Max MFX"
	di in text "Effects"   _col(15) "{c |}  coef.     se           coef.     se"
	di in text "{hline 14}{c +}{hline 41}
	local i = 1
	foreach var of local vars{
		local varn: subinstr local var "." ""
		quietly count if `unique`varn''
		if r(N) > 2 {
			local bmarg     = el(`marg',`i',1)
			local semarg    = el(`marg',`i',2)
			local bmaxmarg  = el(`maxmarg',`i',1)
			local semaxmarg = el(`maxmarg',`i',2)
			output_line `var' `bmarg' `semarg' `bmaxmarg' `semaxmarg'
		}
		local `i++'
	}
	di in text "{hline 14}{c BT}{hline 41}

	di _newline(1)

	di in text "E(" e(depvar) "|x) = " as result %6.0g invlogit(`xb')

	di _newline(1)

	matrix list `summary', noheader noblank format(%6.0g)

	matrix rownames `minmax'  = `vars'
	matrix rownames `sd'      = `vars'
	matrix rownames `unit'    = `vars'
	matrix rownames `marg'    = `vars'
	matrix rownames `maxmarg' = `vars'
	matrix colnames `minmax'  = "coef" "se"
	matrix colnames `sd'      = "coef" "se"
	matrix colnames `unit'    = "coef" "se"
	matrix colnames `marg'    = "coef" "se"
	matrix colnames `maxmarg' = "coef" "se"


	return matrix summary = `summary'
	return matrix maxmarg = `maxmarg'
	return matrix marg    = `marg'
	return matrix unit    = `unit'
	return matrix sd      = `sd'
	return matrix minmax  = `minmax'

	return scalar pr      = invlogit(`xb')
end

program output_line
	args var b1 se1 b2 se2 b3 se3
	noi di in text %13s abbrev("`var'",12) " {c |}" /*
	   */ as result /*
	   */          %6.0g `b1' _col(25) %6.0g `se1' /*
	   */ _col(39) %6.0g `b2' _col(49) %6.0g `se2' /*
	   */ _col(62) %6.0g `b3' _col(72) %6.0g `se3'
end