*! version 1 25Jan2016
program define _reg_distdichoMSK, rclass byable(recall)

	syntax varname [, /*
	*/ tail(string)  cp(real 0) dist(string)  modvar(real 0) dim(real 0)]
	*qui net install st0207, from(http://www.stata-journal.com/software/sj10-4)

	*no. observations
	local dicho = "`varlist'"

	quietly tabulate `dicho' if e(sample), matcell(m)

	tempvar nobs

	matrix `nobs' = m
	forvalues x=1/`dim'{

		local obs`x' = `nobs'[`x',1]
	}

	*margins of dicho. variable
	quietly margins `1'

	*marginal means
	tempvar mmean
	matrix `mmean' = r(b)

	forvalues x=1/`dim'{
		local mean`x' = `mmean'[1,`x']
	}

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

	*value labels
	quietly levelsof `dicho', local(levels)  sep(,)
	mat l = ( `r(levels)')
	qui net install st0207, from(http://www.stata-journal.com/software/sj10-4)
		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)
display ""
display "Alpha:  " `alpha'

	*distributional method

	forvalues x=2/`dim'{
	di as text ""
	local le =l[1,`x']
	di as text "{it}Comparisons of category value `le' and the reference category{sf}"
	di as text ""
		*distdichoi `obs`x'' `mean`x'' `modvar' `obs1' `mean1' `modvar' `cp' `tail'
		*****************************************
	local R=1
	local n1=`obs1'
	local m1=`mean1'
 local s1=`modvar'
 local n2=`obs`x''
 local m2=`mean`x''
 local s2=`modvar'
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"
di "above the cut-point " `cp'
}

else{
local prop1= `p1'
local prop2= `p2'
}
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 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= (`prop2' - `prop1')

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=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')




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'
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 `1'

                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 "`le' "  _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) c(level) _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
	}
end