# 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.