/* ex3_4.do

   example 3.4

   data required: lplaudio_b.dta
      Leisenring neonatal audiology data

   last update:  27 July 2005

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

cap pro drop doit
pro def doit
    log using ${semt_log}ex3_4,replace

    di
    di "{txt} Example 3.4:    Marginal generalized linear models for TPF and FPF
    di "{txt}                       fit to the neonatal audiology data
    di

    use ${semt_data}lplaudio_b,clear

    /* consider only test A from this dataset: */
    keep if test==2
    keep d y sev age loc cid

    gen xage = age-35
    gen xsev = (sev-30)/10

    di "************"
    di
    di "{txt} model 1: log(TPF) = beta_1 + beta_2*X_age + beta_3*X_loc + beta_4*X_sev
    di
    glm y xage loc xsev if d==1, family(binomial) link(log) cluster(cid) nolog

    di
    di "{txt} exponentiated coefficient estimates: "
    glm, eform noheader

    di
    di "************"
    di
    di "{txt} model 2: log(FPF) = alpha_1 + alpha_2*X_age + alpha_3*X_loc
    glm y xage loc if d==0, family(binomial) link(log) cluster(cid) nolog

    di
    di "{txt}  exponentiated coefficient estimates:
    glm, eform noheader
    di

    * call the subprog to generate table 3.10 from these models, w/o age covariate:
    tab3_10
    qui log close
end
*************************************
cap pro drop tab3_10
prog def tab3_10
    * program to generate table 3.10

    preserve
    local alpha = .10   /* for 90% CI's */
    scalar z = invnormal(1.0 - `alpha'/2)

    * fit the log(TPF) and log(FPF) models w/o age covariates, and store results

    qui glm y loc xsev if d==1, family(binomial) link(log) cluster(cid)
    estimates store model_logtpf

    qui glm y loc if d==0, family(binomial) link(log) cluster(cid)
    estimates store model_logfpf

    qui {
        drop in 1/l
        * create a dataset with the location and severity values in table 3.10
        set obs 2
        replace loc = _n-1
        expand 3
        bys loc: replace sev = (_n+2)*10
        replace xsev = (sev-30)/10

        * retrieve predicted values and se's from the model fits:

        est for model_logtpf: predict logtpf, xb
        est for model_logtpf: predict se_logtpf, stdp
        est for model_logfpf: predict logfpf, xb
        est for model_logfpf: predict se_logfpf, stdp

        * generate fitted tpf and fpf values and 90% CI's:

        gen tpf = exp(logtpf)
        gen fpf = exp(logfpf)

        gen tpf_lb = exp(logtpf - z*se_logtpf)
        gen tpf_ub = exp(logtpf + z*se_logtpf)
        gen fpf_lb = exp(logfpf - z*se_logfpf)
        gen fpf_ub = exp(logfpf + z*se_logfpf)
    }

    di
    di "{txt} Table 3.10  Fitted TPF & FPF values (90% confidence intervals)"
    di "{hline 61}"
    di "{txt} Location  Severity           TPF                 FPF
    di "{hline 61}"
    forvalues i = 1/6 {
        local locval = loc[`i']
        local location "`:label (loc) `locval''"
        local tpf "`:di %4.2f `=tpf[`i']' " (" %4.2f `=tpf_lb[`i']' ", " %4.2f `=tpf_ub[`i']' ")" ' "
        local fpf "`:di %4.2f `=fpf[`i']' " (" %4.2f `=fpf_lb[`i']' ", " %4.2f `=fpf_ub[`i']' ")" ' "
        di "{txt}  `location'" _col(12) %4.0f sev[`i'] " dB  {res}   `tpf'  `fpf' "                                      "
    }
end
***********************************************
doit