/* fig5_4.do

   generate and save plot coordinate data for graphics pkg import

   Example 5.2

   programs required
       emroc.ado  (v 1.3.4 or later)
       aucbs.ado  (v 1.3.3 or later)

   last update:  30 Apr 2003

*/
version 7
set more off
cap semt_profile
cap log close
**************************************
cap pro drop doit
pro def doit, rclass
    args y d tpr fpr kwd
    tempvar cutind se selind
    emroc2 `y' `d', conf gense(`tpr') gensp(`fpr') se(`se') nogr kwidth(`kwd')
    replace `fpr' = 1 - `fpr'

    *******  get roc(t) and se ******

    gen byte `cutind' = (fpr <= .2)
    qui bys `cutind' (`fpr' `tpr'): gen byte `selind' = (_n == _N & `cutind' == 1)
    sum `tpr' if `selind'==1, meanonly
    return scalar roct = r(mean)
    sum `se' if `selind'==1,meanonly
    * pause
    return scalar se = r(mean)

    **************
    bys `fpr' `tpr': drop if _n>1

    * restricted to unique fpr,tpr pairs:

    * add record for origin (0,0)
    gen byte expn = _n==1
    expand 2 if expn
    sort expn
    by expn : replace `fpr' = 0 if expn==1 & _n==1
    by expn : replace `tpr' = 0 if expn==1 & _n==1
    sort `fpr' `tpr'

    * empirical ROC pt data for plot

    keep `fpr' `tpr'
    order `fpr' `tpr'
    sort `fpr' `tpr'
end
*********************************************
cap pro drop doit2
pro def doit2
    clear
    use ${semt_data}wiedat2b,clear
    log using ${semt_log}fig5_4, replace
    qui log off

    tempfile f1
    save `f1'
    local alpha = .10
    tempname roct1 roct2 se1 se2
    doit y1 d tpr1 fpr1 20
    scalar `roct1' = r(roct)
    scalar `se1' = r(se)

    save ${semt_log}fig5_4a,replace

    use `f1',clear

    doit y2 d tpr2 fpr2 20
    scalar `roct2' = r(roct)
    scalar `se2' = r(se)
    sca l
    d
    save ${semt_log}fig5_4b,replace

    /* logit based CI's for roc(t) at t=.20

      editing needed beyond this point (pasted from fig5_2.do)

      use emroc2 to get se's ?,   aucbs to get roc(20)?

    */

    scalar z = invnorm(1.0 - `alpha'/2)

    foreach n of numlist 1/2 {
        scalar lba_`n' = `roct`n'' - z*`se`n''
        scalar uba_`n' = `roct`n'' + z*`se`n''
    }

    /* logit based CI's */

    foreach n of numlist 1/2 {
        scalar logitroc`n' = log( (`roct`n'')/(1.0 - `roct`n''))
        scalar lbb_`n' = logitroc`n' - z * (`se`n''/ (`roct`n''*(1 - `roct`n'')))
        scalar lbb_`n' = exp(lbb_`n')/(1.0 + exp(lbb_`n'))
        scalar ubb_`n' = logitroc`n' + z *(`se`n''/ (`roct`n''*(1 - `roct`n'')))
        scalar ubb_`n' = exp(ubb_`n')/(1.0 + exp(ubb_`n'))
    }
    qui log on
    di
    di "   roc           CI                logit-based CI  "
    di
    foreach n of numlist 1/2 {
        di "  " %5.4f `roct`n'' "    (" %5.4f lba_`n' ", " %5.4f uba_`n' ")      (" %5.4f lbb_`n' ", " %5.4f ubb_`n' ")"
    }
    di
    di
    qui log close
end
********************
doit2