Bayesian Shrinkage Estimation – Questions About Paper

bayesiangibbsmarkov-chain-montecarlomultinomial-distributionsampling

I am reading the paper Bayesian Shrinkage Estimation of the Relative Abundance of mRNA Transcripts Using SAGE, and I am trying to work out the calculations for the complete conditionals for the Gibbs scheme presented in Section 3.5 (the prior structure is described in Section 3.1).

I am confused about how to calculate the full joint probability model and how to derive $p(\lambda \mid X, …)$, $p(\pi^* \mid X, …)$, $p(q_A \mid X, …)$, and $p(q_S \mid X, …)$. Is the below expression for the full joint correct?

\begin{align*}p(\pi^*, q_A, q_S, \lambda, X) &\propto {(\pi^*)}^{\sum_{j=1}^{|A|} \lambda_j} (1-\pi^*)^{\sum_{j=|A| + 1}^{N} \lambda_j} \prod_{j=1}^{|A|} q_{A_j}^{\lambda_j} \prod_{j=|A|+1}^N q_{S_j}^{\lambda_j} \cdot \prod_{j = 1}^{|A|} q_{A_j}^{\theta_A – 1} \prod_{j = |A| + 1}^{N} q_{S_j}^{\theta_S – 1} \\ &\cdot \frac{1}{B(a_{\pi^*}, b_{\pi^*})} (\pi^*)^{a_{\pi^*} – 1}(1 – \pi^*)^{b_{\pi^*} – 1} \cdot P^{\lambda} (1 – P)^{1 – \lambda}\end{align*}

I'm not sure where to go from here to calculate $p(\lambda \mid X, …)$, so I'd appreciate any help. Thanks in advance!

Edit: I was thinking of defining $A_1 = \{j : \lambda_j = 1\}$ and $A_0 = \{j : \lambda_j = 0\}$ with $M_1 = |A_1|$ and $M_0 = N – M_1$. Then we assume WLOG that $\lambda_j = 1$ for $j = 1, …, M_1$ and $\lambda_j = 0$ for $j = M_1 + 1, …, M_1 + M_0$. Then would the joint be

\begin{align*}p(\pi^*, q_A, q_S, \lambda, X) &\propto {(\pi^*)}^{\sum_{j=1}^{M_1} \lambda_j} (1-\pi^*)^{\sum_{j=M_1 + 1}^{N} \lambda_j} \prod_{j=1}^{M_1} q_{A_j}^{\lambda_j} \prod_{j=M_1+1}^N q_{S_j}^{\lambda_j} \cdot \prod_{j = 1}^{M_1} q_{A_j}^{\theta_A – 1} \prod_{j = M_1 + 1}^{N} q_{S_j}^{\theta_S – 1} \\ &\cdot \frac{1}{B(a_{\pi^*}, b_{\pi^*})} (\pi^*)^{a_{\pi^*} – 1}(1 – \pi^*)^{b_{\pi^*} – 1} \cdot P^{\lambda} (1 – P)^{1 – \lambda}?\end{align*}

Edit 2: I am still stuck on calculating the conditional posterior for $\lambda$ and $\pi^*$. Would it make sense to evaluate the joint at $\lambda_j = 1$ and $\lambda_j = 0$? I realize that there is a factor to consider in $p(q_A) = Dir(q_A \mid a_\lambda,\ldots,a_\lambda)$ for $\lambda_j = \lambda$.

Best Answer

