Frequentist Perspective

Classically, the approach to this problem is taught from the frequentist perspective. Uncertainty in the values of the parameters that generated the data is represented by probability via the notion of repeated sampling: under the given probability model with the given parameter values, what is the relative frequency with which these same data would be observed, if the experiment that generated the data were repeated again and again? Values of the parameters that have a higher probability of having generated the observed data are thought to be more likely than other values.

As an example, consider a coin with unknown probability of heads \(\theta\). We toss the coin once, and observe random outcome (data) \(X = 1\) if the toss is heads and \(X = 0\) if not. For simplicity, suppose we know we have chosen one of two possible coins, either having \(\theta = 0.7\) or \(\theta = 0.3\). How do we use the observed data to infer which of these two coins we threw?

For any \(0 < \theta < 1\), the probability distribution of the single coin toss is given by \[ P(X = x) = \theta^{x}(1-\theta)^{1-x} \] This just says that \(P(X = 1) = \theta\) and \(P(X = 0) = 1-\theta\).

Let’s say we throw the coin once, and observe \(X = 1\). If \(\theta = 0.7\) the probabilty of observing this result is \(P(X = 1|\theta = 0.7) = 0.7\). That is if \(\theta = 0.7\), we would expect roughly \(70\%\) of repetitions of this experiment to yield the same results as we observed in our data. If \(\theta = 0.3\) on the other hand, \(P(X = 1|\theta = 0.3) = 0.3\); only \(30\%\) of the repetitions of this experiment would yield the observed data if \(\theta = 0.3\). Because \(\theta = 0.7\) would yield the observed data more frequently than \(\theta = 0.3\), we say that \(\theta = 0.7\) is more likely to have generated the observed data than \(\theta=0.3\), and our inference favours \(\theta =0.7\).

For a more comprehensive background on frequentist/likelihood inference, check out the tutorial on maximum likelihood.

Bayesian Perspective

One criticism of the above approach is that is depends not only on the observed data, but also on infinitely many other possible datasets that are not observed. This is an artifact of the manner in which probability is used to represent uncertainty. In contrast, Bayesian statistics represents uncertainty about the value of a parameter directly using probability distributions.

In particular, a prior distribution is placed on the parameter, representing the probable values of that parameter before data is observed. Having observed the data, the prior is updated via Bayes’ Rule, yielding the posterior distribution of the parameter, given the data. The choice of prior distribution is based either on subject-matter knowledge or mathematical convenience, and is a subjective choice on the part of the analyst.

To see how this works, suppose that we think about \(4/5\) coins in our pocket are the \(\theta = 0.3\) coins, and only \(1/5\) are the \(\theta = 0.7\) coins. Our prior distribution on \(\theta\) is then \[ P(\theta = q) = 0.2^{I(q = 0.7)}0.8^{I(q = 0.3)} \] where \(I(q = 0.7) = 1\) if \(q = 0.7\) and \(0\) otherwise. This, like the probability distribution of the actual result of the coin toss, just encodes our notion that \(P(\theta = 0.3) = 0.8\) and \(P(\theta = 0.7) = 0.2\).

So without knowing the result of the coin toss, we think there is a \(20\%\) chance that \(\theta = 0.7\). We know from above that if we observe heads on the coin toss, we have observed a result that would occur about \(70\%\) of the time if \(\theta = 0.7\). In probability terms, we have a marginal distribution for \(\theta\), and a conditional distribution for \(X|\theta\).

These two ideas are combined by computing the conditional distribution of \(\theta|X\), known as the posterior distribution for \(\theta\) having observed \(X\). This is obtained (explaining the name) via Bayes’ Rule: \[ p(\theta|X) = \frac{p(X|\theta)\times p(\theta)}{p(X)} \] where the marginal distribution of \(X\), or the normalizing constant or marginal likelihood or model evidence (this thing has a lot of names) is given by \[ p(X) = \sum_{\theta}p(X|\theta)\times p(\theta) \] and ensures \(p(\theta|X)\) is a proper probability distribution.

