有两种方法可以解决这个问题。首先，使用相关的过去数据以某种方式“自动”创建先验（或以某种方式将这些相关数据包含到我们的新数据的单个模型中）。这个选项通常被认为是有吸引力的，因为它“对它有一定的客观性”。其次，询问专家（在向他们展示他们可能需要记住的任何相关数据之后）。最后，但可能不太相关，可以选择使用弱信息先验（或试图不提供信息的先验）。

在第一类方法中，Schmidli 等人的（稳健的）元分析预测（MAP）先验。已经被提及并且经常使用 - 特别是在添加了额外弱/无信息混合组件的稳健版本中 - 但有各种变体，替代方案，如自适应功率先验，将单个模型拟合到旧模型和新模型的想法以一种对先前数据冲突和其他类似想法具有鲁棒性的方式获取数据。

在第二类方法中，有很多方法可以从专家那里获得先前的意见，以最大限度地减少人们（包括专家）受到的偏见（= “专家启发”）。一个这样的框架是SHELF，你可以在他们的网页上找到整个课程，并且还有一个R 包。我特别提到了一个，因为我在实践中使用它，但还有其他具有不同风格/哲学的。

这里有一些在实践中设置的先验示例，主要用于临床试验/药物开发（仅仅是因为我在那里最熟悉它——更多示例请参见本书）：对于COPD 中的概念验证研究，用于类风湿性关节炎的概念证明（另一个也用于 RA），用于历史数据中的指数风险，用于临床试验中的治疗效果以及用于预测事件发生率和计数结果的离散参数. 在制药行业，使用先验信息和专家知识在临床开发早期分析研究（例如概念验证研究分析和决定是否继续进行）或后期决策时尤其常见，而验证性研究则很少见这旨在支持监管部门的批准（部分原因是，在进行内部决策时，过度乐观的先验对公司来说更成问题，而监管机构则对选择用于验证性研究的先验进行更严格的审查）。

在这里，只是想添加一些补充材料并演示以下内容：使用 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. 频率回归：这就是频率回归模型（即使用普通最小二乘法 (OLS) 估计参数的回归模型——我们在学校都学过的）。

首先，拟合回归模型：

#fit regression model

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

接下来，查看结果：

#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

可选：可视化结果

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）贝叶斯回归：现在，我们尝试将贝叶斯回归模型拟合到相同的数据：

#load library

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

首先，我们指定身高和体重变量的先验（我为它们选择了一个正态分布 - 在我最初的问题中，我们会通过使用其他生物学家对长颈鹿所做的研究来决定这些先验）：

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

接下来，我们运行贝叶斯回归模型

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

现在，我们可以查看结果：

 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

看模型表现：

#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 

可选：可视化结果

#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) 

参考资料：

• 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

注意：我仍在学习贝叶斯回归——请纠正我可能犯的任何错误（例如，由于我选择先验，贝叶斯回归模型的性能似乎比线性回归模型差得多？当我运行具有默认先验（“弱信息先验”）的贝叶斯回归模型，例如model_2 <- stan_glm(age~., data=my_data, seed=111)- 贝叶斯线性回归的结果与线性回归模型相当。我一定做错了什么？）。

谢谢！