/* ex5_7.do

   bootstrap binormal GLM-ROC parameter estimates (from -dfroc-)
        for both tests - and carry out hypothesis tests on difference statistics

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

   programs required:
       dfroc  vers 1.3.0 or later
       aucbs  vers 1.3.3 or later

   last update:  30 Apr 2003

*/
version 7
set more off
cap semt_profile
cap log close
************************************************************************
cap pro drop doit3
pro def doit3
    /* driver prog */
    args nsamp
    clear
    use ${semt_data}wiedat2b

    log using ${semt_log}ex5_7,replace
    qui log off

    local resfile "${semt_log}ex5_7res"

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

    doit2 `nsamp' `resfile'
    use `resfile', clear

    qui gen adiff = a1 - a2
    qui gen bdiff = b1 - b2
    qui gen ap20diff = a1p20 - a2p20
    qui gen bp20diff = b1p20 - b2p20

    local ta1 = a1[1]
    local tb1 = b1[1]
    local ta2 = a2[1]
    local tb2 = b2[1]
    local tadiff = adiff[1]
    local tbdiff = bdiff[1]

    local ta1p20 = a1p20[1]
    local tb1p20 = b1p20[1]
    local ta2p20 = a2p20[1]
    local tb2p20 = b2p20[1]
    local tap20diff = ap20diff[1]
    local tbp20diff = bp20diff[1]

    char a1[bstrap] `ta1'
    char b1[bstrap] `tb1'
    char a2[bstrap] `ta2'
    char b2[bstrap] `tb2'
    char adiff[bstrap] `tadiff'
    char bdiff[bstrap] `tbdiff'

    char a1p20[bstrap] `ta1p20'
    char b1p20[bstrap] `tb1p20'
    char a2p20[bstrap] `ta2p20'
    char b2p20[bstrap] `tb2p20'
    char ap20diff[bstrap] `tap20diff'
    char bp20diff[bstrap] `tbp20diff'

    qui drop in 1
    qui save `resfile',replace

    /* retrieve coeff estimates derived from original dataset
       from characteristics and put them into a vector */

    matrix Beta = (`ta1', `tb1', `ta2', `tb2', `ta1p20', `tb1p20', `ta2p20', `tb2p20',  /*
           */      `tadiff', `tbdiff', `tap20diff', `tbp20diff')
    matrix colnames Beta = a1 b1 a2 b2 a1p20 b1p20 a2p20 b2p20 adiff bdiff ap20diff bp20diff

    /* obtain Covariance matrix based on bootstrap results */

    matrix accum Cov = _all, deviations noconst
    local nobs = r(N)
    scalar adj = 1/(`nobs'-1)
    matrix Cov = adj*Cov

    /* post coeffiecient vector and Covariance to Stata system area
       for later access by post-estimation commands */

    matrix post Beta Cov,obs(`nobs')

    **************************
    qui log on

    di
    di " example 5.7
    di
    di " a) bootstrap results for binormal GLM-ROC parameter estimates (from -dfroc-)
    di "    for both tests (CA 19-9 and CA 125)
    di "     - carry out hypothesis tests on difference statistics
    di
    bstat
    di
    bstat, level(90)
    di
    matrix mlout
    di
    di
    di " 1) Ho: (adiff, bdiff) = (0,0)  (per eq 5.17) "
    di
    test adiff bdiff
    di
    di " 2) Ho: (ap20diff, bp20diff) = (0,0) "
    di
    test ap20diff bp20diff
    di
    di
    di "*************************************************************
    di " b) partial AUC (t<=20) based on empirical ROC, for both tests
    di
    qui log off
    use ${semt_data}wiedat2b
    replace y1 = ln(y1)
    replace y2 = ln(y2)
    qui log on
    aucbs y2 y1 d, cc nsamp(`nsamp') partial(20) res(${semt_log}resex5_7b) replace
    di
    di "  test of Ho: difference in partial AUC = 0
    di "     based on bootstrap distribution of difference
    di
    qui use ${semt_log}resex5_7b,clear
    qui count if aucdelta~=.
    local ngood = r(N)
    qui count if aucdelta < 0
    local p = 2 * r(N)/`ngood'
    di
    di "     number of valid bootstrap samples:   " %4.0f `ngood'
    di "  number of samples with aucdelta < 0 :   " %4.0f `r(N)'
    di "                                    p = " %6.0g `p'
    di
    qui log close
end
*****************************************************
cap pro drop doit2
pro def doit2
    version 7
    * arguments:   nsamp = # bootstrap samples
    *              resfile : bootstrap result file
    args nsamp 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' a1 b1 a2 b2 a1p20 b1p20 a2p20 b2p20 using `resfile',replace

    * post results from oberved data:

    restore
    doit
    post `pfile' (r(a1)) (r(b1)) (r(a2)) (r(b2)) /*
           */ (r(a1p20)) (r(b1p20)) (r(a2p20)) (r(b2p20))

    * bootstrap samples:

    forval i = 1/`nsamp' {
        getsamp `caseobs' `contobs'
        doit
        post `pfile' (r(a1)) (r(b1)) (r(a2)) (r(b2))  /*
               */ (r(a1p20)) (r(b1p20)) (r(a2p20)) (r(b2p20))
    }
    postclose `pfile'
end
***************************************************************************
cap pro drop doit
pro def doit, rclass
    /* get binormal ROC parameter estimates from -dfroc- for both tests */
    qui dfroc y1 d, nog nobs
    return sca a1 = r(a)
    return sca b1 = r(b)
    qui dfroc y2 d, nog nobs
    return sca a2 = r(a)
    return sca b2 = r(b)
    qui dfroc y1 d, nog nobs par(20)
    return sca a1p20 = r(a)
    return sca b1p20 = r(b)
    qui dfroc y2 d, nog nobs par(20)
    return sca a2p20 = r(a)
    return sca b2p20 = r(b)
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
**************************************************

doit3 5000