In our example, the prior probability of \(\theta = 0.7\) is only \(20\%\). But we flip the coin and observe \(X = 1\). We can see how this observation updates our belief about the likely values of \(\theta\) by computing the posterior distribution of \(\theta\) given the observed data: \[ \begin{aligned} &p(\theta|X) = \frac{\theta^{x}(1-\theta)^{1-x}\times 0.2^{I(\theta = 0.7)}0.8^{I(\theta = 0.3)}}{\sum_{\theta = 0.3,0.7}\theta^{x}(1-\theta)^{1-x}\times 0.2^{I(\theta = 0.7)}0.8^{I(\theta = 0.3)}} \\ \implies& P(\theta = 0.7 | X = 1) = \frac{0.7 \times 0.2}{0.7\times0.2 + 0.3\times0.8} \\ &= 0.368 \end{aligned} \] Before observing heads, we would have thought the \(\theta = 0.7\) coin to be very unlikely, but because the observed data favours \(\theta = 0.7\) more strongly than \(\theta = 0.3\), after observing these data we feel that \(\theta = 0.7\) is more likely than before.

Flipping More Coins

Suppose now that we flip \(n\) coins, obtaining a dataset \(X = (X_{1},\ldots,X_{n})\) of heads or tails, represented by 0’s and 1’s. If we’re still considering only two candidate values \(\theta = 0.7\) or \(\theta = 0.3\), we may still ask the question “which value is more likely to have generated the observed data?”. We again form the likelihood function for each value of \(\theta\), the relative frequency with which each value of \(\theta\) would have generated the observed sample. Assuming the tosses are statistically independent: \[ p(X|\theta) = \theta^{\sum_{i=1}^{n}X_{i}} \times (1 - \theta)^{n - \sum_{i=1}^{n}X_{i}} \] where \(\sum_{i=1}^{n}X_{i}\) is just the number of heads observed in the sample. We see that any two samples that have the same number of heads will lead to the same inferences about \(\theta\) in this manner.

Suppose we throw the coin \(10\) times and observe \(6\) heads. The likelihood function for each candidate value of \(\theta\) is \[ \begin{aligned} p(X|\theta = 0.7) &= 0.7^{6} \times 0.3^{4} = 0.000953 \\ p(X|\theta = 0.3) &= 0.3^{6} \times 0.7^{4} = 0.000175 \\ \end{aligned} \] It is much more likely to observe \(6\) heads when \(\theta = 0.7\) than when \(\theta = 0.3\).

In the Bayesian setting, with our prior distribution on \(\theta\) from above, we would form the posterior distribution as follows: \[ p(\theta|X) = \frac{\theta^{\sum_{i=1}^{n}x_{i}}(1-\theta)^{n - \sum_{i=1}^{n}x_{i}}\times 0.2^{I(\theta = 0.7)}0.8^{I(\theta = 0.3)}}{\sum_{\theta = 0.3,0.7}\theta^{\sum_{i=1}^{n}x_{i}}(1-\theta)^{n - \sum_{i=1}^{n}x_{i}}\times 0.2^{I(\theta = 0.7)}0.8^{I(\theta = 0.3)}} \] Computing this for our observed data of \(\sum_{i=1}^{10}x_{i} = 6\) yields \[ \begin{aligned} p(\theta = 0.7 |X) = \frac{0.7^{6}0.3^{4}\times 0.2}{0.7^{6}0.3^{4}\times 0.2 + 0.3^{6}0.7^4\times0.8} = 0.576 \\ p(\theta = 0.3 |X) = \frac{0.3^{6}0.7^{4}\times 0.2}{0.3^{6}0.7^{4}\times 0.8 + 0.7^{6}0.3^4\times0.2} = 0.424 \\ \end{aligned} \]

We can see that the data “updates” our prior belief that \(\theta = 0.3\) was more probable than \(\theta = 0.7\), because the observed data was more likely to have occurred if \(\theta = 0.7\) than if \(\theta = 0.3\).

It is helpful to visualize the prior and posterior, for the observed data. Because both prior and posterior only allow two values, we can do this using a simple bar chart:

visualize_binomial_priorposterior <- function(sumx,n) {
  prior <- function(theta) {
    if (theta == .3) {
      return(.8)
    }
    else if (theta == .7) {
      return(.2)
    }
    0
  }
  likelihood <- function(theta) theta^sumx * (1-theta)^(n - sumx)
  marginal_likelihood <- prior(.7) * likelihood(.7) + prior(.3) * likelihood(.3)
  posterior <- function(theta) likelihood(theta) * prior(theta) / marginal_likelihood
  
  # Plot of the prior and posterior distributions for these observed data
  data_frame(
    theta = c(.3,.7,.3,.7),
    value = c(prior(.3),prior(.7),posterior(.3),posterior(.7)),
    type = c("Prior","Prior","Posterior","Posterior")
  ) %>%
    ggplot(aes(x = theta,y = value,fill = type)) +
    theme_classic() + 
    geom_bar(stat = "identity",position = "dodge",colour = "black") +
    labs(title = "Prior and Posterior for theta",
         subtitle = str_c("Observed data: ",sumx," flips in ",n," throws"),
         x = "Theta, probability of heads",
         y = "Prior/Posterior Probability",
         fill = "") +
    scale_x_continuous(breaks = c(0.30,0.70),labels = c("0.30","0.70")) +
    scale_y_continuous(labels = scales::percent_format()) +
    scale_fill_brewer(palette = "Reds")
  
}

