/* ex5_2.do

   Example 5.2

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

   data required: wiedat2b.dta

   last update:  04 Apr 2003
*/
version 7
set more off
cap semt_profile
cap log close
*****************************************************
cap pro drop doit
pro def doit
    log using ${semt_log}ex5_2,replace
    di
    di "example 5.2"
    di
    use ${semt_data}wiedat2b,clear
    tempname auc1 auc2 se1 se2 logitauc1 logitauc2 alpha z vauc_del auc_del t

    di
    di " a) AUC's and 95% CI's based on logit(AUC) "

    scalar `alpha' = .05

    /* logit based CI's */

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

    forvalues n = 1/2 {
        /* get AUC and se estimates from -emroc- */
        qui emroc y`n' d
        sca `auc`n'' = r(auc)
        sca `se`n'' = sqrt(r(aucvar))
        * these can alternatively be obtained from -roctab- :
        * qui roctab d y`n'
        * sca `auc`n'' = r(area)
        * sca `se`n''  = r(se)
        scalar `logitauc`n'' = log( (`auc`n'')/(1.0 - `auc`n''))
        scalar lbb_`n' = `logitauc`n'' - `z' * (`se`n''/ (`auc`n''*(1 - `auc`n'')))
        scalar lbb_`n' = exp(lbb_`n')/(1.0 + exp(lbb_`n'))
        scalar ubb_`n' = `logitauc`n'' + `z' *(`se`n''/ (`auc`n''*(1 - `auc`n'')))
        scalar ubb_`n' = exp(ubb_`n')/(1.0 + exp(ubb_`n'))
    }

    di
    di "               AUC       logit-based CI  "
    di
    di "    CA 19-9  " %5.4f `auc1' "     (" %5.4f lbb_1 ", " %5.4f ubb_1 ")"
    di "    CA 125   " %5.4f `auc2' "     (" %5.4f lbb_2 ", " %5.4f ubb_2 ")"

    di
    di
    di " b) AUC difference, with variance based on paired data: "

    sca `auc_del' = `auc2' - `auc1'
    vardelta y1 y2 d
    sca `vauc_del' = r(vaucdel)
    sca `t' = `auc_del'/sqrt(`vauc_del')
    local p = 2 * (1-normprob(abs(`t')))
    di
    di "    AUC difference: " %5.3f `auc_del' " ;    se(delta): " %9.0g sqrt(`vauc_del')
    di "                 T: " %5.3f `t'
    di "                 p: " %9.2e `p'
    di
    di
    di " c) test for ROC(t) difference, at t=.20
    di "    based on bootstrap distribution of the difference:
    di
    aucbs y1 y2 d, nsamp(1000) roc(20) ccsamp
    sca `t' = r(rocdelta)/r(se_rocdl)
    local p = 2 * (1-normprob(abs(`t')))
    di
    di "    Ho: roc1_t = roc2_t "
    di "     T = " %4.3f `t'
    di "     p = " %9.2e `p'
    di
    di
    di " d) test for partial AUC difference, for t <=.20
    di "    based on bootstrap distribution of the difference:
    di
    aucbs y1 y2 d, nsamp(1000) partial(20) ccsamp
    sca `t' = r(aucdelta)/r(se_aucdl)
    local p = 2 * (1-normprob(abs(`t')))

    di
    di "    Ho: pAUC_1(20) = pAUC_2(20) "
    di "     T = " %4.3f `t'
    di "     p = " %9.2e `p'
    di

    qui log close

end
*****************************************************
cap pro drop vardelta
pro def vardelta, rclass
    /* program to calculate and return variance of AUC difference for paired data
        see ch.5 eq. 5.11 */
    args y1 y2 d
    tempname vard
    tempvar v1 v2 vd1 vd2 delta_v delta_vd
    qui count if `d'==0
    local ndb = r(N)
    qui count if `d'==1
    local nd  = r(N)
    forvalues n = 1/2 {
        gsort `y`n'' -`d'
        qui gen double `vd`n'' = 1 - (sum(`d'==0)/`ndb')
        gsort -`y`n'' -`d'
        qui gen double `v`n'' = sum(`d'==1)/`nd'
    }
    gen `delta_v' = `v1' - `v2'
    gen `delta_vd' = `vd1' - `vd2'

    qui sum `delta_vd' if `d'==1
    sca `vard' = r(Var)
    qui sum `delta_v' if `d'==0
    return scalar vaucdel = `vard'/`nd' + r(Var)/`ndb'
end
*****************************************************

doit