cap program drop midas
*! version 1.00 August 15, 2007
*! Ben A. Dwamena: bdwamena@umich.edu

program define midas, rclass byable(recall) sortpreserve
version 10
 #delimit;
 syntax varlist(min=4) [if] [in] , [ ID(varname) YEAR(varname)  EStimator(string)
 NIP(integer 15) LEVEL(integer 95) SCHEME (string)
 QUALITab QUALIBar QLAB   RESults(string) TABle(string)
 BIVBOX  CHIPlot GALB(string) CUM INF FUNnel PUBBias MAXbias QQplot(string)
 FORest(string) FORData  HETfor MScale(real 0.45) TEXTScale(real 0.50)
 ROCPlane SROC1 SROC2
 FAGAN PDDAM(string) PRIOR(real 1.0) LRMatrix
 PLOTtype TESTlab(string asis) HSIZE(integer 6) VSIZE(integer 8)
 COVars ZCF(real 0.5)  *];
#delimit cr


qui {
preserve
marksample touse, novarlist
keep if `touse'
}

tokenize `varlist'
global tp `1'
global fp `2'
global fn `3'
global tn `4'
macro shift 4
global varlist2 `*'


/* Check syntax */

if `level' < 10 | `level' > 99 {
di as error "level() must be between 10 and 99"
	exit 198
}


if ( "`pddam'" != "" | "`fagan'" != "" | "`forest'" != "" ///
| "`rocplane'" != "" | "`sroc1'" != "" | "`sroc2'" != "" ///
| "`hetfor'" != ""|"`results'" != "" | "`table'" != "" ///
| "`lrmatrix'" != "")  & ("`estimator'" == "" ) {
  	di as error "estimator() must be specified with this option"
	exit 198
}

if ("`qlab'" != "") & ("`qualitab'" == "" | "`qualibar'" == "") {
 di as error "qlab and quality varlist must be used with qualibar or qualitab"
	exit 198

}


/*Data Management*/

qui {
global alph = (100-`level')/200
local numobs = _N
if ("`id'" != "" & "`year'" != "") {
local id `id'
local year `year'
format `id' %30s/* right aligned string */
egen StudyIds = concat(`id' `year'), p(" ")
}
else gen str30 StudyIds = string(_n)

if "`scheme'" != ""{
 set scheme `scheme'
  }
  if "`scheme'" == ""{
  set scheme s2color
  }

if `"`testlab'"' == "" {
local tlab  ""
}
else if `"`testlab'"' != "" {
local tlab `"`testlab'"'
}

nois di ""
nois di ""
nois di ""

nois di as text "{title:META-ANALYTIC INTEGRATION OF DIAGNOSTIC TEST ACCURACY STUDIES}"
/* QUALITY ASSESSMENT */
if "`qualibar'" == "qualibar" {
if "`qlab'" == "qlab" {
noisily quadas $varlist2,  labvar(`qlab') qgraph
}
else if "`qlab'" == "" {
noisily quadas $varlist2, qgraph
}
}

else  if "`qualitab'" == "qualitab" {
if "`qlab'" == "qlab" {
noisily quadas $varlist2, labvar(`qlab') qtable
}
else if "`qlab'" == "" {
noisily quadas $varlist2, qtable
}
}
version 9
/* CALCULATE TOTALS */
tempvar sum sumtp sumfn sumtn sumfp sumtpfn sumtnfp sumsu
egen `sumtp' = total($tp)
egen `sumfn' = total($fn)
egen `sumtn' = total($tn)
egen `sumfp' = total($fp)
gen `sumtpfn' = `sumtp' + `sumfn'
gen `sumtnfp' = `sumtn' + `sumfp'
egen `sum' = rsum($tp $fn $tn $fp)
global prev = `sumtpfn'/(`sumtnfp' + `sumtpfn')
gen `sumsu' = sum(`sum')


/* Study Specific Adjustment for Zeros */
tempvar zc_tp zc_fn zc_fp zc_tn zero zc_sens zc_fpr zc_spec zc_tpr zc_fnr zc_tot
gen `zc_tp' = $tp
gen `zc_fp' = $fp
gen `zc_fn' = $fn
gen `zc_tn' = $tn
gen `zero' = 0
replace `zero' = 1 if $tp == 0 | $fp == 0 | $fn == 0 | $tn == 0
replace `zc_tp' = `zc_tp' + `zcf' if `zero' == 1
replace `zc_fp' = `zc_fp' + `zcf' if `zero' == 1
replace `zc_fn' = `zc_fn' + `zcf' if `zero' == 1
replace `zc_tn' = `zc_tn' + `zcf' if `zero' == 1
gen `zc_sens' = `zc_tp'/(`zc_tp'+`zc_fn')    /* adjusted sensitivity */
gen `zc_tpr' = `zc_sens'             /* adjusted true pos rate */
gen `zc_fnr' = `zc_fn'/(`zc_tp'+`zc_fn')     /* adjusted false neg rate */
gen `zc_spec' = `zc_tn'/(`zc_tn'+`zc_fp')    /* adjusted specificity */
gen `zc_fpr' = `zc_fp'/(`zc_tn'+`zc_fp')     /* adjusted false pos rate */
gen `zc_tot' = `zc_tp'+`zc_fp'+`zc_fn'+`zc_tn'   /* adjusted total */


/* STUDY-SPECIFIC Sensitivity (True Positive Rate)*/
tempvar sens senslo senshi sensse spec speclo spechi specse FPR
gen `sens' = $tp/($tp+$fn)
gen `senslo' = invbinomial($tp+$fn,$tp,$alph)
gen `senshi' = invbinomial($tp+$fn,$tp,1-$alph)
gen `sensse' = (`senshi'-`senslo')/(2*invnorm(1-$alph))

/* STUDY-SPECIFIC Specificity (True Negative Rate) */

gen `spec' = $tn/($tn+$fp)
gen `speclo' = invbinomial($tn+$fp,$tn,$alph)
gen `spechi' = invbinomial($tn+$fp,$tn,1-$alph)
gen `specse' =(`spechi'-`speclo')/(2*invnorm(1-$alph))
gen `FPR' = 1 - `spec'


/* Study Specific Positive Likelihood Ratio And Confidence Interval */
tempvar lrp llrp llrpvar llrpse lrplo lrphi lrpse
gen `lrp' = `zc_sens'/`zc_fpr'
gen `llrp' = ln(`zc_sens'/`zc_fpr')
gen `llrpvar' = (1/`zc_tp')+(1/`zc_fp')-(1/(`zc_tp'+`zc_fn'))-(1/(`zc_fp'+`zc_tn'))
gen `llrpse' = sqrt((1/`zc_tp')+(1/`zc_fp')-(1/(`zc_tp'+`zc_fn'))-(1/(`zc_fp'+`zc_tn')))
gen `lrplo' = exp(`llrp' - invnorm(1-$alph)*`llrpse')
gen `lrphi' = exp(`llrp' + invnorm(1-$alph)*`llrpse')
gen `lrpse' = (`lrphi'-`lrplo')/(2*invnorm(1-$alph))

/* Study Specific Negative Likelihood Ratio And Confidence Interval */
tempvar lrn llrn llrnvar llrnse lrnlo lrnhi lrnse

gen `lrn' = `zc_fnr'/`zc_spec'
gen `llrn' = ln(`zc_fnr'/`zc_spec')
gen `llrnvar' = (1/`zc_fn')+(1/`zc_tn')-(1/(`zc_tp'+`zc_fn'))-(1/(`zc_fp'+`zc_tn'))
gen `llrnse' = sqrt((1/`zc_fn')+(1/`zc_tn')-(1/(`zc_tp'+`zc_fn'))-(1/(`zc_fp'+`zc_tn')))
gen `lrnlo' = exp(`llrn' - invnorm(1-$alph)*`llrnse')
gen `lrnhi' = exp(`llrn' + invnorm(1-$alph)*`llrnse')
gen `lrnse' = (`lrnhi'-`lrnlo')/(2*invnorm(1-$alph))

/* Study Specific Diagnostic Odds Ratio And Confidence Interval */
tempvar dor dorvar dorse dorlo dorhi ldor ldorvar ldorse ldorlo ldorhi
tempname ecf
scalar `ecf' = sqrt(3)/_pi
gen `dor' = (`zc_tp'*`zc_tn')/(`zc_fp'*`zc_fn')
gen `ldor' = ln(`dor')
gen `dorvar' = (1/`zc_fn')+(1/`zc_tn')+(1/`zc_fp')+(1/`zc_tp')
gen `ldorvar' = (1/`zc_fn')+(1/`zc_tn')+(1/`zc_fp')+(1/`zc_tp')
gen `ldorse' = sqrt(`ldorvar')
gen `ldorlo' = `ldor'-invnorm(1-$alph)*`ldorse'
gen `ldorhi' = `ldor'+invnorm(1-$alph)*`ldorse'
gen `dorlo' = exp(`ldor'-invnorm(1-$alph) * `ldorse')
gen `dorhi' = exp(`ldor'+invnorm(1-$alph) * `ldorse')
gen `dorse' = (`dorhi'-`dorlo')/(2*invnorm(1-$alph))
replace `ldorse' = `ldorse' * `ecf'
replace `ldorlo' = `ldorlo' * `ecf'
replace `ldorhi' = `ldorhi' * `ecf'