Define the index set $\mathcal{I}:=\{j: \lambda_j=1\}$. Then, the posterior becomes $$ \begin{align} p(\pi^*,\lambda,q_A,q_S\mid X)&\propto (\pi^*)^{|\mathcal{I}|}(1-\pi^*)^{N-|\mathcal{I}|}\prod_{j\in \mathcal{I}} q_{A_j}^{\lambda_j+\theta_A-1} \prod_{j \in \mathcal{I}^c} q_{S_j}^{1-\lambda_j +\theta_S-1} \\ &\times (\pi^*)^{a_\pi-1}(1-\pi^*)^{b_\pi-1} \,\prod_{j=1}^N P^{\lambda_j} (1-P)^{1-\lambda_j}, \end{align} $$ where $|\mathcal{I}|$ indicates the cardinality of $\mathcal{I}$. The prior for $\lambda_j$ doesn't make sense without the subscript $j$ so I changed that part accordingly. When deriving a full conditional, ignore the components that don't have the parameter of interest. For example, $$ p(\pi^*\mid -) \propto (\pi^*)^{|\mathcal{I}| + a_\pi-1}(1-\pi^*)^{N-|\mathcal{I}|+b_\pi-1} \sim \mathrm{Be}(a_\pi + |\mathcal{I}|, b_\pi + N - |\mathcal{I}|), $$ and $$ p(\lambda_j\mid - ) \propto (q_{A_j} P)^{\lambda_j}\times \{q_{S_j}(1-P)\}^{(1-\lambda_j)} \sim \mathrm{Bernoulli}(p^*), $$ where $p^* = \frac{q_{A_j}P}{q_{A_j}P+ q_{S_j}(1-P)}$, since the components aren't guaranteed to sum to one.

Related Solutions

Bayesian – Extracting Marginal Likelihood from the Gibbs Output

There is a slight programming mistake in the prior

log_prior <- function(w, mu, va) {
    lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w))
    + sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE))
    + sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va))
}

as it should be instead

log_prior <- function(w, mu, va) {
    lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w)) +
      sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE)) +
      sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va))
}

Rerunning the code this way leads to

> chib
[1] -228.194

which is not the value produced in Chib (1995) for that case! However, in Neal's (1999) reanalysis of the problem, he mentions that

According to one anonymous JASA referee, the figure of -224.138 for the log of the marginal likelihood for the three component model with unequal variances that was given in Chib's paper is a "typo" wtih the correct figure being -228.608.

So this solves the discrepancy issue.

Bayesian Sampling – How to Use Mixture Dirichlet MCMC Model Effectively

The problem is very interesting. Instead of directly imposing a Dirichlet prior on the multinomial probabilities, it reparameterizes $\pi$ in terms of $\pi^*$ that indicates the global probability (which the authors call the "probability of observing some prevalent outcome"), and that global probability is divvied up into pieces and shared by each category. The "difficulty" (if you will) of this parameterization is that the Dirichlet parameters $(q_1,\ldots,q_{M_1})$ and $(r_1,\ldots,r_{M_0})$ are transdimensional throughout the MCMC chain.

Enough with my subjective addendum. First things first, let's write down the likelihood and the priors. To summarize:

(Sampling model) $y \sim \mathrm{Multinomial}_{N-1}(n; \pi_1,\ldots,\pi_N)$, where $$ \pi_j = \begin{cases}\pi^* q_j & \text{if }z_j=1\\ (1-\pi^*)r_j &\text{if }z_j=0 \end{cases} $$ which gives $$ \begin{align} &\binom{n}{y_1 \cdots y_N} \prod_{j=1}^N (\pi^* q_j)^{z_jy_j}((1-\pi^*)r_j)^{(1-z_j)y_j} \\ &= \binom{n}{y_1 \cdots y_N} \prod_{j\in A_1} (\pi^* \tilde{q}_j)^{y_j} \prod_{j \in A_0} ((1-\pi^*)\tilde{r}_j)^{y_j} \end{align} $$ The two expressions are the same and later, the ability to switch to whichever is more convenient comes in quite handy. The first expression puts the dependency on $z_j$ explicit in the exponenets, whereas the second one makes it implicit in the definition of the sets $A_0$ and $A_1$ as well as $\tilde{q}$ and $\tilde{r}$.
(Prior for $z_j$) $z_j \overset{\text{iid}}{\sim} \mathrm{Bernoulli}(\rho)$
(Prior for $\pi^*$) $\pi^* \sim \mathrm{Be}(a^*,b^*)$
(Prior for $\tilde{q}_1,\ldots,\tilde{q}_{M_1}$) $\mathrm{Dirichlet}_{M_1-1}(a_1,\ldots,a_1)$
(Prior for $\tilde{r}_1,\ldots,\tilde{r}_{M_0}$) $\mathrm{Dirichlet}_{M_0-1}(a_0,\ldots,a_0)$

