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

Question

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.

Answer 1

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).

Answer 2

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.