Multivariate Regression


I have a scenario where I have multiple correlated outcomes and I think use of Bayesian multivariate regression (in the seemingly unrelated regression framework) is appropriate. I’ve implemented the following code after encountering it in the guide. My question, though, is about the underlying mechanics of what this code is achieving, especially the ‘model’ component. I’d appreciate if someone could walk me through that.

data {
  int<lower=1> K;
  int<lower=1> J;
  int<lower=0> N;
  vector[J] x[N];
  vector[K] y[N];
parameters {
  matrix[K, J] beta;
  cov_matrix[K] Sigma;
model {
  vector[K] mu[N];
  for (n in 1:N)
    mu[n] = beta * x[n];
  y ~ multi_normal(mu, Sigma);
1 Like

The Stan model is basically an elaborate way of writing \log p(x | \theta) p(\theta) where \theta are the parameters you want to sample and x is your data. Once you have that abstraction, it’s MCMC time.

For animations of algorithms like the ones in Stan, check out the NUTS samplers here:

Michael Betancourt’s HMC paper here is the most up to date documentation on the Stan algorithms (other than the code itself).

I like Radford Neal’s paper a lot too:

The Stan manual has some info on the samplers, but the above references are probably better:


And for info on how the model is like a log density statement, look here: and here: