/* ex9_1.do

  example 9.1 & figure 9.1
     meta-analysis: the summary ROC curve

  data required: fim.dta
     data from 59 studies of the Pap smear screen for cervical cancer

  last update:  14 Aug 2003

*/
version 7
set more off
cap semt_profile
cap log close
clear

log using ${semt_log}ex9_1, replace

* Example 9.1 :
*    meta-analysis: summary ROC curve for 59 studies of the
*                   Pap smear screen for cervical cancer

use ${semt_data}fim

***************************************************************
/*  data "correction":
     record 47 in in irwig.dta and in
     Fahey, Irwig, and Macaskill '95 AJE reference appendix
     is changed here from tpr = .098 to .698
     we have assumed this to be an entry/typographical error
     the change is consistent with the published figure 1 of this reference
     and the estimates in table 2       */

replace tpr = .698 if rec==47

gen d = ln(tpr/(1 - tpr)) - ln(fpr/(1 - fpr))
gen s = ln(tpr/(1 - tpr)) + ln(fpr/(1 - fpr))

regress d s

sca a = _b[_cons]
sca b = _b[s]

gen double roct = a/(1 - b) + ((b + 1)/(1 - b)) * ln(fpr/(1 - fpr))
replace    roct = exp(roct)/(1 + exp(roct))

************
* weighted least squares:

gen double ww = 1/(nd * tpr * (1 - tpr)) + 1/(ndb * fpr * (1 - fpr))
replace    ww = 1/ww

regress d s [aweight=ww]

sca aw = _b[_cons]
sca bw = _b[s]

gen double roctw = aw/(1 - bw) + ((bw + 1)/(1 - bw)) * ln(fpr/(1 - fpr))
replace    roctw = exp(roctw)/(1 + exp(roctw))
qui log close
*****************************

gen byte exrec = _n==1
expand 3 if exrec
gen forty5 = .
local var1lst "roct roctw fpr forty5"

foreach var of local var1lst {
    bys exrec: replace `var' = 0 if _n==1 & exrec==1
    bys exrec: replace `var' = 1 if _n==2 & exrec==1
}
sort fpr
#delimit ;

gr roct roctw tpr forty5 fpr,
    con(ss.l) s(iio.) noaxis xline(0(.2)1) yline(0(.2)1)
    xlab(0(.2)1) ylab(0(.2)1) l1title(TPF) b2t(FPF)
    ;

* save plot coordinate data for import to graphics pkg:
keep roct roctw tpr fpr
order roct roctw tpr fpr
sort fpr
save ${semt_log}fig9_1,replace