*! version 8.1 20Apr2016
program define reg_distdicho, rclass byable(recall)

	syntax varname (min=1 max=1 numeric) [, /*
	*/ tail(string)   dist(string)   nrep(real 2000)] cp(real)

	*tokenize `varname'

	if (e(cmd)=="regress"){
	*Std.Dev. of model fitted
	local std_mod = e(rmse)

	}
	else if (e(cmd)=="xtmixed"|| e(cmd)=="mixed"){
		*Std.Dev. of model fitted
		matrix b = e(b)
		matrix ln_std = b[1, "lnsig_e:"]
		local dim = e(redim)
		local rvar = 0
		forvalues x=1/`dim'{
		matrix ln_rv`x' = b[1, "lns1_1_`x':"]
		local rvar =`rvar'+ exp(ln_rv`x'[1,1])^2
		}
		local std_mod = (exp(ln_std[1,1])^2+`rvar')^0.5

	}
	else if (e(cmd)=="xtreg"){
		*Std.Dev. of model fitted
		local std_mod = (e(sigma_e)^2 + e(sigma_u)^2)^0.5
	}
	else if (e(cmd)==.){
		error 301
	}
	else{
		di as error "cannot perform reg_distdicho for used function"
	}

	*margins of dicho. variable
	quietly  margins `1'

	*dicho. or categ. variable?
	tempvar mmean
	matrix `mmean' = r(b)
	local count = colsof(`mmean')

	if (`count'>=3){

		_reg_distdichoM `0' modvar(`std_mod') dim(`count')
		exit
	}

	*marginal means
	local mean1 = `mmean'[1,1]
	local mean2 = `mmean'[1,2]

	*no. observations
	local dicho = "`varlist'"
	quietly tabulate `dicho' if e(sample), matcell(m)
	tempvar nobs
	matrix `nobs' = m
	local obs1 = `nobs'[1,1]
	local obs2 = `nobs'[2,1]

	*remarks
	di as text ""
	di as text "Comparisons of proportions based on marginal effects of "
	di as input e(cmdline)
	di as text ""

	*distributional method
	if "`dist'"=="sk" {


		*sk_distdichoi `obs2' `mean2' `std_mod' `obs1' `mean1' `std_mod' `cp' `tail' `alpha' 1
**********************************
	local n1=`obs1'

 local m1=`mean1'
 local s1=`std_mod'
 local n2=`obs2'
 local m2=`mean2'
 local s2=`std_mod'
tempvar res
	quietly _predict `res', re
	qui skewnreg `res'
*alpha
local alpha = e(alpha)
*muz
local muz = sqrt(2/_pi)*`alpha'/sqrt(1+(`alpha')^2)
*equal variances, Ratio=1
local sd1 = sqrt(((`n1'-1)*`s1'^2+(`n2'-1)*`s2'^2)/(`n1'+`n2'-2))
local sd2 = `sd1'
*omega
local w1 = `sd1'/sqrt(1-(`muz')^2)
local w2 = `w1'
*alpha'
local alpha21 =`m1'-`w1'*`muz'
local alpha22 =`m2'-`w2'*`muz'
*call psn
psn `cp' `alpha21' `w1' `alpha'
local p1 =r(psn)
psn `cp' `alpha22' `w2' `alpha'
local p2 =r(psn)

if "`tail'"=="upper"{
local prop1= 1-`p1'
local prop2= 1-`p2'
display ""
display "Distributional estimates for the comparison of proportions above the cut-point " `cp'
}

else{
local prop1= `p1'
local prop2= `p2'

display ""
display "Distributional estimates for the comparison of proportions below the cut-point " `cp'

}

display ""
display "Alpha:  " `alpha'


local propdiff= -(`prop1' - `prop2')

