Manipulating Binomial Distribution

Question

Recently, I've been reading Yudi Pawitan's book, In All Likelihood.

In the book, there's a section on profile likelihood; the methods explored in this section are subsequently applied to some data on heart attack prevalence amongst two distinct groups: people taking aspirin and people subjected to a placebo; the groups are modelled as $\text{Bin}(n_{a}, \theta_{a})$ and $\text{Bin}(n_{p}, \theta_{p})$ respectively. Since $n_{a}$ and $n_{p}$ are large while the event rates are small, $X_{a}$ and $X_{p}$ are approximately Poisson with parameters $n_{a}\theta_{a}$ and $n_{p}\theta_{p}$ respectively.

You can see the example here (Pages 87 and 88).

At the end, it is shown that the profile likelihood for a parameter of interest $\theta$ is given by:

$$ L(\theta, \theta_{p}) = \text{constant} \times e^{-\theta_{p}(n_{a}\theta+n_{p})} \theta^{x_{a}} \theta_{p}^{x_{a}+x_{p}} (1)$$

where $\hat{\theta}_{p} = \frac{x_{a}+x_{p}}{n_{a}\theta+n_{p}}$, the MLE for $\theta_{p}$, should be substituted into the above model for $\theta_{p}$.

and

$n_{a}$ : number of people in aspirin group, $n_{p}$ : number of people placebo group, $x_{a}$ : number of heart attacks amongst aspirin group, $x_{p}$ : number of heart attacks amongst placebo group.

also

$\theta = \frac{\theta_{a}}{\theta_{p}}$, $\theta_{a}$ : probability of heart attack in aspirin group and $\theta_{a}$ : probability of heart attack in placebo group.

The above likelihood is then further expressed as binomial, taking the following form:

$$ L(\theta) = \text{constant} \times \left(\frac{n_{a}\theta}{n_{a}\theta+n_{p}}\right)^{x_{a}} \left(1-\frac{n_{a}\theta}{n_{a}\theta+n_{p}}\right)^{x_p} (2)$$

I have tried to understand how to get from the expression given in $(1)$ to the expression given in $(2)$, but when I make the substitution for $\hat{\theta}_{p}$, I get the following:

$$ e^{-(x_{a}+x_{p})}\theta^{x_{a}}(x_{a}+x_{p})^{x_{a}+x_{p}}\left(\frac{1}{n_{a}\theta+np}\right)^{x_{a}+x_{p}} $$

where $e^{-(x_{a}+x_{p})} = \text{constant}$ as it doesn't depend on $\theta$.

I believe this is along the right lines, but I'm not sure where to go next.

Perhaps somebody could explain how to get from expression $(1)$ to expression $(2)$.

Answer 1

First, the expression (1) is the likelihood, not the profile likelihood. It becomes profile likelihood for $\theta$ (eliminating the nuisance parameter $\theta_p$) only after you substitute in (1) $\hat{\theta}_p $ for $\theta_p$, the conditional maximum likelihood estimator for $\theta_p$, which is a function of $\theta$. The likelihood (1) is based on the Poisson approximation, and is in the book by Pawitan only used to develop $\hat{\theta}_p $, since the exact binomial likelihood does not lead to an explicit solution for the conditional mle.

Then, (2) is based on the exact binomial likelihood, but using the formula for $\hat{\theta}_p $ based on the Poisson approximation. This mixed use is not commented in the book, so maybe look somewhat strange.

---

You did not ask for this, but nevertheless. The way to do this computations today is with R and the glm function. I will show solution with both the binomial and Poisson likelihood.

Binomial, we need to use the log link function to get the parametrization used in the book.

library(tidyverse)  
pawdf <- tibble(x=c(139, 239), n=c(11037, 11034), 
           T=factor(c("active", "placebo"), 
           levels=c("placebo", "active")))
    
mod.bin <- glm( cbind(x, n-x)  ~ T, 
    family=binomial(link="log"), data=pawdf)
summary(mod.bin)
    
mod.bin.prof <- profile(mod.bin, which="Tactive")
    
exp(confint(mod.bin.prof))
                   2.5 %    97.5 %
    (Intercept)       NA        NA
    Tactive     0.471508 0.7141749    
    # very close to Pawitan book

Poisson:

mod.po <- glm( x  ~ T + offset(log(n)), 
    data=pawdf, family=poisson(link="log"))
summary(mod.po)
    
mod.po.prof <- profile(mod.po,  which="Tactive")
    
exp(confint(mod.po.prof))
                   2.5 %    97.5 %
    (Intercept)       NA        NA
    Tactive     0.470737 0.7153891

Answer 2

As @kjetilbhalvorsen explains, Equation (1) is the likelihood $L(\theta,\theta_p)$:

$$ \begin{aligned} L(\theta,\theta_p) \propto \exp\big\{-\theta_p(n_a\theta+n_p)\big\}\theta^{x_a}\theta_p^{x_a+x_p} \end{aligned} $$

while Equation (2) is the profile likelihood $L(\theta)$.

To derive the profile likelihood $L(\theta)$ from the likelihood $L(\theta,\theta_p)$ we want to eliminate the nuisance parameter $\theta_p$ by maximimizing $L(\theta,\theta_p)$ for each value of the parameter of interest $\theta$.

One way to find the MLE for $\theta_p$ for each fixed $\theta$ is to notice that, as a function of $\theta_p$, $L(\theta)$ is proportional to the probability density function (pdf) of the Gamma distribution:

$$ \begin{aligned} f(x) \propto x^{\alpha - 1}\exp\left\{-\beta x\right\} \end{aligned} $$ where $\theta_p$ plays the role of $x$, $(x_a+x_p+1)$ is the shape $\alpha$ and $(n_a\theta+n_p)$ is the rate $\beta$.

This observation is handy because we don't have to do any derivations for find the MLE $\widehat{\theta}_p(\theta)$. We know (from Wikipedia) that the density is maximized at the mode $(\alpha-1)/\beta$. And conveniently the mode is positive as long as $\alpha\geq1$. This holds here because the shape $\alpha$ is given by $(x_a+x_p+1)$.

Finally we plug in the MLE $\widehat{\theta}_p(\theta) = (x_a+x_p) / (n_a\theta+n_p)$ into the likelihood $L(\theta,\theta_p)$ and simplify to derive the profile likelihood $L(\theta)$.

$$ \begin{aligned} L(\theta) &\propto \exp\big\{-(x_a+x_b)\big\}\theta^{x_a}\left(\frac{x_a+x_b}{n_a\theta+n_p}\right)^{x_a+x_p} \\ &\propto \big(n_a\theta\big)^{x_a}\left(\frac{n_p}{n_a\theta+n_p}\right)^{x_a+x_p} \\ &\propto \left(\frac{n_a\theta}{n_a\theta+n_p}\right)^{x_a}\left(1 - \frac{n_a\theta}{n_a\theta+n_p}\right)^{x_p} \end{aligned} $$