R code for creating plot showing the effect of different priors

priors.R

library(tidyverse)
library(viridis)

set.seed(92414)


### Draw samples from prior distributions --------------------------------------
# Non-Information: N(-2, 10)
# Weakly Informative: N(-2, 2)
# Informative: N(-2, 1)

# Prior distributions (for plotting)
plot_data <- data_frame(x = seq(-10, 10, by = 0.01)) %>%
  mutate(
    `Non-Informative` = dnorm(x, mean = -2, sd = 10),
    `Weakly Informative` = dnorm(x, mean = -2, sd = 2),
    Informative = dnorm(x, mean = -2, sd = 1)
  ) %>%
  gather(key = "prior_type", value = "prior", -x) %>%
  mutate(prior_type = factor(prior_type, levels = c("Non-Informative",
                                                    "Weakly Informative",
                                                    "Informative")))

# Prior samples
ni_prior <- data_frame(prior_type = "Non-Informative",
                       prior = rnorm(n = 10000, mean = -2, sd = 10))
wi_prior <- data_frame(prior_type = "Weakly Informative",
                       prior = rnorm(n = 10000, mean = -2, sd = 2))
in_prior <- data_frame(prior_type = "Informative",
                       prior = rnorm(n = 10000, mean = -2, sd = 1))

prior_samples = bind_rows(ni_prior, wi_prior, in_prior) %>%
  mutate(prior_type = factor(prior_type, levels = c("Non-Informative",
                                                    "Weakly Informative",
                                                    "Informative")))


### Calculate likelihood -------------------------------------------------------
# Likelihood for sample data
sample_data <- rnorm(n = 1000, mean = 5, sd = 2.5)
plot_data <- plot_data %>%
  mutate(likelihood = dnorm(x = mean(sample_data), mean = x, sd = 2.5))

# Likelihood for prior draws
prior_samples <- mutate(prior_samples,
                        prior_likelihood = dnorm(x = mean(sample_data),
                                                 mean = prior, sd = 2.5))


### Calculate posteriors -------------------------------------------------------
posterior_samples <- prior_samples %>%
  group_by(prior_type) %>%
  sample_n(size = 10000, replace = TRUE, weight = prior_likelihood)


### Create plot ----------------------------------------------------------------
ggplot(plot_data, aes(x = x)) +
  facet_wrap(~ prior_type, nrow = 1) +
  geom_line(aes(y = prior, color = "Prior"), linetype = "dashed") +
  geom_line(aes(y = likelihood, color = "Likelihood"), linetype = "dashed") +
  geom_density(data = posterior_samples,
               aes(x = prior, color = "Posterior", fill = "Posterior"),
               alpha = 0.5, size = 1, show.legend = FALSE) +
  scale_color_manual(values = c("Prior" = viridis(3)[2],
                                "Likelihood" = viridis(3)[1],
                                "Posterior" = viridis(3)[3]),
                     breaks = c("Prior", "Likelihood", "Posterior"),
                     guide = guide_legend(override.aes = list(
                       linetype = "solid", alpha = 1, size = 5
                     ))) +
  scale_fill_manual(values = c("Posterior" = viridis(3)[3])) +
  coord_cartesian(xlim = c(-10, 10)) +
  labs(x = NULL, y = "Density", color = NULL) +
  theme_bw() +
  theme(legend.position = "bottom")