Crossvalidation in hierarchical bayesian models (HBMs)

Question

I am trying to find a way to cross-validate Hierarchical Bayesian Models used for predicting and modelling abundance in Species Distribution Models. For this purpose, I have tried posterior predictive checks proposed in Kéry & Schaub 2011, chapter 12.3.1 with the following modified $\chi^2$ discrepancy measure:

$$\sum_i{{(Y_i - EY_i)^2}\over{var Y_i + 0.5}}$$

(the 0.5 is added to prevent division by zero - this is slighlty modified Gelman 2004 (2nd edition, pg. 175) $\chi^2$ discrepancy).

As I want to test the predictive power of the model - at sites for which we don't have data - I split the data into fitting set and validation set and I did the predictive posterior checks in two variants:

• on the fitted set
• on the validation (i.e. cross-validation) set

Very simplified, the JAGS model including posterior predictive checks looks like this:

for (i in 1:sites) {
    log(lambda[i]) <- lambda_int + inprod(b_lambda[], envCov[i,]) + site_eff[i] 
    site_eff[i] ~ dnorm(0, 1/(site_eff_sd*site_eff_sd))

    N[i] ~ dpois(lambda[i])

    for (j in 1:visits) {
        p[i,j] <- ... some predictor ...

        # posterior predictive checks, a la Kery, BPA chapter 12.3.1
        # done for all the sites (fitted/crossvalidatory)
        Y.new[i, j] ~ dbinom(p[i,j], N[i])
        E[i,j] <-     pow(    Y[i,j] - p[i,j] * N[i], 2) / (p[i,j] * N[i] + 0.5)
        E.new[i,j] <- pow(Y.new[i,j] - p[i,j] * N[i], 2) / (p[i,j] * N[i] + 0.5)
    }
}

# sites determined for fitting the model (other sites are for validation)
for (i in 1:fit_sites) {
    for (j in 1:visits) {
        Y[i, j] ~ dbinom(p[i,j], N[i])
    }
}

# classic PPC: on fitted data
fit <- sum(E[1:fit_sites,])
fit.new <- sum(E.new[1:fit_sites,])

# crossvalidatory PPC: on independent data
fit.cv <- sum(E[(fit_sites+1):sites,])
fit.cv.new <- sum(E.new[(fit_sites+1):sites,])

where Y[i,j] is the observed count of the species at site i on occasion j, N[i] is the actual latent abundance of the species at site i. First two thirds of the sites are used to fit the data, remaining third of the sites is reserved for crossvalidation only.

Questions:

1. Is this approach OK? Kéry doesn't mention crossvalidation, and Gelman 2004 has a note in chapter 8.2 which I do not understand (especially the italic part):

If a model is to be evaluated by its expected predictive fit to new data, a reasonable approach would seem to be cross-validation - fitting the model to part of a dataset and evaluating it on the rest - but it is not clear how to interpret this in terms of posterior inference. Conditional on a model being correct, we can directly assess it's predictive accuracy by comparing its fit to simulated replicated data sets $Y.new$ - but DIC and cross-validation attempt to do better and capture predictive error including model misfit.

2. Is the modification with + 0.5 OK for this purpose? I mean, couldn't it make any mess?

3. How to create easily interpretable measure to assess quality of the model? I mean something that would tell me 0% to 100% - from absolutely poor to absolutely perfect model. Kery is proposing two measures:

• bayesian p-value mean(simulations$fit.new > simulations$fit),
• lack-of-fit ratio mean(simulations$fit / simulations$fit.new)

But how to interpret those? Is e.g. a lack-of-fit ratio of 2 good or bad? I would like some measure which has clear scale like $r$ or $AUC$ where you can e.g. say that 1 is a perfect model, > 0.75 is good, < 0.5 is poor...

Here are few examples of results and plots, where I plot fit or fit.cv on the X axis and fit.new and fit.cv.new on the Y axis.

Example 1 - comparison of two models:

linear model PPC - fitted data:
          p-value lack-of-fit ratio 
        0.9276667         0.9442094 

linear model PPC - independent data:
          p-value lack-of-fit ratio 
         0.000000          1.820381 

superpredictor PPC - fitted data:
          p-value lack-of-fit ratio 
        0.9810000         0.9213055 

superpredictor PPC - independent data:
          p-value lack-of-fit ratio 
         0.000000          1.742552 

Example 2 - comparison of two models:

linear model PPC - fitted data:
          p-value lack-of-fit ratio 
        0.7705333         0.9245741 

linear model PPC - independent data:
          p-value lack-of-fit ratio 
         0.002200          2.008853 

superpredictor PPC -  fitted data:
          p-value lack-of-fit ratio 
        0.7315333         0.9389283 

