Hypothesis test of 2 proportions, with $np < 5$

Question

We are frequently conducting one-tailed hypothesis tests for 2 proportions ($H_0: p_1-p_2=0;\, H_1: p_1-p_2 > 0$). However, $p_2$ is relatively small in terms of $n, x$ and in some cases we find than the condition of $np>5$ is not met, so we can't use the normal approximation and use the normal model for these hypothesis tests.

Any idea or best practice on how to deal with such cases?

Answer 1

One possibility is to take a Bayesian approach using a Beta-Binomial model using Mote Carlo simulations. Specifically, we put independent $\mathrm{Beta}(\alpha_1, \beta_1)$ and $\mathrm{Beta}(\alpha_2, \beta_2)$ priors on $p_1$ and $p_2$, respectively. The posterior for $p_1$ and $p_2$ are independent Beta distributions. The posterior for $p_1$ is a $\mathrm{Beta}(\alpha_1 + x, n-x + \beta_1)$ distribution where $x$ denotes the number of sucesses and $n$ the sample size. By choosing $\alpha_1 = \alpha_2 = \beta_1 = \beta_2 = 1$, we put a uniform distribution on $p_1$ and $p_2$.

Computationally we could proceed as follows:

1. Generate $N$ samples from a $\mathrm{Beta}(\alpha_1 + x_1, n_1 - x_1 + \beta_1)$ and $\mathrm{Beta}(x_2 + \alpha_2, n_2 - x_2 +\beta_2)$ distribution (the posteriors).
2. Calculate the difference of the posteriors.
3. Summarize the posterior (using quantiles to calculate a credible interval, for example). One could also calculate a "Bayesian p-value" by counting the number of samples of the difference that exceed 0.

Here is an R-function that does this:

bayes.prop <- function(x, n, alpha1 = 1, beta1 = 1, alpha2 = 1, beta2 = 1, nsim = 1000, sig.level = 0.95) {

  p1 <- rbeta(nsim, x[1] + alpha1, n[1] - x[1] + beta1)
  p2 <- rbeta(nsim, x[2] + alpha2, n[2] - x[2] + beta2)
  rd <- p1 - p2


  quants <- quantile(rd, c((1-sig.level)/2, (1 + sig.level)/2))

  p.above <- sum(rd > 0)/length(rd)
  p.below <- sum(rd < 0)/length(rd)

  p.two.side <- 2*min(p.above, p.below)

  list(quants = quants, p.val.two = p.two.side, p.greater = p.above, p.less = p.below)

}

For a more sophisticated implementation, check out Rasmus Bååth's blog on this topic.