R语言和STAN,JAGS:用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

简介: R语言和STAN,JAGS:用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

全文链接:http://tecdat.cn/?p=21978


本文将介绍如何在R中用rstan和rjags做贝叶斯回归分析,R中有不少包可以用来做贝叶斯回归分析,比如最早的(同时也是参考文献和例子最多的)R2WinBUGS包。这个包会调用WinBUGS软件来拟合模型,后来的JAGS软件也使用与之类似的算法来做贝叶斯分析。然而JAGS的自由度更大,扩展性也更好。近来,STAN和它对应的R包rstan一起进入了人们的视线。STAN使用的算法与WinBUGS和JAGS不同,它改用了一种更强大的算法使它能完成WinBUGS无法胜任的任务。同时Stan在计算上也更为快捷,能节约时间。


例子


设Yi为地区i=1,…,ni=1,…,n从2012年到2016年支持率增加的百分比。我们的模型

式中,Xji是地区i的第j个协变量。所有变量均中心化并标准化。我们选择σ2∼InvGamma(0.01,0.01)和α∼Normal(0100)作为误差方差和截距先验分布,并比较不同先验的回归系数。


加载并标准化选举数据

# 加载数据
 load("elec.RData")
 Y    <- Y[!is.na(Y+rowSums(X))]
 X    <- X[!is.na(Y+rowSums(X)),]
 n    <- length(Y)
 p    <- ncol(X)


## [1] 3111


p


## [1] 15


X    <- scale(X)
# 将模型拟合到大小为100的训练集,并对剩余的观测值进行预测
 test  <- order(runif(n))>100
 table(test)


## test
## FALSE  TRUE 
##   100  3011


Yo    <- Y[!test]    # 观测数据
 Xo    <- X[!test,]
 Yp    <- Y[test]     # 为预测预留的地区
 Xp    <- X[test,]


选举数据的探索性分析


