/* ex7_3.do

   example 7.3

   programs required:
      binscrn1.ado

   data required: est1.dta
      Weiner - exercise stress test data from CASS study

   last update:  31 Oct 2005
*/
version 9
set more off
cap semt_profile
cap log close
*************************************
cap pro drop doit
pro def doit
    log using ${semt_log}ex7_3, replace
    use ${semt_data}est1, clear
    di
    di "**** example 7.3 ****"
    di

    tempname tpf_est tpf_cph log_or
    tab y1 d, col
    qui binscrn1 y1 d
    sca `tpf_est' = r(tp)

    di
    di "      TPF (EST) : " %4.1f 100*`tpf_est'
    di
    tab y2 d, col nokey
    qui binscrn1 y2 d
    sca `tpf_cph' = r(tp)

    di
    di "      TPF (CPH) : " %4.1f 100*`tpf_cph'
    di

    sca `log_or' = ln( (`tpf_est'/(1 - `tpf_est'))/(`tpf_cph'/(1 - `tpf_cph')) )

    di
    di "  log odds ratio : " %4.2f `log_or'
    di
    di "****************
    di "Assume now that only 30% of those with negative EST (randomly sampled) "
    di "  received the gold standard angiogram for assessment of coronary artery "
    di "  disease outcome "
    di

    set seed 62461    /* for reproducibility of "random" sample */
    gen rn = uniform() if y1==0
    bys y1 (rn): gen byte validate = ( y1==1 | (y1 == 0 & _n/_N <= .30) )
    qui replace d = . if  validate == 0  /* d available only if validate == 1 */

    tab y1 d, col m nokey
    qui binscrn1 y1 d
    sca `tpf_est' = r(tp)

    di
    di "      TPF (EST) : " %4.1f 100*`tpf_est'
    di
    tab y2 d, col m nokey
    qui binscrn1 y2 d
    sca `tpf_cph' = r(tp)

    di
    di "      TPF (CPH) : " %4.1f 100*`tpf_cph'
    di

    sca `log_or' = ln( (`tpf_est'/(1 - `tpf_est'))/(`tpf_cph'/(1 - `tpf_cph')) )

    di
    di "  log odds ratio : " %4.2f `log_or'
    di

    di "****************************
    di "log odds ratio from logistic regression model for the selected data
    di "  with weight = 1/.30 for observations corresponding to negative EST.
    di

    qui gen     pvy1y2 =  1 if y1==1
    qui replace pvy1y2 = .3 if y1==0

    qui gen cid = _n
    qui reshape long y,i(cid) j(test)
    qui gen xtest = 2 - test

    qui qui gen pw = 1/pvy1y2

    glm   y xtest if d==1 & v==1 [pweight=pw],family(binomial) link(logit) cluster(cid) nolog
    di
    di "     ==> log odds ratio = " %4.2f _b[xtest]
    di

    di
    di " and without adjustment (weighting):
    di

    glm   y xtest if d==1 & v==1             ,family(binomial) link(logit) cluster(cid) nolog
    di
    di "     ==> log odds ratio = " %4.2f _b[xtest]
    di

    qui log close
end
**********************************************

doit