/* ex3_8.do

   Example 3.8

   programs required
       binscrn1.ado  (v 1.1.2 or later)
       lrreg.ado, lrreg_ll.ado (v 1.0.2 or later)

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

   last update:  10 October 2005

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

cap pro drop doit
program define doit
    log using ${semt_log}ex3_8, replace
    di
    di "{txt}                       ****  Example 3.8  ****
    di
    di "{txt} Regression models for DLR's - comparing DLR's of two tests
    di
    di "{txt} *****************"
    quietly use ${semt_data}est1,clear
    preserve

    tempname z dlrp_lb dlrp_ub
    local alpha = .10  // (for 90% CI)
    sca `z' = invnormal(1 - `alpha'/2)

    * use returned results from program -binscrn1- to calculate CI's for DLR+

    * a) for CPH
    quietly binscrn1 y2 d
    scalar `dlrp_lb' = exp( ln(r(dlrpos)) - `z'*r(selndlrp) )
    scalar `dlrp_ub' = exp( ln(r(dlrpos)) + `z'*r(selndlrp) )
    di
    di "{txt} Empirical DLR+'s and 90% CI's
    di
    di "{txt}    DLR+(CPH): " as res %6.2f `r(dlrpos)' ///
                 "  (" %4.2f `dlrp_lb' ", "  %4.2f `dlrp_ub' ")"
    * b) for EST
    quietly binscrn1 y1 d
    scalar `dlrp_lb' = exp( ln(r(dlrpos)) - `z'*r(selndlrp) )
    scalar `dlrp_ub' = exp( ln(r(dlrpos)) + `z'*r(selndlrp) )
    di "{txt}    DLR+(EST): " as res %6.2f `r(dlrpos)' ///
                 "  (" %4.2f `dlrp_lb' ", "  %4.2f `dlrp_ub' ")"
    di
    di "{txt} *****************"
    di
    di "{txt} Regression modeling of DLR's as a function of test type
    di "{txt}      (CPH relative to EST)
    di

    * use program -lrreg- to do ml regression modeling of DLR's

    * first, expand to 2 records per subject, 1 for each test type

    quietly {
        gen int subj = _n
        expand 2
        sort subj
        by subj: gen int xtest = (_n==2)
        by subj: gen int    y = y1 if _n==1
        by subj: replace    y = y2 if _n==2
    }

    eq lrpos: xtest
    eq lrneg: xtest
    lrreg y d lrpos lrneg, cluster(subj)

    di
    di
    di "{txt}  rDLR+ (CPH relative to EST)
    lincom [lrpos]xtest, rr

    di
    di
    di "{txt}  rDLR- (CPH relative to EST)
    lincom [lrneg]xtest, rr

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

doit