/* ex6_11.do

   example 6.11
      a) model including disease-specific covariate
      b) model excluding disease-specific covariate

   ROC GLM model for DPOAE test data

   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
clear
set mem 10m
*************************************
* program to fit an ROC-GLM model on a single bootstrap sample:
capture program drop glmfit
program define glmfit
    version 7
    if "`1'"=="?" {
        global S_1 "alpha_1 beta_1 beta_2 beta_3 alpha_0"
        exit
    }
    args postname mod
    drop id
    egen int ndb = sum(d==0), by(cvlv)
    ***********
    /* STEPS 1 & 2 : calculate the covariate specific reference distn, t,
         (the survivor fxn S_db_z(Y))
         This is, equivalently, the fpr, corresponding to the non-dis Y_db_i,
         and yields the placement values for the diseased Y_d_i,
    */
    gsort -y, gen(nysort)
    qui bys cvlv (nysort): gen sumdb = sum(d==0)
    qui bys cvlv y (sumdb): gen t = sumdb[_N]
    qui replace t = t/ndb
    gen plval = t if d==1
    drop nysort sumdb
    ***********
    /* STEP 3a : restructure data to include nt (distinct t (fpr) of T_z)
       obs for each diseased obs.  (nd x nt obs), for each covariate level
    */

    bys d cvlv t: drop if d==0 & _n > 1  /* restrict non-diseased obs to distinct fpr's */
    drop if d==0 & t == 1
    egen int nt = sum(d==0), by(cvlv)

    gen uid = _n  /* unique id for each record */
    quietly bys d cvlv (y): gen id_t =_n if d==0 /* id_t: unique id each t within cvlv */

    expand nt if d==1
    sort uid

    quietly by uid: replace id_t = _n if d==1  /* assigns a set of d's to each fpr */

    bys cvlv id_t (d) : replace t = t[1]

    drop if d==0

    ***********
    /* STEP 3b : calculate the binary placement value indicators */

    gen byte u_it = (plval <= t)

    ***********
    /* STEP 4 : calculate h_s(t) */

    gen phiinvt = invnorm(t)

    *********
    /* STEP 5 : fit the GLM model */

    if "`mod'" == "a" {
        probit u_it phiinvt xf xl xd
        post `1' (_b[phiinvt]) (_b[xf]) (_b[xl]) (_b[xd]) (_b[_cons])
    }
    else {
        probit u_it phiinvt xf xl
        post `1' (_b[phiinvt]) (_b[xf]) (_b[xl]) (1) (_b[_cons])
    }
end

*************************************************************
/* SET UP DATA */

use ${semt_data}dp2,clear

replace y = -y
gen xd = ((amt-20)*d)/10 if d==1 /* cov specific to diseased obs */

#delimit ;
gen cvlv= 1 * ((f==1001)*(l==55))
        + 2 * ((f==1001)*(l==60))
        + 3 * ((f==1001)*(l==65))
        + 4 * ((f==1416)*(l==55))
        + 5 * ((f==1416)*(l==60))
        + 6 * ((f==1416)*(l==65))
        + 7 * ((f==2002)*(l==55))
        + 8 * ((f==2002)*(l==60))
        + 9 * ((f==2002)*(l==65));

/* cvlv = unique levels of cov X common to d & db */;

#delimit cr

gen  xf = f/100
gen  xl = l/10
drop f l amt

compress _all
set seed 156385

foreach model in a b {
    log using ${semt_log}ex6_11`model', replace
    set rmsg on
    bstrap glmfit, cluster(id) reps(5000) args(`model') saving(${semt_log}ex6_11`model') replace
    qui log close
}
set rmsg off