/* Study Specific Logit Transform of Sensitivity (TPR) and CI */
tempvar lsens lsensvar lsensse lsenslo lsenshi
gen `lsens' = logit(`zc_sens')
gen `lsensvar' = 1/(`zc_sens'*(1-`zc_sens')*(`zc_tp'+`zc_fn'))
gen `lsensse' = sqrt(`lsensvar')
gen `lsenslo' = `lsens' - invnormal(1-$alph) * `lsensse'
gen `lsenshi' = `lsens' + invnormal(1-$alph) * `lsensse'

/* Study Specific Logit Transform of Specificity and CI */
tempvar lspec lspecvar lspecse lspeclo lspechi
gen `lspec' = logit(`zc_spec')
gen `lspecvar' = 1/(`zc_spec'*(1-`zc_spec')*(`zc_tn'+`zc_fp'))
gen `lspecse' = sqrt(`lspecvar')
gen `lspeclo' = `lspec' - invnormal(1-$alph) * `lspecse'
gen `lspechi' = `lspec' + invnormal(1-$alph) * `lspecse'


/* Study Specific Logit Transform of 1 - Specificity (FPR) and CI */
tempvar lfpr lfprvar lfprsd lfprlo lfprhi
gen `lfpr' = logit(`zc_fpr')
gen `lfprvar' = 1/(`zc_fpr'*(1-`zc_fpr')*(`zc_tn'+`zc_fp'))
gen `lfprsd' = sqrt(`lfprvar')
gen `lfprlo' = `lfpr' - invnormal(1-$alph) * `lfprsd'
gen `lfprhi' = `lfpr' + invnormal(1-$alph) * `lfprsd'
nois di ""


/*GALBRAITH PLOT FOR INVESTIGATING HETEROGENEITY AND SMALL STUDY BIAS*/


if "`galb'" != ""{
if "`plottype'" != "" {
local plottype "Galbraith Plot"
}
else if "`plottype'" == "" {
local plottype " "
}

if "`galb'" == "ldor" {
nois midagalb `ldor' `ldorse'
}

else if "`galb'" == "lrp" {
nois midagalb `llrp' `llrpse'
}

else if "`galb'" == "lrn" {
nois midagalb `llrn' `llrnse'
}
else if "`galb'" == "tpr" {
nois midagalb `lsens' `lsensse'
}
else if "`galb'" == "tnr" {
nois midagalb `lspec' `lspecse'
}
}


/* BIVARIATE BOX PLOT */
if "`bivbox'"=="bivbox" {
tempvar boxvar1 boxvar2
gen `boxvar1' = `lsens'
label var `boxvar1' "LOGIT_SENS"
gen `boxvar2' = `lspec'
label var `boxvar2' "LOGIT_SPEC"
nois di as text "BIVARIATE BOXPLOT"
nois bvbox `boxvar1' `boxvar2'
}


/* CHI PLOT */

if "`chiplot'" == "chiplot" {
tempvar cvar1 cvar2
gen `cvar1' = `lsens'
label var `cvar1' "LOGIT_SENS"
gen `cvar2' = `lspec'
label var `cvar2' "LOGIT_SPEC"
nois di as text "CHIPLOT"

nois midachi `cvar1' `cvar2'
}

/* CUMULATIVE META-ANALYSIS*/
	if "`cum'" == "cum" {
	tempvar var1 var1lo var1se var1hi var2 var2se var2lo var2hi
     gen `var1' = `sens'
     gen `var1lo' = `senslo'
     gen `var1se' = `sensse'
     gen `var1hi' = `senshi'
     gen `var2' = `spec'
     gen `var2lo' = `speclo'
     gen `var2hi' = `spechi'
     gen `var2se' = `specse'


	tempvar cumvar1 cumvar1se cumvar2 cumvar2se
     local obs=_N
     gen `cumvar1'=`var1' in 1
     gen `cumvar1se' = `var1se' in 1
     gen `cumvar2' = `var2' in 1
     gen `cumvar2se' = `var2se' in 1

	local i 2

     while `i'<=`obs' {
     metan `var1' `var1lo' `var1hi' in 1/`i', random nograph
     replace `cumvar1' = r(ES) in `i'
     replace `cumvar1se' = r(seES) in `i'
     metan `var2' `var2lo' `var2hi' in 1/`i', random nograph
     replace `cumvar2' = r(ES) in `i'
     replace `cumvar2se' = r(seES) in `i'
     local i=`i'+1
      }

 	tempname  obs1 stvar1 stvar1se stvar2 stvar2 stvar2se
	gen `stvar1' = `cumvar1'
	gen `stvar1se' = `cumvar1se'
	gen `stvar2' = `cumvar2'
	gen `stvar2se' = `cumvar2se'
	gen `obs1' = _n
	local xmax = `obs1'[_N]

	set graphics off

	#delimit;
	serrbar `stvar1' `stvar1se' `obs1', scale(2) ylab(, format(%9.2f) angle(hor) labs(*.50))
	xtitle("StudyNo", size(*.75)) ytitle("Sensitivity", size(*0.75)) name(cumplot1, replace)
	xlab(1(1)`xmax',  labs(*.50) format(%9.0f) angle(hor));
	#delimit cr

	#delimit;
	serrbar `stvar2' `stvar2se' `obs1', scale(2) ylab(, format(%9.2f) angle(hor) labs(*.50))
	xtitle("StudyNo", size(*.75)) ytitle("Specificity", size(*0.75)) name(cumplot2, replace)
	xlab(1(1)`xmax', labs(*.50) format(%9.0f) angle(hor)) ;
	#delimit cr

	set graphics on

	nois graph combine cumplot1 cumplot2, cols(1) subtitle("Cumulative Analysis By Publication Year", size(*0.75))
	}


/* INFLUENCE ANALYSIS */
   	if  "`inf'" =="inf" {
   	tempvar var1 var1lo var1hi var2 var2lo var2hi
     gen `var1' = `sens'
     gen `var1lo' = `senslo'
     gen `var1hi' = `senshi'
     gen `var2' = `spec'
     gen `var2lo' = `speclo'
     gen `var2hi' = `spechi'

   	tempvar so
   	qui gen `so' = _n
   	sort `so'

* Meta-analysis estimate omiting one study each step

   tempvar istvar1 istvar1se istvar2 istvar2se
      local n = _N
      gen `istvar1' = .
      gen `istvar1se' = .
      gen `istvar2' = .
      gen `istvar2se' = .

    local i = 1
    tempvar s
    gen `s' = _n
    while (`i' <= `n') {
    sort `so'
    metan `var1' `var1lo' `var1hi' if `s' != `i', random nograph
     replace `istvar1' = r(ES) in `i'
     replace `istvar1se' = r(seES) in `i'

     metan `var2' `var2lo' `var2hi' if `s' != `i', random nograph
     replace `istvar2' = r(ES) in `i'
     replace `istvar2se' = r(seES) in `i'

     local i=`i'+1
   }


	tempvar obs1 studyvar1 studyvar2 studyvar1se studyvar2se

	gen `studyvar2' = `istvar2'
	gen `studyvar2se' = `istvar2se'
	gen `studyvar1' = `istvar1'
	gen `studyvar1se' = `istvar1se'

	tempvar  obs1
	gen `obs1' = _n
	local xmax = `obs1'[_N]
	set graphics off

	#delimit;
	serrbar `studyvar1' `studyvar1se' `obs1', scale(2) ylab(, format(%9.2f) angle(hor) labs(*.75))
	xsc(r(1 `xmax')) xlab(1(1)`xmax', labs(*.50) format(%9.0f) angle(hor))
	ytitle("Sensitivity", size(*0.75)) xtitle("StudyNo", size(*.75)) name(infplot1, replace);
	#delimit cr


	#delimit;
	serrbar `studyvar2' `studyvar2se' `obs1', scale(2) ylab(, format(%9.2f) angle(hor) labs(*.75))
	xsc(r(1 `xmax')) xlab(1(1)`xmax', labs(*.50) format(%9.0f) angle(hor))
	ytitle("Specificity", size(*0.75)) xtitle("StudyNo", size(*.75)) name(infplot2, replace);
	#delimit cr

	set graphics on

	nois graph combine infplot1 infplot2, cols(1) title("Influence Analysis", size(*0.75))

	}


/* PUBLICATION BIAS */

 if "`pubbias'" != "" | "`maxbias'" != "" | "`funnel'" != "" {
local ptitle: di "Log Odds Ratio versus 1/sqrt(Effective Sample Size)(Deeks)"
tempvar  n1 n2 nt ESS zero lthetai thetai sethetai xb yb stpred logn pubwgt
gen `n1' = `zc_tp' + `zc_fn'
gen `n2' = `zc_tn' + `zc_fp'
gen `nt' =  `n1' + `n2'
gen `ESS' =(4 * `n1' * `n2')/(`n1' + `n2')
gen `lthetai' = log((`zc_tp' * `zc_tn')/(`zc_fp' * `zc_fn'))
gen `thetai' =.
gen `sethetai' =.
gen `pubwgt' =.
gen `xb' =.
gen `yb' =.

if "`pubbias'" == "pubbias"  {
nois di " "
nois di " "
nois di as text "{title:STATISTICAL TESTS FOR SMALL STUDY EFFECTS/PUBLICATION BIAS}"
replace `thetai'=`lthetai'
replace `sethetai' = sqrt(`ESS')
replace `xb'=1/sqrt(`ESS')
label var `xb' "1/root(ESS)"
replace `yb' = `thetai'
label var `yb' "Diagnostic Odds Ratio"
replace `pubwgt' = `ESS'
nois di " "
nois di " "
nois di " "
sum `xb', detail
local xbmax=r(max)
sum `yb', detail
local ybmax=r(max)
local ybmin=r(min)
mylabels 1 10 100 1000, myscale(log(@)) local(ylab)
regress  `yb' `xb'[weight=`pubwgt'], level(`level')
estimates store stbias
scalar intercept = _b[_cons]
scalar se_intercept = _se[_cons]
scalar rmse = e(rmse)
scalar df = e(df_r)
scalar p = 2*ttail(e(df_r), abs(return(score_bc)/return(score_se)))
nois matrix define vcov = e(V)
nois matrix define b = e(b)
matrix colnames b =     Bias Intercept
matrix rownames vcov =  Bias Intercept
matrix colnames vcov =  Bias Intercept
nois matrix post b vcov, dep(yb) dof(`e(df_r)') obs(`e(N)')
nois ereturn display, level(`level')
estimates restore stbias
predict `stpred'
local stline "(line `stpred' `xb' ,  clpattern(dash) clwidth(vthin))"
#delimit;
nois twoway `stline'(scatter `yb' `xb' , sort ms(O) `options'),
legend(label(2 "Study") label(1 "Regression" "Line") order(2 1) pos(2) col(1) size(*.75))
 xlab(, angle(horizontal) labsize(*`textscale') format(%7.2f))
ylab( `ylab', labsize(*`textscale')) title("`plottype' `tlab'", size(*0.75))
 xsize(`hsize') plotregion(margin(zero)) ytitle("Diagnostic Odds Ratio")
subtitle(`ptitle', size(*.65)) ;
#delimit cr
}

if "`maxbias'" == "maxbias" {
nois di " "
nois di " "
nois di as text "{title: WORST CASE SENSITIVITY ANALYSIS FOR SMALL STUDY EFFECTS/PUBLICATION BIAS}"

replace `sethetai' = sqrt(`ESS')
nois copas `lthetai' `sethetai'
}

if "`funnel'" == "funnel" {

replace `thetai'=`lthetai'
replace `sethetai' = 1/sqrt(`ESS')
replace `pubwgt' = `ESS'
label var `thetai' "Log Odds Ratio "
nois di " "
nois di " "
nois di " "
if "`plottype'" != "" {
local plottype "Funnel Plot"
}
else if "`plottype'" == "" {
local plottype " "
}

mylabels 1 10 100 1000, myscale(log(@)) local(xlab)
twoway (scatter `sethetai' `thetai', legend(off) ms(O)), yscale(rev) ylab(,  angle(horizontal) format(%7.2f) labsize(*`textscale')) title("`plottype' `tlab'", size(*.75)) /*
*/ ytitle("Precision", size(*.75)) xtitle("Log Odds Ratio", size(*.75)) plotregion(margin(zero)) xlab(`xlab', format(%2.0f) labsize(*`textscale'))

  	drop if missing(`sethetai')
   }
}

/*NORMAL QUANTILE PLOT FOR INVESTIGATING NORMALITY ASSUMPTION, HETEROGENEITY AND PUBLICATION SELECTION BIAS*/
if "`qqplot'" != ""{
if "`plottype'" != "" {
local plottype "Normal Quantile Plot"
}
else if "`plottype'" == "" {
local plottype " "
}

if "`qqplot'" == "ldor" {
nois metanorm `ldor' `ldorse'
}

else if "`qqplot'" == "lrp" {
nois metanorm `llrp' `llrpse'
}

else if "`qqplot'" == "lrn" {
nois metanorm `llrn' `llrnse'
}
else if "`qqplot'" == "tpr" {
nois metanorm `lsens' `lsensse'
}
else if "`qqplot'" == "tnr" {
nois metanorm `lspec' `lspecse'
}
}


if ("`qualitab'`qualibar'`qqplot'`pubbias'`bivbox'`chiplot'" == "") ///
& ("`galb'`maxbias'`funnel'`cum'`inf'" == "") {


/*  MODEL SPECIFICATION AND ESTIMATION   */
if "`estimator'" == "g" {
bbrre tp fp fn tn
}
else if "`estimator'" == "x" {
version 10
xtbbrre tp fp fn tn, `nip'
}

*saving transformed estimates
return scalar mtpr = r(mtpr)
return scalar mtprlo = r(mtprlo)
return scalar mtprhi = r(mtprhi)
return scalar mtprse = r(mtprse)

return scalar mtnr = r(mtnr)
return scalar mtnrlo = r(mtnrlo)
return scalar mtnrhi = r(mtnrhi)
return scalar mtnrse = r(mtnrse)

return scalar mldor =  r(mldor)
return scalar mldorlo = r(mldorlo)
return scalar mldorhi = r(mldorhi)
return scalar mldorse = r(mldorse)

return scalar mdor = r(mdor)
return scalar mdorlo = r(mdorlo)
return scalar mdorhi = r(mdorhi)
return scalar mdorse = r(mdorse)

return scalar mlrp = r(mlrp)
return scalar mlrplo = r(mlrplo)
return scalar mlrphi = r(mlrphi)
return scalar mlrpse = r(mlrpse)

return scalar mlrn = r(mlrn)
return scalar mlrnlo = r(mlrnlo)
return scalar mlrnhi = r(mlrnhi)
return scalar mlrnse = r(mlrnse)

return scalar  reffs1 = r(mreffs1)
return scalar  reffs1lo = r(mreffs1lo)
return scalar  reffs1hi = r(mreffs1hi)
return scalar  reffs1se = r(mreffs1se)


return scalar  reffs2 = r(mreffs2)
return scalar  reffs2lo = r(mreffs2lo)
return scalar  reffs2hi = r(mreffs2hi)
return scalar  reffs2se = r(mreffs2se)


return scalar  rho = r(mrho)
return scalar  rholo = r(mrholo)
return scalar  rhohi = r(mrhohi)
return scalar  covar = r(mcovar)

return scalar fsens = r(fsens)
return scalar fspec = r(fspec)
return scalar fldor = r(fldor)
return scalar fdor =  r(fdor)
return scalar flrp =  r(flrp)
return scalar flrn =  r(flrn)

return scalar Islrt = r(Islrt)
return scalar Islrtlo = r(Islrtlo)
return scalar Islrthi = r(Islrthi)



*SUMMARY ESTIMATES
local cov01 = r(covsnsp)
local mtpr = r(mtpr)
local mtprlo = r(mtprlo)
local mtprhi = r(mtprhi)
local mtprse = r(mtprse)

local mtnr = r(mtnr)
local mtnrlo = r(mtnrlo)
local mtnrhi = r(mtnrhi)
local mtnrse = r(mtnrse)

local mldor =  r(mldor)
local mldorlo = r(mldorlo)
local mldorhi = r(mldorhi)
local mldorse = r(mldorse)

local mdor = r(mdor)
local mdorlo = r(mdorlo)
local mdorhi = r(mdorhi)
local mdorse = r(mdorse)

local mlrp = r(mlrp)
local mlrplo = r(mlrplo)
local mlrphi = r(mlrphi)
local mlrpse = r(mlrpse)

local mlrn = r(mlrn)
local mlrnlo = r(mlrnlo)
local mlrnhi = r(mlrnhi)
local mlrnse = r(mlrnse)

local reffs1 = r(mreffs1)
local reffs1lo = r(mreffs1lo)
local reffs1hi = r(mreffs1hi)
local reffs1se = r(mreffs1se)


local reffs2 = r(mreffs2)
local reffs2lo = r(mreffs2lo)
local reffs2hi = r(mreffs2hi)
local reffs2se = r(mreffs2se)


local rho = r(mrho)
local rholo = r(mrholo)
local rhohi = r(mrhohi)
local covar = r(mcovar)

global fsens = r(fsens)
global fspec = r(fspec)
global fldor = r(fldor)
global fdor =  r(fdor)
global flrp =  r(flrp)
global flrn =  r(flrn)

tempname sp sn spse snse lrtchi lrtpchi lrtdf
tempname Islrt Islrtlo Islrthi
scalar `sp' = r(sp)
scalar `spse' = r(spse)
scalar `sn' = r(sn)
scalar `snse' = r(snse)
scalar `lrtchi'  = r(lrtchi)
scalar `lrtpchi' = r(lrtpchi)
scalar `lrtdf' = r(lrtdf)
scalar `Islrt' = r(Islrt)
scalar `Islrtlo' = r(Islrtlo)
scalar `Islrthi' = r(Islrthi)


*Multi-column Forest Plot

if "`hetfor'" == "hetfor" {
set graphics off
qui {
tempvar obs wgt1 wgt2 wgt3 wgt4 wgt5
gen `obs' = _n
local mscale2 = .25 * `mscale'
count
local max1 = r(N)
label value `obs' obs
forval i = 1/`max1'{
local value = `"`value' `i'"'
label define obs `i' "`=StudyIds[`i']'", modify
}
gen line=_n
gen textx=0
local n1=_N+2
local n2=_N+1
local ylabopt "labsize(*.50) tl(*0) labgap(*5)"

local xlab1 "xlab(minmax, format(%5.1f) labsize(*.50))"
local xlab2 "xlab(minmax, format(%5.1f) labsize(*.50))xsc(log)"

gen `wgt1' = 1/(`sensse' *`sensse')
#delimit ;
twoway (rspike `senslo' `senshi' `obs', ylabel(`"`value'"', valuelabel labsize(*.50) tl(*0) angle(360))
hor s(i) lpat(blank)  `xlab1')(scatter `obs' `sens' [aw = `wgt1'], ms(i) msize(*`mscale2') mcolor(gs10)),
legend(off) xtitle("", size(*.5)) yscale(noline) xscale(off fill) plotregion(style(none)) ytitle("", size(*.5))
title("StudyId", size(*.5) pos(1)) fxsize(0) name(mplot, replace);
#delimit cr

#delimit ;
twoway (rspike `senslo' `senshi' `obs', ylabel(`"`value'"', nolabel
`ylabopt' angle(360)) hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab1')
(scatter `obs' `sens' [aw = `wgt1'], ms(O) msize(*`mscale2') mcolor(gs10))
(scatter `obs' `sens', ms(o) msize(*`mscale') mcolor(black)), xline(`mtpr', lpattern(dash))
ytitle("", size(*.5)) legend(off) xtitle("", size(*.5)) title("Sensitivity", size(*.5)
justification(left)) name(mplot1, replace) ;
#delimit cr

gen `wgt2' = 1/(`specse'*`specse')

#delimit ;
twoway (rspike `speclo' `spechi' `obs', ylabel(`"`value'"', nolabel
`ylabopt' angle(360)) hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab1')
(scatter `obs' `spec' [aw = `wgt2'], ms(O) msize(*`mscale2') mcolor(gs10))
(scatter `obs' `spec', ms(o) msize(*`mscale') mcolor(black)), legend(off)
xtitle("", size(*.5)) ytitle("", size(*.5))  title("Specificity", size(*.5)
justification(left)) nodraw xline(`mtnr', lpattern(dash)) name(mplot2, replace) ;
#delimit cr

gen `wgt3' = 1/(`lrpse' *`lrpse')

#delimit ;
twoway (rspike `lrplo' `lrphi' `obs', ylabel(`"`value'"', nolabel `ylabopt' angle(360))
hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab2')(scatter `obs' `lrp' [aw = `wgt3'],
ms(O) msize(*`mscale2') mcolor(gs10))(scatter `obs' `lrp', ms(o) msize(*`mscale') mcolor(black)),
nodraw legend(off) title("Positive Likelihood Ratio", size(*.5)
justification(left)) xline(`mlrp', lpattern(dash)) ytitle("", size(*.5))
xtitle("", size(*.5)) name(mplot3, replace) ;
#delimit cr

gen `wgt4' = 1/(`lrnse' *`lrnse')


#delimit ;
twoway (rspike `lrnlo' `lrnhi' `obs', ylabel(`"`value'"', nolabel `ylabopt' angle(360))
hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab1')(scatter `obs' `lrn' [aw = `wgt4'],
ms(O) msize(*`mscale2') mcolor(gs10))(scatter `obs' `lrn', ms(o) msize(*`mscale') mcolor(black))
, nodraw legend(off) xtitle("", size(*.5)) ytitle("", size(*.5))
title("Negative Likelihood Ratio", size(*.5) justification(left)) xline(`mlrn', lpattern(dash))
name(mplot4, replace);
#delimit cr

gen `wgt5' = 1/(`dorse' *`dorse')

#delimit ;
twoway (rspike `dorlo' `dorhi' `obs', ylabel(`"`value'"', nolabel `ylabopt' angle(360))
hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab2')(scatter `obs' `dor' [aw = `wgt5'],
ms(O) msize(*`mscale2') mcolor(gs10))(scatter `obs' `dor', ms(o) msize(*`mscale') mcolor(black)),
nodraw legend(off) title("Diagnostic Odds Ratio", size(*.5) justification(left))
xtitle("", size(*.5)) ytitle("", size(*.5)) xline(`mdor', lpattern(dash)) name(mplot5, replace) ;
#delimit cr
}
set graphics on
nois graph combine mplot mplot1 mplot2 mplot3 mplot4 mplot5, row(1)

}


/* HETEROGENEITY: Sensitivity */

tempvar devsens Qsens
tempname  Qsensdf  prsens Isqsens Isqsenslo Isqsenshi
gen `devsens' = ((`sens' - $fsens)^2)/(($fsens*(1-$fsens))/($tp+$fn))
egen `Qsens' = total(`devsens')
scalar `Qsensdf' = `numobs' - 1
scalar `prsens' = chi2tail(`Qsensdf',`Qsens')
homogeni `Qsens' `Qsensdf'
scalar `Isqsens' = r(Isq)
scalar `Isqsenslo' = r(Isqlo)
scalar `Isqsenshi' = r(Isqhi)


/* HETEROGENEITY: Specificity  */

tempvar devspec Qspec
tempname  Qspecdf  prspec Isqspec Isqspeclo Isqspechi
gen `devspec' = ((`spec' - $fspec)^2)/(($fspec*(1-$fspec))/($tn+$fp))
egen `Qspec' = total(`devspec')
scalar `Qspecdf' = `numobs' - 1

scalar `prspec' = chi2tail(`Qspecdf',`Qspec')
homogeni `Qspec' `Qspecdf'
scalar `Isqspec' = r(Isq)
scalar `Isqspeclo' = r(Isqlo)
scalar `Isqspechi' = r(Isqhi)


/* HETEROGENEITY: Positive Likelihood Ratio */

tempvar Qlrp devlrp
tempname prlrp Qlrpdf Isqlrp Isqlrplo Isqlrphi
gen `devlrp' = ((`llrp' - ln($flrp))^2)/`llrpvar'
egen `Qlrp' = total(`devlrp')
scalar `Qlrpdf' = `numobs' - 1
scalar `prlrp' = chi2tail(`Qlrpdf',`Qlrp')
homogeni `Qlrp' `Qlrpdf'
scalar `Isqlrp' = r(Isqlo)
scalar `Isqlrplo' = r(Isqlo)
scalar `Isqlrphi' = r(Isqhi)


/* HETEROGENEITY:  Negative Likelihood Ratio */

tempvar Qlrn devlrn
tempname prlrn Qlrndf Isqlrn Isqlrnlo Isqlrnhi

gen `devlrn' = ((`llrn' - ln($flrn))^2)/`llrnvar'
egen `Qlrn' = total(`devlrn')
scalar `Qlrndf' = `numobs' - 1
scalar `prlrn' = chi2tail(`Qlrndf',`Qlrn')
homogeni `Qlrn' `Qlrndf'
scalar `Isqlrn' = r(Isq)
scalar `Isqlrnlo' = r(Isqlo)
scalar `Isqlrnhi' = r(Isqhi)


/*  HETEROGENEITY: Diagnostic Odds Ratio */

tempvar Qldor devldor
tempname prldor Qldordf Isqldor Isqldorlo Isqldorhi
gen `devldor' = (`ldor' - $fldor)^2/`ldorvar'
egen `Qldor' = total(`devldor')
scalar `Qldordf' = `numobs' - 1
scalar `prldor' = chi2tail(`Qldordf',`Qldor')

homogeni `Qldor' `Qldordf'
scalar `Isqldor' = r(Isq)
scalar `Isqldorlo' = r(Isqlo)
scalar `Isqldorhi' = r(Isqhi)


tempvar Qdor devdor
tempname prdor Qdordf Isqdor Isqdorlo Isqdorhi
gen `devdor' = exp(((`ldor' - ln($fdor))^2)/`ldorvar')
egen `Qdor' = total(`devdor')
scalar `Qdordf' = `numobs' - 1
scalar `prdor' = chi2tail(`Qdordf',`Qdor')

homogeni `Qdor' `Qdordf'
scalar `Isqdor' = r(Isq)
scalar `Isqdorlo' = r(Isqlo)
scalar `Isqdorhi' = r(Isqhi)

tempvar var1 var2 var1se var2se var1lo var2lo var1hi var2hi
tempname Isq1 Isq1lo Isq1hi Isq2 Isq2lo Isq2hi
tempname Qvar1 Qvar1df pr1 Qvar2  Qvar2df pr2

     gen `var1' = .
     gen `var1se' = .
     gen `var1lo' = .
     gen `var1hi' = .
     gen `var2' = .
     gen `var2se' = .
     gen `var2lo' = .
     gen `var2hi' = .

if ("`forest'" == "dss" | "`table'" == "dss" | "`rocplane'" == "rocplane" )  {

     replace `var1' = `sens'
     replace `var1se' = `sensse'
     replace `var1lo' = `senslo'
     replace `var1hi' = `senshi'
     replace `var2' = `spec'
     replace `var2se' = `specse'
     replace `var2lo' = `speclo'
     replace `var2hi' = `spechi'

     local mvar1 = `mtpr'
     local mvar1lo = `mtprlo'
     local mvar1hi = `mtprhi'
     local mvar2 = `mtnr'
     local mvar2lo = `mtnrlo'
     local mvar2hi = `mtnrhi'

     scalar `Isq1' = `Isqsens'
     scalar `Isq1lo' = `Isqsenslo'
     scalar `Isq1hi' = `Isqsenshi'
     scalar `Isq2' = `Isqspec'
     scalar `Isq2lo' = `Isqspeclo'
     scalar `Isq2hi' = `Isqspechi'
	scalar `Qvar1' = `Qsens'
	scalar `Qvar1df' = `Qsensdf'
	scalar `pr1' = `prsens'
	scalar `Qvar2' = `Qspec'
	scalar `Qvar2df' = `Qspecdf'
	scalar `pr2' = `prspec'

local note1a: di " "%4.2f `mvar1' " [" %4.2f `mvar1lo' " - " %4.2f `mvar1hi' "]"
local note1b: di "Q ="%6.2f `Qvar1' ", df = " %3.2f `Qvar1df' ", p = "%5.2f `pr1' "
local note1c: di "I2 = "%3.2f `Isq1' " [" %3.2f `Isq1lo' " - " %3.2f `Isq1hi' "]"
local note2a: di " "%4.2f `mvar2' " [" %4.2f `mvar2lo' " - " %4.2f `mvar2hi' "]"
local note2b: di "Q ="%6.2f `Qvar2' ", df = " %3.2f `Qvar2df' ", p = "%5.2f `pr2' "
local note2c: di "I2 = "%3.2f `Isq2' " [" %3.2f `Isq2lo' " - " %3.2f `Isq2hi' "]"



     local gtitle1 "SENSITIVITY"
     local gtitle2 "SPECIFICITY"
     }


