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
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 = *study_size + *study_size intercept = *study_size*2 x_vars = np.array([treatment, intercept]).transpose() y_var = np.array(*control_events + *(study_size-control_events) + *treatment_events + *(study_size-treatment_events)) fit_model = sm.GLM(y_var, x_vars, family=sm.families.Binomial()).fit() return fit_model.pvalues # 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?