superpredictor PPC -  independent data:
          p-value lack-of-fit ratio 
         0.000000          2.985811 

Example 3 - comparison of two models:

linear model PPC -  fitted data:
          p-value lack-of-fit ratio 
        0.6679667         0.9722895 

linear model PPC -  independent data:
          p-value lack-of-fit ratio 
      0.005516667       1.341490485 

superpredictor PPC - fitted data:
          p-value lack-of-fit ratio 
        0.5848000         0.9869094 

superpredictor PPC - independent data:
          p-value lack-of-fit ratio 
         0.000000          2.012136 

In models 1,2 I can tell that the superpredictor is better than the linear model, but how to get any idea of how good/poor both of the models are on a scale (absolutely poor) - (absolutely perfect)?

Would e.g. something like comparison with null model do the work? I have tried it:

Answer 1

I work on similar models using the same book so I encountered the same issues - thanks for interesting questions.

1. Can posterior predictive checks be performed using both fitted and out-of-sample data?

Yes - at least some authors propose to do both of these in the statistical literature (e.g., http://www.math.ntnu.no/~daniesi/INLA/Papers/festschrift.pdf). But they refer more specifically to CPO (conditional predictive ordinate) and PIT (probability integral transform)

However, Bayesian p-values and out-of-sample-based predictive probabilities seem to measure slightly different things. The latter ones measure your ability to predict new data, which is a blend of goodness of fit and model complexity (keep in mind that in some contexts, IT measures favouring parsimony such as AIC are asymptotically equivalent to leave-one-out CV).

The Bayesian P-value $\text{Pr}(T(Y^{rep})>T(y)|y,\theta)$ (using the fitted data only) would be instead the equivalent of your $R^2$ , a simple measure of fit, except that rather than using only a point prediction from the fitted model, it uses the whole posterior distribution of $y$ values.

Digging a little through the citation pipeline, the above paper refers to

Dawid, A. P. (1984). Statistical theory: The prequential approach, Journal of the Royal Statistical Society. Series A (General) 147: 278-292.

Pettit, L. I. (1990). The conditional predictive ordinate for the normal distribution, Journal of the Royal Statistical Society: Series B 52: 175-184.

Marshall, E. C. and Spiegelhalter, D. J. (2003). Approximate cross-validatory predictive checks in disease mapping models. Statistics in Medicine 22, 1649–1660.

Many papers including a complicated model seem to have their own specific predictive measure based on out of sample data, thus perhaps that's appropriate for you too. Which brings us to what you're trying to do by using cross-validation: selecting a model balancing goodness of fit with complexity. On that point the following review might be up-to-date

Hooten, M. B., & Hobbs, N. T. (2015). A guide to Bayesian model selection for ecologists. Ecological Monographs, 85(1), 3-28.

I don't know if $\chi^2$ is the best distance/discrepancy measure, not entirely sure either if that makes sense to plot the independently obtained discrepancy values with those from the fitted data (but perhaps it does).

2. 0.5 Correction.

$\text{Var}(Y_i)$ is likely to be above 0.5 if $N>>1$, which in turn depends on $\lambda$. You can always use 0.05, 0.1, 0.2 instead to check.

3. Measure of fit.

You seem to be looking for a $R^2$. But even with such measure, good or bad is context dependent, i.e. $75 \%$ will appear barely believable in some fields and ridiculously low in others.

Here's what I found on interpreting Bayesian p-values http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.310.145&rep=rep1&type=pdf (statisticians seem to debate their distribution btw).

One practical solution would be simulate data according to a model which you know to be close to your situation (similar sample sizes, number of sites, approximate abundances/occupancy etc.), then fit the simulated model to simulated data, and use the resulting Bayesian P-values as measures of "good" (as good as it is going to get) when interpreting those from your real data.

Lack-of-fit ratio: Looks like a comparison of mean squared errors. Isn't it mean(simulations$fit) / mean(simulations$fit.new), rather?

The same simulation approach could work to get sensible reference values. Also, if we were comparing to sets of replicated data, the discrepancy measures should follow two $\chi^2$ distributions which means the ratio follows an F - and gives access to a different kind of P-value. [just throwing in ideas here]

The comparison to the null model makes more sense to me when you compare the predictive abilities of the models using the independent data.

EDIT 10/06

The discrepancy measure of Kéry and Schaub is in fact $\sum_{i=1}^{n}{\frac{(Y_i-E(Y_i))^2}{E(Y_i)+0.5}}$, no variance term below.

EDIT2 10/06

For similar hierarchical models, I also get Bayesian P-values always above 0.5 (no matter how good the fit looks graphically) and very often close to 1, with Lack-of-fit ratios just below 1 (e.g. 0.8 or 0.9)