fisher_exact_test_sensitivity_analysis.py

import operator as op
from functools import reduce
from decimal import Decimal

def FisherSens(totalN, treatedN, totalSuccesses, treatedSuccesses, Gammas):
	## 'FisherSens' 
	## Python code by Abhinav Maurya
	## Original code in R by Devin Caughey at https://rdrr.io/cran/rbounds/src/R/FisherSens.R
	##
	## This function performs a sensitivity analysis for Fisher's Exact Test
	## for count (binary) data. It is derived from Section 4.4 of Paul
	## Rosenbaum’s ``Observational Studies" (2nd Ed., 2002). The test is
	## one-sided; that is, the alternative hypothesis to the null is
	## that the number of "successes" (however defined) greater in among
	## the "treated" (however defined).
	##
	## It takes five arguments:
	##   'totalN': total number of observations
	##   'treatedN': number of observations that received treatment
	##   'totalSuccesses': total number of "successes"--i.e., 1's
	##   'treatedSuccesses': number of successes among the treated
	##   'Gammas': a vector Gammas (bounds on the differential odds of 
	##  treatment) at which to test the significance of the results.
	##
	## It returns a matrix of Gammas and upper and lower bounds on the exact 
	## p-value for Fisher's test.

	n = totalN
	m = treatedN
	c_plus = totalSuccesses
	a = treatedSuccesses

	output = []

	for g in range(len(Gammas)):
		p_plus = float(Upsilon(n = n, m = m, c_plus = c_plus, a = a, Gamma = Gammas[g]))
		p_minus = float(Upsilon(n = n, m = m, c_plus = c_plus, a = a, Gamma = (1 / Gammas[g])))
		output.append([Gammas[g], p_plus, p_minus])

	return output

def Upsilon(n, m, c_plus, a, Gamma):
	## Prob A >= a 
	Gamma = Decimal(Gamma)

	numer = Decimal(0.0)
	for k in range(max(a, (m + c_plus - n)), min(m, c_plus)+1):
		numer += Choose(c_plus, k) * Choose((n-c_plus), m - k) * (Gamma ** k)

	denom = Decimal(0.0)
	for k in range(max(0, (m + c_plus - n)), min(m, c_plus)+1):
		denom += Choose(c_plus, k) * Choose((n - c_plus), (m - k)) * (Gamma ** k)

	return numer/denom

def Choose(n, r):
    r = min(r, n-r)
    if r == 0: return 1
    numer = Decimal(reduce(op.mul, range(n, n-r, -1)))
    denom = Decimal(reduce(op.mul, range(1, r+1)))
    return numer/denom

print(FisherSens(3000,2000,200,200,[1,2,3,4,5,6,7,8,9,10]))