program define choi_lr_testi, rclass
version 14

* This is the immediate version of choi_lr_test

* Calculate statistics for 2x2 tables contitioned on the marginal
* success rate. These include the
*    Conditional maximum likelihood estimate of odds ratio (OR).
*    Conditional likelihood ratio for the null hypothesis that OR = 1
*    LiChoi's likelihood ratio chi-squared test for 2x2 tables
*    The 1/6.8259358 likelihood support interval (LSI) for the OR
*    The 1/k LSI for the OR using some other value of k

* See Choi et al. 2015 PLoS ONE. 10(4): e0121263.

* INPUT:
*    a = number of exposed cases
*    b = number of unexposed cases
*    c = number of exposed controls
*    d = number of unexposed controls
*    k = a value used to calculate the 1/k LSI for the OR.
*        The 1/6.8259358 LSI is calculated by 
*        default, which equals the frequentist 95% confidence 
*        interval for normally distrubuted statistics.

* OUTPUT:
*    or_cond  = Maximum likelihood estimate (MLE) of the odds ratio for
*               exposure among cases vs. controls conditioned
*               on the total number of exposed subjects.
*    clr      = Conditional likelihood ratio for the null hypothesis that OR = 1 
*    chi2_clr = Choi's LR chi-squared statistic. This is equation 9 of Choi et al. 
*    p_choi   = P value associated with the null hypothesis that OR = 1 
*    [lr6pt8lsi_lb, lr6pt8lsi_ub]
*             = likelihood ratio support interval for OR with k = 6.8259358.
*    [lrklsi_lb, lrklsi_ub]
*             = the likelihood ratio support interval for the user's value of k

* SOME OTHER VARIABLES:
*    psi     = the log of an odds ratio estimate.  Not necessarily the MLE
*    psi_mle = the conditional MLE of psi
*    y1      = number of successes on treatment 1
*    y2      = number of successes on treatment 2
*    n1      = number of subjects on treatment 1
*    n2      = number of subjects on treatment 2
*    yplus   = y1 + y2 = the marginal (total) number of successes
*    max_lnf = probability mass function of y1|yplus at psi_mle
*  
* NOTATION
*    Choi et al. discuss a 2x2 table with n1 and n2 subjects on treatments 1 and
*    and 2, respectively with y1 and y2 successes observed on these two treatments.
*    This program is written using their notation. Hence,
*    a + b = n1 = total number of cases = total number of subjects on treatment 1
*    c + d = n2 = total number of controls = total number of subjects on treatment 2
*    a = y1 = number of exposed cases = number of successes on treatment 1
*    c = y2 = number of exposed controls = number of successes on treatment 2

* CALLING SYNTAX
*    choi_lr_testi #a #b #c #d, k(#k) 

    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 { 
        display as error "negative numbers invalid"
	exit 498
    }
	
* Identify tables with too many zeros

    if `a'*`d'==0 & `b'*`c'==0 {
        display as error "Too many empty cells to estimate OR"
        exit 498
    }
 
    syntax [, K(numlist max=1) Woolf TB COrnfield Exact  Level(cilevel) * ]

    local levopt="level(`level')"

    preserve
       
* Confirm that needed scalars do not exist in the calling program. If they 
* do, save them with another name so that they can be restored when the 
* program ends.

    foreach v in n1 n2 y1 y2 dummy psi p_choi chi2_clr clr max_lnf psi_mle k {
        capture confirm scalar `v'
        if _rc == 0  {
            scalar ______TMP`v' = `v'
         }
    }      
    scalar y1 = `a'
    scalar y2 = `c'
    scalar n2 = `c' + `d'
    scalar n1 = `a' + `b'
