program define choi_lr_hyperg_prob, rclass version 14 * This program calculates the probability mass function defined by * equation 7 of Choi et al. 2015. * This is the probability of observing y1 successes * given n1 and n2 patients on treatments 1 and 2, respectively, and * is conditioned on the marginal number of successes yplus = y1 + y2. * Input: * n1 = patients on treatment 1 * n2 = patients on treatment 2 * y1 = successes on treatment 1 * y2 = successes on treatment 2 * psi = log odds ratio for success in Rx 1 vs Rx 2 * Output: * lnf = log(f), where f = pmf = equation 7 in Choi et al. 2015 * Equation 7 is of the form f = M/C where M and C can both overflow * for large tables. We manipulate this formula to avoid overflow errors. * Let k = log(M). Then log(f) = k - log(C) = -log(C*exp(-k)) where * C*exp(-k) can be written in such a way to avoid overflow errors * when n1 and/or n2 are large. For details see CalculatingEq7.pdf in * \\biostat1158a\Users\dupontwd.VANDERBILT\Documents\Dropbox\OFFICE_desktop * \ChoiFisherTest\Stata Program\ado\withZeros args n1 n2 y1 y2 psi local yplus = `y1' + `y2' local start = max(0,`yplus' - `n2') local end = min(`n1', `yplus') * Calculate k = log(M), where f =M/C in equation 7 local k = lnfactorial(`n1') - lnfactorial(`y1') - lnfactorial(`n1'-`y1') /// + lnfactorial(`n2') - lnfactorial(`yplus'-`y1') /// - lnfactorial(`n2' - `yplus' +`y1') + `psi'*`y1' local emkc = 0 forvalues u = `start'/`end' { local emkc = `emkc' + exp(lnfactorial(`n1') -lnfactorial(`u') /// -lnfactorial(`n1'-`u') + lnfactorial(`n2') /// -lnfactorial(`yplus'-`u')-lnfactorial(`n2'-`yplus'+`u') /// + `psi'*`u' - `k') } * emkc = C*exp(-k) where C is defined in equation 7 local lnf = -log(`emkc') // = 0 -log(C*exp(-k)) = log(f) return scalar lnf = `lnf' end // **************************************************************