In cbrown5/BayeSens: Prior sensitivity analysis for Bayesian MCMC models

This package provides a few helpful functions for measuring the sensitivity of Bayesian posteriors to priors. In particular these tools are designed to be used with posteriors approximated through MCMC algorithms.

The package also aims to provide a helpful introductory guide to some simple Bayes sensitivity analysis. I recommend reading the references for each function that are in the help files for more information.

Currently all functions take numeric vectors as inputs. Each vector represents random draws from a distribution. In the future I plan to provide specific methods for RJags and other common MCMC algorithms.

Hellinger distance - non-parametric

The Hellinger distance is a standard metric for measuring the difference between two distributions.

It ranges 0-1 where values of zero indicate the distributions are identical and values towards 1 indicate a very large difference. The Hellinger distnace will = 1 if distribution A has positive density everywhere that distribution B has zero density.

As an example, say we have random draws from two normal distributions. To generate some example data:

r1 <- rnorm(10000, 1, 1)
r2 <- rnorm(10000, 1, 3)

In practice r1 and r2 could be draws from a prior and posterior respectively.

To calculate the Hellinger distance between these distributions:

library(BayeSens)
hout <- hellinger(r1, r2, nbreaks = 100)
hout

Which provides two estimates. Hellinger distance is approximated in two ways:

(1) by binning the random variates and calculating the Hellinger distance for discrete distributions and

(2) by creating a continuous approximation of the distributions using (density) and then using numerical integration to calculate the Hellinger distance.

Method (2) - continuous integration - should in genernal be more accurate however, it may give poor approximations for multi-modal distributions. Continuous integration may return NaN if the distributions are near identical.

We can also plot the results, to check the distributions are fitted well by the non-parametric density curves:

plot(hout)

It is recommended to visually check distribution fits, particularly if the number of random variates is small. In general these methods will be inaccurate if analysis is performed on too few samples, e.g. <10 000. >100 000 would be ideal.

Hellinger distance - parametric

You can also fit specific distributions to your random draws and estimate the Hellinger distance. Currently normal-normal and beta-beta comparisons are supported.

For instance:

houtp <- hdistpara(r1, r2, densfun = 'normal')

houtp
plot(houtp)

Posterior shrinkage

Another method for prior sensitivity analysis is to calculate posterior shrinkage, that is how much a posterior estimate has shrunk towards its maximum likelihood value. Posterior shrinkage is the estimate of $\alpha$ in:

$\theta_D = (1 - \alpha) \theta_p + \alpha \theta_{MLE}$

The function simply takes three estimates of a parameter, e.g.

postshrink(thetaprior = 0.1, thetapost = 0.2, thetaMLE = 0.22)

Values close to 1 indicate the prior has little influence on the posterior, whereas values close to zero indicate the prior has a large influence on the posterior.

Some care must be taken in selecting parameter estimates to use in the shrink equation and also in estimating the MLE (which is not alwaysstraightfoward). For complex models the MLE may be estimated using data cloning, see below.

Shrink values can occaisionally be >1 or <0. Values <0 occur when the posterior parameter estimate has moved in the opposite direction from the prior than the MLE.
Values >1 occur when the MLE is closer to the prior estimate than the posterior. Values not in 0-1 can occur if (1) your MLE estimate inaccurate, (2) your posterior is multi-model, or (3) your posterior estimate is constrained by other parameters

If (1) then try other methods for obtaining an MLE or increase replication if using data cloning.

If (2) or (3) posterior shrink may be inappropriate for your model, because the posterior shrink cannot be characterised by a simple univariate measure.

See: Berger JO (1985) Statistical Decision Theory and Bayesian Analysis, Second Edition, Springer, New York.

In practice, the MLE can be hard to obtain. One method available for complex models is data cloning, whereby you clone your dataframe K times and then refit your model using the MCMC algorithm. If K is large enough, the posterior means will be esimates of the MLE (which are uninfluenced by the priors).

The function dataclone is just a helper function that replicates the rows of a dataframe K times. Use it like this:

dataclone(dat, K)

Data cloning can be used to obtain maximum likelihood estimates for parameters using Bayesian MCMC algorithms.

It may be useful for instance when estimating posterior shrinkage with postshrink. It is recommended that multiple values of K are run to find a large enough value that gives stable esimates.

For more info see this paper which provides a pretty accessible introductory guide to data cloning and MCMC:

Lele SR, Dennis B, Lutscher F. Data cloning: easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods. Ecology letters. 2007 Jul 1;10(7):551-63.
Click here for an open access versio?n.

cbrown5/BayeSens documentation built on April 26, 2020, 12:40 a.m.