local sediff = sqrt(((`w1')^2/`n1')*(1-(`muz')^2)*((2*(-exp(-(`cp'-`alpha21')^2/(2*(`w1')^2))/   	/*
			*/ sqrt(2*_pi*(`w1')^2))*normal(`alpha'*(`cp'-`alpha21')/`w1'))^2+(2*(-exp(-(`cp' 	/*
			*/ -`alpha22')^2/(2*(`w1')^2))/sqrt(2*_pi*(`w1')^2))*normal(`alpha'*(`cp'-`alpha22')/ 	/*
			*/ `w1'))^2))
local selogrr = sqrt(((`w1')^2/`n1')*(1-(`muz')^2)*((2*(-exp(-(`cp'-`alpha21')^2/(2*(`w1')^2))/	/*
			*/	sqrt(2*_pi*(`w1')^2))*normal(`alpha'*(`cp'-`alpha21')/`w1'))^2/`prop1'^2+(2*	/*
			*/ (-exp(-(`cp'-`alpha22')^2/(2*(`w1')^2))/sqrt(2*_pi*(`w1')^2))*normal(`alpha'		/*
			*/ *(`cp'-`alpha22')/`w1')/`prop2')^2))
local selogor = sqrt(((`w1')^2/`n1')*(1-(`muz')^2)*((2*(-exp(-(`cp'-`alpha21')^2/(2*(`w1')^2))/		/*
			*/	sqrt(2*_pi*(`w1')^2))*normal(`alpha'*(`cp'-`alpha21')/`w1'))^2/(`prop1'*			/*
			*/ (1-`prop1'))^2+(2*(-exp(-(`cp'-`alpha22')^2/(2*(`w1')^2))/sqrt(2*_pi*(`w1')^2))		/*
			*/  *normal(`alpha'*(`cp'-`alpha22')/`w1')/(`prop2'*(1-`prop2')))^2))


local lev=100-(100-c(level))/2
local ciinf = -(`prop1' - `prop2')-invnorm(`lev'/100)*`sediff'
local cisup = -(`prop1' - `prop2')+invnorm(`lev'/100)*`sediff'
local distrr `prop2'/`prop1'


local ciinfrr = exp(log(`distrr')-invnorm(`lev'/100)*`selogrr')
local cisuprr = exp(log(`distrr')+invnorm(`lev'/100)*`selogrr')

local distor `prop2'*(1-`prop1')/((1-`prop2')*`prop1')

local ciinfor = exp(log(`distor')-invnorm(`lev'/100)*`selogor')
local cisupor = exp(log(`distor')+invnorm(`lev'/100)*`selogor')


*Bootstrap confidence interval
if "`bootci'" != "" {

