/* ex3_7.do

   Example 3.7

   data required: est1.dta
       Exercise Stress test and Chest pain history data from the CASS study
       paired data

   last update:  06 October 2005

*/
version 9
set more off
cap semt_profile
cap log close
pause on

cap pro drop doit
program define doit
    log using ${semt_log}ex3_7, replace
    di
    di "{txt}                       ****  Example 3.7  ****
    di
    di "{txt} Regression modeling of predictive values"
    di
    di "{txt}       What is the incremental improvement in the PPV for coronary
    di "{txt}  artery disease when adding a positive exercise stress test (EST)
    di "{txt}  to positive Chest pain history (CPH)? "
    di
    di "{txt} *****************"
    qui use ${semt_data}est1,clear

    *  expand data to 2 records per subject for subjects with positive CPH test

    quietly {
        preserve
        qui gen int subject = _n
        expand = 2 if y2==1
        bys subj: gen byte x_l = _n==2
        gen xl_y1not = x_l*(1 - y1)
        gen xl_y1 = x_l*y1
    }

    di
    di "{txt} composite model to estimate P[D=1|+CPH] and P[D=1|+CPH,EST]:
    di

    * fitting algorithm hangs up in non-concave region w/o
    *    specifying -difficult- option

    glm d xl_y1not xl_y1 if y2==1, family(binomial) link(log) cluster(subject) difficult nolog

    tempname alpha0 alpha0_lb alpha0_ub beta1 beta1_lb beta1_ub beta2 beta2_lb beta2_ub z

    sca `z' = invnormal(1 - .10/2)

    * obtain model estimates and CI's :

    sca `alpha0' = _b[_cons]
    sca `alpha0_lb' = `alpha0' -  `z'*_se[_cons]
    sca `alpha0_ub' = `alpha0' +  `z'*_se[_cons]

    sca `beta1'  = _b[xl_y1not]
    sca `beta1_lb' = `beta1' -  `z'*_se[xl_y1not]
    sca `beta1_ub' = `beta1' +  `z'*_se[xl_y1not]

    sca `beta2'  = _b[xl_y1]
    sca `beta2_lb' = `beta2' -  `z'*_se[xl_y1]
    sca `beta2_ub' = `beta2' +  `z'*_se[xl_y1]

    * display exponentiated estimates and CI's:

    di
    di
    di "{txt} Estimates and 90% Confidence Intervals
    di
    di "{txt}                        PPV of + CPH alone (e^alpha_0): " as res %5.2f exp(`alpha0') ///
            "  (" %4.2f exp(`alpha0_lb') ", " %4.2f exp(`alpha0_ub') ")"
    di "{txt}  PPV increment, addition of + EST to + CPH (e^beta_2): " as res %5.2f exp(`beta2') ///
            "  (" %4.2f exp(`beta2_lb') ", " %4.2f exp(`beta2_ub') ")"
    di "{txt}  PPV decrement, addition of - EST to + CPH (e^beta_1): " as res %5.2f exp(`beta1') ///
            "  (" %4.2f exp(`beta1_lb') ", " %4.2f exp(`beta1_ub') ")"
    di
    qui log close
end

doit