# Metropolis-Hastings MCMC for Bayesian Regression in R

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?

• 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. Apr 21 '15 at 20:32
• 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. Apr 21 '15 at 21:57
• 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. 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
• 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? Apr 22 '15 at 22:26
• @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).