# Bayesian Multiple regression with Stan

**URL:** <https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729>\
**Category:** Modeling\
**Created:** [September 25, 2024, 4:24pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729 "2024-09-25T16:24:35Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lazhar-Belkhiri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lazhar-belkhiri/32/19057_2.png) [@Lazhar-Belkhiri](https://discourse.mc-stan.org/u/Lazhar-Belkhiri)\
**Post date:** [September 25, 2024, 4:24pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/1 "2024-09-25T16:24:35Z")

</div>

Could you please give me your sussing and corrections of this Stan code  
When I check this code I received this error:

```no-highlight
Error in stanc(filename, allow_undefined = TRUE) : 0
Semantic error in 'string', line 36, column 8 to column 46:
   -------------------------------------------------
    34: }
    35: // Likelihood
    36: y[n] ~ multi_normal_cholesky(mu, L_e);
                 ^
    37: }
    38: }
   -------------------------------------------------

Ill-typed arguments supplied to function 'multi_normal_cholesky':
(real, vector, matrix)
Available signatures:
(vector, array[] row_vector, matrix) => real
  The first argument must be vector but got real
(vector, row_vector, matrix) => real
  The first argument must be vector but got real
(vector, array[] vector, matrix) => real
  The first argument must be vector but got real
(vector, vector, matrix) => real
  The first argument must be vector but got real
(row_vector, array[] row_vector, matrix) => real
  The first argument must be row_vector but got real
(Additional signatures omitted)

```

```stan
data {
    int<lower=0> N; // Total number of valid observations
    int<lower=0> K; // Number of predictors (tree-rings)
    int<lower=0> M; // Number of unique streamflow stations
    matrix[N, K] X; // Predictor matrix (only valid rows)
    vector[N] y; // Flattened response vector (non-missing responses)
    int<lower=1, upper=M> index[N]; // Index for streamflow stations (1 to M)
}

transformed data {
    real<lower=0> df_e; // Degrees of freedom for the inverse Wishart distribution
    df_e = M + 1; // Degrees of freedom for the inverse Wishart distribution
}

parameters {
    vector[M] alphas; // Mean for each streamflow station
    matrix[K, M] betas; // Coefficients for predictors, varying by station
    cov_matrix[M] sigmas_e; // Covariance matrix for residuals
}

transformed parameters {
    cholesky_factor_cov[M] L_e; // Cholesky factor of the covariance matrix
    L_e = cholesky_decompose(sigmas_e); // Cholesky decomposition of the covariance matrix
}

model {
    // Priors
    alphas ~ normal(0, 100);
    to_vector(betas) ~ normal(0, 100);
    sigmas_e ~ inv_wishart(df_e, diag_matrix(rep_vector(1, M)));

    // Likelihood
    for (n in 1:N) {
        vector[M] mu; // Initialize mean vector for the current observation
        mu = alphas; // Start with alphas

        // Compute the contribution from the predictors
        for (m in 1:M) {
            mu[m] += dot_product(X[n], betas[, m]); // Each entry of mu
        }

        // Likelihood
        y[n] ~ multi_normal_cholesky(mu, L_e);
    }
}

```

[edit: escaped code]

---

<div class="post-metadata">

**Author:** ![simonbrauer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simonbrauer/32/5431_2.png) [@simonbrauer](https://discourse.mc-stan.org/u/simonbrauer)\
**Post date:** [September 25, 2024, 5:33pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/2 "2024-09-25T17:33:59Z")

</div>

Welcome to the Discourse! The issue you’re having is that you’re modeling a single value (`y[n]`) using a distribution that needs a vector of values.

> [@Lazhar-Belkhiri](#):
>
> vector[N] y; // Flattened response vector (non-missing responses)

> [@Lazhar-Belkhiri](#):
>
> `y[n] ~ multi_normal_cholesky(mu, L_e);`

Maybe you meant to pass `y` as a matrix?

At least, that’s what the error is telling you; you may have other bugs that are not immediately apparently.

---

<div class="post-metadata">

**Author:** ![Lazhar-Belkhiri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lazhar-belkhiri/32/19057_2.png) [@Lazhar-Belkhiri](https://discourse.mc-stan.org/u/Lazhar-Belkhiri)\
**Post date:** [September 25, 2024, 8:05pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/3 "2024-09-25T20:05:42Z")

</div>

Dear Simonbrauer  
Thank you so much for your quick reply  
Would you please explain more your suggestion?  
Best regards

---

<div class="post-metadata">

**Author:** ![simonbrauer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simonbrauer/32/5431_2.png) [@simonbrauer](https://discourse.mc-stan.org/u/simonbrauer)\
**Post date:** [September 25, 2024, 11:12pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/4 "2024-09-25T23:12:44Z")

</div>

Sure. The first part of the error is saying that you tried to call the function `multi_normal_cholesky` with a three arguments: a (single value) real, a vector, and a matrix. However, there is no implementation of that function that takes those three kinds of arguments.

> [@Lazhar-Belkhiri](#):
>
> Ill-typed arguments supplied to function ‘multi\_normal\_cholesky’:  
> (real, vector, matrix)

The second part of the error is listing out all of the possible type signatures that are implement. For example,

> [@Lazhar-Belkhiri](#):
>
> (vector, vector, matrix) =\> real

means that you can call `multi_normal_cholesky` with a vector, another vector, and a matrix. This makes sense, since `multi_normal_cholesky` is a multivariate distribution, and so the first argument (`y` in your case) needs to be the same length as the second argument (`mu`).

There are multiple ways of doing this. But I’ll just copy your code and illustrate what I had in mind. Note how `y` is now an array of vectors, so `y[n]` refers to a vector of length `M`.

```stan
data {
  int<lower=0> N; // Total number of valid observations
  int<lower=0> K; // Number of predictors (tree-rings)
  int<lower=0> M; // Number of unique streamflow stations
  matrix[N, K] X; // Predictor matrix (only valid rows)
  array[N] vector[M] y; // Flattened response vector (non-missing responses)
  int<lower=1, upper=M> index[N]; // Index for streamflow stations (1 to M)
}
transformed data {
  real<lower=0> df_e; // Degrees of freedom for the inverse Wishart distribution
  df_e = M + 1; // Degrees of freedom for the inverse Wishart distribution
}
parameters {
  vector[M] alphas; // Mean for each streamflow station
  matrix[K, M] betas; // Coefficients for predictors, varying by station
  cov_matrix[M] sigmas_e; // Covariance matrix for residuals
}
transformed parameters {
  cholesky_factor_cov[M] L_e; // Cholesky factor of the covariance matrix
  L_e = cholesky_decompose(sigmas_e); // Cholesky decomposition of the covariance matrix
}
model {
  // Priors
  alphas ~ normal(0, 100);
  to_vector(betas) ~ normal(0, 100);
  sigmas_e ~ inv_wishart(df_e, diag_matrix(rep_vector(1, M)));
  
  // Likelihood
  for (n in 1:N) {
      vector[M] mu; // Initialize mean vector for the current observation
      mu = alphas; // Start with alphas
  
      // Compute the contribution from the predictors
      for (m in 1:M) {
          mu[m] += dot_product(X[n], betas[, m]); // Each entry of mu
      }
  
      // Likelihood
      y[n] ~ multi_normal_cholesky(mu, L_e);
  }
}

```

As a side-note, you can surround your code in ```stan ``` to make it highlight syntax. For example.

```stan  
data{  
vector[N] Y;  
}  
```

becomes

```stan
data{
   vector[N] Y;
}

```

---

<div class="post-metadata">

**Author:** ![Lazhar-Belkhiri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lazhar-belkhiri/32/19057_2.png) [@Lazhar-Belkhiri](https://discourse.mc-stan.org/u/Lazhar-Belkhiri)\
**Post date:** [September 26, 2024, 10:42am UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/5 "2024-09-26T10:42:19Z")

</div>

Thank you so much for your help

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 26, 2024, 9:59pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/6 "2024-09-26T21:59:50Z")

</div>

Aside from the coding issue, which @simonbrauer correctly diagnosed, there’s a computational issue in this model.

> [@Lazhar-Belkhiri](#):
>
> ```no-highlight
> parameters {
> cov_matrix[M] sigmas_e; // Covariance matrix for residuals
> }
> transformed parameters {
> cholesky_factor_cov[M] L_e = cholesky_decompose(sigmas_e); 
> }
> model {
> y[n] ~ multi_normal_cholesky(mu, L_e);
> }
> 
> ```

You’re better off here just using regular `multi_normal` and not doing the Cholesky factorization yourself. The multi-normal does it already.

Given that you’re using a unit matrix in the inverse Wishart prior (i.e., not something informative), I would suggest instead decomposing the prior into an LKJ prior on the correlation matrix (shrink to unit correlation) plus an independent prior on the scales. Here’s a description of why:

[http://www.stat.columbia.edu/~gelman/research/unpublished/Visualization.pdf](http://www.stat.columbia.edu/~gelman/research/unpublished/Visualization.pdf)

Aside from the applied modeling motivation, it’s a whole lot faster. The Cholesky factor operations are quadratic whereas the full covariance operations are cubic and with this approach you don’t ever need to create the full covariance matrix. There’s a description in the _User’s Guide_ chapter on regression of how to code it in Stan. We should really introduce Cholesky-parameterized inverse Wisharts but we don’t use them a lot ourselves given @bgoodri’s paper linked above.

> [@simonbrauer](#):
>
> ```no-highlight
> matrix[N, K] X;
> ...
> matrix[K, M] betas;
> ...
> for (m in 1:M) {
> mu[m] += dot_product(X[n], betas[, m]); // Each entry of mu
> }
> 
> ```

this can probably be sped up with some of our linear algebra and compound operations like `dot_product_rows`. At the very least, if you’re using in this form, it’ll be much more efficient to declare things this way because no copying will be required to extract a row of `X` and column of `betas`.

```stan
array[N] row_vector[K] X;
array[M] vector[K] betas_t;
for (m in 1:M) {
  mu[m] += X[n] * betas_t[m];
}

```

I used `betas_t` because it’s the transpose of your original betas. `X` keeps the same shape, just a different declaration.

```no-highlight

```

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.mc-stan.org/u/spinkney)\
**Post date:** [September 27, 2024, 1:23am UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/7 "2024-09-27T01:23:33Z")

</div>

> [@Bob\_Carpenter](#):
>
> We should really introduce Cholesky-parameterized inverse Wisharts

We do have this. [Covariance Matrix Distributions](https://mc-stan.org/docs/functions-reference/covariance_matrix_distributions.html#inverse-wishart-cholesky-distribution)

In fact, you helped with the docs and review!

---

<div class="post-metadata">

**Author:** ![Aida](https://avatars.discourse-cdn.com/v4/letter/a/77aa72/32.png) [@Aida](https://discourse.mc-stan.org/u/Aida)\
**Post date:** [September 27, 2024, 1:52pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/8 "2024-09-27T13:52:58Z")

</div>

The error is due to mu being a real instead of a vector[M]. Change it to vector[M] and make sure L\_e is correctly set as a Cholesky factor. That should fix the issue!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 27, 2024, 9:36pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/9 "2024-09-27T21:36:51Z")

</div>

> [@spinkney](#):
>
> We do have this.

I knew it was a plan, but I hadn’t realized it’d gotten all the way through to being released! This is super cool. I’ve wanted to have this for years!

---

<div class="post-metadata">

**Author:** ![Lazhar-Belkhiri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lazhar-belkhiri/32/19057_2.png) [@Lazhar-Belkhiri](https://discourse.mc-stan.org/u/Lazhar-Belkhiri)\
**Post date:** [September 30, 2024, 3:18pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/10 "2024-09-30T15:18:12Z")

</div>

I convert Y to a data frame, remove NAs, and flatten  
In this example, I updated Stan Model with Cholesky Decomposition  
But still not working I received this errors:  
“SAMPLING FOR MODEL ‘anon\_model’ NOW (CHAIN 1).  
Chain 1: Unrecoverable error evaluating the log probability at the initial value.  
Chain 1: Exception: multi\_normal\_cholesky\_lpdf: Size of random variable (123) and rows of covariance parameter (4) must match in size (in ‘anon\_model’, line 43, column 4 to column 39)  
[1] “Error : Exception: multi\_normal\_cholesky\_lpdf: Size of random variable (123) and rows of covariance parameter (4) must match in size (in ‘anon\_model’, line 43, column 4 to column 39)””  
###########################################################################  
Here the different steps:  
Y ← readRDS(‘data/data.ann.stramflow.rds’) %\>% .[-c(1,5:36)]%\>% as.matrix()  
X ← readRDS(‘data/data.treerings.rds’) %\>% .[-c(1,4:27)] %\>%  
mutate\_all(~ifelse(is.na(.x), mean(.x, na.rm = TRUE), .x)) %\>% as.matrix()  
nt ← nrow(Y)  
ns ← ncol(Y)  
np ← ncol(X)

# Convert Y to a data frame, remove NAs, and flatten

y\_na ← Y %\>%  
as.data.frame() %\>%  
mutate(row = row\_number()) %\>%  
pivot\_longer(-row, names\_to = “column”, values\_to = “value”) %\>%  
filter(!is.na(value)) %\>%  
mutate(column = match(column, colnames(Y))) %\>% # Convert column names to column numbers  
set\_colnames(c(‘year’,‘site’,‘flow’))

nb = nrow(y\_na)

# Data list for Stan

stan\_data ← list(nb = nb, ns = ns, nt = nt, np = np, x = X, y = y\_na$flow, site = y\_na$site, year = y\_na$year, Wp\_e = diag(ns))

stan\_model\_code ← ‘/Users/samir/Documents/Paleo\_hydrology\_Algeria/Bayesian\_model\_rcode/stan\_np.stan’  
stan\_model ← stan(file = stan\_model\_code, data = stan\_data, #init = init\_list,  
iter = 1000, chains = 1,  
cores = 8, control = list(adapt\_delta = .9, max\_treedepth = 15))  
###########################################################################  
Here the Stan code:  
data {  
int\<lower=0\> nb; // Total number of observations (non-missing)  
int\<lower=0\> ns; // Number of streamflow stations  
int\<lower=0\> nt; // Number of years  
int\<lower=0\> np; // Number of predictors (tree-rings)

```
matrix[nt, np] x; // Predictor matrix (no missing values)
vector[nb] y; // Response vector (non-missing values only)

int<lower=1, upper=ns> site[nb]; // Index for streamflow stations (1 to ns)
int<lower=1, upper=nt> year[nb]; // Index for years (1 to nt)

matrix[ns, ns] Wp_e; // Scale matrix for the inverse-Wishart prior

```

}

transformed data {  
real\<lower=0\> df\_e; // Degrees of freedom for the inverse-Wishart distribution  
df\_e = ns + 1;  
}

parameters {  
vector[ns] alphas; // Mean for each streamflow station  
matrix[np, ns] betas; // Coefficients for predictors, varying by station

```
cov_matrix[ns] sigmas_e; // Covariance matrix for residuals between stations

```

}

transformed parameters {  
vector[nb] mu; // Predicted means for non-missing observations

```
for (n in 1:nb) {
    // Calculate the predicted mean for each non-missing observation
    mu[n] = alphas[site[n]] + dot_product(x[year[n], ], betas[, site[n]]);
}

```

}

model {  
matrix[ns, ns] L\_e; // Cholesky factor of the covariance matrix

```
// Priors
alphas ~ normal(0, 100); // Priors for station means
to_vector(betas) ~ normal(0, 100); // Priors for station-specific coefficients
sigmas_e ~ inv_wishart(df_e, Wp_e); // Prior for the covariance matrix

// Cholesky decomposition of the covariance matrix
L_e = cholesky_decompose(sigmas_e);

// Likelihood
y ~ multi_normal_cholesky(mu, L_e); // Ensure mu is correct size and matches the covariance

```

}

---

<div class="post-metadata">

**Author:** ![simonbrauer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simonbrauer/32/5431_2.png) [@simonbrauer](https://discourse.mc-stan.org/u/simonbrauer)\
**Post date:** [October 1, 2024, 1:41pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/11 "2024-10-01T13:41:40Z")

</div>

> [@Lazhar-Belkhiri](#):
>
> Here the Stan code:

I’m copying your model below but formatted for better legibility using the method I described above.

> [@simonbrauer](#):
>
> ```stan  
> data{  
> vector[N] Y;  
> }  
> ```

```stan
data {
int<lower=0> nb; // Total number of observations (non-missing)
int<lower=0> ns; // Number of streamflow stations
int<lower=0> nt; // Number of years
int<lower=0> np; // Number of predictors (tree-rings)

matrix[nt, np] x; // Predictor matrix (no missing values)
vector[nb] y; // Response vector (non-missing values only)

int<lower=1, upper=ns> site[nb]; // Index for streamflow stations (1 to ns)
int<lower=1, upper=nt> year[nb]; // Index for years (1 to nt)

matrix[ns, ns] Wp_e; // Scale matrix for the inverse-Wishart prior
}

transformed data {
real<lower=0> df_e; // Degrees of freedom for the inverse-Wishart distribution
df_e = ns + 1;
}

parameters {
vector[ns] alphas; // Mean for each streamflow station
matrix[np, ns] betas; // Coefficients for predictors, varying by station

cov_matrix[ns] sigmas_e; // Covariance matrix for residuals between stations
}

transformed parameters {
vector[nb] mu; // Predicted means for non-missing observations

for (n in 1:nb) {
    // Calculate the predicted mean for each non-missing observation
    mu[n] = alphas[site[n]] + dot_product(x[year[n], ], betas[, site[n]]);
}
}

model {
matrix[ns, ns] L_e; // Cholesky factor of the covariance matrix

// Priors
alphas ~ normal(0, 100); // Priors for station means
to_vector(betas) ~ normal(0, 100); // Priors for station-specific coefficients
sigmas_e ~ inv_wishart(df_e, Wp_e); // Prior for the covariance matrix

// Cholesky decomposition of the covariance matrix
L_e = cholesky_decompose(sigmas_e);

// Likelihood
y ~ multi_normal_cholesky(mu, L_e); // Ensure mu is correct size and matches the covariance
}

```

> [@Lazhar-Belkhiri](#):
>
> Chain 1: Exception: multi\_normal\_cholesky\_lpdf: Size of random variable (123) and rows of covariance parameter (4) must match in size (in ‘anon\_model’, line 43, column 4 to column 39)

> [@Lazhar-Belkhiri](#):
>
> `vector[nb] y; // Response vector (non-missing values only)`

> [@Lazhar-Belkhiri](#):
>
> matrix[ns, ns] L\_e; // Cholesky factor of the covariance matrix

> [@Lazhar-Belkhiri](#):
>
> `y ~ multi_normal_cholesky(mu, L_e); // Ensure mu is correct size and matches the covariance`

This error is telling you that `y` has 123 elements (`nb`) while `L_e` only has 4 (`ns`). I think your approach of flattening your data into one long vector is making this more difficult. I would instead pass `y` as an array of vectors or a matrix.

> [@simonbrauer](#):
>
> `array[N] vector[M] y;`

> [@simonbrauer](#):
>
> ```no-highlight
> for (n in 1:N) {
> vector[M] mu; // Initialize mean vector for the current observation
> mu = alphas; // Start with alphas
>   
> // Compute the contribution from the predictors
> for (m in 1:M) {
> mu[m] += dot_product(X[n], betas[, m]); // Each entry of mu
> }
>   
> // Likelihood
> y[n] ~ multi_normal_cholesky(mu, L_e);
> }
> 
> ```

---

<div class="post-metadata">

**Author:** ![Lazhar-Belkhiri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lazhar-belkhiri/32/19057_2.png) [@Lazhar-Belkhiri](https://discourse.mc-stan.org/u/Lazhar-Belkhiri)\
**Post date:** [October 1, 2024, 10:29pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/12 "2024-10-01T22:29:54Z")

</div>

I used Y as one long vector with length nb=123 because I removed the missing data (NAs)

---

<div class="post-metadata">

**Author:** ![simonbrauer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simonbrauer/32/5431_2.png) [@simonbrauer](https://discourse.mc-stan.org/u/simonbrauer)\
**Post date:** [October 2, 2024, 7:51pm UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/13 "2024-10-02T19:51:34Z")

</div>

You can estimate these models with missing data, but it requires slicing the mean vector and covariance matrix at each step, which can get complex. I’m not even sure if that would work for the Cholesky decomposition version (someone can clarify here), so you might have to switch to ` y[n] ~ multi_normal(mu, Sigma)`

With that in mind, I’d recommend that you get your model running with complete data and then figuring out the slicing later. So restructure `Y` either as a matrix (`matrix[N,M] Y`) or an array of vectors (`array[N] vector[M] Y`) as described above. You may still have other errors to resolve (I haven’t looked closely at the rest of your model), but that should get you past this initial error.

---

<div class="post-metadata">

**Author:** ![Lazhar-Belkhiri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lazhar-belkhiri/32/19057_2.png) [@Lazhar-Belkhiri](https://discourse.mc-stan.org/u/Lazhar-Belkhiri)\
**Post date:** [October 3, 2024, 11:18am UTC](https://discourse.mc-stan.org/t/bayesian-multiple-regression-with-stan/36729/14 "2024-10-03T11:18:26Z")

</div>

> [@simonbrauer](#):
>
> equires slicing the mean vector

Thank you so much