Hi Andrew

Thank you

I did this code in multivariate linear regression but I couldn’t get the results and when I check the Stan code, the code is syntactically correct. I don’t know why doesn’t work. The code is:

#### R code

dat <- within(list(), {

y <- as.matrix(dplyr::select(df, F1, F2, F3, F4))

X <- model.matrix(~ 0 + Total_Fines + TOC + Total_solids, data = df)

N <- nrow(y)

K <- ncol(y)

P <- ncol(X)

alpha_loc <- rep(0, K)

alpha_scale <- rep(10, K)

beta_loc <- matrix(0, K, P)

beta_scale <- matrix(2.5, K, P)

Sigma_corr_shape <- 2

Sigma_scale_scale <- 5

})

### Stan code

data {

// multivariate outcome

int N;

int K;

vector[K] y[N];

// covariates

int P;

vector[P] X[N];

// prior

vector[K] alpha_loc;

vector[K] alpha_scale;

vector[P] beta_loc[K];

vector[P] beta_scale[K];

real Sigma_corr_shape;

real Sigma_scale_scale;

}

parameters {

// regression intercept

vector[K] alpha;

// regression coefficients

vector[P] beta[K];

// Cholesky factor of the correlation matrix

cholesky_factor_corr[K] Sigma_corr_L;

vector[K] Sigma_scale;

// student-T degrees of freedom

}

transformed parameters {

vector[K] mu[N];

matrix[K, K] Sigma;

// covariance matrix

Sigma = crossprod(diag_pre_multiply(Sigma_scale, Sigma_corr_L));

for (i in 1:N) {

for (k in 1:K) {

mu[i, k] = alpha[k] + dot_product(X[i], beta[k]);

}

}

}

model {

for (k in 1:K) {

alpha[k] ~ normal(alpha_loc[k], alpha_scale[k]);

beta[k] ~ normal(beta_loc[k], beta_scale[k]);

}

Sigma_scale ~ cauchy(0., Sigma_scale_scale);

Sigma_corr_L ~ lkj_corr_cholesky(Sigma_corr_shape);

y ~ multi_normal_cholesky(mu, Sigma);

}