tempvar fact
gen `fact'=`1'
			if "`tail'" == "upper"{
			quietly bootstrap /*
			*/	pdiff=(-(r(prop2)-r(prop1))) 	/*
			*/	rr=((1-r(prop2))/(1-r(prop1)))		/*
			*/	or=(1/(r(prop2)*(1-r(prop1))/((1-r(prop2))*r(prop1))))/*
			*/, rep(`nrep') nodots strata(`fact'): hlp_boot `res' `fact' , p1(`muz') p2(`alpha') p3(`sd1') cpoint(`cp')
		}
		else{
				quietly bootstrap /*
			*/	pdiff=(r(prop2)-r(prop1)) 	/*
			*/	rr=(r(prop2)/r(prop1))		/*
			*/	or=(r(prop2)*(1-r(prop1))/((1-r(prop2))*r(prop1)))/*
			*/, rep(`nrep') nodots strata(`fact'): hlp_boot `res' `fact', p1(`muz') p2(`alpha') p3(`sd1') cpoint(`cp')
			}


	mat b_ci = e(ci_bc)

	*modify ci
	local ciinf = b_ci[1,1]
	local cisup = b_ci[2,1]

	local ciinfrr = b_ci[1,2]
	local cisuprr = b_ci[2,2]

	local ciinfor = b_ci[1,3]
	local cisupor = b_ci[2,3]

	display "bootstrap confidence interval calculated with `nrep' bootstrap repetitions"
}

local alselogrr= sqrt((exp(`selogrr'^2)-1)*exp(2*log(`distrr')+`selogrr'^2))
local alselogor= sqrt((exp(`selogor'^2)-1)*exp(2*log(`distor')+`selogor'^2))

return scalar prop1=`prop1'
return scalar prop2=`prop2'
return scalar propdiff=`propdiff'
return scalar sediff=`sediff'

return scalar ciinf=`ciinf'
return scalar cisup=`cisup'
return scalar distrr=`distrr'
return scalar selogrr=`alselogrr'

return scalar ciinfrr=`ciinfrr'
return scalar cisuprr=`cisuprr'
return scalar distor=`distor'
return scalar selogor=`alselogor'

return scalar ciinfor=`ciinfor'
return scalar cisupor=`cisupor'

return scalar alpha=`alpha'


if "`twovar'"==""{
display ""

          quietly    summarize `1'
                *if r(N) == 0 { noisily error 2000 }

                local val1 = r(min)
                local val2 = r(max)

				local lab1 = `val1'
                local lab2 = `val2'


                *local bylab : value label `2'

                if "`bylab'"!="" {

                        local lab1 : label `bylab' `lab1'

						if `"`lab1'"'=="" {
                                local lab1 = `val1'
                        }
                        local lab2 : label `bylab' `lab2'
                        if `"`lab2'"'=="" {
                                local lab2 = `val2'
                        }
                }

                local lab1 = substr(`"`lab1'"',1,8)
                local lab2 = substr(`"`lab2'"',1,8)


di as text "{hline 11}{c TT}{hline 50}
di as text "Group      {c      |}   Obs" _col(25) "Mean      Std dev.    Dist. prop."
di as text "{hline 11}{c +}{hline 50}
di as text "`lab1'"  _col(12) as txt "{c      |}" in y _col(11) %7.0f `n1' _col(23) %9.0g `m1' _col(35) %9.0g `s1'  _col(48) %9.0g `prop1'
di as text "`lab2' "  _col(12) as txt "{c      |}" in y _col(11) %7.0f `n2' _col(23) %9.0g `m2' _col(35) %9.0g `s2'  _col(48) %9.0g `prop2'
di as text "{hline 11}{c +}{hline 50}

di as text "Stat       {c      |}   Estimate  " _col(17) "Std error     [" _col(28) `lev' _col(31)"% Conf. Interval]    "
di as text "{hline 11}{c +}{hline 50}
di as text "Diff. prop {c |}" in y _col(14) %9.0g `propdiff' _col(25) %9.0g `sediff' _col(38) %9.0g `ciinf'  _col(50) %9.0g `cisup'
di as text "Risk ratio {c |}" in y _col(14) %9.0g `distrr' _col(25) %9.0g `alselogrr' _col(38) %9.0g `ciinfrr'  _col(50) %9.0g `cisuprr'
di as text "Odds ratio {c |}" in y _col(14) %9.0g `distor' _col(25) %9.0g `alselogor' _col(38) %9.0g `ciinfor'  _col(50) %9.0g `cisupor'
di as text "{hline 11}{c BT}{hline 50}

}




*************************************
		}

		if "`dist'"=="" {


		*distdichoi `obs2' `mean2' `std_mod' `obs1' `mean1' `std_mod' `cp' `tail' `r' 1

				*****************************************

		local R=1
		local n1=`obs1'

 local m1=`mean1'
 local s1=`std_mod'
 local n2=`obs2'
 local m2=`mean2'
 local s2=`std_mod'
 local sd1=sqrt(((`n1'-1)*`s1'^2+(`n2'-1)*`R'*`s2'^2)/(`n1'+`n2'-2))
local sd2=sqrt(((`n1'-1)*`R'^(-1)*`s1'^2+(`n2'-1)*`s2'^2)/(`n1'+`n2'-2))
 if "`tail'"=="upper"{

	local prop1= 1-normal((`cp'-`m1')/`sd1')
	local prop2= 1-normal((`cp'-`m2')/`sd2')
	display "Distributional estimates for the comparison of proportions above the cut-point " `cp'
}

else {
local prop1= normal((`cp'-`m1')/`sd1')
local prop2=normal((`cp'-`m2')/`sd2')
display "Distributional estimates for the comparison of proportions below the cut-point " `cp'

}


local propdiff= -(`prop1' - `prop2')

local ttestse=r(se)
local ttestp=r(p)


local sediff = sqrt(exp(-(`cp'-`m1')^2/`sd1'^2)/(2*_pi*`n1')+exp(-(`cp'-`m2')^2/`sd2'^2)/(2*_pi*`n2'))
local selogrr = sqrt(exp(-(`cp'-`m1')^2/`sd1'^2)/(2*_pi*`prop1'^2*`n1')+exp(-(`cp'-`m2')^2/`sd2'^2)/(2*_pi*`prop2'^2*`n2'))
local selogor = sqrt(exp(-(`cp'-`m1')^2/`sd1'^2)/(2*_pi*`prop1'^2*(1-`prop1')^2*`n1')+exp(-(`cp'-`m2')^2/`sd2'^2)/(2*_pi*`prop2'^2*(1-`prop2')^2*`n2'))
local lev=c(level)
local ciinf = -(`prop1' - `prop2')-invnorm(`lev'/100)*`sediff'
local cisup = -(`prop1' - `prop2')+invnorm(`lev'/100)*`sediff'
local distrr `prop2'/`prop1'


local ciinfrr = exp(log(`distrr')-invnorm(`lev'/100)*`selogrr')
local cisuprr = exp(log(`distrr')+invnorm(`lev'/100)*`selogrr')

local distor `prop2'*(1-`prop1')/((1-`prop2')*`prop1')

local ciinfor = exp(log(`distor')-invnorm(`lev'/100)*`selogor')
local cisupor = exp(log(`distor')+invnorm(`lev'/100)*`selogor')

*Bootstrap confidence interval
if "`bootci'" != "" {

tempvar res
quietly _predict `res', re
tempvar fact
gen `fact'=`1'

		if "`tail'" == "upper"{
				quietly bootstrap pdiff=(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2)))- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd2')^2)/(r(N_1)+r(N_2)-2)))) /*
						*/  rr=(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2)))/ /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2))))/*
						*/  or=(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2)))*(1- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2))))/ /*
						*/ 		((1-normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2))))*	/*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2)))))	/*
						*/, rep(`nrep') nodots strata(`fact'):  `res', by( `fact' )

		}
		else{
		quietly bootstrap pdiff=(-(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2)))- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2))))) /*
						*/  rr=((1-normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2))))/ /*
						*/		(1-normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2)))))/*
						*/  or=(1/(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2)))*(1- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2))))/ /*
						*/ 		((1-normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*`R'*(`sd1')^2)/(r(N_1)+r(N_2)-2))))*	/*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2))))))	/*
						*/, rep(`nrep') nodots strata(`fact'): ttest `res', by( `fact' )
	*strata(`3')
		}





	mat b_ci = e(ci_bc)

	*modify ci
	local ciinf = b_ci[1,1]
	local cisup = b_ci[2,1]

	local ciinfrr = b_ci[1,2]
	local cisuprr = b_ci[2,2]

	local ciinfor = b_ci[1,3]
	local cisupor = b_ci[2,3]

	display "bootstrap confidence interval calculated with `nrep' bootstrap repetitions"
}


