Comparing Bernoulli means across subpopulations in which the number of observed successes may be zero

I've got a binary characteristic and a population $S$ with size $n$ and $P[X] = p$ such that $p$ may be small and $n$ is extremely large. Within this population are subpopulations of various sizes $S_0, S_1, \dots, S_k \subset S$.

I'd like to be able to select each subpopulation in which $p_i < p$ with some concept of statistical significance.

My first inclination is to observe that the standard error on each $p_i$ is $SE_i = \sqrt{\frac{\hat{p_i}(1-\hat{p_i})}{n}}$ and to compare upper bounds on confidence intervals. $\{S_i \; | \; \hat{p_i} + 3 \cdot SE_i < p\}$, for example. But when $\hat{p_i} = 0$, then $SE_i = 0$, and this upper bound is 0 even for the smallest subpopulations (like those where $n_i = 1$).

Is there any way to express uncertainty in $p_i$ when $\hat{p_i} = 0$? Maybe through use of $p$ as a prior?

Edit: It looks like the Jeffreys interval as described in Brown et al. is about what I'm after, though I'm not as-of-yet sure how to apply it.

• You could also try a zero inflated poisson or negative binomial distribution. – probabilityislogic Apr 5 '14 at 23:04
• I really don't understand what "with some concept of statistical significance." is meant to suggest here. It doesn't seem to fit with my understanding of the phrase "statistical significance" – Glen_b Apr 7 '14 at 0:08
• Though I want to compare $p_i$s to my $p$, all I have available to me are $\hat{p_i}$s. I'm looking for a way to take into account the statistical likelihood that a given, true $p_i$ is less than $p$. – BeyondTheZero Apr 7 '14 at 18:42

I coded a small simulation of a Bayesian solution, see below for an explanation. I hope you don't mind Python:

import numpy as np
#beta distribution
from scipy.stats import beta

N_populations = 20
true_p = 0.05*np.random.rand(N_populations) #max p_i = 0.05, and all uniform over [0,0.05]
N_subpopulation_sizes = np.random.poisson(500, N_populations) #S_i have different population sizes. This can be anything, really.

#generate some fake observations
observed_data = []
for i in range(N_populations):
observed_data.append(np.random.binomial(1, true_p[i], size=N_subpopulation_sizes[i]))

#estimate of global p,
global_success = sum(map( sum, observed_data ))
global_total = N_subpopulation_sizes.sum()

print "Total successes observed: ", global_success
print "Total observed: ", global_total
print "Estimate of global p: ", 1.*global_success/global_total
print

#cool, so which subpopulations are *less* than the global average?
global_rv = beta.rvs(1 + global_success, 1 + global_total - global_success , size=10000)
statistical_threshold = 0.05

for i in range(N_populations):
sp_success = observed_data[i].sum()
sp_total = N_subpopulation_sizes[i]
subpopulation_rv = beta.rvs(1 + sp_success, 1 + sp_total - sp_success, size=10000 )

prob = (global_rv > subpopulation_rv).mean()
if prob > 1 - statistical_threshold:
print "subpopulation %d statistically smaller than global population. %.3f"%(i,prob)


example output:

subpopulation 2 statistically smaller than global population. sig=0.972
subpopulation 4 statistically smaller than global population. sig=1.000
subpopulation 9 statistically smaller than global population. sig=0.990
subpopulation 10 statistically smaller than global population. sig=1.000
subpopulation 14 statistically smaller than global population. sig=0.979
subpopulation 15 statistically smaller than global population. sig=0.987


Explantation

Essentially is it comparing the posterior of the global $p$ with the posteriors of each $p_i$ by calculating $P(p > p_i)$ and considering it significant if this probability is greater than 0.95.

This solves your small sample problem (i.e. it will express uncertainty as the posterior will be very spread out) and the zero-estimate problem (i.e. even if $\hat{p} = 0$, the posterior of $p_i$ will still put value on non-zero possibilities).

• This looks to me like the same method that the Jeffreys interval uses, but with a prior that takes into account information about the probability of success for the global population. Thanks for that. Would you have any guidance as to how strongly to weight that prior? Here, I'm assuming, that's what the size=10000 parameter is doing. – BeyondTheZero Apr 7 '14 at 18:50
• All my priors are Uniform(0,1) (i.e. Beta(1,1) ). The 10000 is the number of draws from the posterior I perform, so I can do monte carlo estimates of the $P(p>p_i)$ – Cam.Davidson.Pilon Apr 8 '14 at 1:45
• Just as an addendum, obscureanalytics.com/2012/07/04/… has some background on setting a prior that could be used here based on observations about $p$ from the global population. – BeyondTheZero Apr 8 '14 at 18:31

There are several methods to obtain confidence bounds (or intervals) for a binomial parameter. They are described for example in Ch. 1.4 in "Categorical Data Analysis" by Alan Agresti. As the author explains, the performance of the asyptotic methods (Wald, Score or Likelihood Ratio) for a small number of successes (such as zero) may be questionable. Perhaps exact (Ch. 1.4.4) methods (Clopper-Pearson most famously) that rely on the binomial distribution are worth a try. They suffer from discreteness problems and may be extremely conservative. The refined method by Blaker fares better in this respect. It is implemented in R for example in the BlakerCI package.

• Those two techniques might both end up working out -- I had no idea there were so many methods for calculating confidence intervals. In particular, thanks for the pointer to BlakerCI. Though an ideal method might somehow use the greater population's background probability $p$, it's good to have solutions I know I won't screw up with crummy choices for beta distribution parameters. – BeyondTheZero Apr 6 '14 at 0:32