I am computing the WAIC (widely applicable or Wataname-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][1]). 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?

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 this 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][2]) 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)
    }

  [1]: https://jfiksel.github.io/2017-05-24-waic_aft_models_jags/
  [2]: http://www.cell.com/trends/ecology-evolution/abstract/S0169-5347(15)00240-2