*
* Do conventional case-control analysis
*
    cci `a' `b' `c' `d', `woolf' `tb' `cornfield' `exact' `levopt' `options'
    local log_or_unconditional = log(r(or))
    return scalar p        = r(p)        //two-sided p-value from the 2x2 chi^2 statistic without continuity correction
    return scalar p1_exact = r(p1_exact) //chi-squared or one-sided exact significance
    return scalar p_exact  = r(p_exact)  //two-sided exact significance
    return scalar or       = r(or)       //Unconditioned MLE of the odds ratio
    return scalar lb_or    = r(lb_or)    //lower bound of CI for or
    return scalar ub_or    = r(ub_or)    //upper bound of CI for or
    return scalar afe      = r(afe)      //attributable (prev.) fraction among exposed
    return scalar lb_afe   = r(lb_afe)   //lower bound of CI for afe
    return scalar ub_afe   = r(ub_afe)   //upper bound of CI for afe
    return scalar afp      = r(afp)      //attributable fraction for the population
    return scalar chi2     = r(chi2)     //2x2 chi^2 statistic without continuity correction
    display " "
    clear
    quietly set obs 1
    gen  y1 = y1
    gen  y2 = y2
    

* Calculate the MLE of psi = log(or_cond)

* See http://www.stata.com/features/overview/maximum-likelihood-estimation/
* for an introduction to Stata's ml program. The example given in this URL
* was modified to maximize the likelihood function given in equation 7 of 
* Choi et al. 2015.

* Identify if estimated OR = 0, infinity, of some positive number
    if `a'*`d'>0 & `b'*`c'>0  { // all cells are positive

        quietly ml model lf choi_lr_hypergeom  (y1=)
        quietly ml init `log_or_unconditional', copy
        quietly ml maximize
  
* The preceding command returns _b[_cons] which is the MLE of psi, 
* and e(ll) which is the maximum value of the log likelihood function

        scalar psi_mle = _b[_cons]
        scalar max_lnf = e(ll)

* Calculate the probability of the observed results under the null hypothesis 
* that psi = 0.

        choi_lr_hyperg_prob n1 n2 y1 y2 0
  
* r(lnf) equals the log of the probability mass function when psi = 0

        scalar clr = exp(max_lnf - r(lnf))
        scalar chi2_clr = -2*(r(lnf)-max_lnf)
        scalar p_choi = chi2tail(1, chi2_clr)
        display as result " "
        display as result "Likelihood ratio statistics conditioned on the marginal exposure rate"
        display as result " "
        display as result "Maximum likelihood estimate (MLE) of odds ratio      = "exp(psi_mle)
        display as result "Likelihood ratio under the null hypothesis           = "clr
        display as result "Likelihood ratio chi-squared statistic               = "chi2_clr 
        display as result "P value for the null hypothesis that the OR equals 1 = "p_choi

        return scalar or_cond   = exp(_b[_cons])
        return scalar clr   = clr
        return scalar chi2_clr = chi2_clr
        return scalar p_choi    =  p_choi
        global k_global = 6.8259358
    
* Calculate the 1/6.8259358 LSI for the OR

        choi_lr_support_interval 
        return scalar lr6pt8lsi_lb =LRsupport_lb
        return scalar lr6pt8lsi_ub =LRsupport_ub

* Calculate the 1/k LSI for the OR
        if "`k'"~="" {
	    global k_global = `k'
	    choi_lr_support_interval
	    return scalar lrklsi_lb =LRsupport_lb
	    return scalar lrklsi_ub =LRsupport_ub
        }
    }
    else if `a'*`d'>0 & `b'*`c'==0  { // Estimated OR = infinity
        
* Equation (7) of Choi et al 2015 approaches 1 as psi approaches infinity
* A proof of this is given in Equn7whenAcellEquals0.pdf, which is stored in
* my dropbox folder \OFFICE_desktop\ChoiFisherTest\Stata Program\ado\withZeros
    
        scalar max_lnf = 0
        choi_lr_hyperg_prob n1 n2 y1 y2 0
        scalar clr = exp(max_lnf - r(lnf))
        scalar chi2_clr = -2*(r(lnf)-max_lnf)
        scalar p_choi = chi2tail(1, chi2_clr)
        display as result " "
        display as result "Likelihood ratio statistics conditioned on the marginal exposure rate"
        display as result " "
        display as result "Maximum likelihood estimate (MLE) of odds ratio      = infinity"
        display as result "Likelihood ratio under the null hypothesis           = "clr
        display as result "Likelihood ratio chi-squared statistic               = "chi2_clr 
        display as result "P value for the null hypothesis that the OR equals 1 = "p_choi
        
        return scalar or_cond   = .
        return scalar clr   = clr
        return scalar chi2_clr = chi2_clr
        return scalar p_choi    =  p_choi
        global k_global = 6.8259358
    