if ("`forest'" == "dlr" | "`table'" == "dlr" ) {

     replace `var1' = `lrp'
     replace `var1se' = `lrpse'
     replace `var1lo' = `lrplo'
     replace `var1hi' = `lrphi'
     replace `var2' = `lrn'
     replace `var2se' = `lrnse'
     replace `var2lo' = `lrnlo'
     replace `var2hi' = `lrnhi'

     local mvar1 = `mlrp'
     local mvar1lo = `mlrplo'
     local mvar1hi = `mlrphi'
     local mvar2 = `mlrn'
     local mvar2lo = `mlrnlo'
     local mvar2hi = `mlrnhi'

     scalar `Isq1' = `Isqlrp'
     scalar `Isq1lo' = `Isqlrplo'
     scalar `Isq1hi' = `Isqlrphi'
     scalar `Isq2' = `Isqlrn'
     scalar `Isq2lo' = `Isqlrnlo'
     scalar `Isq2hi' = `Isqlrnhi'

          scalar `Qvar1' = `Qlrp'
	scalar `Qvar1df' = `Qlrpdf'
	scalar `pr1' = `prlrp'
	scalar `Qvar2' = `Qlrn'
	scalar `Qvar2df' = `Qlrndf'
	scalar `pr2' = `prlrn'

	local note1a: di " "%4.2f `mvar1' " [" %4.2f `mvar1lo' " - " %4.2f `mvar1hi' "]"
local note1b: di "Q ="%6.2f `Qvar1' ", df = " %3.2f `Qvar1df' ", p = "%5.2f `pr1' "
local note1c: di "I2 = "%3.2f `Isq1' " [" %3.2f `Isq1lo' " - " %3.2f `Isq1hi' "]"
local note2a: di " "%4.2f `mvar2' " [" %4.2f `mvar2lo' " - " %4.2f `mvar2hi' "]"
local note2b: di "Q ="%6.2f `Qvar2' ", df = " %3.2f `Qvar2df' ", p = "%5.2f `pr2' "
local note2c: di "I2 = "%3.2f `Isq2' " [" %3.2f `Isq2lo' " - " %3.2f `Isq2hi' "]"


     local gtitle1 "DLR POSITIVE"
     local gtitle2 "DLR NEGATIVE"
     }

