Variance of marginal posterior distribution

Question

Suppose $Y_1,\dots,Y_n\mid\mu,\sigma^2 \sim \text{ iid } N(\mu,\sigma^2)$ and suppose the priors $\mu \mid \sigma^2 \sim N(\mu_0, \sigma^2 / \kappa_0)$ and $1/\sigma^2 \sim \text{gamma}(\nu_0/2, \nu_0 \sigma_0^2 / 2)$ are placed on the unknown parameters. Then \begin{align*} p(\sigma^2 \mid y_1,\dots,y_n) &\propto p(\sigma^2)p(y_1,\dots,y_n\mid\sigma^2) \\ &= p(\sigma^2) \int p(y_1,\dots,y_n \mid \mu,\sigma^2)p(\mu\mid\sigma^2)d\mu \end{align*} It turns out (see Hoff) with these priors that the marginal posterior distribution on $\sigma^2$ is $$1/\sigma^2\mid y_1,\dots,y_n \sim \text{gamma}(\nu_n/2,\nu_n\sigma_n^2/2)$$ where $$\nu_n = \nu_0 + n, \text{ and } \sigma_n^2 = \frac{1}{\nu_n}\left( \nu_0\sigma_0^2 + (n-1)s^2 + \frac{\kappa_0 n}{\kappa_0 + n} ( \bar{y}-\mu_0 )^2 \right)$$ where $\bar{y}$ and $s^2$ are the sample mean and unbiased sample variance. Thus we have that $$ \text{Var}\left[\frac{1}{\sigma^2} \mid y_1,\dots,y_n\right] = (\nu_n/2) / (\nu_n\sigma_n^2/2)^2 = \frac{2}{\nu_n\sigma_n^4}. $$

Now suppose we fix all hyperparameters except $\kappa_0$ and analyze the effect of increasing the variance of $\mu\mid\sigma^2$ (by decreasing $\kappa_0$). Since $f(\kappa_0)=\kappa_0 n / (\kappa_0 + n)$ is monotonically increasing as a function of $\kappa_0$, we see that \begin{align*} \kappa_{0,a} < \kappa_{0,b} &\Rightarrow \kappa_{0,a} n / (\kappa_{0,a} + n) < \kappa_{0,b} n / (\kappa_{0,b} + n) \\ &\Rightarrow \sigma_{n,a}^2 < \sigma_{n,b}^2 \\ &\Rightarrow \sigma_{n,a}^4 < \sigma_{n,b}^4 \\ &\Rightarrow \frac{2}{\nu_n\sigma_{n,b}^4} < \frac{2}{\nu_n\sigma_{n,a}^4}, \end{align*} i.e., increasing the variance of $\mu\mid\sigma^2$ results in an increase in the marginal posterior variance $\text{Var}\left[1/\sigma^2\mid y_1,\dots,y_n\right]$.

Intuitively, this result makes a lot of sense because if $p(\mu \mid \sigma^2)$ is more spread out, the averaging (integral) over $p(\mu \mid \sigma^2)$ above makes $p(\sigma^2 \mid y_1,\dots,y_n)$ more spread out.

This makes me think that this result should be able to be proved in general, i.e., given a sampling model $f_{Y\mid\theta,\phi}$ depending on two unknown parameters $\theta$ and $\phi$, and priors $f_{\theta\mid\phi}$ and $f_{\phi}$, an increase in $\text{Var}\left[ \theta\mid\phi\right]$ should result in an increase in $\text{Var}\left[ \phi \mid y_1,\dots,y_n \right]$. Does anyone know of a general result like this?

Answer 1

While this is probably true for most relevant examples an easy counterexample exist in that $\theta$ and $\phi$ can be independent in both prior and likelihood.

I have constructed a realistic example in R with a linear model where $\phi$ is the intercept and $\theta$ the slope with respect to a mean centered, i.e. orthogonal to the intercept, variable $x$.

I use two different priors: $$ \phi \sim \mathcal{N}(0, 0.1)\mathrel{\unicode{x2AEB}} \theta \sim \mathcal{N}(0, 0.1) $$ $$ \phi \sim \mathcal{N}(0, 0.1)\mathrel{\unicode{x2AEB}} \theta \sim \mathcal{N}(0, 2) $$ $\sigma$ uses default priors of arm::bayesglm. $\sqrt{\text{Var}[\phi]}$ can be found as standard error in the summary.

library(arm)

n <- 100
x <- scale(runif(n))
y <- rnorm(n)

m1 <- bayesglm(y~x, prior.scale = 0.1,  prior.scale.for.intercept=0.1)
summary(m1)
m2 <- bayesglm(y~x, prior.scale = 2,  prior.scale.for.intercept=0.1)
summary(m2)

result:

> summary(m1)

Call:
bayesglm(formula = y ~ x, prior.scale = 0.1, prior.scale.for.intercept = 0.1)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.13195    0.09295  -1.420    0.159
x            0.14130    0.09371   1.508    0.135
> summary(m2)

Call:
bayesglm(formula = y ~ x, prior.scale = 2, prior.scale.for.intercept = 0.1)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.13193    0.09286  -1.421   0.1586  
x            0.18414    0.10688   1.723   0.0881 .

In fact the $\sqrt{\text{Var}[\phi]}$ became slightly smaller, which might mean that they are not perfectly independent, but definitely a counterexample. It could also just be an issue with the implementation. Simulation shows this result to always be the case with small decreases from $10^{-5}$ up to $10^{-4}$:

sim_result <- t(replicate(1000, {
  x <- scale(runif(n))
  y <- rnorm(n)
  m1 <- bayesglm(y~x, prior.scale = 0.1,  prior.scale.for.intercept=0.1)
  m2 <- bayesglm(y~x, prior.scale = 2,  prior.scale.for.intercept=0.1)
  return(c(tight_prior = summary(m1)$coefficients[1, 2], wide_prior = summary(m2)$coefficients[1, 2]))
}))

hist(sim_result[, "tight_prior"] -sim_result[, "wide_prior"])

#> summary(sim_result[, "tight_prior"] -sim_result[, "wide_prior"])
#     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#2.360e-05 3.687e-05 4.586e-05 5.016e-05 6.141e-05 9.624e-05 

Answer 2

Here is a counter-example.

Staying in your example of the standard gaussian model and its gamma-normal conjugate prior:

1. let's do a change of variable to $\sigma^2$ instead of $1/\sigma^2$.

2. We still have the property that decreasing $\kappa_0$ increases the prior conditional variance of $\mu$.µ

3. The posterior is an inverse gamma, with the same parameters that you gave:

   $$ \begin{align} \alpha &= \frac{\nu_n}{2} \\ \beta &= \frac{\nu_n}{2} \sigma^2_n \\ \sigma^2_n &= \frac{1}{\nu_n}\left( \nu_0\sigma_0^2 + (n-1)s^2 + \frac{\kappa_0 n}{\kappa_0 + n} ( \bar{y}-\mu_0 )^2 \right) \end{align} $$

4. In particular, $\beta$ decreases with $\kappa_0$, $\alpha$ is fixed.

5. But the variance of $\sigma^2$ is $\beta^2 / (\alpha-1)^2 (\alpha-2)$ and is thus also decreasing.

The key thing to observe is that the posterior of $1 / \sigma^2$ is a gamma distribution and that $\kappa_0$ only affects the beta parameter. That's only a scale parameter. You make the Gamma distribution tighter by modifying $\alpha$ which is only affected by $n$ and $\nu_0$.