I am new to stan and made 3 different stan models to fit to the same dataset of 1380 observations respectively and ranked them using loo::loo_compare(), but I am struggling with some diagnostics.

For two models I get quite small percentage of Pareto k values that are not good (bad or very bad). I have read other posts on the matter and some documentation (LOO package glossary — loo-glossary • loo), but it seems very context/model specific from what I have seen/understood.

Can the models with bad/very bad Pareto k values still be used and compared to other models, or do they absolutely need to be changed?

All three models are basically the same model with increasing complexity (adding some parameters)

Model informations:

Model 1 (has 2 parameters)

model 1 loo outpout:

Model 2 (has 157 parameters)

model 2 loo outpout:

Model 3 (has 311 parameters)

model 3 loo output:

Here is the stan code for model 3 as example:

```
data {
// Size integer
int<lower = 1> n;
int<lower = 1> n_unique; \\ 154 unique
// Vector data
vector[n] bf;
vector[n] bp;
vector[n] ap;
int unique_id[n];
}
transformed data {
vector[n] log_bf = log(bf);
vector[n] log_bp = log(bp);
vector[n] log_ap = log(ap);
}
parameters {
real mu_alp; \\ 1 parameter
vector[n_unique] alp; \\ 154 unique alp for 154 parameters
real<lower = 0> sd_alp; \\ 1 parameter
real<lower = 0> sigma; \\ 1 parameter
vector[n_unique] h_j; \\ 154 unique h_j for 154 parameters
\\ total of 311 parameters
}
model {
vector[n] mu;
vector[n] unique_factor;
// Priors:
mu_alp ~ normal(-5,2);
alp ~ normal(mu_alp, sd_alp);
sd_alp ~ exponential(3);
sigma ~ exponential(5);
h_j ~ normal(-4,2);
// Computing
unique_factor = alp[unique_id] + log_ap;
mu = unique_factor + log_bp - log1p_exp(h_j[unique_id] + unique_factor + log_bp);
log_bf ~ normal(mu, sigma);
}
generated quantities {
vector[n] unique_factor;
vector[n] log_bf_hat;
vector[n] log_lik;
vector[n] y_rep;
real<lower = 0, upper = 1> Rsq_3;
unique_factor = alp[unique_id] + log_ap;
log_bf_hat = unique_factor + log_bp - log1p_exp(h_j[unique_id] + unique_factor + log_bp);
for (i in 1:n) {
log_lik[i] = normal_lpdf(log_bf[i] | log_bf_hat[i], sigma);
y_rep[i] = normal_rng(log_bf_hat[i], sigma);
}
Rsq_3 = variance(log_bf_hat) / (variance(log_bf_hat) + square(sigma));
}
```