I am looking for a teaching example of a multivariate (not bivariate) implementation of Metropolis-Hastings for MCMC in R. I know several packages implement the algorithm more generally, but the code is difficult to follow and typically includes all sorts of other things besides this particular example. How do I implement Metropolis-Hastings for Bayesian multivariate regression like in this bivariate example?

  • 4
    $\begingroup$ I'm not sure what you mean by multivariate / bivariate - the MCMC in the blog post by me to which you link samples from three parameters, so it's "trivariate"? If you want one more dimension, simply add one more parameter to the regression. $\endgroup$ – Florian Hartig Apr 21 '15 at 20:32
  • $\begingroup$ As I understand it, the model has in the example has a Y variable and an x1 variable in the example. It estimates the intercept, a beta coefficient for the single x1 variable, and the standard deviation for the x1 variable. If instead I had Y, x1, x2, and x3 I am unsure how to draw B1 conditional on B2 and B3 as well as how to estimate SD for each of these additional variables. $\endgroup$ – Michael Apr 21 '15 at 21:57
  • 1
    $\begingroup$ No, it estimates intercept, slope, and the standard deviation of the data generating model. The latter is the variation of the data around the regression mean, i.e. the residual variance, which is part of the likelihood. For each of this parameters, you get an uncertainty via the posterior. If you want to extend this example, see my answer below. $\endgroup$ – Florian Hartig Apr 22 '15 at 8:16

If you want to extend the example you link to to a multivariate regression, take the code as it is and:

  1. Add one more predictor in the code chunk generating the data
  2. Add one more parameter in likelihood, as in pred = a1*x1 + a2*x2 + b
  3. Add the additional parameter in the prior specification
  4. Adjust the MCMC and plots to deal with 4 instead of 3 parameters
  • $\begingroup$ So in each iteration the likelihood function is p(y_actual | y_hat) such that it does not matter how many X covariates there are, you simply get a predicted value y_hat_i for iteration i and use that to estimate the likelihood? $\endgroup$ – Michael Apr 22 '15 at 22:26
  • $\begingroup$ @Michael check the two lines of code pred = a*x + b; singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T) as you see there is one x, but you can use possibly any number of covariates here to create pred i.e. the mean of y's distribution. You do not compare actual y's but it's densities (unless you use ABC and some loss function, but it's a different topic). $\endgroup$ – Tim Apr 23 '15 at 6:15
  • $\begingroup$ @Tim Sorry for the very late reply. There is something confusing to me. I thought that the regression coefficients are correlated/dependent on each other so we should NOT add them individually in the prior. I thought that we should represent the coefficients as multivariate distribution with co-variance matrix. Am I right? $\endgroup$ – floyd Aug 25 at 12:45
  • 1
    $\begingroup$ @floyd if you have prior information, you can use such prior. Not assuming correlation a priori doesn't prohibit correlation a posteriori. $\endgroup$ – Tim Aug 25 at 13:55

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.