/* fig5_5.do

   from example 5.4

   data required: wiedat2b.dta
                  (Wieand Pancreatic cancer data)

   programs required:
       emroc  vers 1.3.5 or later

   Semi-parametric ROC's for Wieand data
     (along with empirical ROC's)

   see fig5_5b.do for bootstrap estimation of confidence limits for ROC(.2)
     on this figure

   last update:  30 Apr 2003

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

cap pro drop doit
pro def doit
    version 7
    /* get semi-parametric ROC estimates for Wieand data */
    clear
    use ${semt_data}wiedat2b
    log using ${semt_log}fig5_5, replace
    qui log off

    replace y1 = ln(y1)
    replace y2 = ln(y2)

    qui log on

    foreach k of numlist 1 2 {
        qui count if y`k' ~= .
        local ntot = r(N)
        qui sum y`k' if d==0
        scalar mu`k'_db = r(mean)
        scalar s`k'_db = r(sd)
        qui sum y`k' if d==1
        scalar mu`k'_d = r(mean)
        scalar s`k'_d = r(sd)

        * r: standardized residuals w.r.t sample mean, conditional on disease status

        qui egen mn`k' = mean(y`k'), by(d)
        qui egen sd`k' = sd(y`k'), by(d)
        qui gen r`k' = (y`k' - mn`k')/sd`k'

        * t: survivor fxn (i.e 1-cdf) of residuals (w/o respect to disease status)
        *  thus r is (S_hat_0)_inv in section 5.3.2

        gsort -r`k'
        qui gen int sum`k' = sum(r`k' ~= .)
        qui bys r`k' (sum`k'): gen t`k' = sum`k'[_N]
        qui replace t`k' = t`k'/`ntot'

        * log vbls up to this point

        if `k' ==1 {
            di "scalars corresponding to test measure y`k' "
            di
            sca list
            di
            l y1 r1 sum1 t1
            di
        }

        * post-9/20 msg from MP:

        gen rocterm`k' = ((mu`k'_db - mu`k'_d)/s`k'_d ) + (s`k'_db/s`k'_d) * r`k'
        sort r`k'

        forvalues i = 1/`ntot' {
            qui count if r`k' >= rocterm`k'[`i']
            qui replace sum`k' = r(N) if _n == `i'
        }
        qui gen roc`k' = sum`k'/`ntot'

        if `k' ==1 {
            di
            di "vbl list after calculation of rocterm and roc: "
            di
            l y1 r1 t1 sum1 rocterm1 roc1
            di
        }

    }
    qui log close
end
*********************************************************
cap pro drop doit2
pro def doit2
    /* need to sort each of 4 pairs of plotting vbls here
       and generate sort indicator for merging */
    syntax varlist, Tf(string)
    preserve
    keep `varlist'
    order `varlist'
    qui bys `varlist': drop if _n>1
    sort `varlist'
    gen int mergind = _n
    sort mergind
    save `tf'
end
***********************************************************

doit

* empirical ROC estimates:

emroc y1 d, nogr gense(eroc1) gensp(spec)
gen et1 = 1 - spec
drop spec

emroc y2 d, nogr gense(eroc2) gensp(spec)
gen et2 = 1 - spec
drop spec

unab vlist: t* roc1 roc2 er* et*
keep `vlist'

* add record for origin (0,0)

gen byte expn = _n==1
expand 2 if expn
sort expn

di "vlist: `vlist'"
d `vlist'

foreach var of local vlist {
    by expn : replace `var' = 0 if expn==1 & _n==1
}

tempfile f1 f2 f3 f4

doit2 t1 roc1, tf(`f1')
doit2 t2 roc2, tf(`f2')
doit2 et1 eroc1, tf(`f3')
doit2 et2 eroc2, tf(`f4')

use `f1',clear
merge mergind using `f2'
drop _merge
sort mergind

merge mergind using `f3'
drop _merge
sort mergind

merge mergind using `f4'
drop _merge mergind

save ${semt_log}fig5_5,replace