*! version 6.0.0  13sep1999  (STB-52 sg120)
prog def rocfit, eclass
version 6
	if replay() {
		if `"`e(cmd)'"' != "rocfit" { error 301 } 
		syntax [, Level(integer $S_level) *]
	}
	else {
		syntax varlist(numeric min=2 max=2) [if] [in]  /*
	 	*/ [ fweight] [, Level(integer $S_level)   /*
	 	*/ FROM(string) MLMethod(string) /*
		*/ MLOpt(string) noLOG  * ]
 
		marksample touse
		tempvar sample
		gen int `sample'=`touse'
		tokenize `varlist'
		local D = `"`1'"'
		local C = `"`2'"'
		if "`from'"!="" { local iniopt init(`from') }
		if "`mlmethod'" == "" { local mlmetho  "lf" }
		if "`offset'" !="" { local offopt offset(`offset') }
		mlopts options, `options'

		cap assert `D'==0 | `D'==1 if `touse'
		if _rc~=0 {
		 	noi di in red "true status variable `D' must be 0 or 1"
			exit 198
 		}
		qui summ `D' if `touse', meanonly
		if r(min) == r(max) {
			di in red "Outcome does not vary"
			exit 198
		}
		if "`log'"!="" { 
			local nlog="*" 
		}

		/* BEGIN-create freq weighted summary data */
		 if `"`weight'"' =="" {
			tempvar wv
			qui gen int `wv' = 1 if `touse'
			local weight="fweight"
		}
		else {
			tempvar wv
			qui gen double `wv' `exp'
		}
		tempvar newwt
		sort `touse' `D' `C'
		qui by `touse' `D' `C': gen `newwt'=sum(`wv') if `touse'
		qui by `touse' `D' `C': replace `newwt'=. if _n~=_N
		qui replace `newwt'=. if `newwt'==0
		qui replace `touse'=0 if `newwt'==.
		qui replace `wv'=`newwt'
		local wt=`"[fweight=`wv']"'
		tempvar MN
		sort `touse' `D' `C'
		qui by `touse' `D' : gen long `MN' = sum(`wv')
		qui by `touse' `D' : replace `MN' = `MN'[_N]
		qui replace `MN'=. if `touse'==0
		drop `newwt'
		*preserve
		*qui keep if `touse'
		/* END -create freq weighted summary data */

		/* BEGIN-get initial values */
		tempvar p  
		qui logistic `D' `C' `wt' if `touse' 
		tempvar b
		mat `b'=e(b)
		local bc=colsof(`b')
		if `bc'>1 {
			local coef=_b[`C']
		}
		else {
			logit  d m1, nocoef nolog
			di in blue "fitted ROC is 45 degree line." /*
			*/  " Use roctab command"			
			est clear
			exit	
		}

		
		qui predict double `p', p
		tempvar prob sens spec 
		qui sum `p', meanonly
		cap assert reldif(`p', r(min))<1e-12 if `p'~=.
		if _rc==0 { 
			_rocsen if `touse', class(`C') genprob(`prob') /*
			*/ gensens(`sens') genspec(`spec')
		}
		else {
			qui  lsens if `touse',  nograph genprob(`prob') /*
			*/ gensens(`sens') genspec(`spec')
		}
		qui replace `sens'=. if `prob'==0 | `prob'==1
		qui replace `spec'=. if `prob'==0 | `prob'==1
		qui replace `prob'=. if `sens'==. | `spec'==.
		sort `prob'
		qui replace `prob'=. if `sens'==1 & `spec'==0 & _n==1
		tempvar g
		qui egen `g'=group(`C') if `touse'
		qui sum `g' , meanonly 
		local max=r(max)
		local i= 1
		while `i'<`max' {
			local zv "`zv' /z`i'"			
			local nzv "`nzv' z`i'"
			local zname "`zname' Z:_cut`i'"
			local i=`i'+1
		}
		*noi di in red "`zv'"
		tempvar  zfpf ztpf ztnf
		qui gen double `ztpf'=invnorm(`sens') if `touse'
		qui gen double `zfpf'=invnorm(1-`spec') if `touse'
		qui gen double `ztnf'=invnorm(`spec') if `touse'
			
		if `max'>2 {
			qui regress `ztpf' `zfpf' 
			local a0 = _b[_cons]
			local b0 = _b[`zfpf']
		}
		else {
			local a0 = 1
			local b0 = 1 
		}
		drop `ztpf' `zfpf' `sens' `spec'
		tempname bc0
		qui mat `bc0' = ( `a0', `b0')   
		qui mat colnames `bc0' = a:_cons b:_cons

		tempname ZV

		sort `prob'
		local i= 1
		while `i'<`max' {
			tempname Z
			if `ztnf'[`i'] != . {
				cap mat `Z'=(`ztnf'[`i']+`i'*1e-6)
			}
			cap mat colnames `Z' = z`i':_cons
			if `i'==1 {
				 cap mat `ZV'= `Z'
			}
			else {
			 	 cap mat `ZV'=`ZV' , `Z'
			}
			if _rc~=0 {
				local i=`max'+1
				local this 1
			}
			else { local i=`i'+1 }
		}
		drop `ztnf'
		if "`this'"~="" {
			sort `prob'
			local i= 1
			while `i'<`max' {
				tempname Z
				qui mat `Z'=(invnorm(`prob'[`i']))
				qui mat colnames `Z' = z`i':_cons
				qui mat `bc0'=`bc0' , `Z'
				local i=`i'+1
			}
		}
		else {
			qui mat `bc0'=`bc0' , `ZV'
		}
		if `coef'<0 {
			local tel=colsof(`bc0')
			tempname telm 
			qui mat `telm'=(-`bc0'[1,1], `bc0'[1,2])
			qui mat colnames `telm' = a:_cons b:_cons
			local j 1
			local i=`tel'
			while `i' >2 {
				tempname telc
				mat `telc' =  -`bc0'[1,`i']
				qui mat colnames `telc' = z`j':_cons
				local i=`i'-1
				local j=`j'+1
				mat `telm'= `telm', `telc'
			}
			mat `bc0'= `telm'
		}	
	 	/* noi mat l `bc0' */
		`nlog' di _n in gr "Fitting binormal model:"
		global MD `D'
		global MC `g'
		global ZV `nzv'
		tempvar  p
		qui gen double `p'=.
		global MP `p'
		ml model `mlmethod' rocf_lf /*
			*/  (a: `D' `C'  =  )  /*
			*/  /b `zv' `wt' if `touse', `cont' /*
			*/ `robust' `cluopt' `scopt' init(`bc0') `mlopt' /*
			*/ missing collin nopreserve  /*
			*/ max search(off) `log' `options' 
		global MD
		global MC
		global ZV 
		tempname B
		qui mat `B'=e(b)
		local a = `B'[1,1]
		local b = `B'[1,2]
		mat colnames `B' = :intercep :slope `zname'
		tempname x2 df
		tempvar num
		qui gen double `num'= /*
		*/ (`MN'/$MP)*(($MP - `wv'/`MN')^2) if `touse'
		qui replace `num'=sum(`num')
		scalar `x2'=`num'[_N]	
		scalar `df'=`max'-3	
		est scalar chi2_gf=`x2' 
		est scalar df_gf=`df' 
		est scalar p_gf=chiprob(`df',`x2') 
		est local chi2type="GOF"
		est repost b=`B' , rename esample(`sample')
		if "`e(wtype)'" != "" {
			est local wexp `"`exp'"'
		}
		ROCInd `level'
		est local title "Binormal model" 
		est local predict rocfit_p
		est local cmd rocfit 
		est scalar area =`s(area)'
		est scalar se_area =`s(se_area)'
		est scalar deltam =`s(deltam)'
		est scalar se_delm =`s(se_delm)'
		est scalar de =`s(de)'
		est scalar se_de =`s(se_de)'
		est scalar da =`s(da)'
		est scalar se_da =`s(se_da)'
		est scalar neg=`MN'
		est scalar pos=e(N)-`MN'
		global MP 
	}
	di _n in gr `"`e(title)'"' _col(51) in gr "Number of obs" /*
	*/ _col(67) "= " in ye %10.0g e(N)
	di in gr "Goodness of fit chi2(" in ye "`e(df_gof)'" /*
	*/ in gr ") =  " in ye %10.2f e(chi2_gof)
	di in gr "Prob > chi2             =  " in ye %10.4f e(p_gof)
	di in gr "Log likelihood          =  " in ye %10.0g e(ll) 
	DISTab1 `level' 
	DISTab2 `level'
end

program define ROCInd, sclass
	/* program to calculated indices */
	args level
	tempname B V Z a b va vb vab 
	mat `B'=e(b)
	mat `V'=e(V)
	scalar `a' = `B'[1,1]
	scalar `b' = `B'[1,2]
	scalar `va'=`V'[1,1]
	scalar `vb'=`V'[2,2]
	scalar `vab'=`V'[1,2]
	scalar `Z'=`a'/sqrt(1+`b'^2)
	local alpha=invnorm(.5+`level'/200)
	/* AREA */
	tempname ga gb ase area delta de da
	scalar `area'=normprob(`Z')
	scalar `ga'=normd(`Z')*(1/sqrt(1+`b'^2))
	scalar `gb'=normd(`Z')*(-`a'*`b')/((1+`b'^2)^(3/2))
	scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
	scalar `ase'=sqrt(`ase')
	sret local se_area = `ase' 
	sret local area = `area' 
	/* DELTA M */
	scalar `delta' = `a'/`b'
	scalar `ga'=1/`b'
	scalar `gb'=-`a'/(`b'^2)
	scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
	scalar `ase'=sqrt(`ase')
	sret local se_delm = `ase' 
	sret local deltam = `delta' 
	/* DE Line  */
	scalar `de' = 2*`a'/(`b'+1)
	scalar `ga'=2/(`b'+1)
	scalar `gb'=-2*`a'/((`b'+1)^2)
	scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
	scalar `ase'=sqrt(`ase')
	sret local se_de = `ase' 
	sret local de = `de' 
	/* DA Line  */
	scalar `da' = sqrt(2)*invnorm(`area')
	scalar `ga'=sqrt(2)/sqrt(1+`b'^2)
	scalar `gb'=-sqrt(2)*`a'*`b'/sqrt((1+`b'^2)^3)
	scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
	scalar `ase'=sqrt(`ase')
	sret local se_da = `ase' 
	sret local da = `da' 
end
prog def DISTab1
	args level max
	tempname B V
	qui mat `B'=e(b)
	qui mat `V'=e(V)
	local max=colsof(`B')
	di  in gr  _dup(78) "-"
	di  in gr "          |     Coef.   Std. Err."  /*
	*/ "       z     P>|z|       [`level'% Conf. Interval]" 
	di in gr "----------+" _dup(67) "-"
	local alpha=invnorm(.5+`level'/200)
	di in gr "intercept |" in ye %10.6f `B'[1,1] /*
	*/ _col(23) %10.6f sqrt(`V'[1,1]) _c
 	local ub=`B'[1,1] + `alpha'*sqrt(`V'[1,1])
 	local lb=`B'[1,1] - `alpha'*sqrt(`V'[1,1])
	local z=`B'[1,1]/sqrt(`V'[1,1])
	local p=1-normprob(abs(`z'))
	di %11.3f `z' "  "  %6.3f `p' _col(25) %10.6f `lb' _col(37) %10.6f `ub'
	di in gr "slope (*) |" in ye %10.6f `B'[1,2] /*
	*/ _col(23) %10.6f sqrt(`V'[2,2]) _c
 	local ub=`B'[1,2] + `alpha'*sqrt(`V'[2,2])
 	local lb=`B'[1,2] - `alpha'*sqrt(`V'[2,2])
	local z=(`B'[1,2]-1)/sqrt(`V'[2,2])
	local p=1-normprob(abs(`z'))
	di %11.3f `z' "  "  %6.3f `p' _col(25) %10.6f `lb' _col(37) %10.6f `ub'
	di in gr "----------+" _dup(67) "-"
	local i 3
	while `i'<=`max' {
		di in gr "    _cut"`i'-2 " |" in ye %10.6f `B'[1,`i'] /*
		*/ _col(23) %10.6f sqrt(`V'[`i',`i']) _c
	 	local ub=`B'[1,`i'] + `alpha'*sqrt(`V'[`i',`i'])
	 	local lb=`B'[1,`i'] - `alpha'*sqrt(`V'[`i',`i'])
		local z=`B'[1,`i']/sqrt(`V'[`i',`i'])
		local p=1-normprob(abs(`z'))
		di %11.3f `z' "  " %6.3f `p' _col(25) %10.6f `lb' /*
		*/ _col(37) %10.6f `ub'
		local i=`i'+1
	}
	di  in gr  _dup(78) "="
end
prog def DISTab2
	args level
	local alpha=invnorm(.5+`level'/200)
	local areau=`e(area)'+ `alpha'*`e(se_area)'
	local areal=`e(area)'- `alpha'*`e(se_area)'
	local deltau=`e(deltam)'+ `alpha'*`e(se_delm)'
	local deltal=`e(deltam)'- `alpha'*`e(se_delm)'
	local deu=`e(de)'+ `alpha'*`e(se_de)'
	local del=`e(de)'- `alpha'*`e(se_de)'
	local dau=`e(da)'+ `alpha'*`e(se_da)'
	local dal=`e(da)'- `alpha'*`e(se_da)'
	di in gr  _col(11) "|                    Indices from binormal fit" 
	di  in gr "Index     |  Estimate   Std. Err." /*
	*/ "                         [`level'% Conf. Interval]" 
	di in gr "----------+" _dup(67) "-"
	di in gr "ROC area  |" in ye %10.6f `e(area)' /*
	*/ _col(23) %10.6f `e(se_area)' /* 
	*/ _col(57) %10.6f `areal' _col(69) %10.6f `areau'
	di in gr "delta(m)  |" in ye %10.6f `e(deltam)' /*
	*/ _col(23) %10.6f `e(se_delm)' /* 
	*/ _col(57) %10.6f `deltal' _col(69) %10.6f `deltau'
	di in gr "d(e)      |" in ye %10.6f `e(de)' /*
	*/ _col(23) %10.6f `e(se_de)' /* 
	*/ _col(57) %10.6f `del' _col(69) %10.6f `deu'
	di in gr "d(a)      |" in ye %10.6f `e(da)' /*
	*/ _col(23) %10.6f `e(se_da)' /* 
	*/ _col(57) %10.6f `dal' _col(69) %10.6f `dau'
	di  in gr  _dup(78) "-"
	di  in gr  "(*) z test for slope==1"
end