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 // ********************************************************