/* fig5_5b.do

   from example 5.4

   bootstrap estimation of confidence limits for semiparametric ROC(.2)
      (Wieand data)

   see fig5_5.do for generation of sample semiparametric and empirical ROC estimates
      for this figure

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

   last update:  30 Apr 2003
*/
version 7
cap log close
cap semt_profile
set more off
***************************************************************************
cap pro drop doit
pro def doit, rclass
    /* get semi-parametric ROC estimates for Wieand data */
    version 7
    * args y1 y2 d t
    args t
    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'

        qui gen r`k' = (((y`k' - mu`k'_db)/s`k'_db ) * (d==0) ) +  /*
                 */    (((y`k' - mu`k'_d) /s`k'_d )  * (d==1) )

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

        * 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'
    }
    qui sum roc1 if t1 <= `t',  meanonly
    return sca roct1 = r(max)
    qui sum roc2 if t2 <= `t',  meanonly
    return sca roct2 = r(max)

    drop t1 t2 roc1 roc2 r1 r2 sum1 sum2 rocterm1 rocterm2
end
***********************************************************
cap pro drop getsamp
pro def getsamp
    version 7
    args caseobs contobs
    tempfile cassamp

    qui use `caseobs',clear
    bsample
    qui save `cassamp'
    qui use `contobs',clear
    bsample
    qui append using `cassamp'
end
**************************************************
cap pro drop doit2
pro def doit2
    version 7
    * arguments:   nsamp = # bootstrap samples
    *              t = fpr corresponding to desired roc(t) estimate
    *              resfile : bootstrap result file
    args nsamp t resfile
    tempfile caseobs cassamp contobs

    preserve
    qui keep if d==1
    qui save `caseobs'
    restore, preserve
    qui keep if d==0
    qui save `contobs'

    tempname pfile
    postfile `pfile' roc1 roc2 using `resfile',replace

    * post results from oberved data:

    restore
    doit .2
    post `pfile' (r(roct1)) (r(roct2))

    * bootstrap samples:

    forval i = 1/`nsamp' {
        getsamp `caseobs' `contobs'
        doit .2
        post `pfile' (r(roct1)) (r(roct2))
    }

    postclose `pfile'
end
************************************************************************
cap pro drop doit3
pro def doit3
    /* driver */
    version 7
    args nsamp
    clear
    use ${semt_data}wiedat2b
    log using ${semt_log}fig5_5b, replace
    qui log off
    local resfile "${semt_log}fig5_5bres"

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

    * set trace on
    doit2 `nsamp' .2 `resfile'
    use `resfile', clear

    gen rocdelta = roc1 - roc2
    local x1 = roc1[1]
    local x2 = roc2[1]
    local x3 = rocdelta[1]
    char roc1[bstrap] `x1'
    char roc2[bstrap] `x2'
    char rocdelta[bstrap] `x3'

    qui drop in 1
    qui save `resfile',replace

    count if rocdelta > 0
    local p = 2 * (1 - (r(N)/_N))

    qui log on
    di
    bstat
    di
    bstat, level(90)
    di
    di
    di " 2-sided test for Ho: delta(roc) = 0      p = " in y `p'
    di
    qui log close
end
*****************************************************
doit3 5000