Confidence Interval – How Confidence Interval for Parameter p of a Bernoulli Trial Varies with p Value

bernoulli-distributionconfidence intervalproportion;sample-size

This is a very basic question (I'm currently studying undergrad level statistics), but I was hoping for some clarification regarding an assertion I read in a newspaper article earlier today. The author asserts that

evidence about the risks (of childbirth) has been hard to come by and difficult to interpret. This is partly because the overall risks of maternal and neonatal death are now very small (about five per 100,000 women die in childbirth and four per 1,000 babies), so large numbers of mums are needed to assess relative risks.

Now, intuitively, this seems uncontroversial – if an event occurs rarely, then you'd expect more trials to be needed to achieve a good estimate of the likelihood of the event, when compared to one that occurs more frequently.

However, I'm having trouble seeing how the mathematics bares this out. If we take the example in the article, we can model the probability of the death of a mother during childbirth as a Bernoulli trial with an estimate of the parameter $p$ given by $\hat{p}=\frac{5}{1000}$.

From this, we can construct a 95% confidence interval for the true value of p:

$$ p^{\pm} = \hat{p} \pm 1.96 \frac{p(1-p)}{\sqrt{n}} $$

However, if we take the limit as $p \to 0$ we get

$$ \begin{aligned}
\lim_{p \to 0} p^{\pm} = \hat{p} \pm 1.96 \frac {0(1-0)}{\sqrt{n}}\\
= \hat{p} \pm 1.96 \frac {0}{\sqrt{n}}\\
= \hat{p} \pm 0
\end{aligned} $$

Therefore, it appears that our estimate of $p$ in fact becomes more accurate as $p$ gets smaller, regardless of our value of $n$. This seems pretty counterintuituve to me – am I going wrong somewhere, and if so, where?

Many thanks,

Tim

Best Answer

The method you use, a normal approximation, is an archaicism and should never be taught or even offered as an option in software. It has very poor coverage properties, particularly for small proportions as in your example.

There are many alternative approaches to calculating these intervals, with varying assumptions and coverage characteristics. Some are very ad hoc in design and so are hard to prefer for pedagogic purposes. My preference is the method of Wilson, sometimes called Wilson's scores intervals. It approximates a conditional interval and has excellent frequentist properties.

See this answer for a little more detail: Discrete functions: Confidence interval coverage?

See this question for a formal statement of the meaning of different types of CI for binomial proportions: Statement of result for binomial confidence intervals

This one for confidence interval coverage: Clarification on interpreting confidence intervals?

Related Solutions

Solved – Confidence interval for Bernoulli sampling

If the average, $\hat{p}$, is not near $1$ or $0$, and sample size $n$ is sufficiently large (i.e. $n\hat{p}>5$ and $n(1-\hat{p})>5$, the confidence interval can be estimated by a normal distribution and the confidence interval constructed thus:

$$\hat{p}\pm z_{1-\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$$
If $\hat{p} = 0$ and $n>30$, the $95\%$ confidence interval is approximately $[0,\frac{3}{n}]$ (Javanovic and Levy, 1997); the opposite holds for $\hat{p}=1$. The reference also discusses using using $n+1$ and $n+b$ (the later to incorporate prior information).
Else Wikipedia provides a good overview and points to Agresti and Couli (1998) and Ross (2003) for details about the use of estimates other than the normal approximation, the Wilson score, Clopper-Pearson, or Agresti-Coull intervals. These can be more accurate when above assumptions about $n$ and $\hat{p}$ are not met.

R provides functions binconf {Hmisc} and binom.confint {binom} which can be used in the following manner:

set.seed(0)
p <- runif(1,0,1)
X <- sample(c(0,1), size = 100, replace = TRUE, prob = c(1-p, p))
library(Hmisc)
binconf(sum(X), length(X), alpha = 0.05, method = 'all')
library(binom)
binom.confint(sum(X), length(X), conf.level = 0.95, method = 'all')

Agresti, Alan; Coull, Brent A. (1998). "Approximate is better than 'exact' for interval estimation of binomial proportions". The American Statistician 52: 119–126.

Jovanovic, B. D. and P. S. Levy, 1997. A Look at the Rule of Three. The American Statistician Vol. 51, No. 2, pp. 137-139

Ross, T. D. (2003). "Accurate confidence intervals for binomial proportion and Poisson rate estimation". Computers in Biology and Medicine 33: 509–531.

Solved – Compute a confidence interval for Bernoulli distribution

Your expected value is $3000\times0.02=60$, the variance is $3000\times 0.02\times(1-0.02)=58.8$. Simulating 100,000 trials and plotting a histogram, I would say you can use the normal approximation unless you have very strong requirements on accuracy (in which case the R help page says that "qbinom() uses the Cornish-Fisher Expansion to include a skewness correction to a normal approximation, followed by a search", which may help).

nn <- 3000
n.sim <- 100000
foo <- rbinom(n=n.sim,size=nn,prob=.02)
hist(foo,breaks=seq(min(foo)-.5,max(foo)+.5,by=1))

enter image description here

Best Answer

Related Solutions

Solved – Confidence interval for Bernoulli sampling

Solved – Compute a confidence interval for Bernoulli distribution

Related Question