/* fig5_2.do

  example 5.1:  fig 5.2 a and b
    use (tpr,fpr) for empirical ROC curve from  fig5_1.do
    for these figures

  this file generates confidence bounds and intervals

  programs required
       emroc.ado  (v 1.3.4 or later)

  data required: orchratio2.dta

  last update 11 July 05

*/
version 7
set more off
cap semt_profile
cap log close

use g128 d using ${semt_data}orchratio2,clear
rename g128 y
la var y "gene 128"
compress
sort d y
assert y ~=.

log using ${semt_log}fig5_2, replace

centile y if d==0, c(25 50 75)

local cut1 = r(c_1)
local cut2 = r(c_2)
local cut3 = r(c_3)

log off

gen byte y1 = (y >= `cut1')
gen byte y2 = (y >= `cut2')
gen byte y3 = (y >= `cut3')

log on

foreach yvar in y1 y2 y3 { binscrn1 `yvar' d, level(90) f(%5.3g) }

log off

* plot b
* set trace on

emroc y d, gense(tpr) genspe(spec) conf se(se) kwid(.3)
gen fpr = 1-spec

sort fpr tpr

log on

l y d fpr tpr se

log off

/* 5.2 b & c: confidence intervals for roc(t)
   will use estimates from the listing generated above

   ROC(t) for t = .25, .50, .75

*/
*****************************
cap pro drop doit
pro def doit
    local alpha = .10    /* 90% confidence intervals specified */

    local t1 = .25
    local t2 = .50
    local t3 = .75

    local roct1 = .733
    local roct2 = .90
    local roct3 = .933

    local se1 = .104
    local se2 = .0865
    local se3 = .071

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

    foreach n of numlist 1/3 {
        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/3 {
        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'))
    }
    log on
    di
    di "    roc            CI              logit-based CI  "
    di
    foreach n of numlist 1/3 {
        di "  " %5.4f `roct`n'' "    (" %5.4f lba_`n' ", " %5.4f uba_`n' ")      (" %5.4f lbb_`n' ", " %5.4f ubb_`n' ")"
    di
    di
    sca list _all
    di
    }
end
****************************
doit

log close