/* ms2_fig4.do

  do-file to generate figure 4 for SJ manuscript:
    "Accomodating covariates in ROC analysis".
       Janes H, Longton G, and Pepe M.

  use simulated marker & covariate data: sj_ms2_fig4.dta

  uses -rocreg- version 1.0.8 or later

  last update: 4/22/2008

*/
version 10
cap log close
************************
cap pro drop doit
pro def doit
    set scheme sj
    use sj_ms2_fig4, clear // see fig4_init.do

    range x -4 7 500

    // get kernal density estimates for simulated marker variable
    //    by covariate level and disease status

    forvalues z = 0/1 {
        forvalues d = 0/1 {
            kdensity y if z==`z' & d==`d', width(.3) at(x) ///
              gen(xtmp fx_z`z'_d`d') nogr
            la var fx_z`z'_d`d' ""
            assert xtmp == x
            drop xtmp
        }
    }

    log using sjms2_fig4bnos.txt, replace

    // find cutoff c corresponding to fpr = .2 for z=0

    local fpf0 = .2
    bys z d (y): gen posrate = (_N + 1 - _n)/_N
    sum y if d==0 & z==0 & posrate <= `fpf0'
    local c = r(min)

    // get fpf1, tpf0, fpf1 corresponding to cutoff c

    sum posrate if d==1 & z==0 & y>=`c'
    local tpf0 = r(max)
    sum posrate if d==1 & z==1 & y>=`c'
    local tpf1 = r(max)
    sum posrate if d==0 & z==1 & y>=`c'
    local fpf1 = r(max)

    di "    c: " %4.3f `c'
    di
    di " fpf0: " %6.5f `fpf0'
    di " tpf0: " %6.5f `tpf0'
    di
    di " fpf1: " %6.5f `fpf1'
    di " tpf1: " %6.5f `tpf1'

    // get fit for ROC curves, stratify on, and model z

    rocreg d y, adjcov(z) regcov(z) nobs

    qui log off

    matrix B = e(GLMparm)
    matrix list B

    sca alpha_0 = B[1,1]
    sca alpha_1 = B[1,2]
    sca alpha_2 = B[1,3]

    qui log close

    range t 0 1 500

    **************
    local subt1 "Z = 0"
    local subt2 "Z = 1"

    local xax xscale(r(-3 6.5)) xlab(0 5)
    #delimit ;
    local legend1 `"   legend(order( 1 2 )
              ring(0) pos(2) cols(1) size(3.3)
              label(1 "controls") label(2 "cases")
              keygap(2) rowgap(.6)
              symx(8)
              bmargin(t=2 r=0)
              region(style(none)) ) "'
    ;
    #delimit cr
    local legend2 "legend(off)"

    local txt1 ""
    local txt2 ""

    forvalues i = 1/2 {
        local z = `i' - 1
        tempname g`i'
        #delimit ;
        scatter fx_z`z'_d0 x, connect(l) m(none) lwidth(*1.5) lp(shortdash) lc(gs0) ||
        scatter fx_z`z'_d1 x, connect(l) m(none) lwidth(*1.5) lp(solid) lc(gs0)
            `xax'
            xtitle("")
            xline(`c', lw(*1) lc(gs9) )
            plotregion(style(none))
            yscale(noline r(0,0.405)) ylab(none)
            subtitle("`subt`i''",ring(0) pos(10) margin(l=0 t=10) justification(left) size(3.5))
            `txt`i''
            `legend`i''
            scheme(s1mono)
            name(`g`i'')
        ;

        #delimit cr
    }

    #delimit ;

    graph combine `g1' `g2', cols(1) nocopies
      iscale(*1.2) imargin(t=1.5 b=1.5)
      l1("density", size(5) margin(b=8 r=3))
      graphregion(margin(r=5))
      name(ga, replace)
    ;

    twoway function y = normal(alpha_0 + alpha_1*invnormal(x) + alpha_2*0),
             range(1e-09 1) clw(*1.5) clp(solid) lc(gs0) ||
           function y = normal(alpha_0 + alpha_1*invnormal(x) + alpha_2*1),
             range(1e-09 1) clw(*1.5) clp(solid) lc(gs0) ||
           scatteri `tpf0' `fpf0', m(O) mc(gs0) msize(2) ||
           scatteri `tpf1' `fpf1', m(O) mc(gs0) msize(2)
           text( `tpf0' `fpf0' "TPR  , FPR ", place(3) justification(left) margin(l=4 ) size(2.8))
           text( `tpf0' `fpf0' "0", place(3) justification(left)
                  margin(l=9.3 t=2 ) size(2))
           text( `tpf0' `fpf0' "0", place(3) justification(left)
                  margin(l=18 t=2 ) size(2))
           text( `tpf1' `fpf1' "TPR  , FPR ", place(3) justification(left) margin(l=4 ) size(2.8))
           text( `tpf1' `fpf1' "1", place(3) justification(left)
                  margin(l=9.3 t=2 ) size(2))
           text( `tpf1' `fpf1' "1", place(3) justification(left)
                  margin(l=18 t=2 ) size(2))
           text( .65 .08 "Z = 0", place(3) justification(left) size(3.2))
           text( .54 .305 "Z = 1", place(3) justification(left) size(3.2))
           text( .25 .65 "covariate-specific" "ROC curves", place(c)
                  justification(center) size(4) linegap(1))
         xlab(0 1, grid) ylab(0 1, angle(horizontal) )
         ytitle(TPR, size(4) justification(center) orientation(horizontal))
         yscale(titlegap(2))
         xtitle(FPR, size(4) justification(center) margin(b=0))
         xscale(r(0 1) titlegap(1.5))
         legend(off)
         aspectratio(1)
         plotregion(margin(l=2 r=2 t=2 b=2) )
         graphregion(margin(l=5 t=10))
         name(gb, replace)
     ;
    graph combine ga gb, rows(1) nocopies
        graphregion(margin(t=17 b=17))
        name(fg1, replace)
    ;
    #delimit cr
end
****************
doit