*! version 6 20Jan2016

program distdicho, rclass byable(recall)

syntax  varlist(min=2 max=2 numeric)[if][in][, /*
*/ TWovar tail(string)  varr(real 1) UNEqual CORRection BOOTci nrep(real 2000)] cp(real)

tokenize `varlist'




if "`correction'" != "" {
	local R = 0
}
else {
	if "`unequal'" != "" {
	local R = 0
		}
	else{
		local R = `varr'
		}
		}


if "`twovar'"!=""{
marksample touse
quietly summarize `1' if `touse'
local m1 = r(mean)
local s1 = r(sd)
local n1 = r(N)
quietly summarize `2' if `touse'
local m2 = r(mean)
local s2 = r(sd)
local n2 = r(N)
}
else{

marksample touse

preserve
*display `cp'
quietly drop if `2'==.

quietly levelsof `2', local(levels)  sep(,)
mat l = ( `r(levels)')



quietly bysort `2': summarize `1' if `touse'

*local m1 = r(mean)
*local s1 = r(sd)
*local n1 = r(N)
local l1 = l[1,2]

*tempvar inv_gp
*gen `inv_gp'=1-`2'
*quietly bysort `inv_gp': summarize `1'  if `touse'
*local m2 = r(mean)
*local s2 = r(sd)
*local n2 = r(N)
local l2= l[1,1]



}



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

if "`twovar'"== ""{
	if "`unequal'" != "" {
	ttest `1' if `touse', by(`2') unequal
		local m1 = r(mu_2)
		local s1 = r(sd_2)
		local n1 = r(N_2)
		local m2 = r(mu_1)
		local s2 = r(sd_1)
		local n2 = r(N_1)

	}
	else {
		if `varr'!=1 {
		ttest `1' if `touse', by(`2') unequal
		local m1 = r(mu_2)
		local s1 = r(sd_2)
		local n1 = r(N_2)
		local m2 = r(mu_1)
		local s2 = r(sd_1)
		local n2 = r(N_1)

		}
			else {
			ttest `1' if `touse', by(`2')
			local m1 = r(mu_2)
			local s1 = r(sd_2)
			local n1 = r(N_2)
			local m2 = r(mu_1)
			local s2 = r(sd_1)
			local n2 = r(N_1)

			}
			}
			}
else{
	if "`unequal'" != "" {
	ttest `1'==`2' if `touse', unpaired unequal
			local m1 = r(mu_2)
			local s1 = r(sd_2)
			local n1 = r(N_2)
			local m2 = r(mu_1)
			local s2 = r(sd_1)
			local n2 = r(N_1)

	}
	else {
		if `varr'!=1 {
		ttest `1' ==`2' if `touse', unpaired unequal
			local m1 = r(mu_2)
			local s1 = r(sd_2)
			local n1 = r(N_2)
			local m2 = r(mu_1)
			local s2 = r(sd_1)
			local n2 = r(N_1)

		}
			else {
			ttest `1' ==`2' if `touse', unpaired
			local m1 = r(mu_2)
			local s1 = r(sd_2)
			local n1 = r(N_2)
			local m2 = r(mu_1)
			local s2 = r(sd_1)
			local n2 = r(N_1)

			}
			}
			}

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"
di "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"
di "below the cut-point " `cp'

}


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

local ttestse=r(se)
local ttestp=r(p)
if "`R'"!="0"{
display "Standard error computed under the hypothesis that"
di "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{
	if "`correction'"!=""{
		display "Standard error computed with correction for large effects"
		}
		else {
		display "Standard error computed with correction for"
		di "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 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')

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

	*assumption of unknown variance ratio :lower tail
	if "`R'" == "0" {
		if "`tail'" == "upper"{
				quietly bootstrap /*
			*/	pdiff=(-(normal((`cp'-r(mu_2))/r(sd_2))-normal((`cp'-r(mu_1))/r(sd_1)))) /*
			*/  rr=((1-normal((`cp'-r(mu_2))/r(sd_2)))/(1-normal((`cp'-r(mu_1))/r(sd_1)))) /*
			*/  or=(normal((`cp'-r(mu_1))/r(sd_1))*(1-normal((`cp'-r(mu_2))/r(sd_2)))/((1-normal((`cp'-r(mu_1))/r(sd_1)))*normal((`cp'-r(mu_2))/r(sd_2))))	/*
			*/, rep(`nrep') nodots strata(`2') bca: ttest `1', by( `2' )

		}
		else{

			quietly bootstrap /*
			*/	pdiff=(normal((`cp'-r(mu_2))/r(sd_2))-normal((`cp'-r(mu_1))/r(sd_1))) /*
			*/  rr=(normal((`cp'-r(mu_2))/r(sd_2))/normal((`cp'-r(mu_1))/r(sd_1))) /*
			*/  or=(normal((`cp'-r(mu_2))/r(sd_2))*(1-normal((`cp'-r(mu_1))/r(sd_1)))/((1-normal((`cp'-r(mu_2))/r(sd_2)))*normal((`cp'-r(mu_1))/r(sd_1))))	/*
			*/, rep(`nrep') nodots strata(`2') bca: ttest `1', by( `2' )

		}
	}
	else{
		if "`tail'" == "upper"{

			quietly bootstrap pdiff=(-(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*(`sd2')^2+(r(N_1)-1)*1*(`sd1')^2)/(r(N_1)+r(N_2)-2)))- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*1^(-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)*1^(-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)*1*(`sd1')^2)/(r(N_1)+r(N_2)-2)))*(1- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*1^(-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)*1*(`sd1')^2)/(r(N_1)+r(N_2)-2))))*	/*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*1*(`sd2')^2+(r(N_1)-1)*(`sd1')^2)/(r(N_1)+r(N_2)-2))))))	/*
						*/, rep(`nrep') nodots strata(`2') bca: ttest `1', by( `2' )

		}
		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)*(`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(`2') bca: ttest `1', by( `2' )

		}

	}
	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 intervals 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 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'


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

          quietly    summarize `2' 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 `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 `n2' _col(23) %9.0g `m2' _col(35) %9.0g `s2'  _col(48) %9.0g `prop2'
di as text "`lab2' "  _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 "{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}

}
else{

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 res abbrev(`"`1'"',8) _col(11) as txt " {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 res abbrev(`"`2'"',8) _col(11) as txt " {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 "{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}

}
end