When and how to use approximate leave-future-out cross-validation on hierarchical time series Stan model

Question

I am fitting a hierarchical state space AR(1) model in Stan and am struggling use common model evaluation metrics on the model output. Computing the WAIC or using loo_cv in the loo package give warnings as described in this post: Warnings during WAIC computation: how to proceed?. Specifically:

loo::waic(log_lik)
> (98.6%) p_waic estimates greater than 0.4. We recommend trying loo instead

loo::loo(log_lik)
> Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

The suggestion is to implement leave future out cross validation, following the example in this vignette: https://mc-stan.org/loo/articles/loo2-lfo.html However, this example uses a model built in brms, and my model needs to hierarchically incorporate multiple time series of varying lengths, which I do not think can be done in brms.

My questions are:

1. Am I correct that the warnings I get from using loo indicate that I need to use lfo cross validation for my model, or is it possible I have made a mistake in how I am specifying the log likelihood?
2. How would I modify the approximate lfo approach to work with the following Stan model?

data {
    int<lower=0> N;             // total number of observations
    int<lower=1> K;             // number of covariates
    int<lower=1> S;             // number of sites
    int<lower=1,upper=S> ss[N]; // site for each observation
    matrix[N,K] X;              // model matrix with covariates
    vector[N] y;                // productivity (response)
    vector[N] y_sd;             // known sd of observations
    int new_ts[N];              // vector of 0/1 indicating new site years
}

parameters {
    real<lower=0,upper=1> phi;  // ar1 coefficient
    vector[K+1] gamma;          // population level coefficients
    vector[S] beta;             // site level intercepts
    real<lower=0> tau;          // variation in site intercepts
    real<lower=0> sigma;        // standard deviation of process error
    vector[N] mu;               // underlying mean of process
}

transformed parameters {

}

model {
    //priors
    gamma ~ normal(0,5);
    phi ~ beta(1,1);
    tau ~ cauchy(0, 2.5);
    sigma ~ normal(0,1);

    for(s in 1:S){
        beta[s] ~ normal(gamma[1], tau);
    }

    for(n in 1:N){
        if(new_ts[n] == 1){
            // restart the AR process on each new time series
            mu[n] ~ normal(y[n], sigma);
        }
        else{
            mu[n] ~ normal(beta[ss[n]] + X[n,] * gamma[2:K+1] + phi * mu[n-1], sigma);
        }
    }

    // Likelihood
    y ~ normal(mu, y_sd);

}

generated quantities {

    vector[N] y_tilde;
    vector[N] log_lik;

    for(n in 1:N) {
        y_tilde[n] = normal_rng(mu[n], y_sd[n]);
        log_lik[n] = normal_lpdf(y[n] | mu[n], y_sd[n]);
    }

}