Two-parameter Bayesian credibility Interval

Question

How can I determine the confidence/credibility intervals for the posterior estimates of a multi-parameter model?

I can get the confidence interval for each parameter separately. (Currently using bayestestR, but I don't mind using something else)

library(dplyr)
library(ggplot2)
library(bayestestR)

N <- 10000
#Posterior samples (random example)
p1 <- rnorm(N)
p2 <- p1 + rnorm(N)

df_post <- tibble(p1,p2)

describe_posterior(
  df_post,
  centrality = "median",
  test = c("p_direction", "p_significance")
)

##Yields:
# Summary of Posterior Distribution
# 
# Parameter |    Median |        95% CI |     pd |   ps
# -----------------------------------------------------
#   p1      |      0.02 | [-1.85, 2.08] | 50.64% | 0.46
#   p2      | -2.24e-03 | [-2.82, 2.78] | 50.04% | 0.47

and I can generate a plot with the points and 2D contours, which gives me a visual indication of the posterior distribution (though I have no idea what % each contour represents):

ggplot(df_post, aes(x=p1, y=p2)) +
  geom_density_2d(size=1) +
  geom_point(size=0.1)

My question is, How can I actually compute (and/or plot) the two dimensional X% credibility interval?

Answer 1

Here's one base-R-plotting solution, which plots a 95% highest posterior density region based on a 2-D kernel density estimate:

library(emdbook)
library(coda)
HPDregionplot(as.mcmc(df_post))
with(df_post, points(p1, p2, col = adjustcolor("black", alpha.f = 0.2)))

More smoothly/within ggplot:

library(ggplot2); theme_set(theme_bw())
## see function definition below
L <- with(df_post, get_hpd2d_levels(p1, p2))
gg0 <- ggplot(df_post, aes(p1, p2)) + geom_point(alpha=0.1) +
  geom_density_2d(breaks = L, colour="red")
print(gg0)

The highest posterior density region is the classical Bayesian approach. If you wanted to head down a rabbit hole, you could look at some of the less-parametric approaches to computing central sets (bag plots, functional depths, etc.). This would be analogous to the difference between a highest posterior density region and a quantile-based credible interval.

---

##' Get 2D highest posterior density levels corresponding to probability regions
##' @param x x-coordinate of samples
##' @param y y-coordinate
##' @param probs vector of probability levels
##' @param ... arguments for MASS::kde2d
##' @examples
##' dd <- data.frame(x=rnorm(1000), y=rnorm(1000))
##' get_hpd2d_levels(dd$x,dd$y)
##' gg2 <- ggplot(dd, aes(x,y)) + geom_density_2d(breaks=get_hpd2d_levels(dd$x,dd$y), colour="red")
##' print(gg2)
get_hpd2d_levels <- function(x, y, prob=c(0.9,0.95), ...) {
  post1 <- MASS::kde2d(x, y)
  dx <- diff(post1$x[1:2])
  dy <- diff(post1$y[1:2])
  sz <- sort(post1$z)
  c1 <- cumsum(sz) * dx * dy
  ## remove duplicates
  ## dups <- duplicated(sz)
  ## sz <- sz[!dups]
  ## c1 <- c1[!dups]
  levels <- sapply(prob, function(x) {
    approx(c1, sz, xout = 1 - x, ties = mean)$y
  })
  return(levels)
}

Answer 2

@Ben Bolker's answer works well, but here's another suggestion (based on @dash2's comment) that also works:

car::dataEllipse(p1, p2, levels=0.95)