*! version 4.0.5  19dec1998  (STB-52 sbe31)
/* exactcc - A special version of cc that offers adjusted cornfield confidence
 *    limits and exact confidence limits.
 */
program define exactcc, rclass
	version 6
	syntax varlist(min=2 max=2) [if] [in] [fw] [, BY(varname) BD /*
		*/ Binomial(varname) /*
		*/ noCrude Exact Fast EStandard noHom noHEt IStandard /*
		*/ Level(integer $S_level) Pool Standard(varname) TB Woolf ]
	local hom "`hom'`het'"		/* het is old option name */

	tokenize `varlist'
	local cas `1'
	local ex `2'

	local stdflg=("`standar'"!="")+("`istanda'"!="")+("`estanda'"!="")
	if `stdflg'>1 { error 198 }

	marksample use
	tempvar WGT one first decoded
	
	* checks that weight and binomial are not both present
	if "`binomial'"!="" & "`weight'"!="" {
		di in re "weight not allowed with binomial frequency records"
		exit 198
	}
	if "`binomial'"~="" {
		preserve
		qui keep if `use'
		tempvar case mywt 
		bitowt `cas' `binomial' , c(`case') w(`mywt')
		local exp "= `mywt'"
		qui drop `cas'
		tempvar cas
		rename `case' `cas'
		local weight "fweight"
	}
	if `"`weight'"'!="" { 
		qui gen double `WGT' `exp' if `use'
		local w "[fweight=`WGT']"
	}

	if "`by'"!="" { 
		capture replace `use'=0 if `by'==.	/* fail on str ok */
		quietly { 
			sort `use' `by'
			if "`standar'"!="" {
				capture assert `standar'>=0 if `use'
				if _rc {
					noi di in red /*
					*/ "`standar' has negative values"
					exit 498
				}
				capture by `use' `by': assert /*
				*/ `standar'==`standar'[_n-1] if _n>1 & `use'
				if _rc {
					noi di in red /*
					*/ "`standar' not constant within `by'"
					exit 198
				}
			}
			by `use' `by': gen byte `first'=1 if _n==1 & `use'
			sort `first' `by'

			if "`crude'" == "" {
				local optex = ("`woolf'"=="") +         /*
		*/	("`tb'"=="") + ("`exact'"=="") + ("`level'"=="") 
				
				if "`level'" != "" {
					local ml "level(`level')" 
				}
				else 	local ml

				if `optex' == 4 {
					local myopt
				}
				else {
			  	  local myopt ",`woolf' `exact' `ml' `tb'  "
				}
				cc `cas' `ex' `if' `in' `w' `myopt'
				local CRUDE  `r(or)'
				local C0     `r(lb_or)'
				local C1     `r(ub_or)'
				ret add
				
				local CTYPE  "(Cornfield)"
				if "`tb'" != ""    { local CTYPE  "(tb)" }
				if "`woolf'" != "" { local CTYPE  "(Woolf)" }
			}

		}
	}
	else {
		if `stdflg' { error 198 }
	}
	quietly gen byte `one'=1
	_crcvarl `ex'
	local col "col(`s(varlabel)')"

	if "`by'"=="" { 
		quietly { 
			safesum `one' `w' if `cas' & `ex' & `use'
			local a=r(sum)
			safesum `one' `w' if `cas' & `ex'==0 & `use'
			local b=r(sum)
			safesum `one' `w' if `cas'==0 & `ex' & `use'
			local c=r(sum)
			safesum `one' `w' if `cas'==0 & `ex'==0 & `use'
			local d=r(sum)
		}
		exactcci `a' `b' `c' `d', level(`level') `col' `woolf' `tb' `exact'
		ret add
		exit
	}
	if `level'<10 | `level'>99 { local level 95 } 
	local lby : value label `by'
	if "`lby'"!="" { 
		quietly decode `by', gen(`decoded')
		local lby "`decoded'"
	}
	else	local lby `by'
	if `stdflg'==0 { local hdr2 "M-H" }
	else	local hdr2 "   "

	local lbl : var label `by'
	if "`lbl'"=="" { local lbl "`by'" }
	else local lbl = substr("`lbl'",1,16)
	local skip=16-length("`lbl'")
	di _n in gr _skip(`skip') "`lbl' |" _col(26) "OR" /*
	*/ _col(34) "[`level'% Conf. Interval]    `hdr2' Weight" _n /*
	*/ _dup(17) "-" "+" _dup(49) "-"

	local i 1
	local iz=invnorm(1-(1-`level'/100)/2)
	local sumwgt 0
	local sumor 0
	local sumPR 0
	local sumR 0
	local sumPSQR 0
	local sumS 0
	local sumQS 0
	local sumXnum 0
	local sumXden 0
	local or0 . 
	local or1 . 
	local pnum 0
	local pden 0
	local nf
	local df 0

	tempvar irr virr Sf 
	tempname test
	scalar `test'=0
	local bddf 0
	qui gen `irr' = .
	qui gen `virr' = .
	qui gen `Sf'=.
	
	while 1 { 
		if `first'[`i']!=1 { 
			if `stdflg'==0 {
				local or=`sumor'/`sumwgt'
				local se=sqrt(`sumPR'/(2*((`sumR')^2))+/*
					*/ `sumPSQR'/(2*`sumR'*`sumS')+/*
					*/ `sumQS'/(2*((`sumS')^2)))
				if "`tb'"=="" { 
					local or0=exp(ln(`or')-`iz'*`se')
					local or1=exp(ln(`or')+`iz'*`se')
				}
				else { 
					local x=`sumXnum'/sqrt(`sumXden')
					local or0=(`or')^(1-`iz'/`x')
					local or1=(`or')^(1+`iz'/`x')
					if `or0'>`or1' { 
						local h `or0'
						local or0 `or1'
						local or1 `h'
					}
					local 	note "(tb)"
				}
local x2=((`sumXnum')^2)/`sumXden'
				*local x2=((`sumXnum')^2)/`sumXden'
				global S_1 `x2'
				ret scalar chi2 = `x2'
				local hdr2 "   M-H combined"
			}
			else {
				if "`istanda'"!="" { 
					local hdr2 "I. Standardized" 
				}
				else if "`estanda'"!="" {
					local hdr2 "E. Standardized" 
				}
				else local hdr2 "   Standardized"
				local or=`sumwgt'/`sumor' /* sic */
				local se=sqrt(`sumR')/`sumwgt'
				local or0=exp(ln(`or')-`iz'*`se')
				local or1=exp(ln(`or')+`iz'*`se')
				global S_1
				ret scalar chi2 = .
			}
			di in gr _dup(17) "-" "+" _dup(49) "-" 

			if "`crude'" == "" {
				di in gr "           Crude |  " in ye /*
			*/	%9.0g `CRUDE' "    " %9.0g `C0' "  "  /*
			*/	%9.0g `C1'    in gr _col(69) "`CTYPE'"
			}

			if "`pool'"!="" {
				local POOL = exp(`pnum'/`pden')
				local sd = sqrt(1/`pden')
				local P0 = exp(log(`POOL')-`iz'*`sd')
				local P1 = exp(log(`POOL')+`iz'*`sd')
				
				di in gr " Pooled (direct) |  " in ye /*
			*/	%9.0g `POOL' "    " %9.0g `P0' "  "  /*
			*/	%9.0g `P1'    
			}
			di in gr " `hdr2' |  " in ye /* 
			*/ %9.0g `or' "    " %9.0g `or0' "  " %9.0g `or1' /*
			*/ in gr _col(69) "`note'" 

/* Here is where the test of homogeneity goes */

			if "`hom'"=="" {
				if "`pool'"!="" | `stdflg'==0 {
					di in gr _dup(17)"-" "+" _dup(49) "-"
				}
				if "`pool'"!="" {
					local dfm1=`df'-1
					qui replace `Sf'=sum( /*
					*/ (log(`irr')-log(`POOL'))^2/`virr')/*
					*/ in 1/`df'
					local pchi = `Sf'[`df']
					quietly summ `irr' in 1/`df'
					local cdf = r(N) - 1
					*local pval = chiprob(`dfm1',`pchi')
					local pval = chiprob(`cdf',`pchi')
					if `cdf' > 0 {
						di in gr /*
						*/ "Test of homogeneity " /*
						*/ "(direct)" _col(32) "chi2(" /*
						*/ in ye "`cdf'" in gr ") = " /*
						*/in ye %8.2f `pchi'  /*
						*/ in gr "  Pr>chi2 = " /*
						*/ in ye %6.4f /*
						*/ `pval'
						global S_26 = `cdf'
						global S_27 = `pchi'
						ret scalar df = `cdf'
						ret scalar chi2_p = `pchi'
					}
				}
				if `stdflg'==0 {
					local dfm1=`df'-1
					qui replace `Sf'=sum( /*
					*/ (log(`irr')-log(`or'))^2/`virr')/*
					*/ in 1/`df'
					local pchi = `Sf'[`df']
					quietly summ `irr' in 1/`df'
					local cdf = r(N) - 1
					*local pval = chiprob(`dfm1',`pchi')
					local pval = chiprob(`cdf',`pchi')

					if `cdf' > 0 {
						di in gr /*
						*/"Test of homogeneity "   /*
						*/ "(M-H)" _col(32) "chi2(" /*
						*/in ye "`cdf'" in gr ") = " /*
						*/in ye %8.2f `pchi'  /*
						*/ in gr "  Pr>chi2 = " /*
						*/ in ye %6.4f /*
						*/ `pval'
						global S_26 = `cdf'
						global S_27 = `pchi'
						ret scalar df = `cdf'
						ret scalar chi2_p = `pchi'
					}
					if "`bd'"!="" {
						HOMOBD1 `cas' `ex' `first' /*
						*/ `use' `w', by(`by') or(`or')
						tempname
						scalar `test'= r(homo)
						local bddf=r(bddf)-1
						if `bddf' > 0 {
						di in gr /*
						*/"Test of homogeneity " /*
						*/ "(B-D)" _col(32) "chi2(" /*
						*/in ye `bddf' in gr ") = " /*
						*/in ye %8.2f `test' /*
						*/in gr "  Pr>chi2 = " /*
						*/ in ye %6.4f /*
						*/ chiprob(`bddf', `test')
						ret scalar chi2_bd=`test'
						ret scalar df_bd=`bddf'
						}
					}
				}
				
			}	

			if `stdflg'==0 {
				di _n _col(20) in gr /*
				*/ "Test that combined OR = 1:" _n _col(33) /*
				*/ "Mantel-Haenszel chi2(1) =" /*
				*/ in ye %10.2f `x2' _n in gr /*
				*/ _col(49) "Pr>chi2 =" in ye /*
				*/ %10.4f chiprob(1,`x2')
			}

			ret local p        /* this is the S_2 thing */
			ret local p_exact  /* this is also the S_2 thing */
			ret local rd     
			ret local lb_rd  
			ret local ub_rd  
			ret local rr     
			ret local lb_rr  
			ret local ub_rr  
			ret local afe    
			ret local lb_afe 
			ret local ub_afe 

			ret scalar or     = `or'
			ret scalar lb_or  = `or0'
			ret scalar ub_or  = `or1'

		/* the following because of double saving in S_# and r() */
			global S_1=return(chi2)
			global S_2
			global S_3
			global S_4
			global S_5
			global S_6
			global S_7
			global S_8
			global S_9 `or'
			global S_10 `or0'
			global S_11 `or1'
			global S_12
			global S_13
			global S_14
			global S_24
                        *if "`nf'"=="error" & `stdflg'==0{
                          *di
                          *di in bl "Note: Stratum specific estimates appear"/*
                          **/    " both sides of the null value. "
                          *di in bl "      This problem may require " /*
			  * */ "standardization."
                        *}
			exit
		}
		quietly { 
			safesum `one' `w' /*
			*/ if `cas' & `ex' & `use' & `by'==`by'[`i']
			local a=r(sum)
			safesum `one' `w' /* 
			*/ if `cas' & `ex'==0 & `use' & `by'==`by'[`i']
			local b=r(sum)
			safesum `one' `w' /*
			*/ if `cas'==0 & `ex' & `use' & `by'==`by'[`i']
			local c=r(sum)
			safesum `one' `w' /*
			*/ if `cas'==0 & `ex'==0 & `use' & `by'==`by'[`i']
			local d=r(sum)
		}
		local m1=`a'+`b'
		local m0=`c'+`d'
		local n1=`a'+`c'
		local n0=`b'+`d'
		local t=`n1'+`n0'

		local or=(`a'*`d')/(`b'*`c')
		

		if `stdflg'==0 {
			local wgt=(`b'*`c')/`t'
			*local sumor=`sumor'+`wgt'*`or'
			local sumor=`sumor'+(`a'*`d')/`t'
			local sumwgt=`sumwgt'+`wgt'

			local P=(`a'+`d')/`t'
			local Q=(`b'+`c')/`t'
			local R=(`a'*`d')/`t'
			local S=(`b'*`c')/`t'
			local sumPR=`sumPR'+`P'*`R'
			local sumR=`sumR'+`R'
			local sumPSQR=`sumPSQR'+`P'*`S'+`Q'*`R'	/* sic */
			local sumS=`sumS'+`S'
			local sumQS=`sumQS'+`Q'*`S'

			local sumXnum=`sumXnum'+`a'-(`n1'*`m1')/`t'
			local sumXden=`sumXden'+/*
				*/ (`n1'*`n0'*`m1'*`m0')/((`t')^2*(`t'-1))
		}
		else {
			if "`istanda'"!="" { local wgt `c' }
			else if "`estanda'"!="" { local wgt `d' }
			else local wgt=`standar'[`i']
			local waoc=`wgt'*`a'/`c'
			local sumwgt=`sumwgt'+`waoc'
			local sumor=`sumor'+`wgt'*`b'/`d'
			local sumR=`sumR'+((`waoc')^2)*(1/`a'+1/`b'+1/`c'+1/`d')
		}

		global S_4
		ret local lb_rd		/* -ret drop- */
		if "`fast'"=="" { 
			_crcor `a' `b' `c' `d' `level' `woolf' `tb'
/*
			if "`nf'"=="" {
				if r(nf) < 1.0 { local nf "below" }
				else           { local nf "above" }
			}
			else {
				if r(nf) < 1.0 & "`nf'"=="above" {
					local nf "error"
				}
				else if r(nf) > 1.0 & "`nf'"=="below" {
					local nf "error"
				}
			}
*/
			local or0 = r(lb)
			local or1 = r(ub)
		}
		local lbl=`lby'[`i']
		local lbl=substr("`lbl'",1,16)
		local skip=16-length("`lbl'")
		di in gr _skip(`skip') "`lbl' |  " in ye /*
		*/ %9.0g `or' "    " %9.0g `or0' "  " %9.0g `or1' /*
		*/ "     " %9.0g `wgt' in gr " `r(label)'"

		qui replace `irr'  = (`a'*`d')/(`c'*`b') in `i'
		qui replace `virr' =  1/`a'+1/`b'+1/`c'+1/`d' in `i'
		
		local pnum = `pnum'+log(`irr'[`i'])/`virr'[`i']
		local pden = `pden'+1/`virr'[`i']

		local i=`i'+1
		local df=`df'+1
	}
	/*NOTREACHED*/
end

program define HOMOBD1, rclass
	syntax varlist(min=4 max=4) [fw] , BY(varname) or(string)
	tokenize `varlist'
	local cas `1'
	local ex `2'
	local first `3'
	local use `4'
	tempvar WGT one 
	quietly { 
		if "`weight'"!="" { 
			gen double `WGT' `exp' if `use'
			local w "[fweight=`WGT']"
		}
		sort `first' `by'
		gen byte `one'=1
	}
	tempname test
	scalar `test'=0
	local i 1	
	while 1 { 
		quietly { 
			safesum `one' `w' /*
			*/ if `cas' & `ex' & `use' & `by'==`by'[`i']
			local a=r(sum)
			safesum `one' `w' /* 
			*/ if `cas' & `ex'==0 & `use' & `by'==`by'[`i']
			local b=r(sum)
			safesum `one' `w' /*
			*/ if `cas'==0 & `ex' & `use' & `by'==`by'[`i']
			local c=r(sum)
			safesum `one' `w' /*
			*/ if `cas'==0 & `ex'==0 & `use' & `by'==`by'[`i']
			local d=r(sum)
		}
		
		HOMOBD `a' `b' `c' `d' `or'

		if r(homo)~=. {
			scalar `test'=`test'+ r(homo)
		}
		if `test'~=. {
			local bddf=`bddf'+1
		}
		local i=`i'+1
		if `first'[`i']~=1 {
			ret scalar homo=`test'
			ret scalar bddf=`bddf'
			exit
		}
	}
end

program def HOMOBD, rclass
args A B C D rr 
        local n1=`A'+`B'
        local n0=`C'+`D'
        local m1=`A'+`C'
        tempname a b c d r1 r2 var or test
        scalar `a'=1-`rr'
        scalar `b'=`n0'-`m1'+(`m1'*`rr')+(`n1'*`rr')
        scalar `c'=-`rr'*`n1'*`m1'

        scalar `r1'=(-`b'+sqrt((`b'^2) - (4*`a'*`c') )) / (2*`a')
        scalar `r2'=(-`b'-sqrt((`b'^2) - (4*`a'*`c') )) / (2*`a')

        if `r1'<= `n1' & `r1'<= `m1' & `r1'>=0 {
                scalar `a' = `r1'
        }
        else scalar `a' = `r2'
        scalar `b'=`n1'-`a'
        scalar `c'=`m1'-`a'
        scalar `d'=`n0'-`c'
        scalar `var'=1/( (1/`a')+(1/`b')+(1/`c')+(1/`d') )
        scalar `or'=`a'*`d'/(`b'*`c')
        scalar `test' = ((`A'-`a')^2)/`var'
	ret scalar homo=`test'
end