program define distdichoi, rclass byable(recall)


args n1 m1 s1 n2 m2 s2 cp tail R


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 `R'
                /*local n1 `n1'
                local m1 `m1'
                local s1 `s1'
				local n2 `n2'
                local m2 `m2'
                local s2 `s2'
                local cp `cp'*/

                mac shift 8

                local options "Yname(string)  "
                *local options "`options' Level(cilevel)"
                parse "`*'"
				}
*default values
if "`R'"=="" {
	local R 1
}


ttesti `n1' `m1' `s1' `n2' `m2' `s2'


if "`R'"=="0"{
	display ""
	local sd1=`s1'
	local sd2=`s2'
}
else{
	display ""
	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 ""
display "Distributional estimates for the comparison of proportions below the cut-point " `cp'

}


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

if "`R'"!="0"{
display "Standard error computed under the hypothesis that the ratio of variances is equal to " `R'
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'))

}
else{
	display "Standard error computed with correction for unknown variance ratio"
	if "`tail'"=="upper"{
		if (`cp'-`m1')>0{
		local corr1=(0.05/sqrt(`n1'))*((`cp'-`m1')/`sd1')*((`cp'-`m1')/`sd1'+2.6)
		local corr2=(0.45/sqrt(`n1'))*abs((`cp'-`m1')/`sd1')^(5/2)
		local corr3=`corr2'
		}
		else {
		local corr1=(0.05/sqrt(`n1'))*((`cp'-`m1')/`sd1')*((`cp'-`m1')/`sd1'-2.6)
		local corr2=(0.055/sqrt(`n1'))*((`cp'-`m1')/`sd1')*((`cp'-`m1')/`sd1'-2.6)
		local corr3=(0.45/sqrt(`n1'))*abs((`cp'-`m1')/`sd1')^(5/2)
		}
	}
	else{
		if (`cp'-`m1')>0{
		local corr1=(0.05/sqrt(`n1'))*((`cp'-`m1')/`sd1')*((`cp'-`m1')/`sd1'-2.6)
		local corr2=(0.055/sqrt(`n1'))*((`cp'-`m1')/`sd1')*((`cp'-`m1')/`sd1'-2.6)
		local corr3=(0.45/sqrt(`n1'))*abs((`cp'-`m1')/`sd1')^(5/2)
		}
		else{
		local corr1=-(0.05/sqrt(`n1'))*((`cp'-`m1')/`sd1')*((`cp'-`m1')/`sd1'+2.6)
		local corr2=(0.45/sqrt(`n1'))*abs((`cp'-`m1')/`sd1')^(5/2)
		local corr3=`corr2'
		}
	}

local sediff = sqrt(exp(-(`cp'-`m1')^2/`sd1'^2)/(2*_pi*`n1')+exp(-(`cp'-`m2')^2/`sd2'^2)/(2*_pi*`n2'))+`corr1'
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'))+`corr2'
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'))+`corr2'
}

local ciinf = (`prop1' - `prop2')-1.96*`sediff'
local cisup = (`prop1' - `prop2')+1.96*`sediff'
local distrr `prop1'/`prop2'


local ciinfrr = exp(log(`distrr')-1.96*`selogrr')
local cisuprr = exp(log(`distrr')+1.96*`selogrr')

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

local ciinfor = exp(log(`distor')-1.96*`selogor')
local cisupor = exp(log(`distor')+1.96*`selogor')
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=`selogrr'

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

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


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 "{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 `selogrr' _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 `selogor' _col(38) %9.0g `ciinfor'  _col(50) %9.0g `cisupor'
di as text "{hline 11}{c BT}{hline 50}
end