贝叶斯随机参数模型R语言实现
本文通过使用R,主要通过使用tidyverse包来进行数据清洗处理以及绘图,使用brms包来实现Bayesian回归模型。
1. 贝叶斯数据分析的基本步骤
-
确定与研究问题相关的数据。数据的度量尺度是什么?哪些数据变量是因变量,哪些数据变量应该是解释变量?
-
为相关数据建立一个描述性模型,给予数学形式及其参数描述。
-
指定参数的先验分布。
-
使用贝叶斯推理在参数值之间重新分配可信度。从理论上解释有意义的变量的后验分布(假设模型是对数据的合理描述)。
-
检查后验预测是否以合理的准确度模拟数据(即进行“后验预测检查”)。如果不是,则考虑不同的描述模型。
2. Example
2.1导入数据
# 一个集合了众多数据分析常用包的包,详情https://zhuanlan.zhihu.com/p/80732610
library(bruceR)
##
## ⭐ bruceR: BRoadly Useful Convenient and Efficient R functions
##
## Loaded R packages:
## [Data]: rio / dplyr / tidyr / stringr / forcats / data.table
## [Stat]: psych / emmeans / effectsize / performance
## [Plot]: ggplot2 / ggtext / cowplot / see
##
## Frequently used functions in `bruceR`:
## set.wd() / Describe() / Freq() / Corr() / Alpha() / MEAN()
## MANOVA() / EMMEANS() / model_summary() / theme_bruce()
# bruceR 内置的函数import,可以导入csv、excel等多种格式
df <- import("/Users/cpf/Downloads/salary.txt")
# 观察数据
glimpse(df)
## Rows: 62
## Columns: 6
## $ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ time <int> 3, 6, 3, 8, 9, 6, 16, 10, 2, 5, 5, 6, 7, 11, 18, 6, 9, 7, 7, …
## $ pub <int> 18, 3, 2, 17, 11, 6, 38, 48, 9, 22, 30, 21, 10, 27, 37, 8, 13…
## $ sex <int> 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1…
## $ citation <int> 50, 26, 50, 34, 41, 37, 48, 56, 19, 29, 28, 31, 25, 40, 61, 3…
## $ salary <int> 51876, 54511, 53425, 61863, 52926, 47034, 66432, 61100, 41934…
## id: 就是id; time: 获得博士学位到现在的时间; pub: 出版数量;
# sex: 性别 1是女性,0是男性; citation: 被引用数量; salary: 现在的收入
#去掉第一列id, 数据的基本描述统计, 相关系数描述图1
Describe(df[, -1],plot = TRUE)
## Descriptive Statistics:
## ────────────────────────────────────────────────────────────────────────────
## N Mean SD | Median Min Max Skewness Kurtosis
## ────────────────────────────────────────────────────────────────────────────
## time 62 6.79 4.28 | 6.00 1.00 21.00 1.23 1.29
## pub 62 18.18 14.00 | 13.00 1.00 69.00 1.31 1.61
## sex 62 0.44 0.50 | 0.00 0.00 1.00 0.25 -1.97
## citation 62 40.23 17.17 | 35.00 1.00 90.00 0.65 0.32
## salary 62 54815.76 9706.02 | 53681.00 37939.00 83503.00 0.61 0.25
## ────────────────────────────────────────────────────────────────────────────
# 相关系数描述图2
pairs.panels(df[, -1], ellipses = FALSE)