Gibbs sampler for a multilevel model with no predictors in R

Question

I'm working on multilevel models and want to know how they are estimated in R. For that purpose I'm reading amongst other things "Data Analysis Using Regression and Multilevel/Hierarchical Models" by A. Gelman and J. Hill to be found here. If you want to read in the book see p.395 - 400

They give the R code for Gibbs sampling. For better understanding, I want to code the Gibbs sampler by myself using their template. But there are two questions, I can't find the answers to regarding their code.

The hyperparameters are to be estimated using these formulas:

They describe each step of the iteration on p.399

When they implement the calculation of sigma.y and sigma.a they don't divide by n or J.

sigma.y.update <- function() {
   sigma.y.new <- sqrt(sum((y-a[county])^2)/rchisq(1,n-1))
   return (sigma.y.new)
}

sigma.a.update <- function() {
   sigma.a.new <- sqrt(sum((a-mu.a)^2)/rchisq(1,J-1))
   return (sigma.a.new)
}

What am I missing at this point? I hope you guys can help me on this.

Answer 1

I agree with @StephaneLaurent's comment, that the code is correct, but their description of how it works is not entirely accurate. To quote their exact description:

What they're trying to communicate here, but get the details wrong, is that for each parameter, one first computes an MLE of it, and then samples around that value stochastically. The two problems I see with the details are

1. In the code, they actually use unbiased variance estimators, otherwise the $\chi^2$ degrees of freedom should be $n$ and $J$. That is, the 18.13 and 18.15 formulae should use $\frac{1}{n-1}$ and $\frac{1}{J-1}$, respectively.
2. The sampling formula in (4), for example, should be

$$\sigma_y^2 = \dfrac{(J-1) \hat{\sigma}_y^2}{\chi_{J-1}^2}$$

The intuitive reason for this is that a $\chi^2_{J-1}$ distribution has a mean of $(J-1)$, so in order to keep the sampling unbiased w.r.t the MLE, one needs to scale it to a mean of 1. That is,

$$E[\sigma_y^2] = E[\dfrac{(J-1) \hat{\sigma}_y^2}{\chi_{J-1}^2}] = E[\dfrac{(J-1)}{\chi_{J-1}^2}]\cdot E[ \hat{\sigma}_y^2]=1\cdot \hat{\sigma}_y^2$$

In the code, then, one can let these two factors cancel out, which is why you don't see it. Or another way of thinking of it is that the averaging coefficient is contained in the $\chi^2$ term.