* Calculate the 1/6.8259358 LSI for the OR
        choi_lr_support_interval_big_or
        return scalar lr6pt8lsi_lb =LRsupport_lb
        return scalar lr6pt8lsi_ub =.

* To comfirm this lower bound execute the next three lines
*           local debug_psi = log(LRsupport_lb)		//debug
*           choi_lr_hyperg_prob n1 n2 y1 y2 `debug_psi' //debug
*           di "LR at lb = " 1/exp(r(lnf)) 		//debug

* Calculate the 1/k LSI for the OR
        if "`k'"~="" {
            global k_global = `k'
            choi_lr_support_interval_big_or
            return scalar lrklsi_lb =LRsupport_lb
            return scalar lrklsi_ub =.
* To comfirm this lower bound execute the next three lines    
*               local debug_psi = log(LRsupport_lb)		//debug
*               choi_lr_hyperg_prob n1 n2 y1 y2 `debug_psi' 	//debug
*               di "LR at lb = " 1/exp(r(lnf)) 			//debug
        }

    }
    else if `a'*`d'==0 & `b'*`c'>0  { // Estimated OR = 0
    
    * Equation (7) of Choi et al 2015 approaches 1 as psi approaches minus infinity
    * A proof of this is given in Equn7whenAcellEquals0.pdf, which is stored in
    * my dropbox folder \OFFICE_desktop\ChoiFisherTest\Stata Program\ado\withZeros
    
        scalar max_lnf = 0
        choi_lr_hyperg_prob n1 n2 y1 y2 0
        scalar clr = exp(max_lnf - r(lnf))
        scalar chi2_clr = -2*(r(lnf)-max_lnf)
        scalar p_choi = chi2tail(1, chi2_clr)
        display as result " "
        display as result "Likelihood ratio statistics conditioned on the marginal exposure rate"
        display as result " "
        display as result "Maximum likelihood estimate (MLE) of odds ratio      = 0"
        display as result "Likelihood ratio under the null hypothesis           = "clr
        display as result "Likelihood ratio chi-squared statistic               = "chi2_clr 
        display as result "P value for the null hypothesis that the OR equals 1 = "p_choi
        
        return scalar or_cond   = 0
        return scalar clr   = clr
        return scalar chi2_clr = chi2_clr
        return scalar p_choi    =  p_choi
        global k_global = 6.8259358
    
* Calculate the 1/6.8259358 LSI for the OR
        choi_lr_support_interval_or_eq_0 
        return scalar lr6pt8lsi_lb =0
        return scalar lr6pt8lsi_ub =LRsupport_ub
        
* To comfirm this upper bound execute the next three lines    
*           local debug_psi = log(LRsupport_ub)		//debug
*           choi_lr_hyperg_prob n1 n2 y1 y2 `debug_psi' //debug
*           di "LR at ub = " 1/exp(r(lnf)) 		//debug

* Calculate the 1/k LSI for the OR
        if "`k'"~="" {
            global k_global = `k'
            choi_lr_support_interval_or_eq_0 
            return scalar lrklsi_lb =0
            return scalar lrklsi_ub =LRsupport_ub
* To comfirm this upper bound execute the next three lines    
*               local debug_psi = log(LRsupport_ub)		//debug
*               choi_lr_hyperg_prob n1 n2 y1 y2 `debug_psi' 	//debug
*               di "LR at ub = " 1/exp(r(lnf)) 			//debug
        }
    }
    
*   Restore any scalar values from the calling program that have been changed

    foreach v in n1 n2 y1 y2 dummy psi p_choi chi2_clr clr max_lnf psi_mle k {
        capture confirm scalar ______TMP`v'
        if _rc == 0  {
            scalar `v' = ______TMP`v'
            scalar drop ______TMP`v'
        }
    }

end // ********************************************************