How are Bayesian Priors Decided in Real Life?

Question

I always had the following question: How are Bayesian Priors decided in real life?

I created the following scenario to pose my question: Suppose you are researcher and you are interested in studying if the age of a giraffe can be predicted by the weight and height of a giraffe (e.g. linear regression model : age = b_o + b_1height + b_2weight). You arrive at a national park to record measurements on giraffes - but after only taking measurements on a few giraffes, a terrible storm happens and you have to stop your study. You only had time to measure 15 giraffes:

     weight   height age
1  2998.958 15.26611  53
2  3002.208 18.08711  52
3  3008.171 16.70896  49
4  3002.374 17.37032  55
5  3000.658 18.04860  50
6  3002.688 17.24797  45
7  3004.923 16.45360  47
8  2987.264 16.71712  47
9  3011.332 17.76626  50
10 2983.783 18.10337  42
11 3007.167 18.18355  50
12 3007.049 18.11375  53
13 3002.656 15.49990  42
14 2986.710 16.73089  47
15 2998.286 17.12075  52

Unfortunately, this is not enough information to complete your study. However, you do some research and find these kinds of measurements have been taking on giraffes in the past. For example:

Study 1: A study was done in the 1800's which measured 1000 giraffes and found that that average height of those giraffes was 17ft, the average weight was 2800 lbs and the average age was 35. However this was done in the 1800's and you are doubtful that measurement might not have been as accurate back then, and issues in the environment (e.g. poaching) might have causes giraffes to change in size.

Study 2: A study was done in 2010 were 50 giraffes in zoos across the world and their height was 16ft, weight was 300 lbs and age was 50 years. This study is more recent, but you are skeptical that giraffes in zoos might be different from giraffes in the wild.

Study 3: An expert on giraffe strongly believes that the age, height and weight of giraffes have bell shaped distributions. The expert also believes that giraffes keep growing their whole life, i.e. as age increases, so does weight and height. He does not have any concrete numbers, but he is considered the leading expert.

etc.

Question: In this problem, is it possible to complement the limited measurements you have, along with the prior knowledge available on giraffes (while taking into consideration their reliability)? Is this problem an example of how Bayesian Models (e.g. Bayesian Regression) can be used in real life - or does this problem fundamentally lack enough data to work with?

Suppose you consult several studies where the heights were recorded and manually assess the credibility of these studies (assigning "low weights" to studies deemed not credible, e.g. adjusted_height = credibility_score * average_recorded_height_in_study):

head(my_data)

 average_recorded_height_in_study credibility_score study_number adjusted_height
1                         13.74253         1.0000000            1       13.742525
2                         20.08053         0.3222523            2        6.470999
3                         13.25037         0.5132335            3        6.800532
4                         15.74946         0.2625349            4        4.134783
5                         11.68657         0.5966327            5        6.972592
6                         17.27276         1.0000000            6       17.272759

There are many tools/packages (e.g. using the R programming language) which can attempt to explore this "prior information" and fit distribtuion

library(fitdistrplus)
library(patchwork)
library(ggplot2)

 fg <- fitdist(my_data$adjusted_height, "gamma")
 fln <- fitdist(my_data$adjusted_height, "lnorm")
fg <- fitdist(my_data$adjusted_height, "gamma")
 fw <- fitdist(my_data$adjusted_height, "weibull")

 par(mfrow = c(2, 2))
 plot.legend <- c("Weibull", "lognormal", "gamma")

a <- denscomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")
b <- qqcomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")
c <- cdfcomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")
d <- ppcomp(list(fw, fln, fg), legendtext = plot.legend, plotstyle = "ggplot")

a+b+c+d

The above analysis could be repeated for the other variables in the study. Here, we could see which distribution fit the data the best (e.g. using the - likelihood), and record the parameter estimates for this distribution.

Is this the right idea behind how priors are incorporated into Bayesian Models in the real world? In this example that I have created, can the information from previous studies be analyzed and used to create priors for a Bayesian Linear Regression?

