Confidence interval for the sum of 2 binomially distributed variables

Question

$P_1$ and $P_2$ are uncorrelated, binomially distributed variables with success probabilities $p_1 \neq p_2$. Say I measure:

• $k_1 = 9$ successes out of $n_1 = 10$ trials for $P_1$ and
• $k_2 = 1000$ successes out of $n_2 = 10000$ trials for $P_2$.

I calculate confidence intervals for the estimates of $P_1,P_2$ using Clopper-Pearson with 95% confidence. Then I can put 2 datapoints in my paper:

1. (estimate for $p_1$) $= k_1/n_1 = 90.0\%^{+9.7\%}_{-34.5\%}$
2. (estimate for $p_2$) $= k_2/n_2 = 10.0\%^{+0.6\%}_{-0.6\%}$

Now the question: I also need to show the quantity (estimate for $\mathrm{mean}(p_1+p_2)$) $= (k_1/n_1 + k_2/n_2) / 2 = 50\%$. I understand that $P_1+P_2$ follows the Poisson-binomial distribution with 10 Bernoulli trials having success probability $p_1$ and 10000 Bernoulli trials having success probability $p_2$. How do I calculate a 95% confidence interval for my datapoint $\mathrm{mean}(p_1+p_2)$) $ = 50\%$?

Background info:

• The Clopper-Pearson confidence interval is summed up nicely on page 3 of Måns Thulin, The cost of using exact confidence intervals for a binomial proportion (2013). I understand (to some extent) how it is chosen and how it is calculated.
• $P_1$ and $P_2$ are actually survival probabilities for atoms in optical tweezers. Each measurement of survival is a coin toss: atom in tweezer 1 was either kicked out or it survived, atom in tweezer 2 was either kicked out or it survived. I need to give the tweezer-averaged survival probability, because that is the tweezer-averaged population of a particular energy eigenstate of the atoms.

Answer 1

This special case does not need special methods, maximum likelihood can be used directly. Since the two binomials are independent, estimate $p_1$ based on $P_1$, $p_2$ based on $P_2$. Your interest parameter is then $\bar{p}=\frac{p_1+p_2}{2}$. The maximum likelihood estimator of $\bar{p}$ is then the same function of the individual maximum likelihood estimates, $$ \hat{\bar{p}} =\frac{\hat{p}_1 + \hat{p}_2}{2} $$ with your data, $$ \hat{\bar{p}} = \frac{9/10 + 1000/10000}{2} = \frac12 $$ The variance of this maximum likelihood estimator is found from the usual binomial variance $$ \DeclareMathOperator{\V}{\mathbb{Var}} \V \hat{\bar{p}} = \frac{\frac{p_1(1-p_1)}{n_1} + \frac{p_2(1-p_2)}{n2}}{4} $$ and, using the plug-in estimator this gives $0.00225225$, taking square root gives standard error of $0.04745788$.

From that you can get an approximate confidence interval. Maybe profile likelihood methods can give a better interval (search this site), I will come back to that.

Answer 2

"The most important formula in statistics: $12^2 + 5^2 = 13^2$" (Andrew Gelman: https://statmodeling.stat.columbia.edu/2008/10/25/the_most_import/). What it tells us is that when you have two sources of variation, only the larger one matters, unless they are almost equal.

In your case, one source of variation (the experiment with sample size $10$) is far larger than the other (the experiment with a sample size $1,000$ times greater.) Let's construct a simulation to help us understand what may happen if we just ignore it; while we're at it, we'll compare the results to a confidence interval constructed using the plug-in MLE of the binomial variance and the Normal approximation. We'll generate 10,000 observations from a Binomial(10,0.9) distribution and a Binomial(10000,0.1) distribution for our raw data.

Note that the R function prop.test uses the Wilson score interval with continuity correction, which is generally considered to be much better than the Normal approximation with small to medium sample sizes (https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval).

library(data.table)

p1 <- 0.9
p2 <- 0.1
dt <- data.table(x1 = rbinom(10000, 10, p1),
                 x2 = rbinom(10000, 10000, p2))

dt[, ':='(p1 = x1/10, p2 = x2/10000)]
dt[, ':='(phat = (p1+p2)/2,
          sdhat = sqrt((p1*(1-p1)/10 + p2*(1-p2)/10000))/2)]
dt[, CoverageMLE := (phat - 1.96*sdhat) < 0.5 & 0.5 < (phat + 1.96*sdhat)]

for (i in 1:nrow(dt)) {
    wci <- prop.test(dt$x1[i], 10, conf.level=0.95)$conf.int
    dt[i, lb := wci[1]/2 + p2/2]
    dt[i, ub := wci[2]/2 + p2/2]
}
dt[, CoverageWCI := (lb < 0.5) & (0.5 < ub)]

with coverages:

> dt[, .(mean(CoverageMLE), mean(CoverageWCI))]
       V1     V2
1: 0.6514 0.9891

In this case, ignoring the variability of the estimate of $p_2$ still resulted in a conservative interval.

The key, of course, is using a good confidence interval for $p_1$, since the variability of the estimate of $p_1$ dominates the variability of the estimate of $(p_1+p_2)/2$.

Since you intend to publish, you might get dinged for ignoring the variability of $\hat{p}_2$. You can work around this by simply combining the confidence interval for $p_1$ with the Normal approximation-generated (very small) CI for $p_2$ (acceptable due to the large sample size) in the most obviously conservative way:

dt[, sd_p2 := sqrt(p2*(1-p2)/10000)]
dt[, Coverage2 := (lb - 1.96*sd_p2/2 < 0.5) & (0.5 < ub + 1.96*sd_p2/2)]

which, as it happens, has a coverage in our simulation experiment of 0.9891, identical to that of the CI formed on the basis of $p_1$ only.