I have survey data for women and men following a binomial distribution. Their means are $p_1$ and $p_2$ respectively. I have calculated $(p_1-p_2)/p_2$ and would like to attach a confidence interval to this percentage. How can this be done? Given the number of people surveyed, a normal approx would be fine here.
2 Answers
You can use profile likelihood methods (other answers are surely possible, but I will show that.)
You have binomial counts from two groups, men and women. We write
$$
M \sim \mathcal{Bin}(m, p_m) \\
W \sim \mathcal{Bin}(w, p_w)
$$
and the focus (or interest) parameter is $Q=\frac{p_w-p_m}{p_m}$. Assuming (reasonable) that the counts in the two groups are independent, we find the loglikelihood function is
$$
\ell_0(p_m, p_w)\propto M\log(p_m) + (m-M)\log(1-p_m) + W\log(p_w) + \\ (w-W)\log(1-p_w)
$$
A little algebra shows that $p_w=p_m(1+Q)$, substituting that above we find
$$\ell(Q,p_m)=\ell_0(p_m, p_m(1+Q))
$$
and the profile likelihood for the focus parameter $Q$ is defined by
$$
\ell_Q(Q)= \max_{0\le p_m\le 1}\ell(Q,p_m)
$$ and we can find a confidence interval by the asymptotic theory for profile likelihood, see Profile likelihood confidence interval proof. In this case it might even be possible to do that symbolically, but I will use numerical methods in R
. First, confidence intervals can be read off this plot of the profile likelihood:
What is plotted here is the square root of the profile likelihood-based deviance, so a perfect V shape would indicate a quadratic profile likelihood.
The code used is
set.seed(7*11*13) # My public seed
m <- 50; w <- 40
p_m <- 0.4; p_w <- 0.6; Q <- (p_w-p_m)/p_m
M <- rbinom(1, m, p_m); W <- rbinom(1, w, p_w)
# Loglikelihood function:
mloglik0 <- function(p_m, p_w) -dbinom(M, m, p_m, log=TRUE) -
dbinom(W, w, p_w, log=TRUE)
mloglik1 <- function(Q, p_m) mloglik0(p_m, p_m*(1+Q))
mod <- bbmle::mle2(mloglik1, start=list(p_m=0.4, Q=0))
mod.prof <- bbmle::profile(mod, which=1)
confint(mod.prof)
2.5 % 97.5 %
-0.2176006 0.5851981
Addendum @cdalitz in another answer refers to the R package DescTools
, which implements 11 different methods for the binomial difference (see that answer for references). Among those 3 different profile likelihood methods, one of them called true profile likelihood, the other two modifications using exact tail probability calculations in place of the asymptotic methods used for the true profile likelihood. What I have used here corresponds to the true profile likelihood, as I am doing no exact calculations.
But there is another important difference: This answer is about the interest parameter $Q=\frac{p_w-p_m}{p_m}$, a relative difference, not the difference itself! To facilitate comparisons I redo the calculations for $Q=p_w-p_m$, and compare with results from that R package:
mloglik1_diff <- function(Q, p_m) mloglik0(p_m, p_m + Q)
mod_diff <- bbmle::mle2(mloglik1_diff, start=list(p_m=0.4, Q=0))
mod_diff_prof <- bbmle::profile(mod_diff, which=1)
confint(mod_diff_prof)
2.5 % 97.5 %
-0.1394750 0.2638033
METHODS <- c("ac", "wald", "waldcc", "score", "scorecc", "mn",
"mee", "blj", "ha", "hal", "jp")
CIs <- matrix(NA, length(METHODS), 3)
colnames(CIs) <- c( " est", "lwr.ci", "upr.ci" )
rownames(CIs) <- METHODS
for (method in seq_along(METHODS)) CIs[method, ] <-
DescTools::BinomDiffCI(W, w, M, m, method=METHODS[method])
CIs
est lwr.ci upr.ci
ac 0.065 -0.1381246 0.2608353
wald 0.065 -0.1385663 0.2685663
waldcc 0.065 -0.1610663 0.2910663
score 0.065 -0.1360111 0.2557365
scorecc 0.065 -0.1511457 0.2699042
mn 0.065 -0.1400075 0.2617998
mee 0.065 -0.1388655 0.2607800
blj 0.065 -0.1406425 0.2668979
ha 0.065 -0.1534193 0.2834193
hal 0.065 -0.1378810 0.2607925
jp 0.065 -0.1380281 0.2609785
-
$\begingroup$ To which of the three profile likelihood methods (7,8,9) evaluated by Newcombe ("Interval estimation for the difference between independent proportions: comparison of eleven methods", Stat. Med., 1998) does your method correspond? $\endgroup$– cdalitzCommented Jul 27, 2021 at 13:49
-
$\begingroup$ @cdalitz: I have updated my answer now $\endgroup$ Commented Jul 28, 2021 at 4:16
Unfortunately there is no universally accepted way for computing a confidence interval for a difference in binomial proportions. The R function BinomDiffCI in the package DescTools offers eleven different options, and its help page gives references to publications.
Newcombe recommends the methods by Miettinen and Nurminen (1985) or Mee (1984) (and his own methods, of course) in
R.G. Newcombe: "Interval Estimation for the Difference Between Independent Proportions: Comparison of Eleven Methods." Statistics in Medicine, 17, 873–890, 1998