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:

hyperparameter formulae

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.


migrated from stackoverflow.com Feb 27 at 0:16

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ I have an answer but I voted to close your question because it belongs to Cross Validated (stats.stackechange). Do you know what is the Gibbs sampling ? Gelman & Hill's book do not provide a description of the Gibbs sampling. Each parameter is updated according to its conditional distribution given the other parameters and y. $\endgroup$ – Stéphane Laurent Feb 26 at 12:39
  • $\begingroup$ alright could ask the question in CV, but you could just give an answer on the question and next time ill ask on CV. Sounds fair? I think I got an idea on how the gibbs sampler works. It just doesn't make any sense to me, how they describe the algorithm and code different from that $\endgroup$ – Toby Feb 26 at 12:47
  • 2
    $\begingroup$ It is better to answer with LaTeX, which is available on CV but not here. They exactly code according to the algorithm they describe ("and then draw ..."). Looks like you are confounding sigma and sigma hat. $\endgroup$ – Stéphane Laurent Feb 26 at 13:12
  • $\begingroup$ Sorry, I misread yesterday. You are right, there's an error. The code is correct, the error is in the description of the algorithm. $\endgroup$ – Stéphane Laurent Feb 27 at 8:24
  • 1
    $\begingroup$ There is a family of prior distributions which yields nice full conditional posterior distributions for the Gibbs algorithm. You can see the formulas here. The prior distribution used by Gelman & Hill belongs to this family. $\endgroup$ – Stéphane Laurent Mar 9 at 15:07

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:

enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.