local alselogrr= sqrt((exp(`selogrr'^2)-1)*exp(2*log(`distrr')+`selogrr'^2))
local alselogor= sqrt((exp(`selogor'^2)-1)*exp(2*log(`distor')+`selogor'^2))



tempvar fact
qui gen `fact'=`1'
qui summarize `fact'


*if `touse'
                *if r(N) == 0 { noisily error 2000 }

                local val1 = r(min)
                local val2 = r(max)

				local lab1 = `val1'
                local lab2 = `val2'

              *local bylab : value label `fact'
                if "`bylab'"!="" {

                        local lab1 : label `bylab' `lab1'

						if `"`lab1'"'=="" {
                                local lab1 = `val1'
                        }
                        local lab2 : label `bylab' `lab2'
                        if `"`lab2'"'=="" {
                                local lab2 = `val2'
                        }
                }

                local lab1 = substr(`"`lab1'"',1,8)
                local lab2 = substr(`"`lab2'"',1,8)


di as text "{hline 11}{c TT}{hline 50}
di as text "Group      {c      |}   Obs" _col(25) "Mean      Std dev.    Dist. prop."
di as text "{hline 11}{c +}{hline 50}
di as text "`lab1'"  _col(12) as txt "{c      |}" in y _col(11) %7.0f `n1' _col(23) %9.0g `m1' _col(35) %9.0g `s1'  _col(48) %9.0g `prop1'
di as text "`lab2' "  _col(12) as txt "{c      |}" in y _col(11) %7.0f `n2' _col(23) %9.0g `m2' _col(35) %9.0g `s2'  _col(48) %9.0g `prop2'
di as text "{hline 11}{c +}{hline 50}

di as text "Stat       {c      |}   Estimate  " _col(17) "Std error     [" _col(28) `lev' _col(31)"% Conf. Interval]    "
di as text "{hline 11}{c +}{hline 50}
di as text "Diff. prop {c |}" in y _col(14) %9.0g `propdiff' _col(25) %9.0g `sediff' _col(38) %9.0g `ciinf'  _col(50) %9.0g `cisup'
di as text "Risk ratio {c |}" in y _col(14) %9.0g `distrr' _col(25) %9.0g `alselogrr' _col(38) %9.0g `ciinfrr'  _col(50) %9.0g `cisuprr'
di as text "Odds ratio {c |}" in y _col(14) %9.0g `distor' _col(25) %9.0g `alselogor' _col(38) %9.0g `ciinfor'  _col(50) %9.0g `cisupor'
di as text "{hline 11}{c BT}{hline 50}





		}
	ret add
return clear

return scalar prop1=`prop2'
return scalar prop2=`prop1'
return scalar propdiff=`propdiff'
return scalar sediff=`sediff'

return scalar ciinf=`ciinf'
return scalar cisup=`cisup'
return scalar distrr=`distrr'
return scalar serr=`alselogrr'

return scalar ciinfrr=`ciinfrr'
return scalar cisuprr=`cisuprr'
return scalar distor=`distor'
return scalar seor=`alselogor'

return scalar ciinfor=`ciinfor'
return scalar cisupor=`cisupor'

end

program reg_dist, eclass byable(onecall)
	reg_distdicho `0'
end
exit