if ("`forest'" == "dlor" | "`table'" == "dlor") {

     replace `var1' = `ldor'
     replace `var1se' = `ldorse'
     replace `var1lo' = `ldorlo'
     replace `var1hi' = `ldor'
     replace `var2' = `dor'
     replace `var2se' = `dorse'
     replace `var2lo' = `dorlo'
     replace `var2hi' = `dorhi'

     local mvar1 = `mldor'
     local mvar1lo = `mldorlo'
     local mvar1hi = `mldorhi'
     local mvar2 = `mdor'
     local mvar2lo = `mdorlo'
     local mvar2hi = `mdorhi'

     scalar `Isq1' = `Isqldor'
     scalar `Isq1lo' = `Isqldorlo'
     scalar `Isq1hi' = `Isqldorhi'
     scalar `Isq2' = `Isqdor'
     scalar `Isq2lo' = `Isqdorlo'
     scalar `Isq2hi' = `Isqdorhi'
     scalar `Qvar1' = `Qldor'
     scalar `Qvar1df' = `Qldordf'
     scalar `pr1' = `prldor'
     scalar `Qvar2' = `Qdor'
     scalar `Qvar2df' = `Qdordf'
     scalar `pr2' = `prdor'

local note1a: di " "%4.2f `mvar1' " [" %4.2f `mvar1lo' " - " %4.2f `mvar1hi' "]"
local note1b: di "Q ="%6.2f `Qvar1' ", df = " %3.2f `Qvar1df' ", p = "%5.2f `pr1' "
local note1c: di "I2 = "%3.2f `Isq1' " [" %3.2f `Isq1lo' " - " %3.2f `Isq1hi' "]"
local note2a: di " "%4.2f `mvar2' " [" %4.2f `mvar2lo' " - " %4.2f `mvar2hi' "]"
local note2b: di "Q ="%6.2f `Qvar2' ", df = " %3.2f `Qvar2df' ", p = "%5.2f `pr2' "
local note2c: di "I2 = "%3.2f `Isq2' " [" %3.2f `Isq2lo' " - " %3.2f `Isq2hi' "]"

     local gtitle1 "DIAGNOSTIC SCORE"
     local gtitle2 "ODDS RATIO"
     }



 if "`table'" != "" {
nois di ""
nois di ""
nois di as text "{title: STUDY-SPECIFIC TEST PERFORMANCE ESTIMATES}"
nois di " "
nois di " "
sum `var1', detail
local n1=r(N)
nois di as text "`gtitle1'"
nois di " "
nois di as text"{hline 65}"
nois di in gr _col(2) "Study" _col(20) "{c |}" _col(24) "Estimate" _col(39) "[95%  Conf.  Interval]"
nois di as text"{hline 19}{c +}{hline 47}"
local i = 1
while `i' <= `n1' {
local a1 = StudyIds in `i'
local b1 = `var1' in `i'
local c1 = `var1lo'  in `i'
local d1 = `var1hi'  in `i'
nois di in gr _col(2) %6.2f "`a1'" _col(20) in gr "{c |}" in ye _col(24) %6.2f `b1' _col(39) %6.2f `c1' _col(49) %6.2f `d1'
local i=`i'+1
}
nois di as text"{hline 65}"
nois di in gr _col(2) "Combined" _col(20) in gr "{c |}" in ye _col(24) %6.2f `mvar1' _col(39) %6.2f `mvar1lo' _col(49) %6.2f `mvar1hi'
nois di as text"{hline 19}{c BT}{hline 47}"
nois di ""
nois di as txt "Heterogeneity (Chi-square): Q = "as result %5.2f `Qvar1' ///
as txt ", df = "as result %3.2f `Qvar1df' as txt", p = "as result %5.2f `pr1'
nois di " "
nois di as txt"Inconsistency (I-square): I2 = "as res %3.2f `Isq1' _c
nois di as txt", 95% CI = ["as res %3.2f `Isq1lo'                _c
nois di as txt" - "as res %3.2f `Isq1hi' as txt"]"               _n
nois di " "
nois di ""
nois di as text "`gtitle2'"
nois di " "
nois di as text"{hline 65}"
nois di in gr _col(2) "Study" _col(20) "{c |}" _col(24) "Estimate" _col(39) "[95%  Conf.  Interval]"
nois di as text"{hline 19}{c +}{hline 47}"
local i = 1
while `i' <= `n1' {
local a2 = StudyIds in `i'
local b2 = `var2' in `i'
local c2 = `var2lo'  in `i'
local d2 = `var2hi'  in `i'
nois di in gr _col(2) %6.2f "`a2'" _col(20) in gr "{c |}" in ye _col(24) %6.2f `b2' _col(39) %6.2f `c2' _col(49) %6.2f `d2'
local i=`i'+1
}
nois di as text"{hline 65}"
nois di in gr _col(2) "Combined" _col(20) in gr "{c |}" in ye _col(24) %6.2f `mvar2' _col(39) %6.2f `mvar2lo' _col(49) %6.2f `mvar2hi'
nois di as text"{hline 19}{c BT}{hline 47}"
nois di ""
nois di as txt "Heterogeneity (Chi-square): Q = " as result %5.2f `Qvar2' ///
as txt ", df =" as result %3.2f `Qvar2df' as txt", p ="as result %5.2f `pr2'
nois di " "
nois di as txt "Inconsistency (I-square): I2 = "as res %3.2f `Isq2' _c
nois di as txt ", 95% CI = ["as res %3.2f `Isq2lo'                _c
nois di as txt " - "as res %3.2f `Isq2hi' as txt"]"               _n
nois di " "
nois di ""
}


/*ROC PLANE*/
if "`rocplane'" == "rocplane" {
local msens = `mvar1'
local msenslo = `mvar1lo'
local msenshi= `mvar1hi'
local mspec = 1-`mvar2'
local mspeclo = 1-`mvar2hi'
local mspechi = 1-`mvar2lo'
nois twoway (pci `msens' 0 `msens' 1, lpat(longdash) lwidth(vthin)) /*
*/ (pci  0 `mspec' 1 `mspec' , lpat(shortdash) lwidth(vthin))(scatter `sens' `FPR', /*
*/ sort msymbol(O) mcolor(black)), ytitle("Sensitivity", size(*.90)) /*
*/ xtitle("1-Specificity", size(*.90)) yscale(range(0 1)) /*
*/ yline(`msenslo' `msenshi', lpat(longdash) lwidth(vthin))/*
*/ xline(`mspeclo' `mspechi', lpat(shortdash) lwidth(medthin)) ylabel( 0(.2)1, nogrid angle(horizontal) format(%7.2f)) /*
*/ xscale(range(0 1)) xlabel(0(.2)1, nogrid format(%7.2f)) title(ROC Plane, size(*.5)) /*
*/ legend(order(1 "Sensitivity" "`note1a'" "`note1b'" "`note1c'" 2 "Specificity" "`note2a'" "`note2b'" "`note2c'") /*
*/ pos(2) symxsize(1) forcesize rowgap(1) col(1) size(*.90)) name(ROCplot, replace)
}



/*COMBINED FOREST PLOTS*/
tempname obs obs1 obs2 studyvar1 studyvar2 studyvar1lo studyvar1hi studyvar2lo studyvar2hi
gen `obs' = _n
gen `obs1' = _n
gen `obs2' = _n
local null1: di " "
count
local max = r(N)
local maxx = `max' + 2
label value `obs' obs
forval i = 1/`max'{
local value = `"`value' `i'"'
label define obs `i' "`=StudyIds[`i']'", modify
}

gen `studyvar2' = .
gen `studyvar2lo' = .
gen `studyvar2hi' = .
gen `studyvar1' = .
gen `studyvar1lo' = .
gen `studyvar1hi' = .

local ylabopt "labsize(*`textscale') tl(*0) labgap(*5)"

if "`forest'" == "dss"{
local notef1a: di "Q ="%6.2f `Qvar1' ", df = " %3.2f `Qvar1df' ", p = "%5.2f `pr1' "
local notef1b: di "I2 = "%3.2f `Isq1' " [" %3.2f `Isq1lo' " - " %3.2f `Isq1hi' "]"
local notef2a: di "Q ="%6.2f `Qvar2' ", df = " %3.2f `Qvar2df' ", p = "%5.2f `pr2' "
local notef2b: di "I2 = "%3.2f `Isq2' " [" %3.2f `Isq2lo' " - " %3.2f `Isq2hi' "]"
replace `studyvar2' = `var2'
replace `studyvar2lo' = `var2lo'
replace `studyvar2hi' = `var2hi'
replace `studyvar1' = `var1'
replace `studyvar1lo' = `var1lo'
replace `studyvar1hi' = `var1hi'
local xlab1 "xlab(minmax, format(%5.1f) labsize(*`textscale'))"
local xlab2 "xlab(minmax, format(%5.1f) labsize(*`textscale'))"
tostring `studyvar1lo' `studyvar1' `studyvar1hi', gen(`studyvar1lo'1 `studyvar1'1 `studyvar1hi'1) format(%7.2f) force
replace `studyvar1lo'1 = " [" + `studyvar1lo'1 + " - "
replace `studyvar1hi'1 = `studyvar1hi'1 + "]"
egen studyvar1ci = concat(`studyvar1'1 `studyvar1lo'1 `studyvar1hi'1)
label value `obs1' obs1
forval i = 1/`max'{
local value1 = `"`value' `i'"'
label define obs1 `i' "`=studyvar1ci[`i']'", modify
}
tostring `studyvar2lo' `studyvar2' `studyvar2hi', gen(`studyvar2lo'1 `studyvar2'1 `studyvar2hi'1) format(%7.2f) force
replace `studyvar2lo'1 = " [" + `studyvar2lo'1 + " - "
replace `studyvar2hi'1= `studyvar2hi'1 + "]"
egen studyvar2ci = concat(`studyvar2'1 `studyvar2lo'1 `studyvar2hi'1)
label value `obs2' obs2
forval i = 1/`max'{
local value2 = `"`value' `i'"'
label define obs2 `i' "`=studyvar2ci[`i']'", modify
}
}

else if "`forest'" == "dlr" {
local notef1a: di "Q ="%6.2f `Qvar1' ", df = " %3.2f `Qvar1df' ", p = "%5.2f `pr1' "
local notef1b: di "I2 = "%3.2f `Isq1' " [" %3.2f `Isq1lo' " - " %3.2f `Isq1hi' "]"
local notef2a: di "Q ="%6.2f `Qvar2' ", df = " %3.2f `Qvar2df' ", p = "%5.2f `pr2' "
local notef2b: di "I2 = "%3.2f `Isq2' " [" %3.2f `Isq2lo' " - " %3.2f `Isq2hi' "]"
replace `studyvar2' = `var2'
replace `studyvar2lo' = max(0.01, `var2lo')
replace `studyvar2hi' = min(1.00, `var2hi')
replace `studyvar1' = `var1'
replace `studyvar1lo' = max(0.01, `var1lo')
replace `studyvar1hi' = min(1000, `var1hi')
local xlab1 "xlab(minmax, format(%5.1f) labsize(*`textscale'))xsc(log)"
local xlab2 "xlab(minmax, format(%5.0f) labsize(*`textscale'))"
tostring `studyvar1lo' `studyvar1' `studyvar1hi', gen(`studyvar1lo'1 `studyvar1'1 `studyvar1hi'1) format(%7.2f) force
replace `studyvar1lo'1 = " [" + `studyvar1lo'1 + " - "
replace `studyvar1hi'1 = `studyvar1hi'1 + "]"
egen studyvar1ci = concat(`studyvar1'1 `studyvar1lo'1 `studyvar1hi'1)
label value `obs1' obs1
forval i = 1/`max'{
local value1 = `"`value' `i'"'
label define obs1 `i' "`=studyvar1ci[`i']'", modify
}
tostring `studyvar2lo' `studyvar2' `studyvar2hi', gen(`studyvar2lo'1 `studyvar2'1 `studyvar2hi'1) format(%7.2f) force
replace `studyvar2lo'1 = " [" + `studyvar2lo'1 + " - "
replace `studyvar2hi'1= `studyvar2hi'1 + "]"
egen studyvar2ci = concat(`studyvar2'1 `studyvar2lo'1 `studyvar2hi'1)
label value `obs2' obs2
forval i = 1/`max'{
local value2 = `"`value' `i'"'
label define obs2 `i' "`=studyvar2ci[`i']'", modify

}
}

else if "`forest'" == "dlor" {
local notef1a: di "Q ="%6.2f `Qvar1' ", df = " %3.2f `Qvar1df' ", p = "%5.2f `pr1' "
local notef1b: di "I2 = "%3.2f `Isq1' " [" %3.2f `Isq1lo' " - " %3.2f `Isq1hi' "]"
local notef2a: di "Q ="%6.2f `Qvar2' ", df = " %3.2f `Qvar2df' ", p = "%5.2f `pr2' "
local notef2b: di "I2 = "%3.2f `Isq2' " [" %3.2f `Isq2lo' " - " %3.2f `Isq2hi' "]"
replace `studyvar2' = `var2'
replace `studyvar2lo' = max(0.01, `var2lo')
replace `studyvar2hi' = min(1000, `var2hi')
replace `studyvar1' = `var1'
replace `studyvar1lo' = max(0.01, `var1lo')
replace `studyvar1hi' = `var1hi'
local xlab1 "xlab(minmax, format(%5.1f) labsize(*`textscale')) "
local xlab2 "xlab(minmax, format(%5.0f) labsize(*`textscale')) xsc(log)"
tostring `studyvar1lo' `studyvar1' `studyvar1hi', gen(`studyvar1lo'1 `studyvar1'1 `studyvar1hi'1) format(%7.2f) force
replace `studyvar1lo'1 = " [" + `studyvar1lo'1 + " - "
replace `studyvar1hi'1 = `studyvar1hi'1 + "]"
egen studyvar1ci = concat(`studyvar1'1 `studyvar1lo'1 `studyvar1hi'1)
label value `obs1' obs1
forval i = 1/`max'{
local value1 = `"`value' `i'"'
label define obs1 `i' "`=studyvar1ci[`i']'", modify
}
tostring `studyvar2lo' `studyvar2' `studyvar2hi', gen(`studyvar2lo'1 `studyvar2'1 `studyvar2hi'1) format(%7.2f) force
replace `studyvar2lo'1 = " [" + `studyvar2lo'1 + " - "
replace `studyvar2hi'1= `studyvar2hi'1 + "]"
egen studyvar2ci = concat(`studyvar2'1 `studyvar2lo'1 `studyvar2hi'1)
label value `obs2' obs2
forval i = 1/`max'{
local value2 = `"`value' `i'"'
label define obs2 `i' "`=studyvar2ci[`i']'", modify
}
}

if "`forest'" != "" {
if "`plottype'" != "" {
local plottype "Forest Plot"
}
else if "`plottype'" == "" {
local plottype " "
}

local null  " "

local note1f: di " "%4.2f `mvar1'  "[" %4.2f `mvar1lo' " - " %4.2f `mvar1hi' "]"
local note2f: di " "%4.2f `mvar2'  "[" %4.2f `mvar2lo' " - " %4.2f `mvar2hi' "]"

if ("`fordata'" == "") {
nois twoway (pcarrow `obs' `studyvar1hi' `obs' `studyvar1lo'  if (`studyvar1lo' == 0.01 & "`forest'" == "dlor"),  lwidth(vthin) lpat(solid)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar1lo' `obs' `studyvar1hi'  if (`studyvar1hi' == 1000 & "`forest'" == "dlr"),  lpat(solid) lwidth(vthin)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar1hi' `obs' `studyvar1lo'  if (`studyvar1lo' == 0.01 & "`forest'" == "dlr"),  lwidth(vthin) lpat(solid)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/(rspike `studyvar1lo' `studyvar1hi' `obs', ylabel(`maxx' "STUDY(YEAR)" -1 "COMBINED" -2 "`notef1a'" -3 "`notef1b'" `"`value'"', valuelabel `ylabopt' angle(360)) /*
*/ hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab1')(scatter `obs' `studyvar1', ms(S) msize(*`mscale2') mcolor(gs10))(scatter `obs' `studyvar1', ms(o) msize(*`mscale') mcolor(black) xline(`mvar1', lpattern(-)))/*
*/ (scatteri -1 `mvar1lo' -0.8 `mvar1', clcolor(black) c(l) s(i))(scatteri -1 `mvar1lo' -1.2 `mvar1', clcolor(black) c(l) s(i)) (scatteri -0.8 `mvar1' -1 /*
*/ `mvar1hi', clcolor(black) c(l) s(i))(scatteri -1.2 `mvar1' -1 `mvar1hi', clcolor(black) c(l) s(i)),  /*
*/ legend(off)  xtitle("`gtitle1'", size(*.5)) title("`plottype'" "`tlab'", size(*0.75)) name(forplot1, replace)


nois twoway (pcarrow `obs' `studyvar2hi' `obs' `studyvar2lo'  if (`studyvar2lo' == 0.01 & "`forest'" == "dlor"),  lpat(solid) lwidth(vthin)   barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar2lo' `obs' `studyvar2hi'  if (`studyvar2hi' == 1000 & "`forest'" == "dlor"),  lwidth(vthin) lpat(solid)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar2hi' `obs' `studyvar2lo'  if (`studyvar2lo' == 0.01 & "`forest'" == "dlr"),  lpat(solid) lwidth(vthin)   barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (rspike `studyvar2lo' `studyvar2hi' `obs', ylabel(`maxx' "STUDY(YEAR)" -1 "COMBINED" -2 "`notef2a'" -3 "`notef2b'"  `"`value'"', valuelabel `ylabopt' angle(360)) /*
*/ hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab2')(scatter `obs' `studyvar2', ms(S) msize(*`mscale2') mcolor(gs10))(scatter `obs' `studyvar2', ms(o) msize(*`mscale') mcolor(black) xline(`mvar2', lpattern(-)))/*
*/ (scatteri -1 `mvar2lo' -0.8 `mvar2', clcolor(black) c(l) s(i))(scatteri -1 `mvar2lo' -1.2 `mvar2', clcolor(black) c(l) s(i)) (scatteri -0.8 `mvar2' -1 /*
*/ `mvar2hi', clcolor(black) c(l) s(i))(scatteri -1.2 `mvar2' -1 `mvar2hi', clcolor(black) c(l) s(i)), /*
*/ legend(off) xtitle("`gtitle2'", size(*.5)) title("`plottype'" "`tlab'", size(*0.75)) name(forplot2, replace)
}
else if "`fordata'" == "fordata" {
nois twoway (pcarrow `obs' `studyvar1hi' `obs' `studyvar1lo'  if (`studyvar1lo' == 0.01 & "`forest'" == "dlor"),  lwidth(vthin) lpat(solid)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar1lo' `obs' `studyvar1hi'  if (`studyvar1hi' == 1000 & "`forest'" == "dlr"),  lpat(solid) lwidth(vthin)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar1hi' `obs' `studyvar1lo'  if (`studyvar1lo' == 0.01 & "`forest'" == "dlr"),  lwidth(vthin) lpat(solid)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/(rspike `studyvar1lo' `studyvar1hi' `obs', ylabel(`maxx' "STUDY(YEAR)" -1 "COMBINED" -2 "`notef1a'" -3 "`notef1b'" `"`value'"', valuelabel `ylabopt' angle(360)) /*
*/ hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab1')(scatter `obs' `studyvar1', ms(S) msize(*`mscale2') mcolor(gs10))(scatter `obs' `studyvar1', ms(o) msize(*`mscale') mcolor(black) xline(`mvar1', lpattern(-)))/*
*/ (scatter `obs1' `studyvar1', ms(i) yaxis(2) ylab(`maxx' "`gtitle1' (95% CI)" -3 "`null1'" -2 "`null1'" -1 "`note1f'" `"`value1'"', valuelabel labsize(*`textscale') noticks labgap(*5) angle(360) axis(2)))(scatteri -1 `mvar1lo' -0.8 `mvar1', clcolor(black) c(l) s(i)) /*
*/ (scatteri -1 `mvar1lo' -1.2 `mvar1', clcolor(black) c(l) s(i)) (scatteri -0.8 `mvar1' -1 /*
*/ `mvar1hi', clcolor(black) c(l) s(i))(scatteri -1.2 `mvar1' -1 `mvar1hi', clcolor(black) c(l) s(i)) /*
*/ , legend(off)  yti("", axis(2)) xtitle("`gtitle1'", size(*.5)) title("`plottype'" "`tlab'", size(*0.75)) name(forplot1, replace)


nois twoway (pcarrow `obs' `studyvar2hi' `obs' `studyvar2lo'  if (`studyvar2lo' == 0.01 & "`forest'" == "dlor"),  lpat(solid) lwidth(vthin)   barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar2lo' `obs' `studyvar2hi'  if (`studyvar2hi' == 1000 & "`forest'" == "dlor"),  lwidth(vthin) lpat(solid)  barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (pcarrow `obs' `studyvar2hi' `obs' `studyvar2lo'  if (`studyvar2lo' == 0.01 & "`forest'" == "dlr"),  lpat(solid) lwidth(vthin)   barbsize(0) mlwidth(vthin) ylabel(-1 "" `"`value'"', valuelabel labc(none)))/*
*/ (rspike `studyvar2lo' `studyvar2hi' `obs', ylabel(`maxx' "STUDY(YEAR)" -1 "COMBINED" -2 "`notef2a'" -3 "`notef2b'"  `"`value'"', valuelabel `ylabopt' angle(360)) /*
*/ hor s(i) blpattern(solid) blwidth(vthin) blcolor(black) `xlab2')(scatter `obs' `studyvar2', ms(S) msize(*`mscale2') mcolor(gs10))(scatter `obs' `studyvar2', ms(o) msize(*`mscale') mcolor(black) xline(`mvar2', lpattern(-)))/*
*/ (scatter `obs2' `studyvar2', ms(i) yaxis(2) ylab(`maxx' "`gtitle2' (95% CI)"  -3 "`null1'" -2 "`null1'" -1 "`note2f'" `"`value2'"', valuelabel labsize(*`textscale') noticks  labgap(*5) angle(360) axis(2)))(scatteri -1 `mvar2lo' -0.8 `mvar2', clcolor(black) c(l) s(i)) /*
*/ (scatteri -1 `mvar2lo' -1.2 `mvar2', clcolor(black) c(l) s(i)) (scatteri -0.8 `mvar2' -1 /*
*/ `mvar2hi', clcolor(black) c(l) s(i))(scatteri -1.2 `mvar2' -1 `mvar2hi', clcolor(black) /*
*/ c(l) s(i)), legend(off) xtitle("`gtitle2'", size(*.5))  yti("", axis(2)) title("`plottype'" "`tlab'", size(*0.75)) name(forplot2, replace)

 }

}


/* FAGAN PLOT */
if "`fagan'" == "fagan" {
if `prior' == 1.0 {
	local prev = $prev
}
if `prior' != 1.0 {
local prev = `prior'
}
nois di as text "FAGAN'S(BAYESIAN) NOMOGRAM"
nois di " "
nois di " "

fagani `prev' `mlrp' `mlrn'
}



/*LIKELIHOOD RATIO SCATTERGRAM */

if "`lrmatrix'" == "lrmatrix" {
nois di as text "LIKELIHOOD RATIO SCATTERGRAM "
nois lrmat `lrp' `lrn', sum1(`mlrp' `mlrplo' `mlrphi') sum2(`mlrn' `mlrnlo' `mlrnhi')
}

/*UNCONDITIONAL PREDICTIVE VALUES OR CONDITIONAL PROBABILITY PLOT*/

 if "`pddam'" == "p" {
nois di as text "CONDITIONAL PROBABILITY PLOT"
nois di " "
nois di " "
nois pddami `mtpr' `mtprse' `mtnr' `mtnrse', plot sum1(`mlrp' `mlrplo' `mlrphi') sum2(`mlrn' `mlrnlo' `mlrnhi')
}
else  if "`pddam'" == "r" {
nois di as text "UNCONDITIONAL PREDICTIVE VALUES"
nois di " "
nois di " "
nois pddami `mtpr' `mtprse' `mtnr' `mtnrse', report sum1(`mlrp' `mlrplo' `mlrphi') sum2(`mlrn' `mlrnlo' `mlrnhi')
}



/* BIVARIATE MIXED-EFFECTS METAREGRESSION */

if "`covars'" =="covars" {

noi di in gr "{title: UNIVARIABLE BIVARIATE MIXED-EFFECTS BINARY META-REGRESSION}"
nois midareg $tp $fp $fn $tn $varlist2, `nip' `level'
}

/*SUMMARY ROC CURVE*/
tempvar  x CB1 CB2 CBsens CBspec yroc llyroc ulyroc
tempvar PB1 PB2 PBsn PBsp CPI

/* model-based parameters */

local rhoci = `cov01'/ (`snse' * `spse')
local pred1se = sqrt(`reffs2' + `snse'^2)
local pred2se = sqrt(`reffs1' + `spse'^2)
local rhopred = (`covar' + `cov01')/(`pred1se'*`pred2se')
local NP = 500
range `CPI' 0 `=2* c(pi)' `NP'

/* Parameters for mean operating point and SROC Space */
local mbeta = (max(0.001, `reffs1')/max(0.001, `reffs2'))^.25
local malpha = `sn'*`mbeta' + `sp'/`mbeta'


range `x' 0 1 `NP'
gen double `yroc' = invlogit( ( -logit(`x') /`mbeta' + `malpha') / `mbeta' )
replace `yroc' = 0 if `x' == 1
replace `yroc' = 1 if `x' == 0
integ `yroc' `x',  trapezoid
local AUC = r(integral)
scalar N=r(N_points)
return scalar AUC = `AUC'
local AUClo = min(1.00, (`AUC'+(invnormal(0.975)^2)/(2*N)-invnormal(0.975)*sqrt((`AUC'*(1-`AUC')+/*
*/((invnormal(0.975)^2)/(4*N)))/N))/(1+((invnormal(0.975)^2)/N)))

local AUChi = max(0, (`AUC'+(invnormal(0.975)^2)/(2*N)+invnormal(0.975)*sqrt((`AUC'*(1-`AUC')+/*
*/((invnormal(0.975)^2)/(4*N)))/N))/(1+((invnormal(0.975)^2)/N)))
return scalar AUClo = `AUClo'
return scalar AUChi = `AUChi'

local note: di "AUC = "%3.2f `AUC' " [" %3.2f `AUClo' " - " %3.2f `AUChi' "]"
local snnote: di "SENS = "%3.2f `mtpr' " [" %3.2f `mtprlo' " - " %3.2f `mtprhi' "]"
local spnote: di "SPEC = "%3.2f `mtnr' " [" %3.2f `mtnrlo' " - " %3.2f `mtnrhi' "]"

/* Derivation of parameters for 95% confidence ellipse about mean operating point*/

local kci = sqrt(2*invF(2, `numobs'-2,`level'/100))
gen `CB2' = `sp' + `spse' * `kci' * cos(`CPI')
gen `CB1' = `sn' + `snse' * `kci' * cos(`CPI' + acos(`rhoci'))
gen `CBsens' = invlogit(`CB1')
gen `CBspec' = invlogit(`CB2')

/* Derivation of 95% prediction ellipse*/

gen `PB2' = `sp' + `pred2se' * `kci' * cos(`CPI')
gen `PB1' = `sn' + `pred1se' * `kci' * cos(`CPI' + acos(`rhopred'))
gen `PBsn' = invlogit(`PB1')
gen `PBsp' = invlogit(`PB2')

if "`sroc2'" == "sroc2" {
if "`plottype'" != "" {
local plottype "SROC with Confidence and Predictive Ellipses"
}
else if "`plottype'" == "" {
local plottype " "
}

#delimit;
nois twoway (scatter `sens' `spec', msym(x))
(scatteri `mtpr' `mtnr', msym(D))
(line `yroc' `x', clpat(solid) clc(black))
(line `CBsens' `CBspec', clpat(dash) clc(black) clw(thin))
(line `PBsn' `PBsp', clpat(dot)  clc(black)), plotregion(margin(zero))
xsc(range(0 1) rev) ysc(range(0 1))  xla(0(.5)1, nogrid format(%7.1f))
yla(0(.5)1, nogrid angle(horizontal) format(%7.1f)) xti(Specificity)
yti(Sensitivity)  legend(order(1 "Observed Data" 2 "Summary Operating Point" "`snnote'" "`spnote'"
3 "SROC Curve" "`note'" 4 "`level'% Confidence Ellipse" 5 "`level'% Prediction Ellipse")
pos(2) col(1) size(*.75)) xsize(`hsize') title("`plottype'" "`tlab'", size(*0.75)) ;
#delimit cr
}
if "`sroc1'" == "sroc1" {
if "`plottype'" != "" {
local plottype "SROC without Confidence and Predictive Ellipses"
}
else if "`plottype'" == "" {
local plottype " "
}

#delimit;
nois twoway (scatter `sens' `spec', msym(x))
(scatteri `mtpr' `mtnr', msym(D))
(line `yroc' `x', clpat(solid) clc(black)),
xsc(range(0 1) rev) ysc(range(0 1))  xla(0(.5)1, nogrid format(%7.1f))
yla(0(.5)1, nogrid angle(horizontal) format(%7.1f)) plotregion(margin(zero)) xti(Specificity)
yti(Sensitivity)  legend(order(1 "Observed Data" 2 "Summary Operating Point" "`snnote'" "`spnote'"
3 "SROC Curve" "`note'") pos(2) col(1) size(*.75)) xsize(`hsize') title("`plottype'" "`tlab'", size(*0.75)) ;
#delimit cr
}

nois di ""
nois di ""

if "`results'" == "all" {
nois di "SUMMARY DATA AND PERFORMANCE ESTIMATES"
nois di as text "Bivariate Binomial Mixed Model"
nois di " "
nois di as txt "Number of studies = ", as res `numobs'
nois di " "
nois di as txt "Reference-positive Subjects ="as result %5.0f `sumtpfn'
nois di " "
nois di as txt "Reference-negative Subjects ="as result %5.0f `sumtnfp'
nois di " "
nois di as txt "Pretest Prob of Disease ="as result %5.3f $prev
nois di " "
nois di ""
nois di as txt "Between-study variance(varlogitSEN) ="as result %5.3f `reffs2' _c
nois di as txt ", 95% CI = ["as res %5.3f `reffs2lo'                _c
nois di as txt "-"as res %5.3f `reffs2hi' as txt"]"               _n
nois di " "
nois di as txt "Between-study variance(varlogitSPE)= "as result %5.3f `reffs1' _c
nois di as txt ", 95% CI = ["as res %5.3f `reffs1lo'                _c
nois di as txt "-"as res %5.3f `reffs1hi' as txt"]"               _n
nois di " "
nois di as txt "Correlation (Mixed Model)= " as result %5.3f `rho'
nois di " "
nois di as txt "ROC Area, AUROC = " as res %3.2f `AUC' " [" as res %3.2f `AUClo' " - " as res %3.2f `AUChi' "]"
nois di ""
nois di ""
nois di as txt "Heterogeneity (Chi-square): LRT_Q = " as result %5.3f `lrtchi' ///
as txt ", df =" as result %3.2f `lrtdf' as txt", LRT_p ="as result %5.3f `lrtpchi'
nois di ""
nois di as txt "Inconsistency (I-square): LRT_I2 = "as res %3.2f `Islrt' _c
nois di as txt ", 95% CI = ["as res %3.2f `Islrtlo'                _c
nois di as txt "-"as res %3.2f `Islrthi' as txt"]"               _n
nois di ""
nois di ""
nois di %-28s "Parameter" %8s "Estimate" %16s "`level'% CI"
nois di ""

nois di as text %-28s "Sensitivity" as res %8.3f `mtpr' as text " " "[" as res %8.3f `mtprlo' as text "," as res %8.3f `mtprhi' as text"]"
nois di ""
nois di as text %-28s "Specificity" as res %8.3f `mtnr' as text " " "[" as res %8.3f `mtnrlo' as text "," as res %8.3f `mtnrhi' as text"]"
nois di ""
nois di as text %-28s "Positive Likelihood Ratio" as res %8.3f `mlrp' as text " " "[" as res %8.3f `mlrplo' as text "," as res %8.3f `mlrphi' as text"]"
nois di ""
nois di as text %-28s "Negative Likelihood Ratio" as res %8.3f `mlrn' as text " " "[" as res %8.3f `mlrnlo' as text "," as res %8.3f `mlrnhi' as text"]"
nois di ""
nois di as text %-28s "Diagnostic Score" as res %8.3f `mldor' as text " " "[" as res %8.3f `mldorlo' as text "," as res %8.3f `mldorhi' as text"]"
nois di ""
nois di as text %-28s "Diagnostic Odds Ratio" as res %8.3f `mdor' as text " " "[" as res %8.3f `mdorlo' as text "," as res %8.3f `mdorhi' as text"]"
nois di ""
nois di ""
}
else if "`results'" == "het" {
nois di ""
nois di ""
nois di "GLOBAL HETEROGENEITY STATISTICS"
nois di as txt "Heterogeneity (Chi-square): LRT_Q = " as result %5.3f `lrtchi' ///
as txt ", df =" as result %3.2f `lrtdf' as txt", LRT_p ="as result %5.3f `lrtpchi'
nois di ""
nois di as txt "Inconsistency (I-square): LRT_I2 = "as res %3.2f `Islrt' _c
nois di as txt ", 95% CI = ["as res %3.2f `Islrtlo'                _c
nois di as txt "-"as res %3.2f `Islrthi' as txt"]"               _n
nois di ""
nois di ""
nois di as txt "Between-study variance(varlogitSEN) ="as result %5.3f `reffs2' _c
nois di as txt ", 95% CI = ["as res %5.3f `reffs2lo'                _c
nois di as txt "-"as res %5.3f `reffs2hi' as txt"]"               _n
nois di " "
nois di as txt "Between-study variance(varlogitSPE)= "as result %5.3f `reffs1' _c
nois di as txt ", 95% CI = ["as res %5.3f `reffs1lo'                _c
nois di as txt "-"as res %5.3f `reffs1hi' as txt"]"               _n
nois di " "
}
else if "`results'" == "sum" {
nois di ""
nois di "SUMMARY PERFORMANCE ESTIMATES"
nois di ""
nois di ""
nois di %-28s "Parameter" %8s "Estimate" %16s "`level'% CI"
nois di ""
nois di as text %-28s "Sensitivity" as res %8.3f `mtpr' as text " " "[" as res %8.3f `mtprlo' as text "," as res %8.3f `mtprhi' as text"]"
nois di ""
nois di as text %-28s "Specificity" as res %8.3f `mtnr' as text " " "[" as res %8.3f `mtnrlo' as text "," as res %8.3f `mtnrhi' as text"]"
nois di ""
nois di as text %-28s "Positive Likelihood Ratio" as res %8.3f `mlrp' as text " " "[" as res %8.3f `mlrplo' as text "," as res %8.3f `mlrphi' as text"]"
nois di ""
nois di as text %-28s "Negative Likelihood Ratio" as res %8.3f `mlrn' as text " " "[" as res %8.3f `mlrnlo' as text "," as res %8.3f `mlrnhi' as text"]"
nois di ""
nois di as text %-28s "Diagnostic Score" as res %8.3f `mldor' as text " " "[" as res %8.3f `mldorlo' as text "," as res %8.3f `mldorhi' as text"]"
nois di ""
nois di as text %-28s "Diagnostic Odds Ratio" as res %8.3f `mdor' as text " " "[" as res %8.3f `mdorlo' as text "," as res %8.3f `mdorhi' as text"]"
nois di ""
nois di ""
}
}
}
end


program xtbbrre, rclass sortpreserve byable(recall)
version 10
syntax varlist(min=4 max=4 numeric) [if] [in], [ NIP(integer 7) INDex(string) LEVEL(integer 95) *]

qui {
preserve
marksample touse, novarlist
keep if `touse'
}
tokenize `varlist'
local tp `1'
local fp `2'
local fn `3'
local tn `4'



/* MIXED EFFECTS ESTIMATION */

qui {
local alph = (100-`level')/200
gen study = _n
gen ttruth1 = `tn'
gen ttruth2 = `tp'
gen num1 = `tn'+`fp'
gen num2 = `tp'+`fn'
reshape long num ttruth, i(study) j(dtruth) string
tabulate dtruth, generate(disgrp)
}
xtmelogit (ttruth disgrp1 disgrp2, noc)(study: disgrp1 disgrp2, noc cov(unstr)), bin(num) intp(`nip') var nolr nofet noret nohead

qui {
 nois di " "
 nois di " "
estimates store modr
mat V = e(V)
mat b = e(b)
return local covsnsp = V[1,2]
return scalar mcovar = tanh(_b[atr1_1_1_2:_cons]) * _b[lns1_1_1:_cons] * _b[lns1_1_2:_cons]
nlcom (spbeta: _b[disgrp1])(snbeta: _b[disgrp2])/*
*/(mrho: _b[atr1_1_1_2:_cons])/*
*/(mreffs1: _b[lns1_1_1:_cons])(mreffs2: _b[lns1_1_2:_cons])/*
*/(msens: _b[disgrp2])(mspec: _b[disgrp1])(mldor: _b[disgrp2]+_b[disgrp1]) /*
*/(mdor: _b[disgrp2]+_b[disgrp1]) /*
*/(mlrp: log(invlogit(_b[disgrp2])/(1-invlogit(_b[disgrp1]))))/*
*/(mlrn: log((1-invlogit(_b[disgrp2]))/invlogit(_b[disgrp1]))), post
return scalar mrho = tanh(_b[mrho])
return scalar mrholo = tanh(_b[mrho] - invnorm(1-`alph')*_se[mrho])
return scalar mrhohi = tanh(_b[mrho] + invnorm(1-`alph')*_se[mrho])
ret scalar mreffs1 = exp(_b[mreffs1])^2
ret scalar mreffs1se = exp(_se[mreffs1])^2
ret scalar mreffs1lo = exp(_b[mreffs1] - invnorm(1-`alph') * _se[mreffs1])^2
ret scalar mreffs1hi = exp(_b[mreffs1] + invnorm(1-`alph') * _se[mreffs1])^2
ret scalar mreffs2 = exp(_b[mreffs2])^2
ret scalar mreffs2se = exp(_se[mreffs2])^2
ret scalar mreffs2lo = exp(_b[mreffs2] - invnorm(1-`alph') * _se[mreffs2])^2
ret scalar mreffs2hi = exp(_b[mreffs2] + invnorm(1-`alph') * _se[mreffs2])^2
return scalar mtpr = invlogit(_b[msens])
return scalar mtprlo = invlogit(_b[msens] - invnorm(1-`alph')*_se[msens])
return scalar mtprhi = invlogit(_b[msens] + invnorm(1-`alph')*_se[msens])
return scalar mtprse = (return(mtprhi)-return(mtpr))/invnorm(1-`alph')
return scalar mtnr = invlogit(_b[mspec])
return scalar mtnrlo = invlogit(_b[mspec] - invnorm(1-`alph')*_se[mspec])
return scalar mtnrhi = invlogit(_b[mspec] + invnorm(1-`alph')*_se[mspec])
return scalar mtnrse = (return(mtnrhi)-return(mtnr))/invnorm(1-`alph')
return scalar mldor = _b[mldor]
return scalar mldorlo = _b[mldor] - invnorm(1-`alph')*_se[mldor]
return scalar mldorhi = _b[mldor] + invnorm(1-`alph')*_se[mldor]
return scalar mldorse = (return(mldorhi)-return(mldor))/invnorm(1-`alph')
return scalar mdor = exp(_b[mdor])
return scalar mdorlo = exp(_b[mdor] - invnorm(1-`alph')*_se[mdor])
return scalar mdorhi = exp(_b[mdor] + invnorm(1-`alph')*_se[mdor])
return scalar mdorse = (return(mdorhi)-return(mdor))/invnorm(1-`alph')
return scalar mlrp = exp(_b[mlrp])
return scalar mlrplo = exp(_b[mlrp] - invnorm(1-`alph')*_se[mlrp])
return scalar mlrphi = exp(_b[mlrp] + invnorm(1-`alph')*_se[mlrp])
return scalar mlrpse = (return(mlrphi)-return(mlrp))/invnorm(1-`alph')
return scalar mlrn = exp(_b[mlrn])
return scalar mlrnlo = exp(_b[mlrn] - invnorm(1-`alph')*_se[mlrn])
return scalar mlrnhi = exp(_b[mlrn] + invnorm(1-`alph')*_se[mlrn])
return scalar mlrnse = (return(mlrnhi)-return(mlrn))/invnorm(1-`alph')
return scalar sp = _b[spbeta]
return scalar spse = _se[spbeta]
return scalar splo = _b[snbeta] - invnorm(1-$alph)*_se[spbeta]
return scalar sphi = _b[snbeta] + invnorm(1-$alph)*_se[spbeta]
return scalar sn = _b[snbeta]
return scalar snse = _se[snbeta]
return scalar snlo = _b[snbeta] - invnorm(1-$alph) * _se[snbeta]
return scalar snhi = _b[snbeta] + invnorm(1-$alph) * _se[snbeta]


/* FIXED EFFECTS ESTIMATION */

xtmelogit (ttruth disgrp1 disgrp2, noc)(study: ), bin(num)
estimates store modf
nlcom (fsens: _b[disgrp2])/*
*/(fspec: _b[disgrp1])/*
*/(fldor: (_b[disgrp2]+_b[disgrp1])) /*
*/(fdor: _b[disgrp2]+_b[disgrp1]) /*
*/(flrp: log(invlogit(_b[disgrp2])/(1-invlogit(_b[disgrp1]))))/*
*/(flrn: log((1-invlogit(_b[disgrp2]))/invlogit(_b[disgrp1]))), post

return scalar fsens = invlogit(_b[fsens])
return scalar fspec =  invlogit(_b[fspec])
return scalar fldor = _b[fldor]
return scalar fdor = exp(_b[fdor])
return scalar flrp = exp(_b[flrp])
return scalar flrn = exp(_b[flrn])

/*LRT STATISTICS AND HETEROGENEITY*/

tempname lrtchi lrtpchi lrtdf
lrtest modr modf, stats force
scalar `lrtchi' = r(chi2)
scalar `lrtpchi'= 0.5 * r(p)
scalar `lrtdf' = r(df)
return scalar lrtchi = r(chi2)
return scalar lrtpchi= 0.5 * r(p)
return scalar  lrtdf = r(df)
homogeni `lrtchi' `lrtdf'
return scalar Islrt = r(Isq)
return scalar Islrtlo = r(Isqlo)
return scalar Islrthi = r(Isqhi)
}

end

program bbrre, rclass sortpreserve byable(recall)
version 9
syntax varlist(min=4 max=4 numeric) [if] [in], [ NIP(integer 15) INDex(string) LEVEL(integer 95) *]

qui {
preserve
marksample touse, novarlist
keep if `touse'
}
tokenize `varlist'
local tp `1'
local fp `2'
local fn `3'
local tn `4'

if  !inlist(`nip', 5, 7, 9, 11, 15) {
di as error "nip must be 5, 7, 9, 11 or 15"
error 198
			}


/* MIXED EFFECTS ESTIMATION */

qui {
local alph = (100-`level')/200
gen study = _n
gen ttruth1 = `tn'
gen ttruth2 = `tp'
gen num1 = `tn'+`fp'
gen num2 = `tp'+`fn'
reshape long num ttruth, i(study) j(dtruth) string
tabulate dtruth, generate(disgrp)
eq disgrp1: disgrp1
eq disgrp2: disgrp2
gllamm ttruth disgrp1 disgrp2, nocons i(study) nrf(2) eqs(disgrp1 disgrp2) ///
f(bin) l(logit) denom(num)
mat a = e(b)
local l = e(ll)
local k = e(k)
}
gllamm ttruth disgrp1 disgrp2, nocons i(study) nrf(2) eqs(disgrp1 disgrp2) ///
f(bin) l(logit) denom(num) from(a) lf0(`k' `l') adapt ip(m) nip(`nip')

qui {
 nois di " "
 nois di " "
estimates store modr
mat V = e(V)
mat b = e(b)
scalar cov = V[5,4]
ret scalar mreffs1 =  _b[stu1_1 : disgrp1]^2
ret scalar mreffs1se = 2 * sqrt(return(mreffs1) * _se[stu1_1 : disgrp1]^2)
ret scalar mreffs1lo = return(mreffs1) - invnorm(1-`alph') * return(mreffs1se)
ret scalar mreffs1hi = return(mreffs1) + invnorm(1-`alph') * return(mreffs1se)
ret scalar mreffs2 = _b[stu1_2 : disgrp2]^2 + _b[stu1_2_1 : _cons]^2
return scalar mreffs2se = sqrt( 4*_b[stu1_2 : disgrp2]^2*_se[stu1_2 : disgrp2]^2 /*
*/ + 4* _b[stu1_2_1 : _cons]^2*_se[stu1_2_1 : _cons]^2 /*
 */ + 8* _b[stu1_2 : disgrp2]*_b[stu1_2_1 : _cons]*cov )
ret scalar mreffs2lo = return(mreffs2) - invnorm(1-`alph') * return(mreffs2se)
ret scalar mreffs2hi = return(mreffs2) + invnorm(1-`alph') * return(mreffs2se)
return local covsnsp = V[1,2]
mat chol= e(chol)
mat chol2 = chol*chol'
mat covar= chol2[2,1]
local covar = trace(covar)
return scalar mcovar = trace(covar)
return scalar mrho = `covar'/(sqrt(_b[stu1_1 : disgrp1]^2)* sqrt(_b[stu1_2 : disgrp2]^2 + _b[stu1_2_1 : _cons]^2))

nlcom (spbeta: _b[disgrp1])(snbeta: _b[disgrp2])(msens: _b[disgrp2])/*
*/(mspec: _b[disgrp1])(mldor: _b[disgrp2]+_b[disgrp1]) /*
*/(mdor: _b[disgrp2]+_b[disgrp1]) /*
*/(mlrp: log(invlogit(_b[disgrp2])/(1-invlogit(_b[disgrp1]))))/*
*/(mlrn: log((1-invlogit(_b[disgrp2]))/invlogit(_b[disgrp1]))), post
return scalar mtpr = invlogit(_b[msens])
return scalar mtprlo = invlogit(_b[msens] - invnorm(1-`alph')*_se[msens])
return scalar mtprhi = invlogit(_b[msens] + invnorm(1-`alph')*_se[msens])
return scalar mtprse = (return(mtprhi)-return(mtpr))/invnorm(1-`alph')
return scalar mtnr = invlogit(_b[mspec])
return scalar mtnrlo = invlogit(_b[mspec] - invnorm(1-`alph')*_se[mspec])
return scalar mtnrhi = invlogit(_b[mspec] + invnorm(1-`alph')*_se[mspec])
return scalar mtnrse = (return(mtnrhi)-return(mtnr))/invnorm(1-`alph')
return scalar mldor = _b[mldor]
return scalar mldorlo = _b[mldor] - invnorm(1-`alph')*_se[mldor]
return scalar mldorhi = _b[mldor] + invnorm(1-`alph')*_se[mldor]
return scalar mldorse = (return(mldorhi)-return(mldor))/invnorm(1-`alph')
return scalar mdor = exp(_b[mdor])
return scalar mdorlo = exp(_b[mdor] - invnorm(1-`alph')*_se[mdor])
return scalar mdorhi = exp(_b[mdor] + invnorm(1-`alph')*_se[mdor])
return scalar mdorse = (return(mdorhi)-return(mdor))/invnorm(1-`alph')
return scalar mlrp = exp(_b[mlrp])
return scalar mlrplo = exp(_b[mlrp] - invnorm(1-`alph')*_se[mlrp])
return scalar mlrphi = exp(_b[mlrp] + invnorm(1-`alph')*_se[mlrp])
return scalar mlrpse = (return(mlrphi)-return(mlrp))/invnorm(1-`alph')
return scalar mlrn = exp(_b[mlrn])
return scalar mlrnlo = exp(_b[mlrn] - invnorm(1-`alph')*_se[mlrn])
return scalar mlrnhi = exp(_b[mlrn] + invnorm(1-`alph')*_se[mlrn])
return scalar mlrnse = (return(mlrnhi)-return(mlrn))/invnorm(1-`alph')
return scalar sp = _b[spbeta]
return scalar spse = _se[spbeta]
return scalar splo = _b[snbeta] - invnorm(1-$alph)*_se[spbeta]
return scalar sphi = _b[snbeta] + invnorm(1-$alph)*_se[spbeta]
return scalar sn = _b[snbeta]
return scalar snse = _se[snbeta]
return scalar snlo = _b[snbeta] - invnorm(1-$alph) * _se[snbeta]
return scalar snhi = _b[snbeta] + invnorm(1-$alph) * _se[snbeta]

if "`index'" != "" & "`index'" == "dss" {
return scalar ES1 = invlogit(_b[msens])
return scalar ES1lo = invlogit(_b[msens] - invnorm(1-`alph')*_se[msens])
return scalar ES1hi = invlogit(_b[msens] + invnorm(1-`alph')*_se[msens])
return scalar ES1se = (return(ES1hi)-return(ES1))/invnorm(1-`alph')
return scalar ES2 = invlogit(_b[mspec])
return scalar ES2lo = invlogit(_b[mspec] - invnorm(1-`alph')*_se[mspec])
return scalar ES2hi = invlogit(_b[mspec] + invnorm(1-`alph')*_se[mspec])
return scalar ES2se = (return(ES2hi)-return(ES2))/invnorm(1-`alph')
}
else if "`index'" != "" & "`index'" == "dlor" {

return scalar ES1 = _b[mldor]
return scalar ES1lo = _b[mldor] - invnorm(1-`alph')*_se[mldor]
return scalar ES1hi = _b[mldor] + invnorm(1-`alph')*_se[mldor]
return scalar ES1se = (return(ES1hi)-return(ES1))/invnorm(1-`alph')
return scalar ES2 = exp(_b[mdor])
return scalar ES2lo = exp(_b[mdor] - invnorm(1-`alph')*_se[mdor])
return scalar ES2hi = exp(_b[mdor] + invnorm(1-`alph')*_se[mdor])
return scalar ES2se = (return(ES2hi)-return(ES2))/invnorm(1-`alph')
}
else if "`index'" != "" & "`index'" == "dlr" {
return scalar ES1 = exp(_b[mlrp])
return scalar ES1lo = exp(_b[mlrp] - invnorm(1-`alph')*_se[mlrp])
return scalar ES1hi = exp(_b[mlrp] + invnorm(1-`alph')*_se[mlrp])
return scalar ES1se = (return(ES1hi)-return(ES1))/invnorm(1-`alph')
return scalar ES2 = exp(_b[mlrn])
return scalar ES2lo = exp(_b[mlrn] - invnorm(1-`alph')*_se[mlrn])
return scalar ES2hi = exp(_b[mlrn] + invnorm(1-`alph')*_se[mlrn])
return scalar ES2se = (return(mlrnhi)-return(mlrn))/invnorm(1-`alph')
}


/* FIXED EFFECTS ESTIMATION */

gllamm ttruth disgrp1 disgrp2, nocons i(study)  ///
f(bin) l(logit) denom(num)  init ip(m) nip(`nip') adapt

estimates store modf
nlcom (fsens: _b[disgrp2])/*
*/(fspec: _b[disgrp1])/*
*/(fldor: (_b[disgrp2]+_b[disgrp1])) /*
*/(fdor: _b[disgrp2]+_b[disgrp1]) /*
*/(flrp: log(invlogit(_b[disgrp2])/(1-invlogit(_b[disgrp1]))))/*
*/(flrn: log((1-invlogit(_b[disgrp2]))/invlogit(_b[disgrp1]))), post

return scalar fsens = invlogit(_b[fsens])
return scalar fspec =  invlogit(_b[fspec])
return scalar fldor = _b[fldor]
return scalar fdor = exp(_b[fdor])
return scalar flrp = exp(_b[flrp])
return scalar flrn = exp(_b[flrn])

/*LRT STATISTICS AND HETEROGENEITY*/

tempname lrtchi lrtpchi lrtdf
lrtest modr modf, stats force
scalar `lrtchi' = r(chi2)
scalar `lrtpchi'= 0.5 * r(p)
scalar `lrtdf' = r(df)
return scalar lrtchi = r(chi2)
return scalar lrtpchi= 0.5 * r(p)
return scalar  lrtdf = r(df)
homogeni `lrtchi' `lrtdf'
return scalar Islrt = r(Isq)
return scalar Islrtlo = r(Isqlo)
return scalar Islrthi = r(Isqhi)
}

end


program homogeni, rclass

syntax anything [, Level(int 95) Format(string) ]

tempname Q K df I2 I22 varI2 lb_I2 ub_I2 levelci
tokenize "`anything'"
scalar `Q' = `1'
scalar `df' = `2'
scalar `K' = `df' + 1


if `level' <10 | `level'>99 {
 di in red "level() invalid"
 exit 198
}

scalar `levelci' = `level' * 0.005 + 0.50

if "`format'" == "" {
 local formatI2 = "%4.2f"
 local formatH = "%4.2f"
}
else {
 local formatI2 = "`format'"
 local formatH = "`format'"
}

preserve
tempname varI2 lb_I2 ub_I2
scalar H2 = `Q' / `df'
scalar I2 = max(0, (100*(`Q' -`df')/(`Q' )) )
scalar I22 = max(0, (H2-1)/H2)
if sqrt(H2) < 1 scalar H2 = 1
if `Q' > `K'  {
 scalar SElnH1 = .5*[(log(`Q')-ln(`df')) / ( sqrt(2*`Q') - sqrt(2*`K'-3) )]
}
else {
 scalar SElnH1 = sqrt( ( 1/(2*(`K'-2) )*(1-1/(3*(`K'-2)^2)) )  )
}
scalar `varI2'  = 4*SElnH1^2/exp(4*log(sqrt(H2)))
scalar `lb_I2' = I22-invnorm(`levelci')*sqrt(`varI2')
scalar `ub_I2' = I22+invnorm(`levelci')*sqrt(`varI2')

if  `lb_I2' < 0 {
 scalar  `lb_I2' = 0
}
else scalar `lb_I2' = `lb_I2'
if  `ub_I2' > 1 {
 scalar  `ub_I2' = 1
}
else scalar `ub_I2' = `ub_I2'

return scalar Isq = min(100, 100 * I22)
return scalar Isqlo = max(0, 100 * `lb_I2')
return scalar Isqhi = min(100, 100 * `ub_I2')
return scalar df = `df'
return scalar Q = `Q'

end


program define quadas, sortpreserve
version 10
syntax varlist(min=2) [if] [in] [, LABvar(string) qtable qgraph *]
qui {
preserve
marksample touse, novarlist
keep if `touse'
}

tokenize `varlist'
qui{
tempfile qualires
tempname qualifile
postfile `qualifile' str40 Criterion Yes No Yes_percent No_percent using qualires, replace
foreach var in `varlist' {
count
local totalvar = r(N)
count if `var' == 1
local yesvar = r(N)
count if `var' == 0
local novar = r(N)
local yes_cent = (`yesvar'/`totalvar') * 100
local no_cent = (`novar'/`totalvar') * 100

if "`labvar'" != "" {
local critvar : variable label `var'
post `qualifile' ("`critvar'") (`yesvar') (`novar') (`yes_cent') (`no_cent')
}
else {
post `qualifile' ("`var'") (`yesvar') (`novar') (`yes_cent') (`no_cent')
}
}
postclose `qualifile'
postutil clear
use qualires, clear
summarize Yes
}
local N = r(N)
if "`qtable'"=="qtable" {
di as text "{title: METHODOLOGICAL QUALITY ASSESSMENT}"
di " "
di as text"{hline 83}"
di as text _col(2) "Criterion" _col(44)  "{c |}"  _col(48) "Yes" _col(54) "Yes(%)" _col(64) "No" _col(72) "No(%)"
di as text"{hline 43}{c +}{hline 39}"
local i = 1
while `i' <= `N' {
local a1 = Criterion in `i'
local b1 = Yes in `i'
local c1 = Yes_percent  in `i'
local d1 = No  in `i'
local e1 = No_percent  in `i'
di as text _col(2)  "`a1'" _col(44) in gr "{c |}" as result _col(45) %6.0f `b1' _col(50) %6.0f `c1' _col(60) %6.0f `d1' _col(68) %6.0f `e1'
local i=`i'+1
}
di as text"{hline 43}{c BT}{hline 39}"
}
if "`qgraph'"=="qgraph" {
#delimit;
graph hbar (asis) Yes No, over(Criterion, sort(Total) descending)
nolab  bar (1, fcolor(gs0)) bar(2, fcolor(gs16)) legend(rows(1) position(6))
stack percent title("BAR GRAPH OF QUALITY ASSESSMENT", size(*.75)) lintensity(*.50) ysize(`vsize');
#delimit cr
}
end

program fagani, rclass
version 10
syntax anything, [ *]

tempname prev lrp lrn sens spec
qui{

	tokenize "`anything'"
	scalar `prev' = `1'
	scalar `lrp' = `2'
	scalar `lrn' = `3'
	local prprob = logit(1-`prev')
	local postprob1 = logit(`prev') + log(`lrp')
	local postprob2 = logit(`prev') + log(`lrn')


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'" "'
}

	local priorprob = 100*`prev'
	local postprobpos = 100*invlogit(`postprob1')
	local postprobneg = 100*invlogit(`postprob2')

	local notebb1: di "Prior Prob (%) = " %5.0f `priorprob' "
	local notebb2: di "LR_Positive = " %5.0f `lrp' "
	local notebb3: di "Post_Prob_Pos (%) =" %5.0f `postprobpos' "
	local notebb4: di "LR_Negative = " %5.2f `lrn' "
	local notebb5: di "Post_Prob_Neg (%) = " %5.0f `postprobneg' "

#delimit;
tw (scatteri 0 0, mcolor(none) yaxis(1) ylab(`ylab', angle(0)
tpos(cross)) yscale(reverse axis(1)) ytitle("Pre-test Probability (%)", axis(1)))
(scatteri 0 0, mcolor(none) yaxis(2) ylab(`ylab', angle(0) tpos(cross) axis(2))
 ytitle("Post-test Probability (%)", axis(2)))
(scatteri `lrpts', msymbol(+) mcolor(black) mlabcolor(black) mlabsize(medsmall))
(pci -3.4538776 0 3.4538776 0, recast(pcspike) lcolor(black)
xscale(range(-1 1)) plotregion(margin(zero)) xsize(4) ysize(6)
xscale(off) ylab(, nogrid) text(-4 0 "Likelihood Ratio", place(n)))
(scatteri `prprob' -1, msym(D) yaxis(2))(pcarrowi `prprob' -1 `postprob1' 1,
 yaxis(2) lpat(solid) lwidth(vthin))(pcarrowi `prprob' -1 `postprob2' 1,
 yaxis(2) lpat(dash) lwidth(thin)),subt(Fagan's Nomogram)
 legend(order(5 "`notebb1'" 6 "`notebb2'" "`notebb3'" 7 "`notebb4'"
 "`notebb5'") pos(6) col(1)  rowgap(2) size(*.90)) name(faganplot, replace);
#delimit cr
}
end

program pddami, rclass
version 10
syntax anything, [SUM1(numlist min=3 max=3) SUM2(numlist min=3 max=3) C1(real 0) C2(real 1.0) LEVEL(integer 95) REPort PLOT*]
	tempname contpr contprse contprlo contprhi
	tempname contnr contnrse contnrlo contnrhi
	tempname sigma1 uppv uppvlo uppvhi
	tempname sigma2 unpv unpvlo unpvhi

	qui{
	tokenize "`anything'"
	scalar `contpr' = `1'
	scalar `contprse' = `2'
	scalar `contnr' = `3'
	scalar `contnrse' = `4'

	tokenize `sum1'
          local convar1 `1'
          local convar1lo `2'
          local convar1hi `3'
          tokenize `sum2'
          local convar2 `1'
          local convar2lo `2'
          local convar2hi `3'


       	if `level' < 10 | `level' > 99 {
	di as error "level() must be between 10 and 99"
	exit 198
		}

          if `c1' < 0 | `c2' > 1.0 {
	di as error "c1-c2 must be between 0 and 1"
	exit 198
		}

          if `c1' > `c2' {
	di as error "c1 must always be less than c2 must be between 0 and 1"
	exit 198
		}

          local alph = (100-`level')/200
          tempvar PPP1 PPN1 x1 x2
          twoway__function_gen y= `convar1' * x/(1- x * (1 - `convar1')), r(0 1) x(x) gen(`PPP1' `x1', replace) n(`c(N)')
          twoway__function_gen y= `convar2' * x /(1 - x * (1 - `convar2')), r(0 1) x(x) gen(`PPN1' `x2', replace) n(`c(N)')
          local note1a: di "LR+ =" %4.2f `convar1' " [" %4.2f `convar1lo' " - " %4.2f `convar1hi' "]"
          local note2a: di "LR- =" %4.2f `convar2' " [" %4.2f `convar2lo' " - " %4.2f `convar2hi' "]"


 if "`plot'" == "plot" {
          #delimit;
	nois twoway (line `PPP1' `x1', sort clpat(dash) clwidth(medium) connect(direct ) clcolor(black))
	(line `PPN1' `x2', sort clpat(shortdash_dot) clcolor(black) clwidth(medium) connect(direct ))
 	(function y=x, sort range(0 1)clcolor(black) clpat(solid) clwidth(vthin) connect(direct)),
 	ytitle("Posterior Probability", size(*.90)) yscale(range(0 1))
 	ylabel( 0(.2)1,  angle(horizontal) format(%3.1f)) xtitle("Prior Probability", size(*.90))
 	xscale(range(0 1)) xlabel(0(.2)1, format(%3.1f)) legend(order(1 "Positive Test Result"
 	"`note1a'" 2 "Negative Test Result" "`note2a'") pos(6) colgap(1) row(1) size(*.75))
 	aspect(1) title("Conditional Probability Plot", size(*.75)) plotregion(margin(zero));
	#delimit cr

}

          *estimate unconditional positive or negative predictive values and the 95% confidence intervals
          *assuming Uniform distribution for prior from zero to one.
          *It can be easily adjusted for other uniform priors (by modifying the values of `c1' and `c2')

	scalar `contprlo' = `contpr'-invnorm(1-`alph')*`contprse'
	scalar `contprhi' = `contpr'+invnorm(1-`alph')*`contprse'
	scalar `contnrlo' = `contnr'-invnorm(1-`alph')*`contnrse'
	scalar `contnrhi' = `contnr'+invnorm(1-`alph')*`contnrse'

   	scalar `uppv' = `contpr'/(`contpr'+`contnr'-1)-`contpr'*(1-`contnr')*log((`c2'*`contpr'+(1-`c2')*(1-`contnr'))/(`c1'*`contpr'+(1-`c1')*(1-`contnr')))/((`c2'-`c1')*(`contpr'+`contnr'-1)^2)
   	scalar `sigma1' = (1/(`contpr'+`contnr'-1)^4)*((`contprse')^2*(`contnr'-1-`contpr'*(1-`contnr')^2/((`c2'*`contpr'+(1-`c2')*(1-`contnr'))*(`c1'*`contpr'+(1-`c1')*(1-`contnr')))-((1-`contnr')*(`contnr'-1-`contpr')/((`c2'-`c1')*(`contnr'-1+`contpr')))*log((`c2'*`contpr'+(1-`c2')*(1-`contnr'))/(`c1'*`contpr'+(1-`c1')*(1-`contnr'))))^2+(`contnrse')^2*(-`contpr'-(`contpr'^2)*(1-`contnr')/((`c2'*`contpr'+(1-`c2')*(1-`contnr'))*(`c1'*`contpr'+(1-`c1')*(1-`contnr')))+(`contpr'*(`contpr'-`contnr'+1)/((`c2'-`c1')*(`contpr'+`contnr'-1)))*log((`c2'*`contpr'+(1-`c2')*(1-`contnr'))/(`c1'*`contpr'+(1-`c1')*(1-`contnr'))))^2)
   	scalar `uppvlo' = `uppv'-invnorm(1-`alph')*sqrt(`sigma1')
   	scalar `uppvhi' = `uppv'+invnorm(1-`alph')*sqrt(`sigma1')

   	scalar `unpv' = `contnr'/(`contpr'+`contnr'-1)-`contnr'*(1-`contpr')*log(((1-`c1')*`contnr'+`c1'*(1-`contpr'))/((1-`c2')*`contnr'+`c2'*(1-`contpr')))/((`c2'-`c1')*(`contpr'+`contnr'-1)^2)
	scalar `sigma2' = (1/(`contnr'+`contnr'-1)^4)*((`contnrse')^2*(`contnr'-1-`contnr'*(1-`contnr')^2/((`c2'*`contnr'+(1-`c2')*(1-`contnr'))*(`c1'*`contnr'+(1-`c1')*(1-`contnr')))-((1-`contnr')*(`contnr'-1-`contnr')/((`c2'-`c1')*(`contnr'-1+`contnr')))*log((`c2'*`contnr'+(1-`c2')*(1-`contnr'))/(`c1'*`contnr'+(1-`c1')*(1-`contnr'))))^2+(`contnrse')^2*(-`contnr'-(`contnr'^2)*(1-`contnr')/((`c2'*`contnr'+(1-`c2')*(1-`contnr'))*(`c1'*`contnr'+(1-`c1')*(1-`contnr')))+(`contnr'*(`contnr'-`contnr'+1)/((`c2'-`c1')*(`contnr'+`contnr'-1)))*log((`c2'*`contnr'+(1-`c2')*(1-`contnr'))/(`c1'*`contnr'+(1-`c1')*(1-`contnr'))))^2)
   	scalar `unpvlo' = `unpv'-invnorm(1-`alph')*sqrt(`sigma2')
   	scalar `unpvhi' = `unpv'+invnorm(1-`alph')*sqrt(`sigma2')

   	*returned results
   	return scalar contpr = `contpr'
   	return scalar tprlo = `contprlo'
	return scalar tprhi =  `contprhi'
	return scalar contnr = `contnr'
	return scalar tnrlo = `contnrlo'
	return scalar tnrhi = `contnrhi'
   	return scalar uppv = `uppv'
  	return scalar uppvlo = max(0, `uppvlo')
  	return scalar uppvhi = min(1.00, `uppvhi')
  	return scalar unpv = `unpv'
  	return scalar unpvlo = max(0, `unpvlo')
  	return scalar unpvhi = min(1.00, `unpvhi')
  	return scalar c1 = `c1'
  	return scalar c2 = `c2'

  	if "`report'" == "report" {
  	nois di ""
	nois di as txt "Prior Prevalence = Uniform [" as res %3.2f return(c1) "," as res %3.2f return(c2) as text"]"
          nois di ""
	nois di as txt "Unconditional Positive Predictive Value= "as res %3.2f return(uppv) _c
          nois di as txt ", 95% CI = ["as res %3.2f return(uppvlo)                _c
          nois di as txt "-"as res %3.2f return(uppvhi) as txt"]"               _n
          nois di ""
          nois di as txt "Unconditional Negative Predictive Value= "as res %3.2f return(unpv) _c
          nois di as txt ", 95% CI = ["as res %3.2f return(unpvlo)                _c
          nois di as txt "-"as res %3.2f return(unpvhi) as txt"]"               _n
          }

}
end
program bvbox, sortpreserve
	version 10
	syntax varlist(numeric min=2 max=2) [if] [in][, smooth(str asis) data(str asis) `options' *]

	// observations to use
	marksample touse

	tempvar use diff sum diff1 sum1 work theta radius radius1 sm order
          tempname ymed ymad xmed xmad dmad smad

	// variables set-up
	tokenize `varlist'
	args y x

	quietly {
		// initialize
		gen `work' = .
		gen byte `use' = .
		gen `diff' = .
		gen `sum' = .
		gen `radius' = .
		gen `diff1' = .
		gen `sum1' = .
		gen `radius1' = .
		gen `theta' = .
		gen `order' = .

		rreg `x' `y'
                    predict xpredict
                    rreg `y' `x'
                    predict ypredict
                    local ypred "(line ypredict `x', lpat(dash) lc(black) lw(thin))"
                    local xpred "(line `y' xpredict, lpat(dash) lc(black) lw(thin))"
		tempvar s sx S SX s1 S1  sx1 SX1

		replace `use' = `touse' & `y' < .

		// y <- (y - median y) / MAD y
		mata : median("`y'", "`use'")
		scalar `ymed' = r(p50)
                    replace `work' = abs(`y' - `ymed') if `use'
		mata : median("`work'", "`use'")
		scalar `ymad' = r(p50)
                    clonevar `s' = `y'
		if `"`: variable label `s''"' == "" label var `s' "`y'"
		replace `s' = (`y' - `ymed') / `ymad' if `use'

		// x <- (x - median x) / MAD x
		mata : median("`x'", "`use'")
		scalar `xmed' = r(p50)
		replace `work' = abs(`x' - `xmed') if `use'
		mata : median("`work'", "`use'")
		scalar `xmad' = r(p50)
		gen `sx' = (`x' - `xmed') / `xmad' if `use'

		// (y - x), (y + x) scaled to z / MAD z
		replace `diff' = `s' - `sx' if `use'
		mata : median("`diff'", "`use'")
		replace `work' = abs(`diff' - r(p50))
		mata : median("`work'", "`use'")
		scalar `dmad' = r(p50)
		replace `diff' = `diff' / `dmad'

		replace `sum' = `s' + `sx' if `use'
		mata : median("`sum'", "`use'")
		replace `work' = abs(`sum' - r(p50))
		mata : median("`work'", "`use'")
		scalar `smad' = r(p50)
		replace `sum' = `sum' / `smad'

		// radius = cube root of sum^2 + diff^2
		// theta = arctan of diff / sum
		replace `radius' = sqrt(`sum'^2 + `diff'^2)
		replace `radius' = `radius'^(2/3)
		replace `theta' = atan2(`diff', `sum')

		local sc
		tempvar S C
		gen `S' = sin(`theta')
		gen `C' = cos(`theta')
		local sc `sc' `S' `C'

		capture regress `radius' `sc' if `use'
		if _rc gen `sm' = . if `use'
		else predict `sm' if `use'

		drop `sc'

		// reverse transformation, (x, y) coordinates, scale
		replace `radius' = `sm'^(3/2)
		replace `radius1' = 1.58 * `sm'^(3/2)
		drop `sm'
		replace `diff' = `dmad' * `radius' * sin(`theta')
		replace `sum' = `smad' * `radius' * cos(`theta')
		replace `sx' = `xmed' + `xmad' * (`sum' - `diff')/2
		replace `s' = `ymed' + `ymad' * (`sum' + `diff')/2
		replace `diff1' = `dmad' * `radius1' * sin(`theta')
		replace `sum1' = `smad' * `radius1' * cos(`theta')
		gen `sx1' = `xmed' + `xmad' * (`sum1' - `diff1')/2
		gen `s1' = `ymed' + `ymad' * (`sum1' + `diff1')/2


		// sort order and end points for closing loop
		bysort `use' (`theta') : replace `order' = _n if `use'
		count if !`use'
		local first = 1 + r(N)
		gen `SX' = `sx' in `first'
		gen `S' = `s' in `first'
		replace `SX' = `sx' in l
		replace `S' = `s' in l
		gen `SX1' = `sx1' in `first'
		gen `S1' = `s1' in `first'
		replace `SX1' = `sx1' in l
		replace `S1' = `s1' in l

                    local hshade "fcolor(gs12) nodropbase recast(area)"
                    local fshade "fcolor(gs14) nodropbase recast(area)"
                    local clpcw1 "clpat(solid) clc(black) clw(thin)"
                    local clpcw2 "clpat(dash) clc(black) clw(thin)"


		// construct graph call
		// line plot of smooth
		local l "(line `s' `sx', `clpcw1' `hshade' `smooth')"
		local l1 "(line `s1' `sx1', `clpcw2' `fshade' `smooth')"
		// line plot to close loop
		local p	"(line `S' `SX', `clpcw1' `hshade' `smooth')"
		local p1	"(line `S1' `SX1', `clpcw2' `fshade' `smooth')"

		// scatter of data
		local d "(scatter `y' `x' if `touse', ms(+) `data')"

                    local call "`call'  `l1' `p1' `l' `p' `d' `xpred' `ypred'"


	local off "legend(off)"

	// final graph preparation
	sort `order'
	local yttl : variable label `y'
	if `"`yttl'"' == "" local yttl "`y'"
	local xttl : variable label `x'
	if `"`xttl'"' == "" local xttl "`x'"

	// graph
	#delimit;
    	twoway `call', ti("Bivariate Boxplot", size(*.75)) yti(`"`yttl'"') xti(`"`xttl'"')
	legend(order(`show')) plotregion(margin(zero)) `off' `options';
	#delimit cr
}
end



mata :
void median(string scalar varname, string scalar which)
{
	real scalar L, U
	real matrix X

	st_view(X, ., varname, which)
	X = sort(X,1)
	L = ceil(rows(X) / 2)
	U = floor((rows(X) + 2) / 2)

	st_numscalar("r(p50)", (X[L,1] + X[U,1]) / 2)
}
end



program define midachi
version 10

syntax varlist(num min=2 max=2)[if][in][, *]
marksample touse
tokenize `varlist'
local chivar1 `1'
local chivar2 `2'
qui {
tempvar ry rx Hi Fi Gi Si CHIi Li
 * Gi
 egen `ry' = rank(`chivar1'), field
 gsort -`ry'
 gen `Gi' = (_N - `ry') / (_N - 1)
 * Fi
 egen `rx' = rank(`chivar2'), field
 gsort -`rx'
 gen `Fi' = (_N - `rx') / (_N - 1)
 * Hi
 sort `ry'
  by `ry': replace `ry' = _N
 sort `chivar1'
 tempname xi
 gen `Hi' = 0
 local r1 = 1
 local N = _N
 forvalues i = 1 / `N' {
  if `chivar1'[`i'] == `chivar1'[`i'-1] {
 local r1 = `r1' + 1
 }
    else {
 local r1 = 1
}
    local k = min(`N', `i' + `ry'[`i'] - `r1')
    scalar `xi' = `chivar2'[`i']
    count if `chivar2' <= `xi' & _n != `i' in 1/`k'
    replace `Hi' = r(N) in `i'
 }
 replace `Hi' = `Hi' / (_N - 1)
 * Si, CHIi, Li
  gen `Si'   = sign((`Fi' - .5)*(`Gi' - .5))
 gen `CHIi' = (`Hi' - `Fi'*`Gi') / (`Fi'*(1 - `Fi')*`Gi'*(1 - `Gi'))^.5
 gen `Li'   = 4 * `Si' * max((`Fi' - .5)^2, (`Gi' - .5)^2)
 label var `CHIi' "CHI"
 label var `Li' "LAMDA"


 spearman `chivar1' `chivar2'
 local r = r(rho)
 local note: di "rho = " %3.2f `r' "

*Scatterplot
#delimit;
twoway (scatter `chivar1' `chivar2', ms(o) msize(medium))(lfit `chivar1' `chivar2'),
name(splot, replace) title(ScatterPlot (`note')) plotregion(margin(zero)) legend(off) nodraw `options';
#delimit cr

* cp lines

 local cp = 1.78
 local cphi = `cp' / sqrt(_N)
 local cplo= -`cp' / sqrt(_N)

*chi-plot

#delimit;
twoway (scatter `CHIi' `Li', ms(S) msize(*`mscale')), yline(`cplo' `cphi', lpat(solid) lwidth(vvthin))
yline(0, lpat(dash) lwidth(vvthin)) xline(0, lpat(dash) lwidth(vvthin))
xla(-1(0.5)1) yla(-1(.5)1) title(ChiPlot) plotregion(margin(zero)) nodraw name(cplot, replace) `options';
#delimit cr
nois graph combine splot cplot, xsize(`hsize') ysize(`vsize')
}

end

program define lrmat, rclass byable(recall) sortpreserve
version 10

syntax varlist(min=2 max=2) [if] [in] , SUM1(numlist min=3 max=3) SUM2(numlist min=3 max=3)[ *]


tokenize `varlist'
local var1 `1'
local var2 `2'
tokenize `sum1'
local mvar1 `1'
local mvar1lo `2'
local mvar1hi `3'
tokenize `sum2'
local mvar2 `1'
local mvar2lo `2'
local mvar2hi `3'

qui{
local note11: di "LUQ: Exclusion & Confirmation"
local note11b: di "LRP>10, LRN<0.1"
local note12: di "RUQ: Confirmation Only"
local note12b: di "LRP>10, LRN>0.1"
local note13: di "LLQ: Exclusion Only"
local note13b: di "LRP<10, LRN<0.1"
local note14: di "RLQ: No Exclusion or Confirmation"
local note14b: di "LRP<10, LRN>0.1"

#delimit;
nois twoway (scatter `var1' `var2', sort msymbol(Oh) msize(small) mcolor(black))
(scatteri `mvar1' `mvar2', msymbol(D) msize(large) clcol(black) clwidth (medium))
(scatteri `mvar1' `mvar2lo' `mvar1' `mvar2hi', recast(line) clcol(black) clpat(solid) clwidth (medium))
(scatteri `mvar1lo' `mvar2' `mvar1hi' `mvar2', recast(line) clcol(black) clpat(solid) clwidth (medium)),
xtitle("Negative Likelihood Ratio", size(*.90)) xsc(log) xlab(0.1 "0.1"  1)
xmticks(.02(0.01).09 .2(.1).9) ylab(1 10 100, angle(horizontal))ymticks(2(1)9 20(10)90)
xline(0.1, lpattern(shortdash) lwidth(vthin)) ytitle("Positive Likelihood Ratio", size(*.90))
 ysc(log) yline(10, lpattern(shortdash) lwidth(vthin)) title(Likelihood Ratio Scattergram, size(*.90))
 legend(order(3 "`note11'"  "`note11b'"  "`note12'" "`note12b'" 4 "`note13'"  "`note13b'"  "`note14'"
  "`note14b'" 2 "Summary LRP & LRN for Index Test" "With `level' % Confidence Intervals")
 pos(2) symxsize(0) forcesize rowgap(1) col(1) size(*.75)) xsize(`hsize') plotregion(margin(zero))
name(ScatterMatrix, replace);
#delimit cr
}
end



program midareg, rclass sortpreserve byable(recall)
version 10
syntax varlist(min=4) [if] [in], [ NIP(integer 15) LEVEL(integer 95) *]

qui {
preserve
marksample touse, novarlist
keep if `touse'
}
tokenize `varlist'
local tp `1'
local fp `2'
local fn `3'
local tn `4'
macro shift 4
local varlist2 `*'

if "`varlist2'" =="" {
di as error " covariate varlist must be used with covars"
exit 198
}


/* MIXED EFFECTS ESTIMATION */

qui {
local alph = (100-`level')/200
local ctitle1 "Sensitivity"
local ctitle2 "Specificity"

gen ctruth1 = `tn'
gen ctruth2 = `tp'
gen cnum1 = `tn'+`fp'
gen cnum2 = `tp'+`fn'
gen cstudy = _n
reshape long cnum ctruth, i(cstudy) j(cdtruth) string
tabulate cdtruth, generate(cdisgrp)

tempname covfile
tempfile covresults
postfile `covfile' str30 Parameter coef1 cVAR1 cVAR1lo cVAR1hi z1 p1 coef2 cVAR2 cVAR2lo cVAR2hi /*
*/ z2 p2 LRTCHI LRTPCHI I2 I2lo I2hi using covresults, replace
xtmelogit (ctruth cdisgrp1 cdisgrp2, noc)(cstudy: cdisgrp1 cdisgrp2, noc cov(unstr)), bin(cnum) intp(`nip')
estimates store mod0
foreach var in `varlist2' {
   local varlab "`var'"
   su `var', meanonly
   replace `var'= `var'-r(mean) if r(mean)>1
   forvalues i=1/2 {
     g `var'_`i' = cdisgrp`i'*`var'
   }
nois di " "
nois di " "
nois di in gr "Estimating Covariate Effect Of: " in white " `varlab'"
nois di " "
xtmelogit (ctruth cdisgrp1 cdisgrp2 `var'_1 `var'_2, noc)(cstudy: cdisgrp1 cdisgrp2, noc cov(unstr)), bin(cnum) intp(`nip') nolr nofet noret nohead

estimates store mod`var'
nlcom (csens: _b[cdisgrp2] + _b[`var'_2]) /*
*/(csens0: _b[cdisgrp2]) /*
*/(cspec: _b[cdisgrp1] + _b[`var'_1]) /*
*/(cspec0: _b[cdisgrp1]), post

local coef2 = _b[cspec]
local cvar2= invlogit(_b[cspec])
local cvar2se=_se[cspec]
local cvar2lo=invlogit(_b[cspec]-invnorm(1-$alph) * _se[cspec])
local cvar2hi=invlogit(_b[cspec] + invnorm(1-$alph) * _se[cspec])
local coef1 = _b[csens]
local cvar1= invlogit(_b[csens])
local cvar1se=_se[csens]
local cvar1lo=invlogit(_b[csens]- invnorm(1-$alph) * _se[csens])
local cvar1hi= invlogit(_b[csens] + invnorm(1-$alph) * _se[csens])
local z_cov2=(_b[cspec]-_b[cspec0])/sqrt((_se[cspec]^2) + (_se[cspec0]^2))
if `z_cov2' <=0{
local p_cov2=2*normal(`z_cov2')
}
else {
local p_cov2=2*(1-normal(`z_cov2'))
}
local z_cov1=(_b[csens]-_b[csens0])/sqrt((_se[csens]^2) + (_se[csens0]^2))
if `z_cov1' <=0{
local p_cov1=2*normal(`z_cov1')
}
else {
local p_cov1=2*(1-normal(`z_cov1'))
}
qui lrtest mod0 mod`var', stats force
local LRTchi = r(chi2)
local LRTpchi = r(p)
local LRTdf = r(df)
homogeni `LRTchi' `LRTdf'
scalar I2 = r(Isq)
scalar I2lo = r(Isqlo)
scalar I2hi = r(Isqhi)
post `covfile' ("`var'") (`coef1') (`cvar1') (`cvar1lo') (`cvar1hi') (`z_cov1') (`p_cov1') ///
(`coef2') (`cvar2') (`cvar2lo') (`cvar2hi') (`z_cov2') (`p_cov2') (`LRTchi') (`LRTpchi') ///
(I2) (I2lo) (I2hi)
 }
postclose `covfile'
postutil clear
use covresults, clear
format coef1 z1 p1 cVAR1 cVAR1lo cVAR1hi coef2 z2 p2 cVAR2 cVAR2lo cVAR2hi LRTCHI LRTPCHI I2 I2lo I2hi  %7.2f
format Parameter %-30s
foreach var of varlist cVAR1 cVAR2 I2 {
tostring `var'lo `var' `var'hi, gen(`var'lo1 `var'1 `var'hi1) format(%7.2f) force
replace `var'lo1=" " + "[" + `var'lo1 +" - "
replace `var'hi1= `var'hi1+ "]"
egen `var'_ci= concat(`var'1 `var'lo1 `var'hi1)
format `var'_ci %50s force
}
nois di ""
summarize p1
local N = r(N)
nois di ""
nois di ""
nois di ""
nois di as text "`ctitle1'"
nois di ""
nois di as text"{hline 78}"
nois di as text _col(2) "Parameter" _col(20)  "{c |}"  _col(24) "Estimate(95%CI)" _col(50) "Coef" _col(60) "Z" _col(70) "P>|z|"
nois di as text"{hline 19}{c +}{hline 59}"
local i = 1
while `i' <= `N' {
local a1 = Parameter in `i'
local b1 = cVAR1_ci in `i'
local c1 =  coef1 in `i'
local d1 = z1  in `i'
local e1 = p1  in `i'
nois di as text _col(2)  "`a1'" _col(20) in gr "{c |}" as result _col(22)  "`b1'" _col(48) %6.2f `c1' _col(58) %6.2f `d1' _col(68) %6.2f `e1'
local i=`i'+1
}
nois di as text"{hline 19}{c BT}{hline 59}"
nois di ""
nois di ""
nois di ""
nois di as text "`ctitle2'"
nois di ""
nois di as text"{hline 78}"
nois di as text _col(2) "Parameter" _col(20)  "{c |}"  _col(24) "Estimate(95%CI)" _col(50) "Coef" _col(60) "Z" _col(70) "P>|z|"
nois di as text"{hline 19}{c +}{hline 59}"
local i = 1
while `i' <= `N' {
local a1 = Parameter in `i'
local b1 = cVAR2_ci in `i'
local c1 =  coef2 in `i'
local d1 = z2  in `i'
local e1 = p2  in `i'
nois di as text _col(2)  "`a1'" _col(20) in gr "{c |}" as result _col(22) "`b1'" _col(48) %6.2f `c1' _col(58) %6.2f `d1' _col(68) %6.2f `e1'
local i=`i'+1
}
nois di as text"{hline 19}{c BT}{hline 59}"
nois di ""
nois di ""
nois di ""
nois di "Joint Model"
nois di ""
nois di as text"{hline 78}"
nois di as text _col(2) "Parameter" _col(20)  "{c |}"  _col(24) "I-squared(95%CI)" _col(50) "LRTChi" _col(60) "P value"
nois di as text"{hline 19}{c +}{hline 59}"
local i = 1
while `i' <= `N' {
local a1 = Parameter in `i'
local b1 = I2_ci in `i'
local c1 =  LRTCHI in `i'
local d1 = LRTPCHI  in `i'
nois di as text _col(2)  "`a1'" _col(20) in gr "{c |}" as result _col(22) "`b1'" _col(48) %6.2f `c1' _col(58) %6.2f `d1'
local i=`i'+1
}
nois di as text"{hline 19}{c BT}{hline 59}"
}

end


program define metanorm
version 9.0

// Setup
syntax varlist(numeric min=2 max=2) [if] [in] ///
       [ , Title(str) ID(varname) LEvel(real 95) XTitle(string) YTitle(string) * ]

tokenize `varlist'
tempvar theta setheta

local theta `1'
local setheta `2'

tempvar yq z psubi obs nobs normal normse normallo normalhi
tempname kurtvar skewvar pchi
     qui {
     gen `yq'= `theta'/`setheta'
     label var `yq' "Effect Size"
     sort `yq'
     local xmax=r(max)
     gen `obs'=_n
     gen `nobs'=_N
     gen double `psubi'=(`obs'-0.5)/`nobs'
     gen `z'=invnormal(`psubi')
     label var `z' "Normal Quantile"
     sum `yq', detail
     local sd=r(sd)
     local ymean=r(mean)
     gen double `normal'=`ymean' + `z' * `sd'
     gen double `normse'=(`sd'/((1/sqrt(8 * atan(1)))*exp(-0.5 * `z' * `z')))* sqrt(`psubi' * (1-`psubi')/`nobs')
     gen double `normallo'=`normal'-2*`normse'
     gen double `normalhi'=`normal'+2*`normse'
  	}

     qui sktest `yq'
     scalar `skewvar' =r(P_skew)
     scalar `kurtvar' =r(P_kurt)
     scalar `pchi' =r(P_chi2)
     di as text "TESTING NORMALITY ASSUMPTIONS"
     di ""
     di ""
     di as text "Skewness Test: pvalue = " as result %4.3f `skewvar'
     di ""
     di ""
	di as text "Kurtosis Test: pvalue = " as result %4.3f `kurtvar'
     di ""
     di ""
	di as text "Kurtosis/Skewness Test: pvalue = " as result %4.3f `pchi'

     twoway (scatter `yq' `z', sort ms(O))(line `normal' `normallo' `normalhi' `z', clpat(solid dash dash) clwidth(medium medium medium)),  /*
     */ legend(off) ytitle("Standardized Effect Size", size(medsmall)) plotregion(margin(zero)) /*
     */ylab(,  angle(horizontal) format(%7.0f)) xtitle("Normal Quantile", size(medsmall)) `options'

end

program define midagalb
version 9.0

	syntax varlist(numeric min=2 max=2) [if] [in] ///
       [ , Title(str) ID(varname) LEvel(real 95) XTitle(string) YTitle(string) * ]

	tokenize `varlist'
	tempvar theta setheta

	local theta `1'
	local setheta `2'
qui {
        tempvar x y
        su `theta',detail
        local emax=r(max)
        local emay=r(min)
        gen `x' = 1 / `setheta'
        su `x' , detail
        local maxx = r(max)
        gen `y' = `theta'/`setheta'
        reg `y' `x', noconstant

    local galbropts " yscale(r(-2 2) noline) ylab(-2 2,angle(horizontal) nogrid )"
    local yaxis2 " ysc(r(`emay' `emax') fex)  yti("Unstandardized effect size",  axis(2)) ylab(`emay'(1.0)`emax', angle(horizontal) format(%7.1f) axis(2))"
  }
#delimit;
    twoway(scatter `y' `x', `show' `scatter')(scatteri -2 0 2 0, s(i) recast(line))
    (function fitted = _b[`x'] * x, ra(0 `maxx') `fitted') (function upper = 2 + _b[`x'] * x, ra(0 `maxx')
	clp(dash) clc(green) `upper' )(function lower = -2 + _b[`x'] * x, ra(0 `maxx')
    clp(dash) clc(green) `lower')(scatter `theta' `x', s(i)  `scatter' yaxis(2))
    , legend(off)  yti("standardized effect size",) xti("precision") yline(0, lpat(shortdash))
	`galbropts' `yaxis2' plotregion(margin(zero)) `options' ;
#delimit cr
 end

program define copas, sortpreserve
version 10
syntax varlist(min=2 max=2) [if] [in] [, *]
qui {
preserve
marksample touse, novarlist
keep if `touse'
}
tokenize `varlist'
local y `1'
local yse `2'
qui {
postutil clear
postfile biasfile nstudies maxbias using biasresults, replace
local N = _N
egen sumvar = total(`yse'^-2)
local sumvar = sumvar
egen sumsvar = total(`y' * `yse'^-2)
local sumsvar = sumsvar
egen sumse = total(1/`yse')
local sumse = sumse
scalar upbias = abs(sumsvar/sumvar) - (1.96/sqrt(sumvar))
display upbias
gen maxbias=.
forvalues M = 1/1000 {
gen b`M' = ((_N + `M')/_N) * normalden(invnormal(_N /(_N+ `M'))) * (sumse/sumvar)
replace maxbias = b`M'
post biasfile (`M') (maxbias)
if maxbias >= upbias {
continue, break
}
}
postclose biasfile
postutil clear
use biasresults, clear
count
local nn = r(N)
local prob = `N'/(`nn' + `N')
noi di ""
noi di ""
noi di ""

local note1: di "Estimated Number of Missing Studies: " `nn'
local note2: di "Estimated Publication Probability: " %3.2f `prob'
sum maxbias, detail
local xmin = r(min)
local num = maxbias[_N]
noi di ""
noi di ""
local xmax = `nn' + 5

twoway__function_gen y= ((`N' + x)/`N') * normalden(invnormal(`N' /(`N' + x)))*(`sumse'/`sumvar'), r(1 `xmax') x(x) gen(y x1, replace) n(100)

if `nn' >= 5 {
#delimit;
nois twoway (line y x, lp(solid) lw(medium)lc(black))(scatteri `num' `nn'), xtitle("Number of Missing Studies", size(*.75))
legend(order(2 "Maximum Bias Point" "`note1'" "`note2'") colgap(1) pos(6) col(1) size(*.50)) ylab(, format(%9.3f) angle(hor) labs(*.75) ) xsc(r(`xmin' `xmax')) plotregion(margin(zero))
yline(`num', lpat(dash) lw(thin)) xline(`nn', lpat(dash) lw(thin)) xlab(, format(%9.0f) labs(*.75)) ytitle("Maximum Bias", size(*.75));
#delimit cr
}
else if `nn' < 5 {
nois di "Estimated Number of Missing Studies: " `nn'
nois di "Estimated Publication Probability: " %3.2f `prob'
}
}
end