关于数据挖掘:R语言JAGS贝叶斯回归模型分析博士生延期毕业完成论文时间附代码数据

原文链接：http://tecdat.cn/?p=23652

最近咱们被客户要求撰写对于贝叶斯回归模型的钻研报告，包含一些图形和统计输入。

本文为读者提供了如何进行贝叶斯回归的根本教程。包含实现导入数据文件、摸索汇总统计和回归剖析

在本文中，咱们首先应用软件的默认先验设置。在第二步中，咱们将利用用户指定的先验，对本人的数据应用贝叶斯。

筹备工作

本教程要求:

• 已装置的 JAGS
• 装置 R 软件。
• 假设检验的基本知识
• 相关性和回归的基本知识
• 贝叶斯推理的基本知识
• R 语言编码的基本知识

数据实例

咱们在这个练习中应用的数据是基于一项对于预测博士生实现论文工夫的钻研（Van de Schoot, Yerkes, Mouw and Sonneveld 2013）。钻研人员询问了博士生实现他们的博士论文须要多长时间（n=333）。结果显示，博士学位获得者均匀花了 59.8 个月（5 年 4 个月）来实现他们的博士学位。变量 B3 掂量打算和理论我的项目工夫之间的差别，以月为单位（均匀 =9.97，最小 =-31，最大 =91，sd=14.43）。

对于目前的工作，咱们感兴趣的问题是，博士学位获得者的年龄（M=31.7，SD=6.86）是否与他们我的项目的延期无关。

预计实现工夫和年龄之间的关系是非线性的。这可能是因为在人生的某个阶段（即三十多岁），家庭生存比你在二十多岁时或年长时占用了你更多的工夫。

因而，在咱们的模型中，差距（B3）是因变量，年龄和年龄平方是预测因素。

问题：请写出零假如和备择假如。写下代表这个问题的无效假设和备选假如。你认为哪个假如更有可能？

H0:_年龄与博士我的项目的延期无关。_

H1:  年龄与博士我的项目的延期无关。

H0:_age2 与博士我的项目的延期无关。_

H1:_age2 与博士我的项目的延期无关。_

向下滑动查看后果▼

• *

筹备 – 导入和摸索数据

数据是一个.csv 文件，但你能够应用以下语法间接将其加载到 R 中。

一旦你加载了你的数据，倡议你检查一下你的数据导入是否顺利。因而，首先看看你的数据的汇总统计。你能够应用 describe() 函数。

问题：你所有的数据都被正确地载入了吗？也就是说，所有的数据点都有实质性的意义吗？

describe(data)

描述性统计有意义。

差别。平均值（9.97），SE（0.79）。

年龄。平均值（31.68），SE（0.38）。

age2。平均值（1050.22），SE（35.97）。

向下滑动查看后果▼

**

绘图

在持续剖析数据之前，咱们还能够绘制冀望的关系。

plot(aes(x = age,
             y = diff))

回归

在这个练习中，你将钻研博士生的年龄和 age2 对他们的我的项目工夫延期的影响，这作为后果变量应用回归剖析。

如你所知，贝叶斯推理包含将先验散布与从数据中取得的似然性相结合。指定先验散布是贝叶斯推断中最要害的一点，应该受到高度重视（例如 Van de Schoot 等人，2017）。在本教程中，咱们将首先依赖默认的先验设置。

点击题目查阅往期内容

R 语言实现贝叶斯分位数回归、lasso 和自适应 lasso 贝叶斯分位数回归剖析

左右滑动查看更多

01

02

03

04

要用运行多元回归，首先要指定模型，而后拟合模型，最初取得总结。模型的指定办法如下。

