*! diagt 1.1.2, 5 Oct 2000                       (STB-59: sbe36.1)
*! by PT Seed (paul.seed@kcl.ac.uk)
*! based on diagtest.ado (Aurelio Tobias, STB-56: sbe36)

program define diagt , rclass
version 6.0

* Synatx

	syntax varlist(min=2 max=2 numeric) [if] [in] [fweight] , [prev(real 0) level(real $S_level) *]

	tokenize "`varlist'"
	local diag `1'
	local test `2'	

	preserve
	qui drop if `diag' == . | `test' == .
	qui summ `diag'
	cap assert `diag' == r(min) | `diag' == r(max)
	qui replace `diag' = `diag' == r(max)
	local case = r(max)
	local caselab : label (`diag') `case'
	if "`caselab'" ~= "`case'" { local caselab "`case' (labelled `caselab')" }

	qui summ `test'
	cap assert `test' == r(min) | `test' == r(max)
	qui replace `test' = `test' == r(max)

* Sensitivity
	qui ci `test' if `diag' [`weight'`exp'], bin level(`level')
	local sens = 100*r(mean)
	local sens_lb = 100*r(lb)
	local sens_ub = 100*r(ub)

* Specificity
	qui ci `test' if ~`diag' [`weight'`exp'], bin level(`level')
	local spec = 100*(1-r(mean))
	local spec_lb = 100*(1-r(ub))
	local spec_ub = 100*(1-r(lb))

	if `prev'==0 {
* PPV
		qui ci `diag' if `test'  [`weight'`exp'], bin level(`level')
		local ppv = 100*r(mean)
		local ppv_lb = 100*r(lb)
		local ppv_ub = 100*r(ub)
* NPV
		qui ci `diag' if ~`test'  [`weight'`exp'], bin level(`level')
		local npv = 100*(1-r(mean))
		local npv_lb = 100*(1-r(ub))
		local npv_ub = 100*(1-r(lb))
* Prevalence
		qui ci `diag' [`weight'`exp'], bin level(`level')
		local prev = 100*r(mean)
		local prev_lb = 100*r(lb)
		local prev_ub = 100*r(ub)

		return scalar ppv_lb  = `ppv_lb'
		return scalar ppv_ub = `ppv_ub'
		return scalar npv_lb  = `npv_lb'
		return scalar npv_ub = `npv_ub'
		return scalar prev_lb  = `prev_lb'
		return scalar prev_ub = `prev_ub'

	}

	else {
	local prev_ub "."
	local prev_lb "."
	local ppv = ((`sens'*`prev') / ((`sens'*`prev')+((1-`spec')*(1-`prev'))))*100	  
	local ppv_ub "."
	local ppv_lb "."
	local npv = ((`spec'*(1-`prev')) / (((1-`sens')*`prev')+(`spec'*(1-`prev'))))*100 
	local npv_ub "."
	local npv_lb "."

	return scalar prev_ub = .
	return scalar prev_lb = .
	return scalar ppv = ((`sens'*`prev') / ((`sens'*`prev')+((1-`spec')*(1-`prev'))))*100	  
	return scalar ppv_ub = .
	return scalar ppv_lb = .
	return scalar npv = ((`spec'*(1-`prev')) / (((1-`sens')*`prev')+(`spec'*(1-`prev'))))*100 
	return scalar npv_ub = .
	return scalar npv_lb = .
	}


* Table
	qui replace `test'  = 1 - `test' 
	qui replace `diag'  = 1 - `diag'  

	tempname labtest
	label define `labtest' 0 "Pos." 1 "Neg."
	label values `test' `labtest'

	tempname labdiag
	label define `labdiag' 0 "Abnormal" 1 "Normal"
	label values `diag' `labdiag'

	tabulate `diag' `test' [`weight'`exp'], matcell(T) `options'
	local a=T[1,1]
	local b=T[1,2]
	local c=T[2,1]
	local d=T[2,2]    


	di in gr "True abnormal diagnosis defined as `diag' = `caselab'" 
	di _n in gr _col(55) "[`level'% Conf. Inter.]"

    di in gr "-------------------------------------------------------------------------"  
    di in gr "Sensitivity                     Pr( +| D) " _col(44) in ye %5.2f `sens' "%" _col(56) in ye %5.2f `sens_lb' "%" _col(65) in ye %5.2f `sens_ub' "%"
    di in gr "Specificity                     Pr( -|~D) " _col(44) in ye %5.2f `spec' "%" _col(56) in ye %5.2f `spec_lb' "%" _col(65) in ye %5.2f `spec_ub' "%"
    di in gr "Positive predictive value       Pr( D| +) " _col(44) in ye %5.2f `ppv' "%" _col(56) in ye %5.2f `ppv_lb' "%" _col(65) in ye %5.2f `ppv_ub' "%"
    di in gr "Negative predictive value       Pr(~D| -) " _col(44) in ye %5.2f `npv' "%" _col(56) in ye %5.2f `npv_lb' "%" _col(65) in ye %5.2f `npv_ub' "%"
    di in gr "-------------------------------------------------------------------------"                              
    di in gr "Prevalence                      Pr(D)     " _col(44) in ye %5.2f `prev' "%" _col(56) in ye %5.2f `prev_lb' "%" _col(65) in ye %5.2f `prev_ub' "%"
    di in gr "-------------------------------------------------------------------------"

	return scalar sens = `sens'
	return scalar sens_lb  = `sens_lb'
	return scalar sens_ub = `sens_ub'
	return scalar spec = `spec'
	return scalar spec_lb  = `spec_lb'
	return scalar spec_ub = `spec_ub'
	return scalar ppv = `ppv'
	return scalar npv = `npv'
	return scalar prev = `prev'

end
exit