Plotting is nice as it lets us compare how different observed data, and different experiments (number of throws) affect the prior/posterior balance of belief:

cowplot::plot_grid(
  visualize_binomial_priorposterior(6,6),
  visualize_binomial_priorposterior(6,10),
  visualize_binomial_priorposterior(6,20),
  visualize_binomial_priorposterior(6,50),
  visualize_binomial_priorposterior(0,10),
  visualize_binomial_priorposterior(1,10),
  visualize_binomial_priorposterior(7,10),
  visualize_binomial_priorposterior(10,10),
  ncol=2
)

Estimating \(\theta\)

The case where there are only two possible parameter values is useful for examples, but not common in practice. More realistic is the case where \(0 < \theta < 1\), and we need to use the observed data to estimate \(\theta\).

Intuition and frequentist results dictate that a “good” estimator of \(\theta\) is the sample proportion number of heads, \[ \hat{\theta}_{freq} = \frac{1}{n}\sum_{i=1}^{n}X_{i} \] This makes sense; of course the best guess of the value of the probability of heads based on one sample is just the relative frequency with which \(X\) is heads in that sample. It’s also the maximum likelihood estimator, and the unbiased estimator with minimum variance. Let’s see how this estimator, which is optimal from a frequentist perspective, behaves compared to what we come up with using Bayesian estimation.

# Generate plot data- dbeta evaluated at x and a,b for x on a grid between 0 and 1 and various values of a,b
# Run this code line by line if you have trouble reading it and figuring out what it does
expand.grid(a = c(.5,1,2),b = c(.5,1,2)) %>%
  pmap_dfr(~data_frame(x = seq(0.01,0.99,by=0.01),y = dbeta(x = x,shape1=.x,shape2=.y),a = .x,b = .y)) %>%
  ggplot(aes(x = x,y = y)) +
  theme_classic() +
  facet_wrap(a~b) +
  geom_line(colour = "purple") +
  labs(title = "Beta Distribution, various a and b",
       subtitle = "Top value is a, bottom is b",
       x = "Datapoint x",
       y = "Density f(x;a,b)") +
  theme(text = element_text(size = 22))

data_frame(x = c(0.01,0.99)) %>%
  ggplot(aes(x = x)) +
  theme_classic() + 
  stat_function(fun = dbeta,
                args = list(shape1 = 12,shape2 = 12),
                colour = "blue") +
  labs(title = "Beta Prior for Theta",
       subtitle = "Bayesian Coin Flipping Example",
       x = "Theta",
       y = "Prior Density, p(Theta)") +
  scale_x_continuous(breaks = seq(0,1,by=0.1))

# Prior probability of theta being between 0.3 and 0.7
pbeta(0.7,shape1=12,shape2=12) - pbeta(0.3,shape1=12,shape2=12)

## [1] 0.957104

The Posterior

