Manipulating Binomial Distribution 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)$.
 A: 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

A: 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}
$$
