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

program define diagti , rclass
version 6.0

	tokenize "`*'", parse(" ,")
	local i= 1
	while `i' <= 4 {
		confirm integer number ``i'' 
		local i = `i' + 1
	}
	local a = `1'
	local b = `2'
	local c = `3' 
	local d = `4'

	mac shift 4
	local options "prev(real 0) level(real $S_level) *"
	parse "`*'" 
	
* Sensitivity
	local n1_ = `a' + `b'
	qui cii `n1_' `a' , level(`level')
	local sens = 100*r(mean)
	local sens_lb = 100*r(lb)
	local sens_ub = 100*r(ub)

* Specificity
	local n2_ = `c' + `d'
	qui cii `n2_' `d', level(`level')
	local spec = 100*r(mean)
	local spec_lb = 100*r(lb)
	local spec_ub = 100*r(ub)

	if "`prev'" == "0" {

* PPV
		local n_1 = `a' + `c'
		qui cii `n_1' `a', level(`level')
		local ppv = 100*r(mean)
		local ppv_lb = 100*r(lb)
		local ppv_ub = 100*r(ub)
* NPV
		local n_2 = `b' + `d'
		qui cii `n_2' `d', level(`level')
		local npv = 100*r(mean)
		local npv_lb = 100*r(lb)
		local npv_ub = 100*r(ub)

* Prevalence
		local n = `n_1' + `n_2'
		qui cii `n' `n1_', 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_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
	if _N >0 { preserve }
	clear 
	qui tabi `a' `b' \ `c' `d'
	qui expand pop
	qui drop if pop == 0
	qui drop pop

	rename col test
	label var test "Test result"
	label define test 1 "Pos." 2 "Neg."
	label values test test

	rename row true
	label var true "True disease status"
	label define true 1 "Abnormal" 2 "Normal"
	label values true true

	tabulate true test, matcell(T) `options'

	di
    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