关于数据挖掘: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与博士我的项目的延期无关。

向下滑动查看后果▼

*
*

相干视频

**

拓端

，赞37

筹备--导入和摸索数据

数据是一个.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/s10654-016-0149-3** \_

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/journal.pone.0068839

点击文末 “浏览原文”

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

本文选自《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采样用于回归的贝叶斯预计