cap program drop fagan
*! version 1.00 June 14, 2009
*! Ben A. Dwamena: bdwamena@umich.edu

program define fagan, rclass
version 10
syntax varlist(min=2 max=2) [if] [in], GRPvar(varname) [PRior(real 0.5) LEGENDopts(string asis) YSIZE(integer 6) XSIZE(integer 5) *]


tokenize `varlist'

local lrp `1'
local lrn `2'

foreach p in 0.1 0.2 0.3 0.5 0.7 1 2 3 5 7 10 ///
20 30 40 50 60 70 80 90 93 95 97 98 99  99.3 99.5 99.7 99.8 99.9 {
         local ylab `"`ylab' `=ln(`p' / (100 - `p'))' "`p'" "'
}
foreach lr in 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10 20  ///
        50 100 200 500 1000 {
        local lrpts `"`lrpts' `=-.5*ln(`lr')' 0 "`lr'" "'
}

tempvar id
gen `id'=_n
qui levelsof `id', local(slevels)
local l 5
local prprob = logit(1-`prior')
local note0: di  "PreProb:"%3.0f 100*`prior' "% "

foreach s of local slevels {
local id "`grpvar'[`s']"
local note: di `id'
local lrp`s' = `lrp'[`s']
local ppp`s' = logit(`prior') + log(`lrp'[`s'])
local note1: di "PostProb:"%5.1f 100*invlogit(`ppp`s'') "% "
local note2: di "LRP:"%5.2f `lrp`s'' "

local pppplot `"`pppplot' (pcarrowi `prprob' -1 `ppp`s'' 1, yaxis(2) lpatt(shortdash) barbsize(1) mlw(vvthin))"'
local legend1 `"`legend1' label(`l' "`note'(+)" "`note2'" "`note1'" )"'
local order "`order' `l++'"
}


foreach s of local slevels {
local id "`grpvar'[`s']"
local note: di `id'
local ppn`s' = logit(`prior') + log(`lrn'[`s'])
local lrn`s' = `lrn'[`s']
local note3: di "PostProb:"%5.1f 100*invlogit(`ppn`s'') "% "
local note4: di "LRN:"%5.2f `lrn`s'' "
local ppnplot `"`ppnplot' (pcarrowi `prprob' -1 `ppn`s'' 1, yaxis(2) lpatt(solid) barbsize(1) mlw(vvthin))"'
local legend2 `"`legend2' label(`l' "`note'(-)" "`note4'" "`note3'" )"'
local order "`order' `l++'"
}

if "`legendopts'" != "off" {
	if "`legendopts'" == "" { // default legend options
		local legendopts `"pos(2) size(*.550) symy(0.5) symx(5) cols(1) rowgap(2)"'
	}
	local legendopts `" order(4 "`note0'" `order') `legend1' `legend2' `legendopts'"'
}





#delimit;
tw (scatteri 0 0, mcolor(none) yaxis(1) ylab(`ylab', angle(0)
tpos(cross)) yscale(reverse axis(1)) ytitle("Prior Probability (%)", size(*.75) axis(1)))
(scatteri 0 0, mcolor(none) yaxis(2) ylab(`ylab', angle(0) tpos(cross) axis(2))
 ytitle("Posterior Probability (%)", size(*.75) axis(2)))
(scatteri `lrpts', msymbol(+) mcolor(black) mlabcolor(black) mlabsize(medsmall))
(pci -3.4538776 0 3.4538776 0, msymbol(none) recast(pcspike) lcolor(black)
xscale(range(-1 1)) plotregion(margin(zero)) xscale(off) ylab(, nogrid) text(-4 0 "Likelihood Ratio", place(n)))
`pppplot' `ppnplot', legend(`legendopts') subti(Fagan's nomogram) xsize(`xsize') ysize(`ysize') `options';
#delimit cr

end