*! version 2 25Jan2016
program define sk_distdichoi, rclass byable(recall)


args n1 m1 s1 n2 m2 s2 cp tail alpha


global S_9            /* will be z    */
        parse "`*'", parse(" ,")
        cap confirm number `8'
        if "`4'" != "" & _rc == 0 {
                confirm integer number `n1'
                confirm number `m1'
				confirm number `s1'
                confirm integer number `n2'
                confirm number `m2'
				confirm number `s2'
				confirm number `cp'
				confirm string `tail'
				confirm number `alpha'
				mac shift 8
                local options "Yname(string)  "
                *local options "`options' Level(cilevel)"
                parse "`*'"
				}

ttesti `n1' `m1' `s1' `n2' `m2' `s2'
*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"
di "above the cut-point " `cp'
}

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

display ""
display "Distributional estimates for the comparison of proportions"
di  "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 `prop1'/`prop2'


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

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

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'
return scalar alpha=`alpha'


display ""
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 "x          {c      |}" in y _col(12) %7.0f `n1' _col(23) %9.0g `m1' _col(35) %9.0g `s1'  _col(48) %9.0g `prop1'
di as text "y          {c      |}" in y _col(12) %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    [95% Conf. Interval]    "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}
end