*! version 3 20Jan2016

program sk_distdicho, rclass byable(recall)

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

tokenize `varlist'

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'"== ""{
	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{
			ttest `1' ==`2' if `touse'
			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)
			}




qui net install st0207, from(http://www.stata-journal.com/software/sj10-4)
qui skewnreg `1' if `touse'
*alpha
local alpha = e(alpha)
*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')


*Bootstrap confidence interval
if "`bootci'" != "" {
			if "`tail'" == "upper"{
			quietly bootstrap /*
			*/	pdiff=(-(r(prop2)-r(prop1))) 	/*
			*/	rr=((1-r(prop2))/(1-r(prop1)))		/*
			*/	or=(1/(r(prop2)*(1-r(prop1))/((1-r(prop2))*r(prop1))))/*
			*/, rep(`nrep') nodots strata(`2'):  hlp_boot `1' `2' , p1(`muz') p2(`alpha') p3(`sd1') cpoint(`cp')
		}
		else{
				quietly bootstrap /*
			*/	pdiff=(r(prop2)-r(prop1)) 	/*
			*/	rr=(r(prop2)/r(prop1))		/*
			*/	or=(r(prop2)*(1-r(prop1))/((1-r(prop2))*r(prop1)))/*
			*/, rep(`nrep') nodots strata(`2'): hlp_boot `1' `2', p1(`muz') p2(`alpha') p3(`sd1') cpoint(`cp')
			}


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


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 "`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 "`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 "{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 `n1' _col(23) %9.0g `m1' _col(35) %9.0g `s1'  _col(48) %9.0g `prop1'
di as res abbrev(`"`2'"',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 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