I am actually exploring multivariate latent variable models, and I am wondering how to compute information criteria (loo, waic) for such models. I have an hard time to find references dealing with such questions, except this one, which presents criteria I never heard of.

Specifically, I can’t figure out how to write the likelihood computation in the generated quantity block, and the type and dimensions of output (array of matrices?).

I know that brms can compute such criteria, but the likelihood computation is made after model fitting, so I can’t use it to produce an exemple of correct generated quantity block.

There exist many multivariate latent variable models, could you be more specific? Also if you would show Stan code for your model, it would be easier to say how to write the likelihood computation in the generated quantity block.

generated quantities{
matrix[P,P] cov_lambda;
matrix[N,P] log_lik1;
vector[N*P] log_lik;
cov_lambda = multiply_lower_tri_self_transpose(lambda); //Create the covariance matrix
//Compute the likelihood for each observation
for(n in 1:N){
for(p in 1:P){
log_lik1[n,p] = poisson_lpmf(Y[n,p] | mu[n,p]);
}
}
log_lik = to_vector(log_lik1); //Tranform in a vector usable by loo package
}

I get this output : Computed from 3000 by 336 log-likelihood matrix

Estimate SE
elpd_loo -779.0 47.2
p_loo 284.8 27.0
looic 1558.0 94.3
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 124 36.9% 1
(0.5, 0.7] (ok) 75 22.3% 0
(0.7, 1] (bad) 90 26.8% 0
(1, Inf) (very bad) 47 14.0% 0
See help('pareto-k-diagnostic') for details.
There were 29 warnings (use warnings() to see them)
warnings()
Messages d'avis :
1: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
2: In log(z) : production de NaN
3: In log(z) : production de NaN
4: In log(z) : production de NaN
...

When I did what you suggested, that is

generated quantities{
matrix[P,P] cov_lambda;
matrix[N,P] log_lik1;
vector[N*P] log_lik;
cov_lambda = multiply_lower_tri_self_transpose(lambda); //Create the covariance matrix
//Compute the likelihood for each observation
for(n in 1:N){
for(p in 1:P){
log_lik1[n,p] = poisson_lpmf(Y[n,p] | mu[n]);
}
}
log_lik = to_vector(log_lik1); //Tranform in a vector usable by loo package
}

The loo computation was impossible, and returned NA only.

Sorry, I didn’t notice that mu is a matrix since you had

Y[i,] ~ poisson(mu[i]);

instead of

Y[i,] ~ poisson(mu[i,]);

If Y has 336 observations, and mu has 336 unknowns then this is possible, and implies that your prior is weak (too flexible) or data is overdispersed compared to Poisson.