*! diagtest AT 1.1.0, 6 January 1999    (STB-56: sbe36)
* diagtest AT 1.0.0, 21 December 1998

program define diagtest
    version 5.0
    local varlist "req ex min(2) max(2)"
    local if "opt"
    local in "opt"
    local weight "fweight pweight"
    local options "prev(real 0) `options' level(int 95)*"
    parse "`*'"
    parse "`varlist'", parse(" ")
    local test `1'
    local diag `2'

    preserve
    if ("`if'"!="") {
       keep `if'
    } 
    if ("`in'"!="") {
       keep `in'
    } 

    local z=invnorm((`level'/100)+((1-(`level'/100))/2))

    qui summ `diag'
    local min=_result(5)

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

    global S_1 = (`a'/(`a'+`c'))*100
    global S_10 =( ($S_1/100)-`z'*sqrt( ($S_1/100)*(1-($S_1/100))/`n') )*100
    global S_11 =( ($S_1/100)+`z'*sqrt( ($S_1/100)*(1-($S_1/100))/`n') )*100
    global S_2 = (`d'/(`b'+`d'))*100
    global S_20 =( ($S_2/100)-`z'*sqrt( ($S_2/100)*(1-($S_2/100))/`n') )*100
    global S_21 =( ($S_2/100)+`z'*sqrt( ($S_2/100)*(1-($S_2/100))/`n') )*100
    
    if `prev'==0 {
		global S_3 = (`a'/(`a'+`b'))*100       
    		global S_30 =( ($S_3/100)-`z'*sqrt( ($S_3/100)*(1-($S_3/100))/`n') )*100
		global S_31 =( ($S_3/100)+`z'*sqrt( ($S_3/100)*(1-($S_3/100))/`n') )*100

		global S_4 = (`d'/(`c'+`d'))*100       
    		global S_40 =( ($S_4/100)-`z'*sqrt( ($S_4/100)*(1-($S_4/100))/`n') )*100
    		global S_41 =( ($S_4/100)+`z'*sqrt( ($S_4/100)*(1-($S_4/100))/`n') )*100

		local pD=((`a'+`c')/(`a'+`b'+`c'+`d')) 
                global S_5 = `pD'*100  
                global S_50 =( ($S_5/100)-`z'*sqrt( ($S_5/100)*(1-($S_5/100))/`n') )*100
                global S_51 =( ($S_5/100)+`z'*sqrt( ($S_5/100)*(1-($S_5/100))/`n') )*100
	}
    else {
		local prev=`prev'/100
		local sens=$S_1/100
		local spec=$S_2/100
		global S_3 = ((`sens'*`prev') / ((`sens'*`prev')+((1-`spec')*(1-`prev'))))*100      
		global S_30="."
		global S_31="."
		global S_4 = ((`spec'*(1-`prev')) / (((1-`sens')*`prev')+(`spec'*(1-`prev'))))*100 
		local pD=`prev'
		global S_40="."
		global S_41="."
    		global S_5 = `pD'*100
		global S_50="."
		global S_51="."
	}

    di
    di in gr "True D defined as " in gr "`diag'" in gr " ~= " in gr `min' _col(55) in gr "[" in gr `level' in gr "% Conf. Inter.]"                                                  
    di in gr "-------------------------------------------------------------------------"                              
    di in gr "Sensitivity                     Pr( +| D) " _col(44) in ye %5.2f $S_1 "%" _col(56) in ye %5.2f $S_10 "%" _col(65) in ye %5.2f $S_11 "%"
    di in gr "Specificity                     Pr( -|~D) " _col(44) in ye %5.2f $S_2 "%" _col(56) in ye %5.2f $S_20 "%" _col(65) in ye %5.2f $S_21 "%"
    di in gr "Positive predictive value       Pr( D| +) " _col(44) in ye %5.2f $S_3 "%" _col(56) in ye %5.2f $S_30 "%" _col(65) in ye %5.2f $S_31 "%"
    di in gr "Negative predictive value       Pr(~D| -) " _col(44) in ye %5.2f $S_4 "%" _col(56) in ye %5.2f $S_40 "%" _col(65) in ye %5.2f $S_41 "%"
    di in gr "-------------------------------------------------------------------------"                              
    di in gr "Prevalence                      Pr(D)     " _col(44) in ye %5.2f $S_5 "%" _col(56) in ye %5.2f $S_50 "%" _col(65) in ye %5.2f $S_51 "%"
    di in gr "-------------------------------------------------------------------------"
end