/* ex6_7.do

   example 6.7
   figure 6.6

   example 6.7: figure 6.6

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

   last update:  23 Apr 2003
*/
version 7
set more off
cap semt_profile
cap log close
pause on

tempfile fa fb

log using ${semt_log}ex6_7,replace
use ${semt_data}dp2,clear

replace y=-y
gen xl=l/10
gen xf=f/100
gen xd=(amt-20)/10 if d==1
bys d:sum x* y
*********************************************************************
* a) normal linear model curve (disease severity covariate omitted):

regress y xf xl if d==0
sca a0_db  = _b[_cons]
sca a1_db  = _b[xf]
sca a2_db  = _b[xl]
sca mse_db = e(rmse)
predict rhat0 if d==0,rstandard

regress y xf xl if d==1
sca a0_d  = _b[_cons]
sca a1_d  = _b[xf]
sca a2_d  = _b[xl]
sca mse_d = e(rmse)

sca a0 = (a0_d - a0_db)/mse_d
sca a1 = (a1_d - a1_db)/mse_d
sca a2 = (a2_d - a2_db)/mse_d
sca b  = mse_db/mse_d

sca l

predict rhat1 if d==1,rstandard
gen res = rhat1 if d==1
replace res = rhat0 if d==0
pnorm res,grid
qui log off
preserve /* will need residuals for part (b) */

* plot binormal curves based on parameters from the linear models:

drop _all
range t 0 1 500
local xf = 10.01  /* corresponding to stimulus frequency f = 1001 */

foreach l in 55 60 65 {
    local xl = `l'/10.
    gen roc_a`l' = norm(a0 + a1*`xf' + a2*`xl' + b*invnorm(t))
    replace roc_a`l' = t if t==0 | t==1
}

sort t
gen forty5 = t if t==0 | t==1

sort t roc*
gr roc* t, s(...) c(lll)
keep t roc*
sort t
save `fa'

* pause

**************************************************************
* b) semiparametric linear model curve:

restore

* first, restrict to observations with freq = 1001:
keep if f==1001

preserve

/* z = (1-t) quantile for the residuals:
   so, define k_i = quantile of resid_i; k_i = fraction of r <= r_i
   i.e. k is the cdf of the resids
*/

local xf = 10.01  /* corresponding to stimulus frequency f = 1001 */

tempfile f55 f60 f65

foreach l in 55 60 65 {
    keep if l==`l'
    local xl = `l'/10.
    qui count if res ~= .
    local n = r(N)

    gsort -res
    qui gen int sum_r = sum(res~=.)
    qui bys res (sum_r): gen t = sum_r[_N]
    qui replace t = t/`n'

    gen lincomb = -(a0 + a1*`xf' + a2*`xl') + b*res
    gen roc_b`l' = .
    assert res ~= . & lincomb ~= .
    forvalues i = 1/`n' {
        qui count if res >= lincomb[`i']
        qui replace roc_b`l' = r(N) if _n == `i'
    }
    replace roc_b`l' = roc_b`l'/`n'

    * restrict to unique fpr,tpr pairs:
    bys t roc_b`l': drop if _n>1

    * add record for origin (0,0)
    gen byte expn = _n==1
    expand 2 if expn
    sort expn
    by expn : replace t = 0 if expn==1 & _n==1
    by expn : replace roc_b`l' = 0 if expn==1 & _n==1

    keep t roc_b`l'
    sort t roc_b`l'
    save `f`l''
    restore, preserve
}
use `f55',clear
merge t using `f60'
drop _merge
sort t
merge t using `f65'
drop _merge

sort t
gen forty5 = t if _n==1 | _n==_N
sort t roc*
gr roc* forty5 t, s(....) c(llll)

keep t roc*
sort t
save `fb'

**********************************************
* c) empirical ROC curves for f = 1001:

restore, preserve
forvalues i = 55(5)65 {
    tempfile f`i'

    di "i: `i'"
    ta d

    emroc y d if l ==`i', gensens(roc_c`i') genspec(t)
    pause q <return> to continue . . .
    replace t = 1 - t

    * restricted to unique fpr,tpr pairs:
    bys t roc_c`i': drop if _n>1

    * add record for origin (0,0)
    gen byte expn = _n==1
    expand 2 if expn
    sort expn
    by expn : replace t = 0 if expn==1 & _n==1
    by expn : replace roc_c`i' = 0 if expn==1 & _n==1
    sort t roc_c`i'

    * empirical ROC pt data for plot

    keep  t roc_c`i'
    order t roc_c`i'
    sort  t roc_c`i'

    save `f`i''
    restore, preserve
}
restore

use `f55', clear
merge t using `f60'
drop _merge
sort t roc_c55 roc_c60
merge t using `f65'
drop _merge
sort t roc_c55 roc_c60 roc_c65

merge t using `fa'
sort t roc_c55 roc_c60 roc_c65 roc_a55
drop _merge
merge t using `fb'
drop _merge

sort t roc_c55 roc_c60 roc_c65
gen forty5 = t if _n==1 | _n==_N
set textsize 200
tempfile g55 g60 g65
foreach l in 55 60 65 {
    #delimit ;
    gr roc_a`l' roc_b`l' roc_c`l' forty5 t,
       s(....) c(lsll) border
       xscale(0,1) yscale(0,1)
       yline(.2(.2).8) xline(.2(.2).8)
       xlab(0(.2)1) ylab(0(.2)1)
       t1("stimulus intensity = `l'")
       saving(`g`l'')
    ;
    #delimit cr
    pause q <return> to continue . . .
}
set textsize 125
gr using `g55' `g60' `g65', margin(5) t1("Fig. 6.6") saving(${semt_log}fig6_6,replace)

pause q <return> to continue . . .
save ${semt_log}fig6_6,replace

qui log close