bayes.Rmd

---
title: "Comparison to Bayesian Posterior Distributions"
description: "This vignette discusses how consonance distributions and Bayesian posterior distributions differ in interpretation but how they can also often converge. This is shown with several examples of how to calculate both distributions in a range of scenarios."
output: 
  rmarkdown::html_vignette:
    toc: true
opengraph:
  image: 
    src: "https://upload.wikimedia.org/wikipedia/commons/d/d4/Thomas_Bayes.gif"
  twitter:
    card: summary
    creator: "@dailyzad"    
bibliography: references.bib
link-citations: yes
csl: american-medical-association.csl
vignette: >
  %\VignetteIndexEntry{Comparison to Bayesian Posterior Distributions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  message = TRUE,
  warning = TRUE,
  collapse = TRUE,
  comment = "#>"
)
```


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. 

```{r}
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](https://rpsychologist.com/d3/CI/).

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. However, for the sake of this demonstration, we will be using flat priors.

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 `rstanarm` package.[@goodrichRstanarmBayesianApplied2020] 

```{r echo=TRUE}
library(concurve)
library(rstan)
library(rstanarm)
library(ggplot2)
library(cowplot)
library(bayesplot)
library(scales)
```

We will simulate some data (two variables) from a normal distribution with a location parameter of 0 and scale parameter of 1 (something very simple) and then regress the second variables (GroupB) on the first using the base lm function. We will take the regression coefficient and construct a consonance function for it. 
```{r echo=TRUE}
GroupA <- rnorm(50)
GroupB <- rnorm(50)
RandomData <- data.frame(GroupA, GroupB)
model_freq <- lm(GroupA ~ GroupB, data = RandomData)
```

Now we will do the same using Bayesian methods, but instead of specifying a prior, we will use a flat prior to show the convergence of the posterior with the consonance function.

```{r, results = 'hide', message = FALSE}
rstan_options(auto_write = TRUE)

# Using flat prior
model_bayes <- stan_lm(GroupA ~ GroupB,
  data = RandomData, prior = NULL,
  iter = 5000, warmup = 1000, chains = 4
)
```

Now that we've fit the models, we can graph the functions.

```{r echo=TRUE}
randomframe <- curve_gen(model_freq, "GroupB", steps = 10000)

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

color_scheme_set("teal")

function2 <- mcmc_dens(model_bayes, pars = "GroupB") +
  ggtitle("Posterior Distribution") +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") +
  scale_y_continuous(breaks = c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0))


(breaks1 <- c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0))

(adjustment <- function(x) {
  x / 3
})

(labels <- adjustment(breaks1))

breaks <- labels
labels1 <- labels

(function3 <- mcmc_dens(model_bayes, pars = "GroupB") +
  ggtitle("Posterior Distribution") +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") +
  scale_x_continuous(expand = c(0, 0), breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(expand = c(0, 0), breaks = waiver(), labels = waiver(), n.breaks = 10, limits = c(0, 3.25)) +
  yaxis_text(on = TRUE) +
  yaxis_ticks(on = TRUE) +
  annotate("segment",
    x = 0, xend = 0, y = 0, yend = 3,
    color = "#990000", alpha = 0.4, size = .75, linetype = 3
  ))
```

I made some adjustments above to the bayesplot code so that we could more easily compare the consonance distribution to the posterior distribution. We will be using plot_grid() from cowplot to achieve this. 

```{r echo=TRUE, fig.height=9, fig.width=7}
plot_grid(function1, function3, ncol = 1, align = "v")
```

As you can see, the results end up being very similar. You can likely get very similar results with weakly informative priors normal(0, 100) or with much larger datasets, where the likelihood will end up swamping the prior, though this isn't always the case.

Here's another example, but this time the variables we simulate have different location parameters. 

```{r echo=TRUE}
GroupA <- rnorm(500, mean = 2)
GroupB <- rnorm(500, mean = 1)
RandomData <- data.frame(GroupA, GroupB)
model_freq <- lm(GroupA ~ GroupB, data = RandomData)
```

```{r, results = 'hide', message = FALSE}
# Using flat prior
model_bayes <- stan_lm(GroupA ~ GroupB,
  data = RandomData, prior = NULL,
  iter = 5000, warmup = 1000, chains = 4
)
```


```{r echo=TRUE}
randomframe <- curve_gen(model_freq, "GroupB", steps = 10000)

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

color_scheme_set("teal")

function2 <- mcmc_dens(model_bayes, pars = "GroupB") +
  ggtitle("Posterior Distribution") +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") +
  scale_y_continuous(breaks = c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0))


(breaks1 <- c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0))

(adjustment <- function(x) {
  x / 3
})

(labels <- adjustment(breaks1))

breaks <- labels
labels1 <- labels

(function3 <- mcmc_dens(model_bayes, pars = "GroupB") +
  ggtitle("Posterior Distribution") +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") +
  scale_x_continuous(expand = c(0, 0), breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(expand = c(0, 0), breaks = waiver(), labels = waiver(), n.breaks = 10, limits = c(0, 9)) +
  yaxis_text(on = TRUE) +
  yaxis_ticks(on = TRUE) +
  annotate("segment",
    x = 0, xend = 0, y = 0, yend = 9,
    color = "#990000", alpha = 0.4, size = .75, linetype = 3
  ))
```
```{r echo=TRUE, fig.height=9, fig.width=7}
plot_grid(function1, function3, ncol = 1, align = "v")
```


Here's another dataset, however, here we're not generating random numbers. 


```{r, results = 'hide', message = FALSE}
data(kidiq)

# flat prior

post1 <- stan_lm(kid_score ~ mom_hs,
  data = kidiq, prior = NULL,
  seed = 12345
)
```

```{r}
post2 <- lm(kid_score ~ mom_hs, data = kidiq)

df3 <- curve_gen(post2, "mom_hs")

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

summary(post1)

color_scheme_set("teal")

(function101 <- mcmc_areas(post1, pars = "mom_hs", point_est = "none", prob = 1, prob_outer = 1, area_method = "equal height") +
  ggtitle("Posterior Distribution") +
  labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") +
  yaxis_text(on = TRUE) +
  yaxis_ticks(on = TRUE))
```

```{r echo=TRUE, fig.height=9, fig.width=7}
cowplot::plot_grid(function99, function101, ncol = 1, align = "v")
```

# Cite R Packages

Please remember to cite the packages that you use.

```{r}
citation("concurve")
citation("ggplot2")
citation("rstan")
citation("rstanarm")
citation("cowplot")
citation("scales")
citation("bayesplot")
```

# References

* * *