Hello!

I’ve read this thread, and I think it adresses a problem I encounter when I try to fit multivariate regressions.

Some components of L_Omega, the cholesky factor of my correlation matrix, go to zero and my chains do not explore the posterior space (rhat = NAN).

Is this behaviour due to the cholesky prior, or did I make mistakes in my models?

Here is my model :

require(rstan)

my_model <- stan_model(model_code =

"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(1);

L_sigma ~ cauchy(0, 2.5);

Y ~ multi_normal_cholesky(mu, L_Sigma);

}

generated quantities{

matrix[K, K] L_Sigma;

L_Sigma = diag_pre_multiply(L_sigma, L_Omega);

}")

And here is a process generating a dataset similar to mine :

require(mvtnorm)

## Number of plot sampled

n_plot = 20

## Number of individuals sampled in each plots

n_ind = 40

## Matrix of predictors

X <- rmvnorm(n_plot*n_ind, mean = c(0,0), sigma = matrix(c(1,0.07,0.07,1), ncol = 2))

## Matrix of regression coefficients

Beta <- matrix(c(

0.42, -0.04,

0.65,-0.05,

0.51,-0.03,

-0.70, -0.12,

-0.38, 0.10,

0.11,-0.01,

-0.04,-0.13

), ncol = 2, byrow = T)

## Compute the deterministic mean for each observation

mu <- matrix(NA, nrow = nrow(X), ncol = nrow(Beta))

for(n in 1:nrow(X)) mu[n,] <- Beta %*% X[n,]

## Create correlation matrix

Sigma <- matrix(c(1.00, 0.67, -0.09, -0.20, -0.15, -0.36, 0.25,

0.67, 1.00, 0.23, -0.49, -0.18, -0.13, 0.04,

-0.09, 0.23, 1.00, -0.63, -0.43, 0.43, -0.52,

-0.20, -0.49, -0.63, 1.00,0.47,-0.17, 0.15,

-0.15, -0.18, -0.43, 0.47, 1.00, 0.01, 0.17,

-0.36, -0.13, 0.43, -0.17, 0.01, 1.00, -0.29,

0.25, 0.04, -0.52, 0.15, 0.17, -0.29, 1.00), ncol = 7)

## Compute observations based on multivariate normal distribution

Y <- matrix(NA, nrow = nrow(mu), ncol = ncol(mu))

for(n in 1:nrow(mu)) Y[n,] <- rmvnorm(1, mu[n,], sigma = Sigma)

pairs(cbind(Y,X))

## Create objects to pass to stan

## Dimensions

N <- nrow(Y) ## Number of observations

K <- ncol(Y) ## Number of response variables

J <- ncol(X) ## Number of predictors

Data_stan <- list(Y = Y, X = X,

N = N, K = K, J = J)

## Fit the model

mod_reparam <- sampling(my_model,

data = Data_stan,

chains = 2, cores = 2, iter = 1000)

# Check diagnostics

check_all_diagnostics(mod_reparam)

print(mod_reparam)

I do not have that kind of trouble if I do this analysis directly with the non-reparametred Sigma, that is, without doing Cholesky decomposition.

Thank you!

Lucas