Combining all the above, the posterior is proportion to $$ \begin{align} &\binom{n}{y_1 \cdots y_N} \prod_{j=1}^N (\pi^* q_j)^{z_jy_j}((1-\pi^*)r_j)^{(1-z_j)y_j}\\ &\times \prod_{j=1}^N \rho^{z_j}(1-\rho)^{1-z_j} (\pi^*)^{a^*-1}(1-\pi^*)^{b^*-1} \prod_{i=1}^{M_1}\tilde{q}_i^{a_1-1} \prod_{i=1}^{M_0}\tilde{r}_i^{a_0-1}. \end{align} $$

In part (c), I am stuck on calculating the conditional posterior for $z_j$ and $q$, $r$.

For this, discard irrelevant components in the joint likelihood, since from the perspective of the complete conditional, they're all constants and do not contribute to the distribution. You can normalize the complete conditional density after. That being said, for $z_j$, $$ \begin{align} p(z_j\mid - ) &\propto (\pi^* q_j)^{z_jy_j}((1-\pi^*)r_j)^{(1-z_j)y_j} \rho^{z_j}(1-\rho)^{1-z_j}\\ &= \left\{(\pi^* q_j)^{y_j} \rho \right\}^{z_j} \left\{((1-\pi^*)r_j)^{y_j} \right\}^{1-z_j} \end{align} $$ which is a Bernoulli distribution. The probability should be normalized though.

For $q$, by virtue of being a complete conditional and having information about $z_j$, we know which $y_j$ are to be used. Then, $p(q\mid -)\propto \prod_{j=1}^{M_1}\tilde{q}_j^{y_j+a_1-1}$, which is a Dirichlet distribution. Likewise, $p(r\mid -) \propto \prod_{j=1}^{M_0} \tilde{r}_j^{y_j+a_0-1}$.

I am also stuck on part (g) - how would I compute an appropriate acceptance probability for the Metropolis-Hastings scheme?

In a Metropolis-Hastings scheme, the acceptance probability is given to us (how lucky!). It's always, $$ \alpha = \min\left\{1, \frac{L(\theta')\pi(\theta') q(\theta\mid \theta')}{L(\theta)\pi(\theta) q(\theta'\mid \theta)} \right\}, $$ where $L(\cdot)$ is the likelihood function, $\pi(\cdot)$ is the prior densities, $\theta$ is the vector of parameters as you have it in the chain, and $\theta'$ is the candidate parameters sampled from the proposal distribution, which if accepted will be $\theta \gets \theta'$. For example, the proposal distribution of $z_i$ is given by $P(z_i'=1)=0.5$ in the problem. This is a Bernoulli distribution whose probability mass function is $0.5^{z_i}0.5^{1-z_i} = 0.5$. Since the proposal density does not depend on $z_i$, the Hastings ratio $\frac{q(z_i\mid z_i')}{q(z_i'\mid z_i)}$ reduces to 1, but the proposal density of $q'$ is the complete conditional where $z$ conditioned is the candidate value $z'$ so you won't be able to ignore the Hastings ratio. In an implementation though, you can create a variable accept_rate = 0 outside the loop and increment it each time a proposal is accepted. At the end of the program, divide accept_rate by the number of iterations, and it will give you a Monte Carlo estimate of the acceptance rate.

Best Answer

Related Solutions

Bayesian – Extracting Marginal Likelihood from the Gibbs Output

Bayesian Sampling – How to Use Mixture Dirichlet MCMC Model Effectively

Related Question