Unlike Bayesian posterior distributions, confidence/consonance functions do not have any distributional properties and also lack the interpretation that should be given to Bayesian posterior intervals. For example, a Bayesian 95% posterior interval has the proper interpretation of having a 95% probability of containing the true value.

This does not apply to 95% frequentist intervals, where the 95% refers to the long run coverage of these intervals containing the true parameter if the study were repeated over and over. Thus, either a 95% frequentist interval contains the true parameter or it does not. In the code below, we simulate some data where the true population parameter is 20 and we know this because we’re the deities of this world. A properly behaving statistical procedure with a set alpha of 0.05 will yield at least 95% intervals in the long run that will include this population parameter of 20. Those that do not are marked in red.

sim <- function() {
  fake <- data.frame((x <- rnorm(100, 100, 20)), (y <- rnorm(100, 80, 20)))
  intervals <- t.test(x = x, y = y, data = fake, conf.level = .95)$conf.int[]
}

set.seed(1031)

z <- replicate(100, sim(), simplify = FALSE)

df <- data.frame(do.call(rbind, z))
df$studynumber <- (1:length(z))
intrvl.limit <- c("lower.limit", "upper.limit", "studynumber")
colnames(df) <- intrvl.limit
df$point <- ((df$lower.limit + df$upper.limit) / 2)
df$covered <- (df$lower.limit <= 20 & 20 <= df$upper.limit)
df$coverageprob <- ((as.numeric(table(df$covered)[2]) / nrow(df) * 100))

library(ggplot2)


ggplot(data = df, aes(x = studynumber, y = point, ymin = lower.limit, ymax = upper.limit)) +
  geom_pointrange(mapping = aes(color = covered), size = .40) +
  geom_hline(yintercept = 20, lty = 1, color = "red", alpha = 0.5) +
  coord_flip() +
  labs(
    title = "Simulated 95% Intervals",
    x = "Study Number",
    y = "Estimate",
    subtitle = "Population Parameter is 20"
  ) +
  theme_bw() + # use a white background
  theme(legend.position = "none") +
  annotate(
    geom = "text", x = 102, y = 30,
    label = "Coverage (%) =", size = 2.5, color = "black"
  ) +
  annotate(
    geom = "text", x = 102, y = 35,
    label = df$coverageprob, size = 2.5, color = "black"
  )

Although the code above demonstrates this, one of the best visualization tools to understand this long-run behavior is the D3.js visualization created by Kristoffer Magnusson, which can be viewed here.

However, despite these differences in interpretation, Bayesian and frequentist intervals often end up converging, especially when there are large amounts of data. They also end up converging when a Bayesian posterior distribution is computed with a flat or weakly informative prior. However, there are several problems with using flat priors, such as giving equal weight to all values in the interval including implausible ones. These sorts of priors should generally be avoided.

Here, I demonstrate with a simple example how Bayesian posterior distributions and frequentist confidence functions end up converging in some scenarios. For these first few examples, I’ll be using the brms package.1

library(concurve)
library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.12.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar
library(ggplot2)
library(cowplot)
#> 
#> ********************************************************
#> Note: As of version 1.0.0, cowplot does not change the
#>   default ggplot2 theme anymore. To recover the previous
#>   behavior, execute:
#>   theme_set(theme_cowplot())
#> ********************************************************
library(bayesplot)
#> This is bayesplot version 1.7.1
#> - Online documentation and vignettes at mc-stan.org/bayesplot
#> - bayesplot theme set to bayesplot::theme_default()
#>    * Does _not_ affect other ggplot2 plots
#>    * See ?bayesplot_theme_set for details on theme setting


GroupA <- rnorm(50)
GroupB <- rnorm(50)
RandomData <- data.frame(GroupA, GroupB)
model_freq <- lm(GroupA ~ GroupB, data = RandomData)
# Using default prior
model_bayes <- brm(GroupA ~ GroupB,
  data = RandomData,
  iter = 2000, warmup = 1000, chains = 2, family = gaussian()
)
randomframe <- curve_gen(model_freq, "GroupB")

(function1 <- ggcurve(type = "c", randomframe[[1]], levels = c(0.99), nullvalue = TRUE))

color_scheme_set("teal")

(function2 <- mcmc_areas(model_bayes, pars = "b_GroupB", point_est = "none", prob = 1, prob_outer = 1, area_method = "scaled height") +
  ggtitle("Posterior Distribution") +
  theme_minimal() +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  annotate("segment",
    x = 0, xend = 0, y = 0.95, yend = 2,
    color = "#990000", alpha = 0.4, size = .75, linetype = 3
  ))
#> Warning: `expand_scale()` is deprecated; use `expansion()` instead.
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.

plot_grid(function1, function2)

# Informative Priors

GroupA <- rnorm(500, mean = 2)
GroupB <- rnorm(500, mean = 1)
RandomData <- data.frame(GroupA, GroupB)
model_freq <- lm(GroupA ~ GroupB, data = RandomData)
# Using default prior
model_bayes <- brm(GroupA ~ GroupB,
  data = RandomData, prior = prior("normal(0, 1)", class = "b"),
  iter = 2000, warmup = 1000, chains = 2, family = gaussian()
)
randomframe <- curve_gen(model_freq, "GroupB")

