/* fig5_8.do

   Figure 5.8 (Example 5.6)

     (generate and save coordinate data for graphics pkg import)

   programs required
       emroc.ado  (v 1.3.4 or later)
       dfroc.ado  (v 1.3.0 or later)

   data required: dp2.dta
      Stover/Gorga DPOAE test data

   last update:  14 Aug 2003

  get ROC estimates per
    i)   empirical
    ii)  binormal fit (Dorfman & Alf)  (rocfit)
    iii) ROC-GLM (dfroc) (also get bootstrap se estimates for parameters)

   for one frequency/level combination

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

use ${semt_data}dp2,clear

replace y = -y
tempfile f1 f2 f3 f4 f5 f6 f7 f8 f9

/*   for inspection of all 9 freq/level combinations    ******

set textsize 200
local k = 1

foreach i of numlist 55 60 65 {
    foreach j of numlist 1001 1416 2002 {
        emroc y d if l == `i' & f == `j', t2(`" level = `i' ; freq = `j' "') /*
            */ saving(`f`k'') s(..)
        local k = `k' + 1
        * pause
    }
}
gr using `f1' `f2' `f3' `f4' `f5' `f6' `f7' `f8' `f9', saving(${semt_log}fig5_8a,replace)

*******     */

log using ${semt_log}fig5_8, replace

/*  choose f (stimulus frequency) = 1001; l (stimulus intensity) = 65  */

keep if f==1001 & l==65
preserve

emroc y d, nogr gense(etpr) genspe(efpr)
replace efpr = 1-efpr
qui bys efpr etpr: drop if _n>1

* add record for origin (0,0)
gen byte expn = _n==1
expand 2 if expn
sort expn
by expn : replace efpr = 0 if expn==1 & _n==1
by expn : replace etpr = 0 if expn==1 & _n==1

keep efpr etpr
order efpr etpr
sort efpr etpr

save ${semt_log}fig5_8empir,replace

restore, preserve

/* get binormal parameters from fit with dfroc
   get b.s. se's separately via fig5_8b.do */

dfroc y d, nogr nobs

sca a_df = r(a)
sca b_df = r(b)
sca auc_df = r(auc)

/* get binormal parameters from -rocfit-

    will use cont(10) for now,

    estimates don't seem to be affected much by increasing number of categories
      with cont(#):

                              a     b
                 cont(5)   1.432  .806
                cont(10)   1.434  .818
                cont(20)   1.426  .806

*/

rocfit d y, cont(10)

matrix mat_b = e(b)
scalar a_da = mat_b[1,1]
scalar b_da = mat_b[1,2]

drop _all
set obs 200
range t 0 1
gen roc_da = normprob(a_da + b_da * invnorm(t))
replace roc_da = 0 if t==0  /* otherwise missing/undefined */
replace roc_da = 1 if t==1

gen roc_df = normprob(a_df + b_df * invnorm(t))
replace roc_df = 0 if t==0  /* otherwise missing/undefined */
replace roc_df = 1 if t==1

sort t
save ${semt_log}fig5_8par,replace

log close