PSIS and LOOIC returning different Pareto k values

I am fitting a GEV distribution with a time trend in the location parameter \mu, such that

Y_t \sim \mathrm{GEV}(\mu_t, \sigma, \xi), \\ \mu_t = \mu_0 \big(1 + \Delta (t - t_c) \big), \\ \mu_0 \sim \mathcal N(0, 10), \\ \Delta \sim \mathcal N(0, 10), \\ \sigma \sim \mathrm{Exp}(3), \\ \xi_0 \sim \mathrm{Beta}(4, 4), \\ \xi := \xi_0 - 0.5,

where t_c = \lfloor N / 2 \rfloor, i.e. the index of the observation in the middle of the data set. Shown below is the model implementation in code:

functions {
  real gevt_lpdf(vector y, real mu_0, real sigma, real xi, real delta) {
    int n = rows(y);
    int c = n / 2;
    vector[n] mu;

    for (i in 1:n) { mu[i] = mu_0 * (1 + delta * (i - c)); }

    vector[n] eta = (abs(xi) < 1e-12)
                    ? exp(-(y - mu) / sigma)
                    : fmax(1e-12, pow(1 + xi * (y - mu) / sigma, -1 / xi));
    
    real lp = -n * log(sigma) + sum((1 + xi) * log(eta)) - sum(eta);

    return lp;
  }
}

data {
  int<lower = 0> n;
  vector[n] y;
}

parameters {
  real mu_0;
  real<lower = 0> sigma;
  real<lower = 0, upper = 1> xi_0;
  real delta;
}

transformed parameters {
  real xi = xi_0 - 0.5;
}

model {
  mu_0 ~ normal(0, 10);
  sigma ~ exponential(3);
  xi_0 ~ beta(4, 4);
  delta ~ normal(0, 10);

  y ~ gevt(mu_0, sigma, xi, delta);
}

generated quantities {
  vector[n] log_lik;
  for (i in 1:n)
    log_lik[i] = gevt_lpdf([y[i]]' | mu_0, sigma, xi, delta);
}

Model convergence is good, \widehat R \approx 1 and the effective sample size is larger than the number of draws. But I am having this confusing issue where the loo functions psis and loo seemingly disagree on the Pareto k values.

This is what happens when I run psis on the samples from the above model (timetrend_fit):

> timetrend_loglik <- loo::extract_log_lik(timetrend_fit)
> loo::psis(timetrend_loglik)
Computed from 40000 by 132 log-weights matrix.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

Whereas for the loo function this is the produced output:

> loo::loo(timetrend_loglik)

Computed from 40000 by 132 log-likelihood matrix.

         Estimate   SE
elpd_loo   -494.3 20.7
p_loo        15.6 13.1
looic       988.6 41.4
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     131   99.2%   23774   
   (0.7, 1]   (bad)        0    0.0%   <NA>    
   (1, Inf)   (very bad)   1    0.8%   <NA>    
See help('pareto-k-diagnostic') for details.
Warning message:
Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

There is one problematic observation w.r.t. to the Pareto k value. Is this a bug in the implementation, or is there a difference in how psis and loo estimate the Pareto k diagnostic values? What can I do to fix this?

Welcome to the Stan forum. The loo() and psis() functions use the same method for estimating the Pareto k values, but the functions take different arguments. For loo() you should give it log-likelihood values, like you did. However, the input to psis() should be log-importance-ratios, which in this context would be the negative of the log-likelihood values (when loo calls psis internally it does the negative for you). If you try that does it resolve the issue?