# Warnings during WAIC computation: how to proceed?

I am computing the WAIC (widely applicable or Watanabe-Akaike information criterion) using the waic() function from the 'loo' package in R. When I do so, I see a warning message as follows:

Warning messages: 1: 118 (5.1%) p_waic estimates greater than 0.4. We recommend trying loo() instead.

What should I do about this? I have noticed that several online examples/tutorials ignore this warning (e.g. here--this is a reproducible example). Is 5.1% (in my case) or 3.3% (in the example's case) a dangerously high number? How is this likely to impact my inference when comparing two models? What if the WAIC values are REALLY different in my comparison (like more than 40 WAIC units different); can I trust the inference that one model is to be preferred?

Using loo() instead (for leave-one-out cross-validation), I get a different warning:

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

My model uses Bernoulli error, so the error term cannot be misspecified. In my particular class of model, the linear predictor is also unlikely to be badly misspecified.

WARNING, GORY DETAILS FOLLOW

My model is a latent variable model (see here) to jointly model the distributions of fifteen bird species across ~100 sites. The model includes no covariates, and so can be thought of as a model-based ordination. The occurrence probability of species i at site j is modeled on the probit scale as a species-specific intercept plus the impact of several (1-3) latent factors that behave like covariates in a GLM. Each latent factor takes a particular value (estimated from data) at each site, and each species gets a species-specific factor loading that governs its abundance across sites. The matrix of factor loadings times its transpose is an approximation to the full (15x15) variance-covariance matrix among species. I want to fit the model with two latent factors (because this produces ordinations that are easy to visualize in two dimensions). However, changing the number of latent factors to one or three does not improve the percentage of p_waic estimates greater than 0.4; this leads me to believe that I can't ameliorate the issue by choosing a better specification for the linear predictor.

My reason for performing model selection is not to dredge through a large number of possible covariates, but rather to perform inference on exactly two alternative models that represent meaningfully different biological hypotheses/scenarios. I want to know if the evidence strongly favors one hypothesis over the other.

For those interested, the full model specification in JAGS, including the computation of the log-likelihood, is below:

##### Model specification #####
LV_C <- function() {
## Data Level ##
for(i in 1:n) { # sites
for(j in 1:p){ # species
eta[i,j] <- inprod(lv.coefs[j,2:(num.lv+1)],lvs[i,]) # LV part of the linear predictor
Z[i,j] ~ dnorm(lv.coefs[j,1] + eta[i,j], 1) # This line and the next one implement a probit link
y[i,j] ~ dbern(step(Z[i,j]))
loglikelihood[(i-1)*(p)+j] <- y[i,j]*log(phi(lv.coefs[j,1] + eta[i,j])) + (1 - y[i,j])*log(1 - phi(lv.coefs[j,1] + eta[i,j]))
}
}

## Latent variables ##
for(i in 1:n) { for(k in 1:num.lv) { lvs[i,k] ~ dnorm(0,1) } } # Says what Latent Variable values are

## Process level and priors ##
for(j in 1:p) { lv.coefs[j,1] ~ dnorm(0,0.01) } ## Separate species intercepts

for(i in 1:(num.lv-1)) { for(j in (i+2):(num.lv+1)) { lv.coefs[i,j] <- 0 } } ## Constraints to 0 on upper diagonal
for(i in 1:num.lv) { lv.coefs[i,i+1] ~ dunif(0,20) } ## Sign constraints on diagonal elements
for(i in 2:num.lv) { for(j in 2:i) { lv.coefs[i,j] ~ dunif(-20,20) } } ## Free lower diagonals
for(i in (num.lv+1):p) { for(j in 2:(num.lv+1)) { lv.coefs[i,j] ~ dunif(-20,20) } } ## All other elements
}

#  Set up for Run

jags.data = list(y=ssm2, n=dim(ssm2)[1], p=dim(ssm2)[2],
num.lv=2)

params = c('lvs', 'lv.coefs', 'loglikelihood')

p <- dim(ssm2)[2]
n <- dim(ssm2)[1]
inits <- function(jjj) {
Tau <- rWishart(1,p+1,diag(p))[,,1]
Sigma <- solve(Tau)
Z <- abs(t(mvtnorm::rmvnorm(n,rep(0,p),Sigma)))
Z <- ifelse(as.matrix(ssm2), Z, -1 * Z)

list(Z=Z)
}