Warnings during WAIC computation: how to proceed?

Question

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)
}

Answer 1

I'm happy that you care about these diagnostic messages. I'm responsible for that specific warning message. WAIC doesn't have a good way to diagnose the reliability and that 0.4 threshold is empirically chosen. Pareto k diagnostic in PSIS-LOO is much better. You did not mention the specific k values, but if they are larger than 1, then theory and experiments show that the error can be arbitrary large (see Aki Vehtari, Andrew Gelman and Jonah Gabry (2017). Pareto smoothed importance sampling. https://arxiv.org/abs/1507.02646). WAIC and PSIS-LOO are connected so that if PSIS-LOO fails then WAIC fails even more (Aki Vehtari, Andrew Gelman and Jonah Gabry (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. In Statistics and Computing, 27(5):1413–1432, https://arxiv.org/abs/1507.04544). Thus, if even one p_waic can have arbitrary large error, then the total WAIC can have arbitrary large error and should not be trusted for final results.

Since you have a large difference in WAIC, it's likely that the difference is similar with more accurate/robust approach (e.g. k-fold-CV), but you cannot be sure without computing that more accurate/robust result. In your case it is likely that WAIC and PSIS-LOO fail because you have a flexible model and some of the observation are highly influential, that is, the full data posterior and leave-one-out-posterior are so different that using full posterior as the proposal distribution in importance sampling LOO fails. WAIC uses instead Taylor series approximation, which also fails in case of highly influential observations (which explains why p_waic is used as the ad-hoc diagnostic).