Thanks

Note: Suppose the 15 giraffes you measured happened to be diseased giraffes and their height/weight measurements are not representative of the general population of giraffes - but perhaps the information encoded within the priors represent a wide range of giraffes. Thus, combining your measurements with the prior information could result in a more realistic model that could generalize to a larger population of giraffes (this fact being unknown to you at this time).

Answer 1

There's two big approaches to this problem. Firstly, using relevant past data to somehow "automatically" create a prior (or to somehow include this relevant data into a single model with our new data). This option is often considered attractive, because it "has a certain objectivity to it". Secondly, to ask experts (after showing them any relevant data they may need to have in mind). Finally, but perhaps less relevant, there's the option of using weakly informative priors (or priors that attempt to be uninformative).

In the first class of approaches, the (robust) meta-analytic predictive (MAP) prior of Schmidli et al. was already mentioned and is used quite often - especially in the robust version with an extra weakly/uninformative-mixture component added -, but there are various variants, alternatives like adaptive power priors, ideas to fit a single model over the old and the new data in a fashion robust to prior-data-conflict, and other similar ideas.

In the second class of approaches, there's many ways of getting prior opinions out of experts in ways that minimize the biases that people (including experts) are subject to (="expert elicitation"). One such framework is SHELF, on which you can find a whole course on their webpage and for which there's also a R package. I'm mentioning that one specifically, because I use it in practice, but there are others with different flavors/philosophies.

Here's a few examples of priors being set in practice, mostly drawn for clinical trials/drug development (simply because I'm the most familiar with it there - for more examples see e.g. this book): for a proof of concept study in COPD, for a proof of concept in rheumatoid arthritis (and another one also for RA), for an exponential hazard from historical data, for treatment effects in clinical trials and for predicting event rates and dispersion parameter for count outcomes. In the pharmaceutical industry, using prior information and expert knowledge is especially common for analyzing studies early in clinical development (e.g. analysis of proof of concept studies and deciding whether to proceed) or for decision making later on, while it is rarer for the confirmatory studies that are meant to support regulatory approval (in part, an overoptimistic prior is more a problem for the company when it is for internal decision making, while regulatory authorities put priors chosen for confirmatory studies under much more scrutiny).

Answer 2

OP here, just wanted add some supplementary material and demonstrate the following: a comparison between Frequentist Regression and Bayesian Regression using R

#cool trick to directly bring this data into R

