/*
 *  ======== exactcci ===============================================
 *
 * Special version of cci with additional output.  Cornfield's confidence
 * interval for the odds ratio is calculated both with and without
 * continuity correction.
 *
 * Reference:
 * Dupont, W.D., Plummer, W.D. 1999. SBE31: Exact confidence intervals for odds ratios
 * from case-control studies.  Stata Technical Bulletin 52: 12-16. 
 * 
 *  Programs defined in this file are:
 *
 *  exactcci		The main exactcci command routine
 *  _xcrnfd		Calculates adjusted Cornfield confidence limits for the odds ratio
 *  _crcrnfd		Calculates unadjusted Cornfield confidence limits for the odds ratio
 *  xcnfint		Calculates exact confidence limits for the odds ratio
 *  finda		Support routine for xcnfint
 *  mkdata		Support routine for xcnfint
 */
version 6.0
program drop _all
program define exactcci, rclass
	version 6
	preserve
	local args="`0'"
/*
 *  Break out a, b, c, and d from the argument list.
 */
	gettoken a 0 : 0, parse(" ,")
	gettoken b 0 : 0, parse(" ,")
	gettoken c 0 : 0, parse(" ,")
	gettoken d 0 : 0, parse(" ,")

	confirm integer number `a'
	confirm integer number `b'
	confirm integer number `c'
	confirm integer number `d'

	if `a'<0 | `b'<0 | `c'<0 | `d'<0 { 
		di in red "negative numbers invalid"
		exit 498
	}
/*
 *  Parse the remaining command arguments.
 */
	syntax [, Exact Level(integer $S_level) COL(string) Woolf TB]
/*
 *  Call the standard cci command.
 */
    if `"`woolf'"'~="" | `"`tb'"'~=""{ 
		cci `args'
	}
	else {
		quiet cci `args'
	}
/*
 *  Make sure cci's returns get returned to the user.
 */
 	local ccior=r(or)
	local ccilb=r(lb_or)
	local cciub=r(ub_or)
	local chi2 = r(chi2)
	local p = r(p)
	local p1_ex = r(p1_exact)
	local p_ex = r(p_exact)
	local afe = r(afe)
	local lb_afe = r(lb_afe)
	local ub_afe = r(ub_afe)
	local afp = r(afp)
	ret scalar chi2 = r(chi2)
	ret scalar p = r(p)
	ret scalar p1_exact = r(p1_exact)
	ret scalar p_exact = r(p_exact)
	ret scalar or = r(or)
	ret scalar lb_or = r(lb_or)
	ret scalar ub_or = r(ub_or)
	ret scalar afe = r(afe)
	ret scalar lb_afe = r(lb_afe)
	ret scalar ub_afe = r(ub_afe)
	ret scalar afp = r(afp)
/*
 *  If user specifies Woolf or TB arguments then return here.
 */
    if `"`woolf'"'~="" | `"`tb'"'~="" {exit}

