*! version 1.0.0 28feb2017 program heckroc, rclass version 12 syntax varlist (numeric min=2 max=2) [if] [in] [fweight], SELect(string) /// [ /// collinear Level(cilevel) table /// noci NOGraph NOEMPirical NOREFline cbands vce(string) Robust /// DIFficult TECHnique(string) TOLerance(real 1e-6) /// LTOLerance(real 1e-7) NRTOLerance(real 1e-5) /// * /// ] local maxopts `difficult' tech(`technique') tol(`tolerance') /// ltol(`ltolerance') nrtol(`nrtolerance') _get_gropts, graphopts(`options') getallowed(irocopts erocopts rlopts cbopts) local twowayopts `"`s(graphopts)'"' local rlopts `"`s(rlopts)'"' local irocopts `"`s(irocopts)'"' local erocopts `"`s(erocopts)'"' local cbopts `"`s(cbopts)'"' _check4gropts irocopts, opt(`irocopts') _check4gropts erocopts, opt(`erocopts') _check4gropts rlopts, opt(`rlopts') _check4gropts cbopts, opt(`cbopts') if "`cbands'"!="" & "`ci'"!="" { di as err "{p}option {bf:noci} may not be specified with {bf:cbands}{p_end}" exit 198 } /* Call preserve because we will need to create new variables later for integration*/ preserve /* Standardize */ gettoken y rhs : varlist qui sum `rhs' local mean =r(mean) local sd =r(sd) tempvar std_rhs gen `std_rhs' =(`rhs'-`mean')/`sd' /* Display table*/ if `"`table'"' != "" table `y' `rhs', row col /* Get the empirical AUC for later */ quietly roctab `y' `rhs' local emp_auc= r(area) /* Call heckprob and derive estimates from the output */ cap qui heckprob `y' `std_rhs' `if' `in' `fweight', select(`select') /// `collinear' `maxopts' vce(`vce') `robust' if _rc { if e(converged)==0 { di as err "{p}{bf:heckprobit} failed to converge{p_end}" exit 198 } else { heckprob `y' `std_rhs' `if' `in' `weight',select(`select') /// `collinear' `maxopts' vce(`vce') `robust' } } local rho_ap = sign([`y']_b[`std_rhs'])*([`y']_b[`std_rhs']^2/(1+ /// [`y']_b[`std_rhs']^2))^.5 local p_star = ([`y']_cons)*(1-`rho_ap'^2)^.5 /* Create confidence intervals for rho_ap and p_star */ tempvar heckprob_var matrix `heckprob_var'=e(V) local rho_ap_lb=sign([`y']_b[`std_rhs']-invnormal(1-`level'/200)* /// `heckprob_var'[1,1]^.5)*(([`y']_b[`std_rhs']-invnormal(0.975)* /// `heckprob_var'[1,1]^.5)^2/(1+([`y']_b[`std_rhs']-invnormal(0.975)* /// `heckprob_var'[1,1]^.5)^2))^.5 local rho_ap_ub=sign([`y']_b[`std_rhs']+invnormal(1-`level'/200)* /// `heckprob_var'[1,1]^.5)*(([`y']_b[`std_rhs']+invnormal(0.975)* /// `heckprob_var'[1,1]^.5)^2/(1+([`y']_b[`std_rhs']+invnormal(0.975)* /// `heckprob_var'[1,1]^.5)^2))^.5 if sign(`rho_ap_lb'*`rho_ap_ub'){ if abs(`rho_ap_lb') < abs(`rho_ap_ub') { local p_star_lb=([`y']_cons-invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-(`rho_ap_lb')^2)^.5 local p_star_ub=([`y']_cons+invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-(`rho_ap_ub')^2)^.5 } else { local p_star_lb=([`y']_cons-invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-(`rho_ap_ub')^2)^.5 local p_star_ub=([`y']_cons+invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-(`rho_ap_lb')^2)^.5 } } else{ if abs(`rho_ap_lb') < abs(`rho_ap_ub') { local p_star_lb=([`y']_cons-invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-0)^.5 local p_star_ub=([`y']_cons+invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-(`rho_ap_ub')^2)^.5 } else { local p_star_lb=([`y']_cons-invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-0)^.5 local p_star_ub=([`y']_cons+invnormal(1-`level'/200)* /// `heckprob_var'[2,2]^.5)*(1-(`rho_ap_lb')^2)^.5 } } /* End of creating CI for rho_ap and p_star */ di as text _n "Estimating inferred ROC curve..." /* I am creating the variables x_for_integ and y_for_integ to calcuate the integrals using integ. The values from integ (when 1000 points are used) are very similar to those obtained by using the user-written command "integrate". I am using min(max((...),0),1) out of an abundance of caution. */ tempvar TPR FPR tempvar x_for_integ y_for_integ qui range `x_for_integ' -4 10 1000 qui gen `y_for_integ'=normalden(`x_for_integ')*(1-normal((`p_star'-`rho_ap'*`x_for_integ')/ /// (1-`rho_ap'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `TPR'=min(max((r(integral)/(1-normal(`p_star'))),0),1) qui drop `x_for_integ' qui range `x_for_integ' -10 -4 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(normal((`p_star'-`rho_ap'*`x_for_integ')/ /// (1-`rho_ap'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `FPR'=min(max((r(integral)/normal(`p_star')),0),1) forval c = -3.9(.1)4 { qui drop `x_for_integ' qui range `x_for_integ' `c' 10 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(1-normal((`p_star'-`rho_ap'*`x_for_integ')/ /// (1-`rho_ap'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `TPR'=`TPR' \ min(max((r(integral)/(1-normal(`p_star'))),0),1) qui drop `x_for_integ' qui range `x_for_integ' -10 `c' 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(normal((`p_star'-`rho_ap'*`x_for_integ')/ /// (1-`rho_ap'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `FPR'=`FPR' \ min(max((r(integral)/normal(`p_star')),0),1) } svmat `TPR' svmat `FPR' qui replace `FPR' = 1 - `FPR' /* Calling rocreg to facilitate adding the empirical ROC curve to plot */ qui rocreg `y' `rhs', nobootstrap qui integ `TPR' `FPR' local AUC=r(integral) if "`ci'" == "" { if sign(`p_star_lb'*`p_star_ub'){ if abs(`p_star_lb') < abs(`p_star_ub') { local p_star_for_lb_ROC = abs(`p_star_lb') local p_star_for_ub_ROC = abs(`p_star_ub') } else { local p_star_for_lb_ROC = abs(`p_star_ub') local p_star_for_ub_ROC = abs(`p_star_lb') } } else { if abs(`p_star_lb') < abs(`p_star_ub') { local p_star_for_lb_ROC = 0 local p_star_for_ub_ROC = abs(`p_star_ub') } else { local p_star_for_lb_ROC = 0 local p_star_for_ub_ROC = abs(`p_star_lb') } } /* Creating lower bound for inferred ROC curve */ tempvar TPR_lb FPR_lb qui drop `x_for_integ' qui range `x_for_integ' -4 10 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(1-normal((`p_star_for_lb_ROC'- /// `rho_ap_lb'*`x_for_integ')/(1-`rho_ap_lb'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `TPR_lb'=min(max((r(integral)/(1-normal(`p_star_lb'))),0),1) qui drop `x_for_integ' qui range `x_for_integ' -10 -4 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(normal((`p_star_for_lb_ROC'- /// `rho_ap_lb'*`x_for_integ')/(1-`rho_ap_lb'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `FPR_lb'=min(max((r(integral)/normal(`p_star_for_lb_ROC')),0),1) forval c = -3.9(.1)4 { qui drop `x_for_integ' qui range `x_for_integ' `c' 10 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(1-normal((`p_star_for_lb_ROC'- /// `rho_ap_lb'*`x_for_integ')/(1-`rho_ap_lb'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `TPR_lb'=`TPR_lb' \ min(max((r(integral)/(1-normal(`p_star_for_lb_ROC'))),0),1) qui drop `x_for_integ' qui range `x_for_integ' -10 `c' 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(normal((`p_star_for_lb_ROC'- /// `rho_ap_lb'*`x_for_integ')/(1-`rho_ap_lb'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `FPR_lb'=`FPR_lb' \ min(max((r(integral)/normal(`p_star_for_lb_ROC')),0),1) } svmat `TPR_lb' svmat `FPR_lb' qui replace `FPR_lb' = 1 - `FPR_lb' /* Creating upper bound for inferred ROC curve*/ tempvar TPR_ub FPR_ub qui drop `x_for_integ' qui range `x_for_integ' -4 10 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(1-normal((`p_star_for_ub_ROC'- /// `rho_ap_ub'*`x_for_integ')/(1-`rho_ap_ub'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `TPR_ub'=min(max((r(integral)/(1-normal(`p_star_for_ub_ROC'))),0),1) qui drop `x_for_integ' qui range `x_for_integ' -10 -4 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(normal((`p_star_for_ub_ROC'- /// `rho_ap_ub'*`x_for_integ')/(1-`rho_ap_ub'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `FPR_ub'=min(max((r(integral)/normal(`p_star_for_ub_ROC')),0),1) forval c = -3.9(.1)4 { qui drop `x_for_integ' qui range `x_for_integ' `c' 10 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(1-normal((`p_star_for_ub_ROC'- /// `rho_ap_ub'*`x_for_integ')/(1-`rho_ap_ub'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `TPR_ub'=`TPR_ub' \ min(max((r(integral)/(1-normal(`p_star_for_ub_ROC'))),0),1) qui drop `x_for_integ' qui range `x_for_integ' -10 `c' 1000 qui replace `y_for_integ'=normalden(`x_for_integ')*(normal((`p_star_for_ub_ROC'- /// `rho_ap_ub'*`x_for_integ')/(1-`rho_ap_ub'^2)^.5)) qui integ `y_for_integ' `x_for_integ' matrix `FPR_ub'=`FPR_ub' \ min(max((r(integral)/normal(`p_star_for_ub_ROC')),0),1) } svmat `TPR_ub' svmat `FPR_ub' qui replace `FPR_ub' = 1 - `FPR_ub' /* Report estimates */ qui integ `TPR_lb' `FPR_lb' local AUC_lb=r(integral) qui integ `TPR_ub' `FPR_ub' local AUC_ub=r(integral) di as txt _n "{col 2}Empirical {col 15}Inferred {col 31} Inferred AUC" di as txt "{col 2}ROC area {col 15}ROC area {col 28} `level'% Conf. Interval" di as txt "{hline 50} di as res %8.04f `emp_auc' _col(14) %8.04f `AUC' _col(28) %8.04f /// `AUC_lb' _col(40) %8.04f `AUC_ub' di as txt "{hline 50} } else { di as txt _n "{col 2}Empirical {col 15}Inferred" di as txt "{col 2}ROC area {col 15}ROC area" di as txt "{hline 23} di as res %8.04f `emp_auc' _col(14) %8.04f `AUC' di as txt "{hline 23} } return local cmdline = "`0'" return scalar EmpAUC = `emp_auc' return scalar AUC = `AUC' if "`ci'"=="" { return scalar AUC_lb = `AUC_lb' return scalar AUC_ub = `AUC_ub' } /* Graphs */ if "`nograph'"=="" { local titles ytitle("Sensitivity") xtitle("1 - Specificity") local line1 (line `TPR' `FPR', `irocopts') if `"`norefline'"'=="" { local line2 (line `FPR' `FPR', `rlopts') local legend 1 "Inferred ROC" 2 "Reference" } else local legend 1 "Inferred ROC" if `"`noempirical'"'=="" { local scatter (scatter _roc_ _fpr_, connect(L) sort `erocopts') if `"`norefline'"'=="" local legend `legend' 3 "Empirical ROC" else local legend `legend' 2 "Empirical ROC" } if "`cbands'"!="" { local cbands (line `TPR_lb' `FPR_lb', lpattern(dash) `cbopts') /// (line `TPR_ub' `FPR_ub', lpatter(dash) `cbopts') if `"`norefline'"'=="" & `"`noempirical'"'=="" { local legend `legend' 4 "`level'% Lower bound" 5 /// "`level'% Upper bound" } else if (`"`norefline'"'=="" & `"`noempirical'"'!="") | /// (`"`norefline'"'!="" & `"`noempirical'"'=="") { local legend `legend' 3 "`level'% Lower bound" 4 /// "`level'% Upper bound" } else { local legend `legend' 2 "`level'% Lower bound" 3 /// "`level'% Upper bound" } } graph twoway `line1' `line2' `scatter' `cbands' , `twowayopts' /// legend(order(`legend')) `titles' || `addplot' } /* Call restore to remove variables created for integration */ restore end