# 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?

• How are you varying the data? It's conceivable you can get the information you need with very little work, because once you know the p-value is less than 0.05 for a certain $2\times 2$ table you can determine it is less than 0.05 for many other related $2\times 2$ tables that exhibit larger odds ratios for the same study size. BTW, you intended to write "< 0.05," not "> 0.05," I hope. – whuber Jan 26 at 23:09
• The counts are samples from a beta distribution. I like the idea of ruling out related 2x2 tables, I'll check that out. You're right, I had the wrong sign on the p-value – Kyle Heuton Jan 27 at 3:10
• I do not understand what you mean because values from a Beta distribution by definition are numbers strictly between $0$ and $1,$ so they cannot possibly be counts. – whuber Jan 27 at 14:44
• Samples from the beta distribution determine the event frequency for each simulated study- it's then scaled by the study size – Kyle Heuton Jan 28 at 20:59

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) ) }  • Wouldn’t the$\chi^2\$ test not pin down which simulations have significant p-values? – Dave Jan 26 at 22:38
• Hat tip on the needed clarification, Dave. I've specified that each of the 10K comparisons needs its own 2x2 ChiSquare test. – Robert Alan Greevy Jr PhD Jan 26 at 22:57
• This worked great, and brought the time down from 9 minutes to 40 seconds or so. I wouldn't mind going a little faster. Would I have to calculate the smallest y for every unique x, or is there a better way? – Kyle Heuton Jan 27 at 3:13
• That's an exciting improvement in time. My answer to your question on improving further is too long for a comment, so I'll expand my answer above. – Robert Alan Greevy Jr PhD Jan 27 at 17:07