/* fig5_6.do

   Figure 5.6,  example 5.5

   programs required:
      emroc.ado

   data required:
      tostbegg2.dta  (Liver metastasis data, Tosteson and Begg, 1988)

   last update:  30 Apr 2003

*/
version 7
cap log close
cap semt_profile
set more off
pause on

local alpha = .10   /* for 90% CI */

log using  ${semt_log}fig5_6, replace

use ${semt_data}tostbegg2
* primary cancer type: restrict to colon cancer
keep if type==0

tempname b V est1
preserve

/* plot empirical ROC curve with 90% confidence bands
   and save empirical roc estimates in file for
   later plotting */

emroc y d, conf level(90) gense(tpr) genspe(fpr)
replace fpr = 1-fpr
qui bys fpr tpr: drop if _n>1
order fpr tpr
sort fpr tpr
save ${semt_log}fig5_6empir,replace

pause q <return> to continue . . .

restore

* ml estimation of binormal ROC via Stata's -rocfit- :

rocfit d y

* plot fitted curve:

rocplot, conf level(90)
pause q <return> to continue . . .

matrix mat_b = e(b)
matrix mat_V = e(V)

matrix l mat_b
matrix l mat_V

matr Beta = mat_b[1,1..2]
matr V = mat_V[1..2,1..2]

matr list Beta
matr list V

matr C = [1, invnorm(.2)]

matr CB = (C*Beta')'
matr Var_CB = (C*V)*C'

matr list CB

matr list Var_CB

* conf coeff for 90% CI;
sca z = invnorm( 1.-(`alpha'/2.) )

sca roct = normprob(CB[1,1])
sca lb = normprob( CB[1,1] - z*sqrt(Var_CB[1,1]) )
sca ub = normprob( CB[1,1] + z*sqrt(Var_CB[1,1]) )

sca l

pause q <return> to continue . . .

* check the above matrix operations with scalar operations:

sca ahat = mat_b[1,1]
sca bhat = mat_b[1,2]
sca var_a = mat_V[1,1]
sca var_b = mat_V[2,2]
sca cov_ab = mat_V[2,1]

sca rocterm = ahat + invnorm(.2)*bhat
sca roct2 = normprob(rocterm)
sca se = sqrt( var_a + ((invnorm(.2))^2)*var_b + 2*invnorm(.2)*cov_ab)
sca lb2 = normprob(rocterm - z*se)
sca ub2 = normprob(rocterm + z*se)

sca l

/* generate points for the binormal ROC: */

* tempname a b
* sca `a' = ahat
* sca `b' =  bhat

* tempvar t roc forty5

qui {
   range t 0 1 201
   gen roc3 = norm(ahat + bhat*invnorm(t))
   replace roc3 = t if t==0 | t==1
   sort t
   gen forty5 = t if _n==1 | _n==_N
}

local bstr = string(round(bhat,.01))
local astr = string(round(ahat,.01))

#delimit ;
gr roc3 forty5 t,
  c(sl) s(ii)
  xscale(0,1) yscale(0,1)
  xlabel(0,.2,.4,.6,.8,1) ylabel(0,.2,.4,.6,.8,1)
  xtick(.1,.3,.5,.7,.9) ytick(.1,.3,.5,.7,.9)
  /* xline(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)
     yline(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1) */
  xline(0,1) yline(0,1) noaxis
  b1("FP rate") b2(" ")
  l1("TP rate")
  t1("binormal ROC; a = `astr', b = `bstr' ")
;
#delimit cr

keep roc3 t
sort t roc3
save ${semt_log}fig5_6,replace

qui log close