my_data <- data.frame(read.table(header=TRUE,
row.names = 1,
text="
                         weight   height age
                      1  2998.958 15.26611  53
                      2  3002.208 18.08711  52
                      3  3008.171 16.70896  49
                      4  3002.374 17.37032  55
                      5  3000.658 18.04860  50
                      6  3002.688 17.24797  45
                      7  3004.923 16.45360  47
                      8  2987.264 16.71712  47
                      9  3011.332 17.76626  50
                      10 2983.783 18.10337  42
                      11 3007.167 18.18355  50
                      12 3007.049 18.11375  53
                      13 3002.656 15.49990  42
                      14 2986.710 16.73089  47
                      15 2998.286 17.12075  52
"))

1. Frequentist Regression : This is how a Frequentist Regression Model (i.e. a Regression Model where the parameters are estimated using Ordinary Least Squares (OLS) - what we all learn in school).

First, fit the regression model:

#fit regression model

model_1 <- lm(age ~ weight + height, data = my_data)

Next, view the results:

#view results

 summary(model_1)

Call:
lm(formula = age ~ weight + height, data = my_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2369 -1.8688  0.3864  2.1065  5.6170 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -525.2843   369.9144  -1.420    0.181
weight         0.1875     0.1238   1.515    0.156
height         0.6871     1.0859   0.633    0.539

Residual standard error: 3.796 on 12 degrees of freedom
Multiple R-squared:  0.1954,    Adjusted R-squared:  0.06135 
F-statistic: 1.457 on 2 and 12 DF,  p-value: 0.2712

Optional : Visualize Results

library(scatterplot3d)

s3d <- scatterplot3d(my_data$weight, my_data$height,my_data$age, pch = 19, type = c("p"), color = "darkgrey",
                     main = "Regression Plane", grid = TRUE, box = FALSE,  
                     mar = c(2.5, 2.5, 2, 1.5), angle = 55)



# regression plane
s3d$plane3d(model_1, draw_polygon = TRUE, draw_lines = TRUE, 
            polygon_args = list(col = rgb(.1, .2, .7, .5)))

# overlay positive residuals
wh <- resid(model_1) > 0
s3d$points3d(my_data$height, my_data$weight, my_data$age, pch = 19)

2) Bayesian Regression: Now, we try to fit a Bayesian Regression Model to the same data:

#load library

library(rstanarm)
library(see)
library(bayestestR)
library(performance)

First, we specify priors on the Height and Weight variables (I picked a normal distribution for both of them - in my original question, we would have decided on these priors by using the research done on giraffes by other biologists):

#specify priors
my_prior <- normal(location = c(3000, 17), scale = c(1, 2))

Next, we run the Bayesian Regression Model

#run bayesian regression model
model_2 <- stan_glm(age~., data=my_data, prior = my_prior,    seed=111)

Now, we can view the results:

 summary(model_2)

Model Info:
 function:     stan_glm
 family:       gaussian [identity]
 formula:      age ~ .
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 15
 predictors:   3

Estimates:
              mean       sd         10%        50%        90%     
(Intercept) -9000290.7     3116.3 -9004290.9 -9000230.6 -8996293.9
weight          2999.7        1.0     2998.4     2999.7     3001.1
height            17.0        2.0       14.4       17.0       19.6
sigma           3207.5       65.0     3124.2     3207.2     3291.0

Fit Diagnostics:
           mean    sd      10%     50%     90%  
mean_PPD    55.5   824.4 -1002.3    66.1  1107.1

Look at the model performance:

#model performance 

 performance(model_2)

# Indices of model performance

ELPD     | ELPD_SE |    LOOIC | LOOIC_SE |     WAIC |    R2 | R2 (adj.) |      RMSE |    Sigma
----------------------------------------------------------------------------------------------
-574.459 | 154.366 | 1148.918 |  308.733 | 1160.324 | 0.983 |    -1.000 | 23876.735 | 3207.163
> se <- sqrt(diag(vcov(model_2)))
> se
    (Intercept)      weight      height 
3116.342642    1.038384    2.040471 

Optional: Visualize Results

#MCMC Trace

x <- as.array(model_2, pars = c("(Intercept)", "height", "weight"))
bayesplot::mcmc_trace(x, facet_args = list(nrow = 2))

#Posterior Distributions

plot_title <- ggplot2::ggtitle("Posterior Distributions")
plot(model_2, "hist", "weight", "height") + plot_title

#confidence ellipse
bayesplot::color_scheme_set("green")
plot(model_2, "scatter", pars = c("height", "weight"),
     size = 3, alpha = 0.5) +
    ggplot2::stat_ellipse(level = 0.9) 

References:

• https://rpubs.com/Qsheep/BayesianLinearRegression
• https://www.theoj.org/joss-papers/joss.01541/10.21105.joss.01541.pdf
• https://cran.r-project.org/web/packages/rstanarm/vignettes/priors.html#default-priors-and-scale-adjustments
• https://mc-stan.org/rstanarm/reference/plot.stanreg.html

Note: I am still learning about Bayesian Regression - please feel to correct any mistakes that I might have made (e.g. It seems like the Bayesian Regression Model is performing far worse than the Linear Regression Model due to my choice of priors? When I run the Bayesian Regression Model with the default priors ("weakly informative priors"), e.g. model_2 <- stan_glm(age~., data=my_data, seed=111) - the results of the Bayesian Linear Regression are comparable with the Linear Regression Model. I must be doing something wrong?).

Thank you!