Bayesian model: binomial(s) conditional on Poisson

Question

I am struggling a little with how to conceptualize and specify a particular Bayesian model. I suspect the solution is rather simple but for some reason I am having a hard time thinking about this.

My primary purpose is to model the incidence of a disease (let's denote this as $\lambda$), where the total number of cases of the disease in the population (let's denote this as $X$) is unobserved. What we DO observe is the number of cases of the disease as detected by some surveillance program, let's call this $Y$. Note that $Y$ is just a subset of $X$. We have made some assumptions based on a combination of previous data and our operative understanding of the structural context to relate $Y$ to $X$. Essentially, some subset of $X$ has a probability of seeking care (let's call this $C$); then, some subset of $C$ has a probability of being eligible for the surveillance program (let's call this $E$); finally, some subset of $E$ will have a positive test result that is a function of the sensitivity of the assay. This final step is what we observe ($Y$).

We can write this all out using simplified notation as:

$X|\lambda \sim Poisson(\lambda)$

$C|X,\pi_c \sim Binomial(X,\pi_c)$

$E|C,\pi_e \sim Binomial(C,\pi_e)$

$Y|E,\pi_y \sim Binomial(E, \pi_y)$

Therefore, through some simple algebra, we can show that:

$Y|\lambda,\pi_c,\pi_e,\pi_y \sim Poisson(\lambda\pi_c\pi_e\pi_y)$

Where I am struggling is then how to take the next step and evaluate the posterior, since the parameter is a product of multiple parameters. The simplest way of doing this would be to simply assume a conjugate prior on $\lambda$ and treat each of the $\pi$ parameters as fixed, but in practice we will have some amount of information on each of those parameters that we would want to incorporate into the model (e.g. we may reasonable expect those probabilities to be different among different subpopulations).

But I am struggling with how to move beyond that simplest case. If a different prior distribution (which may be functions of covariates) is assumed for each of the four parameters in that final equation, how do I derive the posterior (i.e. $P(\lambda|Y)$)? What if I am also interested in inference on the $\pi$ terms? Is that even possible within a single model? Or am I just thinking about this all wrong?

I suspect there is a very simple solution here, but for some reason I am feeling completely lost on how to go about doing this, and none of the examples I can find online are helping me.

EDIT: Here is my work trying to figure out the posterior. If by the definition of conditional probability we have:

$P(\lambda|Y,\pi_c,\pi_e,\pi_y)=\frac{P(Y,\pi_c,\pi_e,\pi_y,\lambda)}{P(Y,\pi_c,\pi_e,\pi_y)}$

The numerator can be decomposed using the chain rule as:

$P(\lambda|Y,\pi_c,\pi_e,\pi_y)=\frac{P(Y|\lambda,\pi_c,\pi_e,\pi_y)P(\pi_y|\lambda,\pi_e,\pi_c)P(\pi_e|\lambda,\pi_c)P(\pi_c|\lambda)P(\lambda)}{P(Y,\pi_c,\pi_e,\pi_y)}$

We can then perform a similar operation on the denominator:

$P(\lambda|Y,\pi_c,\pi_e,\pi_y)=\frac{P(Y|\lambda,\pi_c,\pi_e,\pi_y)P(\pi_y|\lambda,\pi_e,\pi_c)P(\pi_e|\lambda,\pi_c)P(\pi_c|\lambda)P(\lambda)}{P(Y|\pi_c,\pi_e,\pi_y)P(\pi_y|\pi_c,\pi_e)P(\pi_e|\pi_c)P(\pi_c)}$

Since we can (reasonably, I think) assume the parameters are independent, this all simplifies to:

$P(\lambda|Y)=\frac{P(Y|\lambda,\pi_c,\pi_e,\pi_y)P(\pi_y)P(\pi_e)P(\pi_c)P(\lambda)}{P(Y|\pi_c,\pi_e,\pi_y)P(\pi_y)P(\pi_e)P(\pi_c)}$

Which further simplifies as:

$P(\lambda|Y)=\frac{P(Y|\lambda,\pi_c,\pi_e,\pi_y)P(\lambda)}{P(Y|\pi_c,\pi_e,\pi_y)}$

But this doesn't look right to me. I don't quite understand what the quantity in the denominator would be (the conditional likelihood of $Y$ after marginalizing over $\lambda$?) or how I would calculate it. (Even the calculation of the likelihood for this model I am vague on, but that's more of an implementation question than a theoretical one).

Note that this would just be the posterior distribution of $\lambda$, if we were interested in a joint posterior so we could make inference over other parameters, we would have something like:

$P(\lambda,\pi_c,\pi_e,\pi_y|Y)=\frac{P(Y|\lambda,\pi_c,\pi_e,\pi_y)P(\pi_y)P(\pi_e)P(\pi_c)P(\lambda)}{P(Y)}$

Which looks a lot cleaner, but then I'm not sure how to draw marginal inference (e.g. on $\lambda$) from such a joint posterior (though I suppose that's where the assumption of independence comes in handy again?). I have a feeling that the answer is sitting right here in front of me, but something about the intuition hasn't yet "clicked" in my head and I'm feeling rather lost...

Answer 1

Looks like you answered most of your own questions already!

Where I am struggling is then how to take the next step and evaluate the posterior, since the parameter is a product of multiple parameters.

Bayes theorem still holds with joint distributions, all you would have to do is to specify a joint distribution for the set of parameters. Typically in a Bayesian model, this could be defined as a product of conditionals and marginals. In your scenario, it might be reasonable to define them as a product of marginals of $\pi_i$'s distributions.

The simplest way of doing this would be to simply assume a conjugate prior on λ and treat each of the π parameters as fixed, but in practice we will have some amount of information on each of those parameters ...

That's right, in that case you would probably use a beta prior for each of these $\pi_i$s. However, your posterior $P_{all} = P(\lambda, \pi_c, \pi_e, \pi_y| Y)$ is not going to look pretty. This is where Markov chain Monte Carlo (MCMC) methods come in. If you're interested, I suggest studying the Metropolis algorithm (or Metropolis-Hastings, slightly different). In short: you sample from the posterior to approximate the distribution (e.g., you would have 10000 points of $(\lambda^{(t)}, \pi_c^{(t)}, \pi_e^{(t)},\pi_y^{(t)})$, for $t=1,2,...,10000$ whose distribution is approximately $P_{all}$.

To obtain the marginal posterior $P(\lambda | Y )$, you would have to integrate over $\pi_c,\pi_e,\pi_y$ from $P_{all}$. That is $P(\lambda | Y) = \int_{[0,1]^3}P(\lambda, \pi_c,\pi_e,\pi_y| Y)d\pi_cd\pi_ed\pi_y$. Suppose you have approximated the posterior $P_{all}$ via MCMC, then deriving the marginal distribution of $P(\lambda | Y)$ simply comes down to ignoring the $\pi_c,\pi_e,\pi_y$ samples (summing over discrete samples).

Here's how you would derive the posterior (and/or marginals)

$$P(\lambda, \pi_c,\pi_e,\pi_y | Y) = \frac{P(Y,\lambda,\pi_c,\pi_e,\pi_y)}{P(Y)} \\ = \frac{ P(Y|\lambda, \pi_c, \pi_e,\pi_y) P(\lambda)P(\pi_c)P(\pi_e)P(\pi_y)}{\int_{\lambda,\pi_c,\pi_e,\pi_y} P(Y|\lambda, \pi_c, \pi_e,\pi_y) P(\lambda)P(\pi_c)P(\pi_e)P(\pi_y)d\lambda d\pi_c d\pi_e d\pi_y}$$

$$P(\lambda| Y) = \int_{\pi_c,\pi_e,\pi_y}P(\lambda, \pi_c, \pi_e,\pi_y | Y) d\pi_c d\pi_e d\pi_y $$

$$P(\pi_c| Y) = \int_{\lambda,\pi_e,\pi_y}P(\lambda, \pi_c, \pi_e,\pi_y | Y) d\lambda d\pi_e d\pi_y $$