• 咱们想要预测的因变量。
• “~”，咱们用它来示意咱们当初给其余感兴趣的变量。(相当于回归方程的 ”=”）。
• 用求和符号 ’+’ 分隔的不同自变量。
• 最初，咱们插入因变量有一个方差，有一个截距。

上面的代码是如何指定回归模型的。

# 1) 指定模型
 '# 回归模型
                    diff ~ age + age2

                    #显示因变量有方差
                    diff ~~ diff

                    #有一个截距
                    diff ~~ 1'

而后，咱们须要应用以下代码来拟合模型。咱们指定 target = “jags “ 来应用 Jags 而不是 Stan 编译器。

fitbayes(model, data, target = "jags", test = "none", seed = c(12,34,56) )

# test="none" 的输出进行了一些后验的计算，咱们当初不须要，放慢了计算过程。# 种子命令只是为了保障在屡次运行采样器时有雷同的精确后果。你不须要设置这个。当应用 Jags 时，你须要设置尽可能多的种子链（默认）。

当初咱们用 summary(fit.bayes) 来看看总结。

显示输入

频率主义模型与贝叶斯分析模型所提供的后果的确不同。

贝叶斯统计推断和_频率_主义统计办法之间的要害区别在于预计的未知参数的性质。在_频率_主义框架中，一个感兴趣的参数被假设为未知的，但却是固定的。也就是说，假如在人口中只有一个实在的人口参数，例如，一个实在的平均值或一个实在的回归系数。在贝叶斯的主观概率观中，所有的未知参数都被视为不确定的，因而要用一个概率分布来形容。每个参数都是未知的，而所有未知的货色都会失去一个散布。

这就是为什么在_频率_推断中，你次要失去的是一个未知但固定的群体参数的点估计。这是一个参数值，思考到数据，它最有可能呈现在人群中。附带的置信区间试图让你进一步理解这个估计值的不确定性。重要的是要意识到，置信区间只是形成一个模拟量。在从人口中抽取的有限多的样本中，构建（95%）置信区间的程序将使其在 95% 的工夫内蕴含实在的人口值。这并没有为你提供任何信息，即人口参数位于你所剖析的十分具体和惟一的样本中的置信区间边界内的可能性有多大。

在贝叶斯剖析中，你推断的要害是感兴趣的参数的后验散布。它满足了概率分布的每一个属性，并量化了人口参数位于某些区域的概率。一方面，你能够通过它的模式来形容后验的特点。这是一个参数值，思考到数据和它的先验概率，它在人群中是最有可能的。另外，你也能够应用后验的平均数或中位数。应用雷同的散布，你能够构建一个 95% 的置信区间，与_频率_主义统计中的置信区间绝对应。除了置信区间之外，贝叶斯的对应区间间接量化了人口值在肯定范畴内的概率。所关注的参数值有 95% 的概率位于 95% 置信区间的边界内。与置信区间不同，这不仅仅是一个模拟量，而是一个扼要直观的概率申明。

问题：解释预计成果、其区间和后验散布

年龄_仿佛是预测博士延期的一个相干因素，后验均匀回归系数为 2.317，95%HPD（可信区间）[1.194 3.417]。另外，age2 仿佛也是预测博士延期的一个相干因素，后验平均值为 -0.022，95% 可信区间为 [-0.033-0.01]。95% 的 HPD 显示，人口中的这些回归系数有 95% 的概率位于相应的区间内，也请看上面的数字中的后验散布。因为 0 不蕴含在可信区间内，咱们能够相当必定存在影响。

向下滑动查看后果▼

问题：每个贝叶斯模型都应用一个先验散布。形容一下回归系数的先验散布的形态。

查看应用了哪些默认的先验。

（Jags）利用一个十分宽的正态分布来得出这个无信息的先验。默认状况下，平均值为 0，标准差为 10（精度为 0.01）。

向下滑动查看后果▼

• *

**

回归 – 用户指定的先验

你也能够手动指定你的先验散布。实践上，你能够应用你喜爱的任何一种散布来指定你的先验常识。然而，如果你的先验散布不遵循与你的似然雷同的参数模式，计算模型可能会很麻烦。共轭先验防止了这个问题，因为它们采纳了你构建的模型的函数模式。对于你的正态线性回归模型，如果你的回归参数的预设是用正态分布来指定的，就能够达到共轭性（残差失去一个反伽马散布，这里忽略不计）。你能够很灵便地指定信息性先验。

让咱们用共轭先验来从新指定下面练习的回归模型。咱们临时不波及截距和残差的预设。对于你的回归参数，你须要指定其正态分布的超参数，即均值和方差。平均值示意你认为哪一个参数值最有可能。方差示意你对此的确定水平。咱们为 β 年龄回归系数和 β 年龄 2 系数尝试了 4 种不同的先验标准。

首先，咱们应用以下先验。

Age ~ N(3,0.4)

Age2 ~ N(0,0.1)

先验指标是在模型制订步骤中设置的。请留神，精度而不是正态分布的方差。精度是方差的倒数，所以方差为 0.1 对应的精度为 10，方差为 0.4 对应的精度为 2.5。

先验参数在代码中出现如下。

  '# 带有 priors 的回归模型
prior("dnorm(3,2.5)")*age + prior("dnorm(0,10)")*age2

当初拟合模型并提供汇总统计。

答复：

summary(fit)

向下滑动查看后果▼

问题：在下表中填写后果。

答复：

下一步，尝试改编代码，应用其余列的先验标准，而后实现该表。

 '# 有 priors 的回归模型
  ~ prior("dnorm(3,0.001)")*age + prior("dnorm(0,0.001)")*age2

summary(prior2)

summary(prior3)

summary(prior4)

向下滑动查看后果▼

问题：比拟不同先验的后果。后果是否与默认模型有可比性？**
**

问题：应用不同的先验，咱们最终的论断是否类似？

要答复这些问题，按以下步骤进行。咱们能够计算出相对偏差来示意这种差别。咱们将只计算两个回归系数的偏差，比拟默认（无信息）模型和应用 N(20,.4) 和 N(20,.1) 先验的模型。

#1）减去 MCMC 链的内容
fitbayes(what = "mcmc")

#2) 绑定不同的链，计算回归系数的平均值（估计值）。colMeans(as.matrix(mcmc.list)

#3) 计算偏差
100*((estimat-estimat)/estimat)

beta[1,2,1] 和 beta[1,3,1] 指数别离代表了 βage 和 βage2。Beta[x,x,x] 是回归系数（依照咱们在模型中指定的程序，所以首先是 age，而后是 age2），alpha[x,x,x] 是截距，psi[x,x,x] 是方差，def[x,x,x] 是间接效应（如果你的模型中有这些）。它们的排列程序与 summary() 输入中的程序雷同。因而，首先是回归系数，而后是截距，而后是协方差，而后是间接效应。

咱们还能够通过绘制咱们运行的五个不同模型的后验和先验来绘制这些差别。在这个例子中，咱们只绘制年龄 βage 的回归系数。

首先咱们提取 5 个不同模型的 MCMC 链，只针对这一个参数（βage=beta[1,2,1]）。

 binrows(posterior1.5, prior1.5)

而后，咱们能够通过应用以下代码绘制不同的后验和前验。

qplot(data = post,)+
  density(size = 1)+
  scale_x_continuous(limit)

当初，依据表格中的信息、偏差预计和图表，你能够答复对于先验因素对后果的影响的两个问题。

答复：

#1）减去 MCMC 链
fit.bayes(what = "mcmc")

#2) 绑定不同的链，计算回归系数的平均值（估计值）。colMeans(as.matrix(mcmc.list)


#3) 计算偏差
100*((imstim-estima)/estimat)

## beta[1,2,1] beta[1,3,1]. 
## 374.55 397.31

咱们看到，这种高信息量先验的影响对两个回归系数的影响别离为 375% 和 397% 左右。这是一个很大的差别，因而咱们必定不会得出相似的论断。

不同的先验，后果会发生变化，但仍具备可比性。只有对年龄应用 N(20,.4)，才会产生真正不同的系数，因为这个先验均值离数据的均值很远，而其方差却相当确定。然而，一般来说，其余的后果是能够比拟的。因为咱们应用了一个大的数据集，先验的影响绝对较小。如果应用一个较小的数据集，先验的影响就会更大。为了查看这一点，你能够所有案例的大概 20% 进行抽样，而后从新进行同样的剖析。后果当然会不同，因为咱们应用的案例少了很多。应用这段代码。

indices   <- sample.int(333, 60)
smalldata <- data[indices,]

咱们做了一个新的数据集，从原始数据集的 333 个观测值中随机抉择了 60 个。用同样的代码反复剖析，只扭转数据集的名称，以察看先验因素对较小数据集的影响。用这个新的数据集运行 priors2 模型

fit.bayes(model = priors2,smalldata)
summary(fit)

向下滑动查看后果▼

参考文献

Benjamin, D. J., Berger, J., Johannesson, M., Nosek, B. A., Wagenmakers, E.,… Johnson, V. (2017, July 22). Redefine statistical significance_. Retrieved from psyarxiv.com/mky9j_

Greenland, S., Senn, S. J., Rothman, K. J., Carlin, J. B., Poole, C., Goodman, S. N. Altman, D. G. (2016). Statistical tests, P values, confidence intervals, and power: a guide to misinterpretations_._ European Journal of Epidemiology 31 (4_). https://doi.org/10.1007/s1065…** _

van de Schoot R, Yerkes MA, Mouw JM, Sonneveld H (2013) What Took Them So Long? Explaining PhD Delays among Doctoral Candidates_._ PLoS ONE 8(7): e68839. https://doi.org/10.1371/journ…

点击文末 “浏览原文”

获取全文残缺代码数据资料。

本文选自《R 语言 JAGS 贝叶斯回归模型剖析博士生延期毕业实现论文工夫》。

点击题目查阅往期内容

R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真剖析
python 贝叶斯随机过程：马尔可夫链 Markov-Chain，MC 和 Metropolis-Hastings，MH 采样算法可视化
Python 贝叶斯推断 Metropolis-Hastings（M-H）MCMC 采样算法的实现
Metropolis Hastings 采样和贝叶斯泊松回归 Poisson 模型
Matlab 用 BUGS 马尔可夫区制转换 Markov switching 随机稳定率模型、序列蒙特卡罗 SMC、M H 采样剖析工夫序列 R 语言 RSTAN MCMC：NUTS 采样算法用 LASSO 构建贝叶斯线性回归模型剖析职业声望数据
R 语言 BUGS 序列蒙特卡罗 SMC、马尔可夫转换随机稳定率 SV 模型、粒子滤波、Metropolis Hasting 采样工夫序列剖析
R 语言 Metropolis Hastings 采样和贝叶斯泊松回归 Poisson 模型
R 语言贝叶斯 MCMC：用 rstan 建设线性回归模型剖析汽车数据和可视化诊断
R 语言贝叶斯 MCMC：GLM 逻辑回归、Rstan 线性回归、Metropolis Hastings 与 Gibbs 采样算法实例
R 语言贝叶斯 Poisson 泊松 - 正态分布模型剖析职业足球比赛进球数
R 语言用 Rcpp 减速 Metropolis-Hastings 抽样预计贝叶斯逻辑回归模型的参数
R 语言逻辑回归、Naive Bayes 贝叶斯、决策树、随机森林算法预测心脏病
R 语言中贝叶斯网络（BN）、动静贝叶斯网络、线性模型剖析错颌畸形数据
R 语言中的 block Gibbs 吉布斯采样贝叶斯多元线性回归
Python 贝叶斯回归剖析住房累赘能力数据集
R 语言实现贝叶斯分位数回归、lasso 和自适应 lasso 贝叶斯分位数回归剖析
Python 用 PyMC3 实现贝叶斯线性回归模型
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真剖析
R 语言和 STAN,JAGS：用 RSTAN,RJAG 建设贝叶斯多元线性回归预测选举数据
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言贝叶斯线性回归和多元线性回归构建工资预测模型
R 语言贝叶斯推断与 MCMC：实现 Metropolis-Hastings 采样算法示例
R 语言 stan 进行基于贝叶斯推断的回归模型
R 语言中 RStan 贝叶斯层次模型剖析示例
R 语言应用 Metropolis-Hastings 采样算法自适应贝叶斯预计与可视化
R 语言随机搜寻变量抉择 SSVS 预计贝叶斯向量自回归（BVAR）模型
WinBUGS 对多元随机稳定率模型：贝叶斯预计与模型比拟
R 语言实现 MCMC 中的 Metropolis–Hastings 算法与吉布斯采样
R 语言贝叶斯推断与 MCMC：实现 Metropolis-Hastings 采样算法示例
R 语言应用 Metropolis-Hastings 采样算法自适应贝叶斯预计与可视化
视频：R 语言中的 Stan 概率编程 MCMC 采样的贝叶斯模型
R 语言 MCMC:Metropolis-Hastings 采样用于回归的贝叶斯预计