Bernoulli distribution p-value shortcut I am trying to evaluate the impact of a hypothetical treatment that has a binary outcome (adverse health event or no event). I have ten thousand of simulated "studies". I want to quickly identify which studies have a significant p-value (< 0.05).
My data looks like this: for each simulated study I have 2 numbers, the number of adverse health events in the control population, and the number of adverse events in the treatment population. Currently, to find the p-value, I construct an X and y matrix for each study, where X consists of an intercept term and a binary treatment feature, and y is 0 for no-event and 1 for an event, and I fit a binomial GLM to each model. This is slow, as there are thousands of simulated studies. Is there a faster way to identify studies with significant p-values? Can it be ascertained directly from the number of events?
My current approach:
import numpy as np
import statsmodels.api as sm
# The study consists of this many patients 
# in each the treatment and control arms
study_size = 500

def get_pvalue(inputs):
    control_events, treatment_events = tuple(inputs)
    treatment = [0]*study_size + [1]*study_size
    intercept = [1]*study_size*2
    x_vars = np.array([treatment, intercept]).transpose()
    y_var = np.array([1]*control_events + [0]*(study_size-control_events) + [1]*treatment_events + [0]*(study_size-treatment_events))
    fit_model = sm.GLM(y_var, x_vars, family=sm.families.Binomial()).fit()
    return fit_model.pvalues[0] 

# these are 10000 length arrays, each is a sample number of events
control_event_samples = [31, 24, 53, ... ]
treatment_event_samples = [26, 25, 48 ...]

p_value_array = np.array([arm_1_event_draws, arm_2_event_draws])
get_pvalue = get_pvalue_func(study_size)
p_values = np.apply_along_axis(get_pvalue, axis=0, arr=p_value_array)

It's not a bad approach, but it is slower than I like, and it seems inefficient to have to construct such large arrays for every sample. Is there an obvious shortcut I am missing?
 A: Fitting 10K GLMs will be slow. A 2x2 Chi-squared test for each comparison of two event counts will be faster and will work well enough provided your event counts don't get too extreme, e.g. your event counts range from [5, 495] for your sample size of 500. You can work with the event counts directly so you won't need to create the vectors of zeroes and ones. 
For info on executing the ChiSquare tests, see  https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html
There are other possible ways to speed things up more. For example, grouping identical event count comparisons together so you only calculate the p-value for each comparison once; or predetermining the smallest event count y needed to achieve p<0.05 for event count pairs (x, y) where y >= x. However, those tricks are leveraging the fixed sample sizes of your set up, and you may want to make those more flexible later on. I suspect simply changing to the ChiSquare test may make your code fast enough for your needs.
In response to comment below, here are some more thoughts on speeding up the code further.
For a fixed n=500 for each of the two groups, you could leverage clusters of (y-x) in checking for p < 0.05. For example, for x in [172, 297], (y-x)=32. And for x >= 495, no y > x exists that gives p < 0.05. These clusters would help if writing a bunch of if statements. However, even running through just 28 cleverly constructed if statements might not be the fastest approach.
I can see creating a dictionary of the (x, y) for all x in [0, 500] as keys. I suspect that would run faster than even a cleverly reduced set of if statements.
Below is the R code I used to get the (x, y). I used a continuity correction in the ChiSquare test. You might get a few different y when not using the continuity correction, but those differences are small and it is debatable as to which should be preferred.
n <- 500
for(x in 0:(n-1)){
  p <- 1 # initialize p-value
  y <- x
  while(p > 0.05 & y < n){
    y <- y + 1
    z <- matrix( c(x, n-x, y, n-y), ncol=2)
    p <- chisq.test(z, correct=T)$p.value
  }
  d <- y-x
  if(p > 0.05){ y <- -9; d <- -9 ; } # forgive the magic number -9 for NA
  #print( round( c(x, y, d, p), 3 ) )
  print( c(x, y, d) )
}

