Prior selection
Suppose we have data \(Y\), which we take to be normally distributed: \(y_i \sim \mathcal{N}(\mu_i, \sigma^2)\) for \(i=1,\ldots,n\). Suppose further that we want to model \(\mu_i\) as a linear function of a covariate/predictor \(x_i\), using an intercept \(\alpha\) and coefficient \(\beta\). That is, \[ \mu_i = \alpha + \beta x_i. \]
Let’s simulate data (with very small sample size) under this model.
create_df <- function(post_draws, prior_draws){
post_draws <- data.frame(post_draws)
post_draws$distribution <- "posterior"
prior_draws <- data.frame(prior_draws)
colnames(prior_draws) <- "alpha"
prior_draws$distribution <- "prior"
dat <- rbind(post_draws, prior_draws)
return(dat)
}set.seed(689934)
alpha <- 1
beta <- -0.25
sigma <- 1
N <- 5
x <- array(runif(N, 0, 2), dim=N)
y <- array(rnorm(N, beta * x + alpha, sigma), dim=N)Note the true values of \(\alpha\) and \(\beta\). Once we observe the data \(\{x_i,y_i\}\), both \(\alpha\) and \(\beta\) are unknown, so we would like to perform inference on them. However, we only have 5 data points. Given so few data points, we can be quite sure that the resulting posteriors will be sensitive to the choice of priors.
Flat priors
One possible choice of prior is a flat prior for \(\alpha\) and \(\beta\). That is, \(\pi(\alpha) \propto 1\) and \(\pi(\beta) \propto 1\). Let’s look at how our posterior beliefs about \(\alpha\) and \(\beta\) (especially \(\alpha\)) act under these priors.
stan_dat <- list(y = y, x=x, N=N)
fit.flat <- stan(file = "lab-04-flat_prior.stan", data = stan_dat, chains = 1, refresh = 0, iter = 2000, warmup = 500, seed=48)
alpha.flat <- as.matrix(fit.flat, pars = "alpha")
beta.flat <- as.matrix(fit.flat, pars = "beta")ggplot(alpha.flat %>% as.data.frame, aes(x = alpha)) +
geom_histogram(bins = 30)
print(fit.flat, pars = c("alpha"))## Inference for Stan model: lab-04-flat_prior.
## 1 chains, each with iter=2000; warmup=500; thin=1;
## post-warmup draws per chain=1500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 98% n_eff Rhat
## alpha 0.78 0.09 1.6 -2.6 0.03 0.8 1.5 4 315 1
##
## Samples were drawn using NUTS(diag_e) at Mon May 10 12:31:09 2021.
## For each parameter, 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).Notice how the posterior for \(\alpha\) is quite diffuse – there is much uncertainty about what \(\alpha\) is. While the true value for \(\alpha\) is 1, what is the posterior mean? What is the 95% credible interval?
- Write down the posterior means of \(\alpha\) and \(\beta\). Give 95% credible intervals for each. Considering the amount of data we have, do the results seem surprising?
By doing inference with a flat/diffuse prior, we might have thought we were using the least prior information possible. However, flat priors may actually bias our estimates for a parameter by allowing the posterior to be pulled towards extreme and unlikely values.
Diving a bit deeper, consider another flat prior: \(\alpha \sim Unif(a, b)\). Under this prior, we are saying that we believe \(a \leq \alpha \leq b\). We exhibit this in the following code, with the prior \(\alpha \sim Unif(-10, 10)\)
stan_dat <- list(y = y, x=x, N=N, lb = -10, ub = 10)
fit.unif <- stan(file = "lab-04-unif_prior.stan", data = stan_dat, chains = 1, refresh = 0, iter = 2000, warmup = 500, seed=48)## Warning: There were 3 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: Examine the pairs() plot to diagnose sampling problemsalpha.unif <- as.matrix(fit.unif, pars = c("alpha"))
beta.unif <- as.matrix(fit.unif, pars = c("beta"))ggplot(alpha.unif %>% as.data.frame, aes(x = alpha)) +
geom_histogram(bins = 30)
print(fit.unif, pars = c("alpha"))## Inference for Stan model: lab-04-unif_prior.
## 1 chains, each with iter=2000; warmup=500; thin=1;
## post-warmup draws per chain=1500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 98% n_eff Rhat
## alpha 0.65 0.1 1.5 -2.8 0.05 0.78 1.4 3.5 214 1
##
## Samples were drawn using NUTS(diag_e) at Mon May 10 12:31:29 2021.
## For each parameter, 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).While the posterior mean under this uniform prior is closer to the true value, the posterior is still very spread out. So we have seen that a diffuse or flat prior is necessarily non-informative, and in these cases is actually extremely informative! Diffuse priors inherently spread probability mass across large regions of parameter space. We often assume that the data will overwhelm the prior, so a diffuse prior will let the data dominate posterior inference. However as showcased here, having only a small amount of observed data may allow the diffuse prior to become informative. Therefore, it would be wise to make the conscious choice to have an informative prior. However, we might specify just how informative we would like our prior beliefs to be.
Weakly informative priors
We often must consider the scale of the parameters we wish to estimate. In applied problems where we know how to interpret the parameters, the scale is easier to identify. We may consider a weakly informative prior to be such that if there is a reasonably large amount of data, the likelihood will dominate the posterior and the prior is not important. This sort of prior ought to rule out unreasonable parameter values, but is not so strong as to rule out possible values which might make sense. As a general rule, it is wise to not use hard constraints unless the bounds represent true constraints (ex. bounding a prior for a variance parameter below by 0). As an example, we might think that \(\alpha\) could be between 0 and 1. Instead of setting the prrior \(\alpha \sim Unif(0,1)\), it would be wise to use a Normal(0.5, 1) prior instead. In the following, we consider some weakly informative priors for the parameters. The data we have are on unit scale, so we consider priors also on the unit scale.
Light-tailed
In this section, we consider a \(N(0,1)\) prior.
stan_dat <- list(y = y, x=x, N=N)
fit.norm <- stan(file = "lab-04-normal_prior.stan", data = stan_dat, chains = 1, refresh = 0, iter = 2000, warmup = 500, seed=49)
alpha.norm<- as.matrix(fit.norm, pars = c("alpha"))ggplot(alpha.norm %>% as.data.frame, aes(x = alpha)) +
geom_histogram(bins = 30)
print(fit.norm, pars = c("alpha"))## Inference for Stan model: lab-04-normal_prior.
## 1 chains, each with iter=2000; warmup=500; thin=1;
## post-warmup draws per chain=1500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 98% n_eff Rhat
## alpha 0.58 0.02 0.57 -0.56 0.22 0.6 0.98 1.6 625 1
##
## Samples were drawn using NUTS(diag_e) at Mon May 10 12:31:49 2021.
## For each parameter, 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).- Compute the posterior means of \(\alpha\) and \(\beta\). Give 95% credible intervals for each. How does the posterior inference under this \(N(0,1)\) prior compare to the diffuse priors above? How informative is this weakly informative prior?
Heavy-tailed
The Cauchy distribution (the Cauchy(0,1) distribution is just the t-distribution with df=1) has heavier/fatter tails than the Normal distribution.
stan_dat <- list(y = y, x=x, N=N)
fit.cauchy <- stan(file = "lab-04-cauchy_prior.stan",data = stan_dat, chains = 1, refresh = 0, iter = 2000, warmup = 500, seed=55)
alpha.cauchy<- as.matrix(fit.cauchy, pars = c("alpha"))ggplot(alpha.cauchy %>% as.data.frame, aes(x = alpha)) +
geom_histogram(bins = 30)
print(fit.cauchy, pars = c("alpha"))## Inference for Stan model: lab-04-cauchy_prior.
## 1 chains, each with iter=2000; warmup=500; thin=1;
## post-warmup draws per chain=1500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 98% n_eff Rhat
## alpha 0.56 0.03 0.56 -0.54 0.21 0.53 0.91 1.7 465 1
##
## Samples were drawn using NUTS(diag_e) at Mon May 10 12:32:10 2021.
## For each parameter, 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).The following plot displays the posteriors for \(\alpha\) under these two priors
plot_dat <- create_df(alpha.norm, alpha.cauchy) %>%
mutate(distribution = if_else(distribution == "posterior", "Normal","Cauchy"))
ggplot(plot_dat, aes(alpha, fill = distribution)) +
geom_histogram(binwidth = 0.25, alpha = 0.7, position = "identity")+
geom_vline(xintercept = alpha) +
scale_fill_brewer()
The Cauchy prior allocates higher probabiity mass to extreme values as compared to the Normal prior, while still concentrating most of the posterior mass for \(\alpha\) within a desired scale.
- Would you say that a Cauchy prior is more or less informative than a Normal prior (assume that their inter-quartile ranges are comparable)?
