/* ex7_4.do

   example 7.4

   data required: psa_dre.dta
      Smith - Prostate Ca data: paired PSA & DRE cohort data

   last update:  24 May 2007
*/
version 9
set more off
cap semt_profile
cap log close
*************************************
cap pro drop doit
pro def doit
    log using ${semt_log}ex7_4, replace
    // use ${semt_data}psa_dre, clear
    use ${semt_data}psa_dre_v2, clear    // data update/correction 5/24/2007
    di
    di "**** example 7.4 ****"
    di
    di "            Table 7.1"
    table psa dre d, by(race)  missing  // col row scol

    tempvar touse
    gen byte `touse' = race==0   //  will restrict the following to observations for white males:
    qui count if `touse'
    local n = r(N)
    di
    di
    di "Detection probabilities (DP) and False referral probabilities (FP)"
    di "  among" %6.0f `n' " white men:"
    di
    assert  ~missing(psa) & ~missing(dre)
    assert ~missing(d) if psa | dre
    qui count if `touse' & d==1
    local nd = r(N)
    qui count if `touse' & d==0
    local ndb = r(N)

    tempname dp_dre fp_dre dp_psa fp_psa

    qui count if `touse' & dre & d==1
    local n_d_pos = r(N)
    qui count if `touse' & dre & d==0
    local n_db_pos = r(N)
    sca `dp_dre' = `n_d_pos'/`n'
    sca `fp_dre' = `n_db_pos'/`n'

    qui count if `touse' & psa & d==1
    local n_d_pos = r(N)
    qui count if `touse' & psa & d==0
    local n_db_pos = r(N)
    sca `dp_psa' = `n_d_pos'/`n'
    sca `fp_psa' = `n_db_pos'/`n'

    di "{txt}  DP(DRE) = {res}" %3.1f `dp_dre'*100 " %"
    di "{txt}  FP(DRE) = {res}" %3.1f `fp_dre'*100 " %"
    di
    di "{txt}  DP(PSA) = {res}" %3.1f `dp_psa'*100 " %"
    di "{txt}  FP(PSA) = {res}" %3.1f `fp_psa'*100 " %"
    di
    di
    di "{txt}Relative detection & false referral probabilities, PSA compared to DRE,
    di "{txt}   with 90% CI's:

    tempname rtpf rfpf ddp ffp varlog_rtpf varlog_rfpf tpf_lb tpf_ub fpf_lb fpf_ub z Z

    sca `rtpf' = `dp_psa'/`dp_dre'
    sca `rfpf' = `fp_psa'/`fp_dre'
    * CI's
    local level = .90
    sca `z' = invnorm(1 - (1-`level')/2)
    qui count if `touse' & psa & dre & d==1
    sca `ddp' = `r(N)'/`n'
    qui count if `touse' & psa & dre & d==0
    sca `ffp' = `r(N)'/`n'
    sca `varlog_rtpf' = (`dp_dre' + `dp_psa' - 2*`ddp')/(`n'*`dp_dre'*`dp_psa')
    sca `varlog_rfpf' = (`fp_dre' + `fp_psa' - 2*`ffp')/(`n'*`fp_dre'*`fp_psa')
    sca `tpf_lb' = exp( ln(`rtpf') - `z'*sqrt(`varlog_rtpf') )
    sca `tpf_ub' = exp( ln(`rtpf') + `z'*sqrt(`varlog_rtpf') )
    sca `fpf_lb' = exp( ln(`rfpf') - `z'*sqrt(`varlog_rfpf') )
    sca `fpf_ub' = exp( ln(`rfpf') + `z'*sqrt(`varlog_rfpf') )

    di
    di "{txt}  rDP(PSA,DRE) = rTPF = {res}" %4.2f `rtpf' " (" %4.3f `tpf_lb' ", " %4.3f `tpf_ub' ")"
    di "{txt}  rFP(PSA,DRE) = rFPF = {res}" %4.2f `rfpf' " (" %4.3f `fpf_lb' ", " %4.3f `fpf_ub' ")"
    di

    sca `Z' = ln(`rtpf')/sqrt(`varlog_rtpf')

    di
    di "{txt} test of Ho: rTPF(PSA,DRE) = 1
    di "
    di "{txt}   Z = ln(rTPF(PSA,DRE))/se = {res}" %4.2f `Z' ",{txt}   p = {res}" %4.0g 2*(1-normal(abs(`Z')))
    di

    di "{txt}  McNemar's test  (among those with ascertainment of disease) "
    ta dre psa if `touse' & d==1
    qui mcc dre psa if `touse' & d==1
    di
    di "{txt}   McNemar's chi2(1) = {res}" %6.0g `r(chi2)' "{txt},   p = {res}" %4.0g chi2tail(1,`r(chi2)')
    di

    di "{txt}**************"
    di "Covariate effects, eg. effect of race on DP & FP for DRE:
    di
    qui gen byte dpi = d & dre
    qui gen byte fpi = ~d & dre

    di "logit P[D=1,Y=1|X_R] = beta_0 + beta_1 * X_R: "
    logit dpi race, nolog
    di
    di
    di "logit P[D=0,Y=1|X_R] = theta_0 + theta_1 * X_R: "
    logit fpi race, nolog
    di
    di
    * table 7.4:
    qui gen byte detcat = dre + d
    qui replace detcat = 0 if dre==0
    la def detcat 0 "Y=0" 1 "Y=1,D=0" 2 "Y=1,D=1"
    la val detcat detcat
    la var detcat "Detection and False Referral"
    di
    di "{txt} Table 7.4"
    tab race detcat,r // nokey
    di

    * expand to 2 records per subject; 1 for each test type
    qui {
       gen patid=_n
       expand=2
       bys patid: gen byte tpsa=_n-1
       gen xtest=tpsa
       gen race_xtest = race*xtest

       gen      udr=(d==1)*(psa==1) if tpsa==1   /* the psa record */
       replace  udr=(d==1)*(dre==1) if tpsa==0   /* the dre record */

       gen      ufr=(d==0)*(psa==1) if tpsa==1   /* the psa record */
       replace  ufr=(d==0)*(dre==1) if tpsa==0   /* the dre record */
    }
    di
    di "{txt}**************"
    /* compare PSA versus DRE in whites */
    di " Comparison of PSA with DRE among white men: "
    di "   loglinear models "
    di
    di "log P[D=1,Y=1|Test] = beta_0 + beta_1 * X_T: "
    glm udr xtest if race==0, family(binomial) link(log) cluster(patid) level(90) nolog
    di
    di
    di "log P[D=0,Y=1|Test] = theta_0 + theta_1 * X_T: "
    glm ufr xtest if race==0, family(binomial) link(log) cluster(patid) level(90) nolog
    di
    di
    di "{txt}**************"
    di " Comparison of PSA with DRE, covariate effect (race): "
    di "   loglinear models "
    di
    di "log P[D=1,Y=1|Race,Test] = alpha_0 + alpha_1*X_R + beta_0*X_test + beta_1*X_R*X_T: "
    glm udr race xtest race_xtest, family(binomial) link(log) cluster(patid) level(90) nolog
    di
    di
    qui log close
end
*************************************
doit