/*
 *  Calculate the Yate's contunity corrected chi2.
 */
	local m1=`a'+`b'
	local m0=`c'+`d'
	local n1=`a'+`c'
	local n0=`b'+`d'
	local T=`n1'+`n0'
	local ae=`n1'*`m1'/`T'
	local be=`n0'*`m1'/`T'
	local ce=`n1'*`m0'/`T'
	local de=`n0'*`m0'/`T'
	local chi_cnt=((abs(`a'*`d'-`b'*`c')-0.5*`T')^2*`T')/(`n1'*`n0'*`m1'*`m0')
	local yatesp = chiprob(1,`chi_cnt')
/*
 *  calulate chi2.
 */
	local xchi2=(abs(`a'-`ae'))^2*(1/`ae'+1/`be'+1/`ce'+1/`de')
	local xp = chiprob(1,`xchi2')
/*
 *  If level is not specified then default to 95% confidence interval.
 */
	if `level'<10 | `level'>99 { 
		local level 95
	}
/*
 * ?? Calculate UNADJUSTED Cornfield confidence limits for the odds ratio.
 */
		_crcrnfd `a' `b' `c' `d' `level'
		local ccilb = r(lb)
		local cciub = r(ub)
		local lb_afe=(`ccilb'-1)/`ccilb'
		local ub_afe=(`cciub'-1)/`cciub'
/*
 *  Calculate ADJUSTED Cornfield confidence limits for the odds ratio.
 */
	_xcrnfd `a' `b' `c' `d' `level'
	local adj_lb = r(lb)
	local adj_ub = r(ub)
	ret scalar adj_lb = r(lb)
	ret scalar adj_ub = r(ub)
	
	local savadlb=`adj_lb'
	local savadub=`adj_ub'
/* 
 *  Calculate starting values for iterative estimate of the
 *  exact confidence interval.
 *
 *  If the lower limit<=0 or the upper limit=. then replace the
 *  offending limit with the corresponding woolf limit.
 */
    local flag_lb=0
    if `adj_lb'<=0 {local flag_lb=1}
    local flag_ub=0
    if `adj_ub' == . {local flag_ub=1}
/*
 *  If called for get and display the exact limits.
 */
    if `"`exact'"'~="" {
/*
 *  Calculate UNADJUSTED Cornfield confidence limits for the odds ratio.
 */
		_crcrnfd `a' `b' `c' `d' `level'
		local unad_lb = r(lb)
		local unad_ub = r(ub)
	
		local savlb=`unad_lb'
		local savub=`unad_ub'
/*
 *  If the lower limit<=0 or the upper limit=. or flag_ub or flag_lb = true
 *  then replace the offending unadjusted limit with the corresponding
 *  woolf limit;  replace the adjusted limit by a multiple of the
 *  unadjusted limit.  This is done to ensure convergence to the exact
 *  limts.
 */
		if `unad_lb'<=0 | `unad_ub'==. | `flag_ub'==1 | `flag_lb'==1 {
			local r=(`a'*`d')/(`b'*`c')
			local iz=invnorm(1-(1-`level'/100)/2)
			local sdlnr=sqrt(1/`a'+1/`b'+1/`c'+1/`d')
			local woolflb = exp(ln(`r')-`iz'*`sdlnr')
			local woolfub = exp(ln(`r')+`iz'*`sdlnr')
			if `unad_lb'<=0 | `flag_lb'==1 {
				local unad_lb=`woolflb'
				local adj_lb=`woolflb'*0.75
			}
			if `unad_ub'==. | `flag_ub'==1 {
				local unad_ub=`woolfub'
				local adj_ub=`woolfub'*1.3
			}
		}
/*
 *  Calculate exact confidence limits for the odds ratio.
 */
		quietly xcnfint `a' `b' `c' `d' `level' `unad_lb' `unad_ub' `adj_lb' `adj_ub' 
		local ex_lb = r(ex_lb)
		local ex_ub = r(ex_ub)
	}
/*
 *  Set up to return data.
 */
	ret scalar ex_lb = r(ex_lb)
	ret scalar ex_ub = r(ex_ub)
/*
 *  Produce the table part of the output.
 */
	 _crc4fld `a' `b' `c' `d' `"`col'"' "" Exposed Unexposed Cases Controls "" yes
/*
 *  Set a text field for part of the output.
 */
	if `ccior'>=1 {
		local hdr " Attr. frac."
		local adjfrlb=(`savadlb'-1)/`savadlb'
		local adjfrub=(`savadub'-1)/`savadub'
		if `"`exact'"'~="" {
			local exfrlb=(`ex_lb'-1)/`ex_lb'
			local exfrub=(`ex_ub'-1)/`ex_ub'
		}
	}
	else {
		local hdr " Prev. frac."
		local adjfrlb=1-`savadub'
		local adjfrub=1-`savadlb'
		if `"`exact'"'~="" {
			local exfrlb=1-`ex_ub'
			local exfrub=1-`ex_lb'
		}
	}

	di in gr _col(18) "|" _col(43) "|"
	di in gr _col(18) "|      Point estimate    |  [`level'% Conf. Interval]" 
	di in gr _col(18) "|" _dup(24) "-" "+" _dup(22) "-"
	di in gr _col(18) "|" _col(43) "|" " Cornfield's limits"

	di in gr "      Odds ratio | " in ye _col(27) %9.0g `ccior' in gr _col(43) "| " in ye %9.0g `savadlb' "   " %9.0g `savadub' in gr `"  Adjusted"'
	di in gr _col(18) "|" _col(43) "| " in ye %9.0g `ccilb' "   " %9.0g `cciub' in gr `"  Unadjusted"'
	if `"`exact'"'~="" {
		di in gr _col(18) "|" _col(43) "|" " Exact limits"
		di in gr _col(18) "|" _col(43) "| " in ye %9.0g `ex_lb' "   " %9.0g `ex_ub'
		di in gr _col(18) "|" _col(43) "|" " Cornfield's limits"
    }

	di in gr `"`hdr' ex. | "' in ye _col(27) %9.0g `afe' in gr _col(43) "| " in ye %9.0g `adjfrlb' "   " %9.0g `adjfrub' in gr `"  Adjusted"'
	di  in gr _col(18) "|" _col(43) "| " in ye %9.0g `lb_afe' "   " %9.0g `ub_afe' in gr `"  Unadjusted"'
	if `"`exact'"'~="" {
		di in gr _col(18) "|" _col(43) "|" " Exact limits"
		di in gr _col(18) "|" _col(43) "| " in ye %9.0g `exfrlb' "   " %9.0g `exfrub'
    }

	di in gr `"`hdr' pop | "' in ye _col(27) %9.0g `afp' in gr _col(43) "|"
		
	di in gr _col(18) "+-----------------------------------------------"
	
	di in gr /*
	*/ _col(19) "           chi2(1) =" in ye %9.2f `xchi2' /*
	*/ in gr "  Pr>chi2 =" in ye %7.4f `xp'
	di in gr /*
	*/"             Yates' adjusted chi2(1) =" in ye %9.2f `chi_cnt' /*
	*/ in gr "  Pr>chi2 =" in ye %7.4f `yatesp'

	if `"`exact'"'~="" {
        di in gr _col(25) "        1-sided Fisher's exact P =" in ye %7.4f `p1_ex'
		di in gr _col(25) "        2-sided Fisher's exact P =" in ye %7.4f `p_ex'
        di in gr _col(25) "2 times 1-sided Fisher's exact P =" in ye %7.4f `p1_ex'*2
	}

end
/*
 *  ======== _xcrnfd =====================================================
 */
program define _xcrnfd, rclass
	version 3.1
	local a `1'
	local b `3' /* sic */
	local c `2' /* sic */
	local d `4'

	tempname iz m1 n1 n2 al alold i

	scalar `iz'=invnorm(1-(1-`5'/100)/2)
	scalar `m1'=`a'+`b'
	scalar `n1'=`a'+`c'
	scalar `n2'=`b'+`d'
	scalar `i' = 0 
	scalar `al'= `a'
	scalar `alold'= .
	while abs(`al'-`alold')>.001 & `al'!=. { 
		scalar `alold' = `al'
		scalar `al'=`a'- 0.5 -`iz'*1/sqrt(1/`al'+1/(`m1'-`al')+/*
			*/ 1/(`n1'-`al')+ 1/(`n2'-`m1'+`al'))
		if `al'==. {
			scalar `i'=`i'+1
			scalar `al'=`a'-`i'
			if (`al'<0 | (`n2'-`m1'+`al')<0) { scalar `al'= . } 
		}
	}
	if `al'==. { scalar `al'= 0 } 
	ret scalar lb = `al'*(`n2'-`m1'+`al')/((`m1'-`al')*(`n1'-`al'))
	scalar `al'= `a'
	scalar `alold'= . 
	scalar `i'= 0 
	while abs(`al'-`alold')>.001 & `al'!=. {
		scalar `alold'= `al'
		scalar `al'=`a'+ 0.5 +`iz'*1/sqrt(1/`al'+1/(`m1'-`al')+/*
			*/ 1/(`n1'-`al')+ 1/(`n2'-`m1'+`al'))
		if `al'==. {
			scalar `i'=`i'+1
			scalar `al'=`a'+`i'
			if (`al'>`n1'|`al'>`m1') { scalar `al' = . } 
		}
	}
	ret scalar ub = `al'*(`n2'-`m1'+`al')/((`m1'-`al')*(`n1'-`al'))
end

/*
 *  ======== _crcrnfd =====================================================
 */
program define _crcrnfd, rclass
	version 3.1
	local a `1'
	local b `3' /* sic */
	local c `2' /* sic */
	local d `4'

	tempname iz m1 n1 n2 al alold i

	scalar `iz'=invnorm(1-(1-`5'/100)/2)
	scalar `m1'=`a'+`b'
	scalar `n1'=`a'+`c'
	scalar `n2'=`b'+`d'
	scalar `i' = 0 
	scalar `al'= `a'
	scalar `alold'= .
	while abs(`al'-`alold')>.001 & `al'!=. { 
		scalar `alold' = `al'
		scalar `al'=`a'-`iz'*1/sqrt(1/`al'+1/(`m1'-`al')+/*
			*/ 1/(`n1'-`al')+ 1/(`n2'-`m1'+`al'))
		if `al'==. {
			scalar `i'=`i'+1
			scalar `al'=`a'-`i'
			if (`al'<0 | (`n2'-`m1'+`al')<0) { scalar `al'= . } 
		}
	}
	if `al'==. { scalar `al'= 0 } 
	ret scalar lb = `al'*(`n2'-`m1'+`al')/((`m1'-`al')*(`n1'-`al'))
	scalar `al'= `a'
	scalar `alold'= . 
	scalar `i'= 0 
	while abs(`al'-`alold')>.001 & `al'!=. {
		scalar `alold'= `al'
		scalar `al'=`a'+`iz'*1/sqrt(1/`al'+1/(`m1'-`al')+/*
			*/ 1/(`n1'-`al')+ 1/(`n2'-`m1'+`al'))
		if `al'==. {
			scalar `i'=`i'+1
			scalar `al'=`a'+`i'
			if (`al'>`n1'|`al'>`m1') { scalar `al' = . } 
		}
	}
	ret scalar ub = `al'*(`n2'-`m1'+`al')/((`m1'-`al')*(`n1'-`al'))
end

/*
 *  ======== xcnfint =====================================================
 */
/*   (xcnfint.do)
 *
 *  Usage:
 *		. xcnfint a b c d level `unadjlb' `unadjub' `adjlb' adjub' 
 *
 *  Calculate the exact level% confidence limits for the contigency table:
 *
 *              exposed|un-exposed
 *             +-------+-------+
 *     cases   |   a   |   b   |
 *             +-------+-------+
 *     controls|   c   |   d   |
 *             +-------+-------+
 *
 *  Returns:
 *		Printed output
 *		r(ex_lb) - exact lower bound
 *		r(ex_ub) - exact upper bound
 */
program define xcnfint, rclass  /* a b c d level */
    version 6
/*
 *  Get the arguments.
 */
	gettoken a 0 : 0, parse(" ,")
	gettoken b 0 : 0, parse(" ,")
	gettoken c 0 : 0, parse(" ,")
	gettoken d 0 : 0, parse(" ,")
	gettoken level 0 : 0, parse(" ,")
	gettoken unadjlb 0 : 0, parse(" ,")
	gettoken unadjub 0 : 0, parse(" ,")
	gettoken adjlb 0 : 0, parse(" ,")
	gettoken adjub 0 : 0, parse(" ,")

	confirm integer number `a'
	confirm integer number `b'
	confirm integer number `c'
	confirm integer number `d'
	confirm number `level'
/*
 *  Calculate the marginals.
 */
    local m1=`a'+`b'
	local m0=`c'+`d'
	local n1=`a'+`c'
	local n0=`b'+`d'
/*
 *  Calculate alpha/2 (needed in a calculation later).
 */
	local alpha=(100-`level')/100
	local ao2=`alpha'/2
/* 
 * ---- Calculate the exact lower limit -----------------------------------
 */
/*
 *  Get the un-adjusted (Stata) Cornfield limits.
 */
/*
 *  Extract the un-adjusted lower limit.
 */
    local psi_im1=`unadjlb'
/*
 *  Get the adjusted Cornfield limits.  psi(i) is the lower limit.
 */
/*
 *  Extract the adjusted lower limit as a starting point.
 */
	local psi_i=`adjlb'
/*
 *  If the limits are close enough then we are done and the Cornfield
 *  is the exact lower limit.
 */
	if `psi_i' ~= . & `psi_im1' ~= . & abs(`psi_i'-`psi_im1') >= 0.005*`psi_i' {
/*
 *  Otherwise, call mkdata0 to create a table of all possible a's 
 *  (given the marginals)
 */
		mkdata0 `n1' `n0' `m1' `m0'
/*
 *  Call mkdata to derive  the cumulative distribution of a conditioned
 *  on the marginal totals. (un-adjusted)
 */
		mkdata `psi_im1'
/*
 *  Scan the in-memory data set for the row where a is
 *  equal to the a in the 2x2 table.  i points to that row.
 */
		finda `a'
		local aindex=r(index)
/*
 *  Get p(i-1).
 */
		local p_im1=1
		if `aindex' ~= 1 {local p_im1=1-cum[`aindex'-1]}
/*
 *  Call mkdata again using psi_i.
 */
		mkdata `psi_i'
/*
 *  Get p(i).
 */
		local p_i=1
		if `aindex' ~= 1 {local p_i=1-cum[`aindex'-1]} 
/*
 *  Loop ...
 */
 		local lcount=0
		while abs(`psi_i'-`psi_im1') >= 0.005*`psi_i' {
/*
 *  Check the number of iterations.  Quit if too many.
 */
			local lcount=`lcount'+1
			if `lcount'==100 {
			    display "Too many iterations:lower"
				exit 19
			}
/*
 *  Calculate the slope.
 */
			local slope=(`p_i'-`p_im1')/(`psi_i'-`psi_im1')
/*
 *  Calculate psi(i+1)
 */
			local psi_ip1=(`ao2'+`slope'*`psi_i'-`p_i')/`slope'
/*
 *  Calculate the cumulative distribution using psi_ip1.
 */
			mkdata `psi_ip1'
/*
 *  Get p(i+1)
 */
 			local p_ip1=1
			if `aindex' ~= 1 {local p_ip1=1-cum[`aindex'-1]} 
/*
 *  Reset for next iteration.
 *	 psi(i-1)=psi(i)
 *	 p(i-1=p(i)
 *	 psi(i)=psi(i+1)
 *	 p(i)=p(i+1)
 */
			local psi_im1=`psi_i'
			local p_im1=`p_i'
			local psi_i=`psi_ip1'
			local p_i=`p_ip1'
		}
	}
	local ex_lb=`psi_i'
	return scalar ex_lb=`psi_i'
	
/* 
 * ---- Calculate the exact upper limit -----------------------------------
 */
/*
 *  Get the un-adjusted (Stata) Cornfield limits.  psi(i-1) is the upper limit.
 */
	local psi_im1=`unadjub'
/*
 *  Get the adjusted Cornfield limits.  psi(i) is the upper limit.
 */
	local psi_i=`adjub'
/*
 *  If the limits are close enough then we are done and the adjusted Cornfield
 *  is the exact upper limit.
 */
	if `psi_i' ~= . & `psi_im1' ~= . & abs(`psi_i'-`psi_im1') >= 0.005*`psi_i' {

		mkdata `psi_im1'
		local p_im1=cum[`aindex']

		mkdata `psi_i'
		local p_i=cum[`aindex']

		local lcount=0
		while abs(`psi_i'-`psi_im1') >= 0.005*`psi_i' {
/*
 *  Check the number of iterations.  Quit if too many.
 */
			local lcount=`lcount'+1
			if `lcount'==100 {
			    display "Too many iterations:upper"
				exit 29
			}

			local slope=(`p_i'-`p_im1')/(`psi_i'-`psi_im1')
			local psi_ip1=(`ao2'+`slope'*`psi_i'-`p_i')/`slope'

			mkdata `psi_ip1'
			
			local p_ip1=cum[`aindex']

			local psi_im1=`psi_i'
			local p_im1=`p_i'
			local psi_i=`psi_ip1'
			local p_i=`p_ip1'
		}
	}
	local ex_ub=`psi_i'
	return scalar ex_ub=`psi_i'
end

/*
 *  ======== finda =====================================================
 */
program define finda, rclass
    version 6
/*
 *  Scan the in-memory data set for the row where a is
 *  equal to the argument to this routine.  Return the index.
 */
    local a=`1'
	local index=1
	while a[`index'] ~= `a' {
		local index=`index'+1
		if `index'>_N {
			display "index out of range"
			exit
		}
	}
	return scalar index = `index'
end

/*
 *  ======== mkdata =====================================================
 */
program define mkdata  /* psi  */
    version 6

	local psi= `1'	/* psi  */
/*
 *  Let g1ao =g1(a,0).  We need to be careful about our
 *  arithmetic to avoid machine overflows in calculating
 *  f(a|psi).
 */
    gen g1a0=a*log(`psi')+fac
    sort g1a0
/*
 *  k1= - max g1(a,0)
 */
    gen k1 = - g1a0[_N]
/*
 *  g2a1 = g2(a,1)
 */
    gen g2a1 = exp(k1 + g1a0)
    gen k2 = sum(g2a1)
/*
 *  Calculate f(a|psi)
 */
    gen fa = g2a1/k2[_N]
/*
 *  Calculate cumulative f(a,K).
 */
    sort a
    replace cum=sum(fa)
    keep a cum fac
/*
 *  We now have 2 relevant variables in the data set:  a and cum.  a contains
 *  all possible values of a and cum contains the cumulative distribution of
 *  a conditioned on the marginal totals.
 */
end
/*
 *  ======== mkdata0 =====================================================
 */
program define mkdata0  /* exposed unexposed cases controls */
    version 6
/* 
 *  Set up for data iterations.
 */
    drop _all
	local n1= `1'   /* exposed */
	local n0= `2'   /* unexposed */
	local m1= `3'   /* cases */
	local m0= `4'   /* controls */
/*
 *  Calculate the range of a.  The possible range of a is
 *  max(0,m1-n0)...min(n1,m1).
 */
    local alow=max(0,`m1'-`n0')
	local ahigh=min(`n1',`m1')
/*
 *  Calculate the number of values in the range of a.
 */
	local nobs=`ahigh'-`alow'+1
/*
 *  Set data set to that size...
 */
	set obs `nobs'
/*
 *  Put our fixed values in the data set.  We can remove them later.
 */
	gen n1=`n1'
	gen n0= `n0'
	gen m1= `m1'
	gen m0= `m0'
/*
 *  Put values of a in the table
 */
	gen a=`alow'
	quietly replace a=a[_n-1]+1 if _n ~=1
/*
 *  Define by log factorial terms.
 */
    gen fac = -lnfact(a)-lnfact(n1-a)-lnfact(m1-a)-lnfact(n0-m1+a)
    sort fac
/*
 *  Initialize cum
 */
	generate cum=.
end