Our prior has density \[ p(\theta;a,b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\theta^{a-1}(1-\theta)^{b-1} \] Our likelihood remains the same as before: \[ p(X|\theta) = \theta^{\sum_{i=1}^{n}x_{i}}(1-\theta)^{n - \sum_{i=1}^{n}x_{i}} \] Bayes’ Rule is still used to compute the posterior from these quantities; however, it looks slightly different now: \[ \begin{aligned} p(\theta|X) &= \frac{p(X|\theta)p(\theta)}{p(X)} \\ &= \frac{p(X|\theta)p(\theta)}{\int_{0}^{1}p(X|\theta)p(\theta)d\theta} \end{aligned} \]

Now, because \(\theta\) is defined on a continuous interval, the marginal likelihood/model evidence/normalizing constant is computed via integrating the joint distributionn of \(X\) and \(\theta\) over the range of \(\theta\).

In this example, the marginal likelihood and the posterior can be computed explicitly as follows: \[ \begin{aligned} p(X) &= \int_{0}^{1}p(X|\theta)p(\theta)d\theta \\ &= \int_{0}^{1} \theta^{\sum_{i=1}^{n}x_{i}}(1-\theta)^{n - \sum_{i=1}^{n}x_{i}} \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\theta^{a-1}(1-\theta)^{b-1} d\theta \\ &= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \times \int_{0}^{1} \theta^{\sum_{i=1}^{n}x_{i} + a - 1}(1-\theta)^{n - \sum_{i=1}^{n}x_{i} + b - 1} d\theta \\ &= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \times \frac{\Gamma(\sum_{i=1}^{n}x_{i} + a)\Gamma(n - \sum_{i=1}^{n}x_{i} + b)}{\Gamma(n + a + b)} \end{aligned} \]

How did we evaluate that integral and get to the last line? We recognized the integrand as being the \(\theta\)-dependent part of a \(Beta(\sum_{i=1}^{n}x_{i} + a,n - \sum_{i=1}^{n}x_{i} + b)\) density; hence it integrates to the reciprocal of the appropriate normalizing constant. This trick is commonly used in examples illustrating Bayesian inference; it shouldn’t be taken from this, though, that this integral is always easy to evaluate like this. It is almost never easy, or even possible, to evaluate this integral in anything beyond these simple examples- more on this later.

With \(p(X)\) available, we can explicitly compute the posterior: \[ \begin{aligned} p(\theta|X) &= \frac{p(X|\theta)p(\theta)}{p(X)} \\ &= \frac{\theta^{\sum_{i=1}^{n}x_{i}}(1-\theta)^{n - \sum_{i=1}^{n}x_{i}} \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\theta^{a-1}(1-\theta)^{b-1}}{\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \times \frac{\Gamma(\sum_{i=1}^{n}x_{i} + a)\Gamma(n - \sum_{i=1}^{n}x_{i} + b)}{\Gamma(n + a + b)}} \\ &= \frac{\Gamma(n + a + b)}{\Gamma(\sum_{i=1}^{n}x_{i} + a)\Gamma(n - \sum_{i=1}^{n}x_{i} + b)} \times \theta^{\sum_{i=1}^{n}x_{i} + a - 1}(1-\theta)^{n - \sum_{i=1}^{n}x_{i} + b - 1} \end{aligned} \] which we recognize as a \(Beta(a + \sum_{i=1}^{n}x_{i},b + n - \sum_{i=1}^{n}x_{i})\) distribution. Priors where the posterior belongs to the same family of distributions as the prior are called conjugate priors, and while they represent the minority of practical applications, they are very useful for examples.

In this scenario, we can interpret the likelihood as directly updating the prior parameters, from \((a,b)\) to \((a + \sum_{i=1}^{n}x_{i},b + n - \sum_{i=1}^{n}x_{i})\). Let’s visualize this for a few different datasets and sample sizes, for our chosen prior:

prior <- function(theta) dbeta(theta,shape1 = 12,shape2 = 12)
posterior <- function(theta,sumx,n) dbeta(theta,shape1 = 12 + sumx,shape2 = 12 + n - sumx)

data_frame(x = c(0.01,0.99)) %>%
  ggplot(aes(x = x)) +
  theme_classic() + 
  stat_function(fun = prior,
                colour = "blue") +
  stat_function(fun = posterior,
                args = list(sumx = 5,n = 10),
                colour = "purple") +
  stat_function(fun = posterior,
                args = list(sumx = 0,n = 10),
                colour = "red") +
  stat_function(fun = posterior,
                args = list(sumx = 10,n = 10),
                colour = "orange") +
  labs(title = "Beta Prior vs Posterior for Theta, 10 coin flips",
       subtitle = "Blue: Prior. Purple: 5 heads. Red: 0 heads. Orange: 10 heads",
       x = "Theta",
       y = "Density") +
  scale_x_continuous(breaks = seq(0,1,by=0.1))

Some interesting points can be seen from this plot:

• When the observed data “matches” the prior, in the sense that we observe a dataset for which the original frequentist estimate of \(\theta\) is very probable under the prior, the posterior becomes more peaked around that value.
• When the observed data are extreme, as in the case of 0 or 10 heads in 10 flips, the frequestist inference would also be extreme. In these cases, we would have estimated \(\hat{\theta}_{freq} = 0\) or \(1\). Because our prior distribution is not extreme, though, the posterior is more moderate, and ends up being peaked at the low/high end of the range of values that were reasonable under the prior.

Check out this shiny app to get a better feel for how the prior, observed data and sample size influence the location and shape of the posterior under this model.

# E.g. for n = 10 and sumx = 5
c(qbeta(0.025,shape1=12 + 5,shape2 = 12 + 5),qbeta(0.975,shape1=12 + 5,shape2 = 12 + 5))

## [1] 0.3354445 0.6645555

c(.5 - sqrt(.5*.5/10)*1.96,.5 + sqrt(.5*.5/10)*1.96)

## [1] 0.1900968 0.8099032

c(.1 - sqrt(.1*.9/10)*1.96,.1 + sqrt(.1*.9/10)*1.96)

## [1] -0.08594193  0.28594193

# Point estimate- find the posterior mode, which is a critical point
# Use R's built in optimization tools, function nlminb() performs constrained minimization
# Pass it a function that returns minus the posterior; minimizing -f(x) is the same as maximizing f(x)
minus_posterior <- function(theta,sumx,n) -1 * posterior(theta,sumx,n)
opt <- nlminb(objective = minus_posterior,start = 0.3,sumx = 1,n = 10,lower = 0.01,upper = 0.99)
opt$par # Return the value at which the maximum occurs, i.e. the posterior mode

## [1] 0.375

# Credible interval
c(qbeta(0.025,shape1=12 + 1,shape2 = 12 + 9),qbeta(0.975,shape1=12 + 1,shape2 = 12 + 9))

## [1] 0.2290662 0.5487546

Conjugate Priors

A conjugate prior, in relation to a specific likelihood, is a prior that when combined with that likelihood gives a posterior with the same functional form as that prior. The Beta/Bernoulli we saw above is an example of this, because we found: \[ \begin{aligned} \mbox{Prior: } & p(\theta) \propto \theta^{a-1}(1-\theta)^{b-1} \\ \mbox{Likelihood: } & \ell(X|\theta) \propto \theta^{\sum_{i=1}^{n}x_{i}}(1-\theta)^{n - \sum_{i=1}^{n}x_{i}} \\ \implies \mbox{ Posterior: } & p(\theta | X) \propto \theta^{a + \sum_{i=1}^{n}x_{i} - 1}(1-\theta)^{b + n - \sum_{i=1}^{n}x_{i} - 1} \end{aligned} \] The prior has the form \(\theta^{a-1}(1-\theta)^{b - 1}\), and the posterior has the form \(\theta^{c-1}(1-\theta)^{d - 1}\), with \(c\) and \(d\) depending on \(a\) and \(b\) as well as the data. The posterior has the same functional form as the prior, with parameters that are updated after the data is observed.

Conjugate priors are great because they are mathematically tractible, and allow us to easily evaluate the impact of the prior distribution on the posterior under different datasets. It is often not possible, though, to find a conjugate prior for a given likelihood in anything but the most simple examples. Here are some common likelihoods and their conjugate priors; as an exercise, verify that each posterior is in the same family as the prior, and find expressions for the updated parameters:

Wikipedia has a great list containing many more examples.

Setting Hyperparameters

When using a conjugate prior (or any prior), once a family of distributions like Beta, Normal, Gamma, etc is chosen, the analyst still needs to set hyperparameter values. We did this above- first we chose a \(Beta(a,b)\) family of distributions, then we went a step further and actually specified \(a = 12\) and \(b = 12\). How did we come up with such wonky looking values of \(a\) and \(b\)? We will discuss two ways here.

x <- rbeta(1000,12,12)
c(mean(x),sd(x))

## [1] 0.5062132 0.1031181

n <- 100
N <- 500
truetheta <- 0.4
set.seed(80545)
flips <- data_frame(x = rbinom(N,n,truetheta) / n)

flips %>%
  ggplot(aes(x = x)) +
  theme_classic() + 
  geom_histogram(bins = 30,colour = "black",fill = "purple") +
  labs(title = "Observed Sample Proportions",
       subtitle = str_c(N," repetitions of flipping a coin ",n," times"),
       x = "Sample Proportion of Heads",
       y = "Number of Times")

# No closed-form formula for (a,b); estimate numerically
minusloglik <- function(params) -sum(dbeta(flips$x,shape1 = params[1],shape2 = params[2],log = TRUE))

# Minimize minus the log-likelihood
# Get starting values using the method of moments
xb <- mean(flips$x)
vb <- var(flips$x)
astart <- xb * ((xb * (1-xb))/vb - 1)
bstart <- (1-xb) * ((xb * (1-xb))/vb - 1)
c(astart,bstart)

## [1] 40.12836 59.95234

opt <- nlminb(start = c(astart,bstart),objective = minusloglik,lower=c(0,0))
opt$par

## [1] 40.18212 60.03573

flips %>%
  ggplot(aes(x = x)) +
  theme_classic() + 
  geom_histogram(aes(y = ..density..),bins = 30,colour = "black",fill = "purple") +
  stat_function(fun = dbeta,
                args = list(shape1 = opt$par[1],shape2 = opt$par[2]),
                colour = "orange") +
  labs(title = "Observed Sample Proportions - Beta Distribution Fit by Maximum Likelihood",
       subtitle = str_c(N," repetitions of flipping a coin ",n," times"),
       x = "Sample Proportion of Heads",
       y = "Number of Times")

# Moment-matching parameters: E(theta) = 0.4, SD(theta) = 0.1
a_momatch <- 0.4 * (0.4*0.6 / 0.1^2 - 1)
b_momatch <- 0.6 * (0.4*0.6 / 0.1^2 - 1)
c(a_momatch,b_momatch)

## [1]  9.2 13.8

plt1 <- data_frame(x = c(0,1)) %>%
  ggplot(aes(x = x)) +
  theme_classic() +
  stat_function(fun = dbeta,
                args = list(shape1 = a_momatch,shape2 = b_momatch),
                colour = "purple") +
  stat_function(fun = dbeta,
                args = list(shape1 = opt$par[1],shape2 = opt$par[2]),
                colour = "blue") +
  labs(title = "Comparison of Prior Distributions",
       subtitle = "Blue: Empirical Bayes. Purple: Moment-Matching",
       x = "Theta",
       y = "Prior Density"
  ) +
  scale_x_continuous(breaks = seq(0,1,0.1))

plt2 <- data_frame(x = c(0,1)) %>%
  ggplot(aes(x = x)) +
  theme_classic() +
  stat_function(fun = dbeta,
                args = list(shape1 = a_momatch + n * sum(flips$x),shape2 = b_momatch + N*n - n*sum(flips$x)),
                colour = "purple") +
  stat_function(fun = dbeta,
                args = list(shape1 = opt$par[1] + n * sum(flips$x),shape2 = opt$par[2] + N*n - n*sum(flips$x)),
                colour = "blue") +
  labs(title = "Comparison of Posterior Distributions",
       subtitle = "Blue: Empirical Bayes. Purple: Moment-Matching",
       x = "Theta",
       y = "Posterior Density"
  ) +
  scale_x_continuous(breaks = seq(0,1,0.1))

cowplot::plot_grid(plt1,plt2,nrow=1)

A Bayesian Formulation

There are a number of ways, many of them equivalent to each other, to address this. Following the theme of this lab, we’ll do it in a Bayesian framework: put a prior on \(\beta\), make it so that prior says that large values of the coefficients are extremely unlikely, and then compute and analyze the posterior for \(\beta\) given the prior and the data.

The likelihood for the data \(y\) in this problem is \[ y_{i}|\beta \sim N(x_{i}^{\prime}\beta,\sigma^{2}) \] From the discussion of conjugate priors above, and from your own completion of the associated exercise, we know that the conjugate prior for this distribution is also normal: \[ \beta \sim N(m,\Sigma) \] with \(m = (m_{0},m_{1},\ldots,m_{p})\) representing the prior mean and \(\Sigma\) representing the prior covariance matrix. Moment-matching is easy here, because the prior hyperparameters correspond directly to the first two moments of the prior distribution, so we don’t have to solve any equations like for the Beta example. We can figure this out as follows:

• We’ll take \(m = 0\), because we simply don’t have any reason to believe that the curve we’re fitting takes any particular shape, prior to looking at the dataset.
• We’ll take the prior covariance matrix to be diagonal, \(\Sigma = s^{2} \times I\) partly because multivariate stats is not a prerequisite for this lab; also partly because this is reasonable. Inducing prior covariance is sometimes a great idea, but it’s too complicated a discussion for this example.
• We’ll choose the prior variance, \(s^{2}\), to enforce our belief that the values of the coefficients in the model that generated the data are actually small. What does “small” mean? This is highly subjective; interpreting these regression coefficients as being “standard deviations away from the predicted mean response due to a unit standard deviation change in the covariates”, which is a mouthful, leads to the notion that a \(\beta_{j}\) of 1 is probably getting pretty big. Hence we will choose \(s^{2} = 0.5^{2}\), which gives a prior standard deviation of \(0.5\), and hence any value more than 2 standard deviations away from the prior mean is unlikely under the prior- i.e. any value greater than 1. When running the code here, though, you should definitely change this and note the effect that this has on the resulting curve, because this is super important.

For simplicity, we will take the data variance \(\sigma^{2}\) to be fixed at its sample value, \(\hat{\sigma}^{2} = 1\). It equals 1 precisely because we divided the data by the sample variance, when we standardized. A fully Bayesian treatment would put a prior on this as well.

Using standard results about joint and conditional normal random variables leads us to the posterior distribution: \[ \begin{aligned} \beta | y &\sim N(m_{post},\Sigma_{post}) \\ m_{post} &= \left( \frac{1}{s^{2}}I + \frac{1}{\sigma^{2}}X^{\prime}X\right)^{-1}\left( \frac{1}{s^{2}}m + \frac{1}{\sigma^{2}}X^{\prime}y\right) \\ \Sigma_{post} &= \left( \frac{1}{s^{2}}I + \frac{1}{\sigma^{2}}X^{\prime}X\right)^{-1} \end{aligned} \] Here \(X\) is the design matrix, with rows consisting of the individual \(x_{i}\). Because the posterior is normal, both the mean and mode are equal to \(m_{post}\)- our Bayesian point estimate of \(\beta\).

From normal regression theory, if we fit a standard (non-Bayesian) linear regression, the resulting point estimate of \(\beta\) would be \[ \hat{\beta}_{normal} = \left( X^{\prime}X \right)^{-1}X^{\prime}y \] with a Gaussian sampling distribution, and variance \[ Var\left( \hat{\beta}_{normal} \right) = \sigma^{2}\left( X^{\prime}X \right)^{-1} \] We can interpret the effects of the prior distribution intuitively as shrinking the frequentist point estimate towards the prior mean, by an amount depending on the data variance and the prior variance.

Let’s compare the Bayesian regression to the normal-theory frequentist regression. For different values of \(p\), we’ll plot the curves corresponding to the Bayesian posterior mean and the frequentist point estimate, and pointwise credible (Bayesian) and confidence (frequentist) intervals. A pointwise interval for \(\beta\) at a given \(x\) is calculated as \[ \hat{\beta} \pm 2\times x^{\prime}Var(\hat{\beta})x \] These are meant to give a rough idea of the range of plausible curves at each given \(x\); they are misleading in that when calculated this way they don’t actually give a measure of the variability associated with the entire curve. Nonetheless they are useful for comparing estimated curves. Along with the plots, we’ll compute the actual model coefficients, and the residual sum of squares associated with the fit for that value of \(p\). We can then compare the two approaches as \(p\) gets high, and make a subjective choice as to the value of \(p\) that we prefer.

All of this takes some work to code up, so here we go! We’re going to need

• A function to estimate the Bayesian linear regression for a given p, returning the posterior mean and variance
• A function to estimate the frequentist linear regression for a given p, returning the point estimate and variance
• A function to take a value of p, fit the above two models, and return a dataframe consisting of the observed responses, their predictions under the model (the value of the curve at that x point), and the lower and upper points of the associated interval, for each model
• A function to plot the above data, and
• A script that fits the models for various values of \(p\), returning a dataframe with the residual sum of squares for the model

# Function to fit the Bayesian linear regression described above
# Returns the posterior mean and variance
fit_bayesian_lm <- function(p,m=rep(0,p+1),s=0.5) {
  # p is the degree of polynomial to fit
  # m is the prior mean, and s the prior standard deviation
  
  m <- cbind(m)
  X <- model.matrix(answer_speed ~ poly(date,degree = p,raw = TRUE),data = contact_std)
  y <- cbind(contact_std$answer_speed)
  
  # Posterior mean and variance
  varpost <- solve((1/s^2)*diag(dim(X)[2]) + t(X) %*% X)
  meanpost <- varpost %*% ((1/s^2)*m + t(X) %*% y)
  
  list(
    estimate = meanpost,
    var = varpost
  )
}

# Function to fit the frequentist linear model
# Returns the point estimate and variance
fit_frequentist_lm <- function(p) {
  mod <- lm(answer_speed ~ poly(date,degree = p,raw = TRUE),data = contact_std)
  list(
    estimate = cbind(coef(mod)),
    var = vcov(mod)
  )
}

# Function to predict new y at given x
# For given x, returns predicted y and lower/upper bounds of a 2-sd interval
predict_with_interval <- function(model_data,x,p) {
  # model_data is the output of one of the above functions
  # x here is a scalar
  # Create a vector of length p+1 containing 1, x, x^2, ... , x^p
  xvec <- c(1)
  for (i in 1:p) {
    xvec <- c(xvec,x^i)
  }
  xvec <- cbind(xvec)
  
  vn <- t(xvec) %*% model_data$var %*% xvec
  mn <- t(xvec) %*% model_data$estimate
  
  c(
    x = x,
    predicted = mn,
    lower = mn - 2*vn,
    upper = mn + 2*vn
  )
}

# Now, a function that will take a value of p, and fit the above two models, 
# returning predictions and pointwise interval bounds for each training
# set point
fit_both_models <- function(p,...) {
  # The ... can be used to pass additional arguments to the fit_bayesian_lm function,
  # like a different prior mean and variance
  bayesian_lm <- fit_bayesian_lm(p,...)
  frequentist_lm <- fit_frequentist_lm(p)
  
  # Get all unique training set points, and map the predict_with_interval function
  # over them, reducing the results into a data_frame. Do for each of frequentist and bayesian.
  # Also, manually add back on the observed values of y
  bayes_predictions <- contact_std %>%
    pull(date) %>%
    map(~predict_with_interval(bayesian_lm,.x,p)) %>%
    reduce(bind_rows) %>%
    bind_cols(dplyr::select(contact_std,observed = answer_speed))
  freq_predictions <- contact_std %>%
    pull(date) %>%
    map(~predict_with_interval(frequentist_lm,.x,p)) %>%
    reduce(bind_rows) %>%
    bind_cols(dplyr::select(contact_std,observed = answer_speed))
  
  # Return a list containing the model parameters and predictions
  list(
    bayesian = list(
      predictions = bayes_predictions,
      model = bayesian_lm
    ),
    frequentist = list(
      predictions = freq_predictions,
      model = frequentist_lm
    )
  )
}

# Fit the models for several values of p
ptotry <- c(1:10)
model_fits <- ptotry %>%
  map(~list(
            p = .x,
            data = fit_both_models(p = .x))
  )

# Plot them
plot_one_p <- function(model) {
  # Stack the data for the frequentist and bayesian together, create a new variable
  # indicating frequentist/bayesian. This allows us to make a single ggplot with
  # two frames, one for each curve
  plt <- model$data$bayesian$predictions %>%
    mutate(type = "Bayesian") %>%
    bind_rows(model$data$frequentist$predictions %>% mutate(type = "Frequentist")) %>%
    ggplot(aes(x = x,colour = type,fill = type)) +
    theme_classic() + 
    facet_grid(~type) +
    geom_ribbon(aes(ymin = lower,ymax = upper),alpha = 0.3) +
    geom_point(aes(y = observed),colour = "black",fill = "grey",pch = 21) +
    geom_line(aes(y = predicted)) + 
    labs(title = "Contact Time Model",
         subtitle = str_c("p = ",model$p),
         x = "Date (Standardized)",
         y = "Contact Time (Standardized)",
         fill = "Type") +
    scale_fill_brewer(palette = "Set1") +
    scale_colour_brewer(palette = "Set1") +
    guides(colour = FALSE) +
    ylim(c(-2.5,2.5))
  
  return(plt)
}

# Plot them all
model_fits %>%
  map(plot_one_p) %>%
  cowplot::plot_grid(plotlist = .,ncol = 1)

As the frequentist curve gets more wiggly, the Bayesian curve stays reasonably stable. What about for higher \(p\)?

ptotry <- c(11:20)
ptotry %>%
  map(~list(
            p = .x,
            data = fit_both_models(p = .x))
  ) %>%
  map(plot_one_p) %>%
  cowplot::plot_grid(plotlist = .,ncol = 1)

While it’s probably not a good idea to fit models with such high degree polynomials, at least the Bayesian implementation guards somewhat against high variability in the fitted curves.

Based on this, we can decide what degree of polynomial to use. Since the quadratic is shaped almost exactly like the curves of higher degrees (up to around 10), we can stick with the simplest reasonable model and choose that one; we also could be more formal about it and split the data into a training and a test set, fit the polynomials on each, and report the test set error.