*Version 3.0 08 Oct 2014


program distdicho, rclass byable(recall)

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

tokenize `varlist'

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


if "`method'"=="var"{

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'
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 "`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')

local ttestse=r(se)
local ttestp=r(p)
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')

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

	*assumption of unknown variance ratio :lower tail
	if "`R'" == "0" {
		if "`tail'" == "lower"{
				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'): 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=((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'): ttest `1', by( `2' )

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

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

		}
		else{
			quietly bootstrap pdiff=(-(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*r(sd_2)^2+(r(N_1)-1)*`R'*r(sd_1)^2)/(r(N_1)+r(N_2)-2)))- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*r(sd_2)^2+(r(N_1)-1)*r(sd_1)^2)/(r(N_1)+r(N_2)-2))))) /*
						*/  rr=((1-normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*r(sd_2)^2+(r(N_1)-1)*`R'*r(sd_1)^2)/(r(N_1)+r(N_2)-2))))/ /*
						*/		(1-normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*r(sd_2)^2+(r(N_1)-1)*r(sd_1)^2)/(r(N_1)+r(N_2)-2)))))/*
						*/  or=(1/(normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*r(sd_2)^2+(r(N_1)-1)*`R'*r(sd_1)^2)/(r(N_1)+r(N_2)-2)))*(1- /*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*r(sd_2)^2+(r(N_1)-1)*r(sd_1)^2)/(r(N_1)+r(N_2)-2))))/ /*
						*/ 		((1-normal((`cp'-r(mu_2))/sqrt(((r(N_2)-1)*r(sd_2)^2+(r(N_1)-1)*`R'*r(sd_1)^2)/(r(N_1)+r(N_2)-2))))*	/*
						*/		normal((`cp'-r(mu_1))/sqrt(((r(N_2)-1)*`R'^(-1)*r(sd_2)^2+(r(N_1)-1)*r(sd_1)^2)/(r(N_1)+r(N_2)-2))))))	/*
						*/, rep(`nrep') nodots strata(`2'): 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 interval calculated"
}


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 "`l1'          {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 "`l2'          {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