(function1 <- ggcurve(type = "c", randomframe[[1]], levels = c(0.99), nullvalue = TRUE))

color_scheme_set("teal")

summary(model_bayes)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: GroupA ~ GroupB 
#>    Data: RandomData (Number of observations: 500) 
#> Samples: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 2000
#> 
#> Population-Level Effects: 
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept     1.94      0.07     1.81     2.06 1.00     2214     1537
#> GroupB        0.08      0.05    -0.02     0.17 1.00     2042     1506
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     1.02      0.03     0.96     1.08 1.00     2129     1544
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

(function2 <- mcmc_areas(model_bayes, pars = "b_GroupB", point_est = "none", prob = 1, prob_outer = 1, area_method = "scaled height") +
  ggtitle("Posterior Distribution") +
  theme_minimal() +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  annotate("segment",
    x = 0, xend = 0, y = 0.95, yend = 2,
    color = "#990000", alpha = 0.4, size = .75, linetype = 3
  ))
#> Warning: `expand_scale()` is deprecated; use `expansion()` instead.
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.

plot_grid(function1, function2)

Here I use a simple experimental design taken from a brms vignette

income_options <- c("below_20", "20_to_40", "40_to_100", "greater_100")
income <- factor(sample(income_options, 100, TRUE),
  levels = income_options, ordered = TRUE
)
mean_ls <- c(30, 60, 70, 75)
ls <- mean_ls[income] + rnorm(100, sd = 7)
dat <- data.frame(income, ls)
dat$income_num <- as.numeric(dat$income)
fit1 <- brm(ls ~ income_num, data = dat)
fit2 <- lm(ls ~ income_num, data = dat)

randomframe <- curve_gen(fit2, "income_num")

(function3 <- ggcurve(type = "c", randomframe[[1]], levels = c(0.99), nullvalue = FALSE))


summary(fit1)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: ls ~ income_num 
#>    Data: dat (Number of observations: 100) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Population-Level Effects: 
#>            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept     21.23      2.44    16.39    25.88 1.00     3944     2768
#> income_num    15.20      0.93    13.38    16.99 1.00     4043     3072
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma    10.22      0.72     8.93    11.75 1.00     4012     2903
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

(function4 <- mcmc_areas(fit1, pars = "b_income_num", point_est = "none", prob = 1, prob_outer = 1, area_method = "scaled height") +
  ggtitle("Posterior Distribution") +
  theme_minimal() +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  annotate("segment",
    x = 15, xend = 15, y = 0.95, yend = 2,
    color = "#990000", alpha = 0.4, size = .75, linetype = 3
  ))
#> Warning: `expand_scale()` is deprecated; use `expansion()` instead.
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.

plot_grid(function3, function4)

Here we use another example from an rstanarm vignette.2

library(rstanarm)
data(kidiq)

post1 <- stan_glm(kid_score ~ mom_hs,
  data = kidiq,
  family = gaussian(link = "identity"),
  seed = 12345
)
post2 <- lm(kid_score ~ mom_hs, data = kidiq)

df3 <- curve_gen(post2, "mom_hs")

function99 <- ggcurve(df3[[1]])

summary(post1)
#> 
#> Model Info:
#>  function:     stan_glm
#>  family:       gaussian [identity]
#>  formula:      kid_score ~ mom_hs
#>  algorithm:    sampling
#>  sample:       4000 (posterior sample size)
#>  priors:       see help('prior_summary')
#>  observations: 434
#>  predictors:   2
#> 
#> Estimates:
#>               mean   sd   10%   50%   90%
#> (Intercept) 77.6    2.1 74.9  77.6  80.2 
#> mom_hs      11.8    2.4  8.7  11.7  14.9 
#> sigma       19.9    0.7 19.0  19.9  20.8 
#> 
#> Fit Diagnostics:
#>            mean   sd   10%   50%   90%
#> mean_PPD 86.8    1.3 85.1  86.8  88.5 
#> 
#> The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
#> 
#> MCMC diagnostics
#>               mcse Rhat n_eff
#> (Intercept)   0.0  1.0  3755 
#> mom_hs        0.0  1.0  3776 
#> sigma         0.0  1.0  3456 
#> mean_PPD      0.0  1.0  3341 
#> log-posterior 0.0  1.0  2069 
#> 
#> For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

(function100 <- mcmc_areas(post1, pars = "mom_hs", point_est = "none", prob = 1, prob_outer = 1, area_method = "scaled height") +
  ggtitle("Posterior Distribution") +
  theme_minimal() +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  annotate("segment",
    x = 15, xend = 15, y = 0.95, yend = 2,
    color = "#990000", alpha = 0.4, size = .75, linetype = 3
  ))
#> Warning: `expand_scale()` is deprecated; use `expansion()` instead.
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.

cowplot::plot_grid(function99, function100)


References


1. Burkner P-C. Advanced bayesian multilevel modeling with the r package brms. The R Journal. 2018;10(1):395-411. doi:10/gfxzpn

2. Goodrich B, Gabry J, Ali I, Brilleman S. Rstanarm: Bayesian Applied Regression Modeling via Stan.; 2020. https://mc-stan.org/rstanarm.