Multivariate outcomes with separate independent variables

This is the example of a multivariate model from the users guide :

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;
  cholesky_factor_corr[K] L_Omega;
  vector<lower=0>[K] L_sigma;
}
model {
  vector[K] mu[N];
  matrix[K, K] L_Sigma;

  for (n in 1:N)
    mu[n] = beta * x[n];

  L_Sigma = diag_pre_multiply(L_sigma, L_Omega);

  to_vector(beta) ~ normal(0, 5);
  L_Omega ~ lkj_corr_cholesky(4);
  L_sigma ~ cauchy(0, 2.5);

  y ~ multi_normal_cholesky(mu, L_Sigma);
}

My model is something like this, but for each of my K outcomes Y, I have a baseline value Yb. So my model looks something like this, where the new bits are **ed

data {
  int<lower=1> K;
  int<lower=1> J;
  int<lower=0> N;
  vector[J] x[N];
  vector[K] y[N];
  **vector[K] yb[N];**
}
parameters {
  matrix[K, J] beta;
  **vector[K] beta_yb**;
  cholesky_factor_corr[K] L_Omega;
  vector<lower=0>[K] L_sigma;
}
model {
  vector[K] mu[N];
  matrix[K, K] L_Sigma;

  for (n in 1:N)
    mu[n] = beta * x[n] + yb[n]*beta_yb;      // ???? 

  L_Sigma = diag_pre_multiply(L_sigma, L_Omega);

  to_vector(beta) ~ normal(0, 5);
  L_Omega ~ lkj_corr_cholesky(4);
  L_sigma ~ cauchy(0, 2.5);

  y ~ multi_normal_cholesky(mu, L_Sigma);
}

I am struggling with what goes before the “???”. I’ve tried every combination I can think of to produce a column vector that looks like yb[n]*beta_yb, but they all give errors, either when parsing or when trying to initialize any chain. I have tried declaring both as every combination of matrix and vector or row_vector that makes sense to combine, but nothing will get past the first iteration of any chain.

I also tried doing this manually,

vector[K] product[N];
for (n in 1:N) {
   for (k in 1:K) {
     product[n,k]=beta_yb[k]*yb[n,k];
   }
}
 for (n in 1:N)
    mu[n] = beta * x[n] + product[n];

with similar (non) results. Does anyone have any thoughts on how to do this?

I believe you are looking for element-wise multiplication, indicated by .*

yb[n] .* beta_yb

That said, your double-loop version looks fine to me. So the problem may lie elsewhere.

Thanks so much for the suggestion. I thought I had tried that but tried it again anyway, and get the same error as always, which I should have included in my original post:

SAMPLING FOR MODEL 'model' NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: add: m1 ((1201775849648561985, 1)6, 1) must match in size (in 'string', line 45, column 4 to column 46)
[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                    
[2] "  Exception: add: m1 ((1201775849648561985, 1)6, 1) must match in size (in 'string', line 45, column 4 to column 46)"
[1] "error occurred during calling the sampler; sampling not done"

Column 46 is where yb[n].*beta_yb starts; if I delete everything starting with the last “+” the model runs fine. So it seems that something is the wrong size.

When I searched for this error message, I learn that it is a dimension problem - it seems like whatever I putting there (either my double loop solution or your suggestion), I’m not getting the expected vector.

In R, I can see that yb has 6 columns, just like y, and has no missings.

I’m beginning to think it’s a data problem. I tried treating yb like x; that is, adding the baseline for all Ys to each of the regressions:

  vector[J] x[N];
  vector[K] yb[N];
...
 matrix[K, J] beta; 
 matrix[K,K] beta_yp;             
...
  vector[K] mu[N];
  for (n in 1:N) {
      mu[n] = beta * x[n] + beta_yp*yp[n];
  }

This gives the same error as above. (“must match in size”).

Yep, dumb data error.