boxplot(X, las = 3


image(1:p, 1:p, main = "预测因子之间的相关性")


点击标题查阅往期内容


R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归


01

02

03

04

rstan中实现


统一先验分布


如果模型没有明确指定先验分布,默认情况下,Stan将在参数的合适范围内发出一个统一的先验分布。注意这个先验可能是不合适的,但是只要数据创建了一个合适的后验值就可以了。

data {
  int<lower=0> n; // 数据项数
  int<lower=0> k; // 预测变量数
  matrix[n,k] X; // 预测变量矩阵
  vector[n] Y; // 结果向量
}
parameters {
  real alpha; // 截距
  vector[k] beta; // 预测变量系数
  real<lower=0> sigma; // 误差


rstan_options(auto_write = TRUE)
#fit <- stan(file = 'mlr.stan', data = dat)


print(fit)


hist(fit, pars = pars)


dens(fit)


traceplot(fit)


rjags中实现


用高斯先验拟合线性回归模型

library(rjags)
model{
#  预测
  for(i in 1:np){
    Yp[i]  ~ dnorm(mup[i],inv.var)
    mup[i] <- alpha + inprod(Xp[i,],beta[])
  # 先验概率
  alpha     ~ dnorm(0, 0.01)
  inv.var   ~ dgamma(0.01, 0.01)
  sigma     <- 1/sqrt(inv.var)


在JAGS中编译模型

# 注意:Yp不发送给JAGS
jags.model(model, 
                    data = list(Yo=Yo,no=no,np=np,p=p,Xo=Xo,Xp=Xp))


coda.samples(model, 
        variable.names=c("beta","sigma","Yp","alpha"),


从后验预测分布(PPD)和JAGS预测分布绘制样本

#提取每个参数的样本
 samps       <- samp[[1]]
 Yp.samps    <- samps[,1:np] 
#计算JAGS预测的后验平均值
 beta.mn  <- colMeans(beta.samps)
# 绘制后验预测分布和JAGS预测
 for(j in 1:5)
    # JAGS预测
    y  <- rnorm(20000,mu,sigma.mn)
    plot(density(y),col=2,xlab="Y",main="PPD")
    # 后验预测分布
    lines(density(Yp.samps[,j]))
    # 真值
    abline(v=Yp[j],col=3,lwd=2)


# 95% 置信区间
 alpha.mn+Xp%*%beta.mn - 1.96*sigma.mn
 alpha.mn+Xp%*%beta.mn + 1.96*sigma.mn


## [1] 0.9452009


# PPD 95% 置信区间
 apply(Yp.samps,2,quantile,0.025)
 apply(Yp.samps,2,quantile,0.975)


## [1] 0.9634673


请注意,PPD密度比JAGS预测密度略宽。这是考虑β和σ中不确定性的影响,它解释了JAGS预测的covarage略低的原因。但是,对于这些数据,JAGS预测的覆盖率仍然可以。

相关文章
|
2月前
|
数据采集 机器学习/深度学习 数据可视化
R语言从数据到决策:R语言在商业分析中的实践
【9月更文挑战第1天】R语言在商业分析中的应用广泛而深入,从数据收集、预处理、分析到预测模型构建和决策支持,R语言都提供了强大的工具和功能。通过学习和掌握R语言在商业分析中的实践应用,我们可以更好地利用数据驱动企业决策,提升企业的竞争力和盈利能力。未来,随着大数据和人工智能技术的不断发展,R语言在商业分析领域的应用将更加广泛和深入,为企业带来更多的机遇和挑战。
|
3月前
|
存储 数据采集 数据处理
R语言数据变换:使用tidyr包进行高效数据整形的探索
【8月更文挑战第29天】`tidyr`包为R语言的数据整形提供了强大的工具。通过`pivot_longer()`、`pivot_wider()`、`separate()`和`unite()`等函数,我们可以轻松地将数据从一种格式转换为另一种格式,以满足不同的分析需求。掌握这些函数的使用,将大大提高我们处理和分析数据的效率。
|
2月前
R语言基于表格文件的数据绘制具有多个系列的柱状图与直方图
【9月更文挑战第9天】在R语言中,利用`ggplot2`包可绘制多系列柱状图与直方图。首先读取数据文件`data.csv`,加载`ggplot2`包后,使用`ggplot`函数指定轴与填充颜色,并通过`geom_bar`或`geom_histogram`绘图。参数如`stat`, `position`, `alpha`等可根据需要调整,实现不同系列的图表展示。
|
2月前
|
数据采集 数据可视化 数据挖掘
R语言在金融数据分析中的深度应用:探索数据背后的市场智慧
【9月更文挑战第1天】R语言在金融数据分析中展现出了强大的功能和广泛的应用前景。通过丰富的数据处理函数、强大的统计分析功能和优秀的可视化效果,R语言能够帮助金融机构深入挖掘数据价值,洞察市场动态。未来,随着金融数据的不断积累和技术的不断进步,R语言在金融数据分析中的应用将更加广泛和深入。
|
3月前
|
数据采集 机器学习/深度学习 数据挖掘
R语言数据清洗:高效处理缺失值与重复数据的策略
【8月更文挑战第29天】处理缺失值和重复数据是数据清洗中的基础而重要的步骤。在R语言中,我们拥有多种工具和方法来有效地应对这些问题。通过识别、删除或插补缺失值,以及删除重复数据,我们可以提高数据集的质量和可靠性,为后续的数据分析和建模工作打下坚实的基础。 需要注意的是,处理缺失值和重复数据时,我们应根据实际情况和数据特性选择合适的方法,并在处理过程中保持谨慎,以避免引入新的偏差或错误。
|
3月前
|
数据采集 存储 数据可视化
R语言时间序列分析:处理与建模时间序列数据的深度探索
【8月更文挑战第31天】R语言作为一款功能强大的数据分析工具,为处理时间序列数据提供了丰富的函数和包。从数据读取、预处理、建模到可视化,R语言都提供了灵活且强大的解决方案。然而,时间序列数据的处理和分析是一个复杂的过程,需要结合具体的应用场景和需求来选择合适的方法和模型。希望本文能为读者在R语言中进行时间序列分析提供一些有益的参考和启示。
|
3月前
|
资源调度 数据挖掘
R语言回归分析:线性回归模型的构建与评估
【8月更文挑战第31天】线性回归模型是统计分析中一种重要且实用的工具,能够帮助我们理解和预测自变量与因变量之间的线性关系。在R语言中,我们可以轻松地构建和评估线性回归模型,从而对数据背后的关系进行深入的探索和分析。
|
3月前
|
数据处理
R语言数据合并:掌握`merge`与`dplyr`中`join`的巧妙技巧
【8月更文挑战第29天】如果你已经在使用`dplyr`进行数据处理,那么推荐使用`dplyr::join`进行数据合并,因为它与`dplyr`的其他函数(如`filter()`、`select()`、`mutate()`等)无缝集成,能够提供更加流畅和一致的数据处理体验。如果你的代码中尚未使用`dplyr`,但想要尝试,那么`dplyr::join`将是一个很好的起点。
|
6月前
|
数据采集 数据可视化
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
|
6月前
|
Web App开发 数据可视化 数据挖掘
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)