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

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

全文下载链接：/?p=21978

本文将介绍如何在R中用rstan和rjags做贝叶斯回归分析，R中有不少包可以用来做贝叶斯回归分析，比如最早的（同时也是参考文献和例子最多的）R2WinBUGS包（点击文末“阅读原文”获取完整代码数据）。

这个包会调用WinBUGS软件来拟合模型，后来的JAGS软件也使用与之类似的算法来做贝叶斯分析。然而JAGS的自由度更大，扩展性也更好。近来，STA和它对应的R包rstan一起进入了人们的视线。STA使用的算法与WinBUGS和JAGS不同，它改用了一种更强大的算法使它能完成WinBUGS无法胜任的任务。同时Stan在计算上也更为快捷，能节约时间。

相关视频

例子

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

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

加载并标准化选举数据代码语言：javascript代码运行次数：0运行复制

# 加载数据



 load("elec.RData")

 Y    <- Y[!(Y+rowSums(X))]
 X    <- X[!(Y+rowSums(X)),]
 n    <- length(Y)
 p    <- ncol(X)

代码语言：javascript代码运行次数：0运行复制

## [1] 111

代码语言：javascript代码运行次数：0运行复制

p

代码语言：javascript代码运行次数：0运行复制

## [1] 15

代码语言：javascript代码运行次数：0运行复制

 X    <- scale(X)

# 将模型拟合到大小为100的训练集，并对剩余的观测值进行预测

 test  <- order(runif(n))>100
 table(test) 

代码语言：javascript代码运行次数：0运行复制

## test
## FALSE  TRUE 
##   100  011 

代码语言：javascript代码运行次数：0运行复制

 Yo    <- Y[!test]    # 观测数据
 Xo    <- X[!test,]

 Yp    <- Y[test]     # 为预测预留的地区
 Xp    <- X[test,]

选举数据的探索性分析

代码语言：javascript代码运行次数：0运行复制

boxplot(X, las = 

代码语言：javascript代码运行次数：0运行复制

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

---

点击标题查阅往期内容

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

01

02

0

04

rstan中实现

统一先验分布

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

代码语言：javascript代码运行次数：0运行复制

 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; // 误差

代码语言：javascript代码运行次数：0运行复制

rstan_opti(auto_write = TRUE)

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

代码语言：javascript代码运行次数：0运行复制

print(fit)

代码语言：javascript代码运行次数：0运行复制

hist(fit, pars = pars)

代码语言：javascript代码运行次数：0运行复制

dens(fit)

代码语言：javascript代码运行次数：0运行复制

traceplot(fit)

rjags中实现

用高斯先验拟合线性回归模型代码语言：javascript代码运行次数：0运行复制

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中编译模型代码语言：javascript代码运行次数：0运行复制

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

代码语言：javascript代码运行次数：0运行复制

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

从后验预测分布(PPD)和JAGS预测分布绘制样本代码语言：javascript代码运行次数：0运行复制

#提取每个参数的样本

 samps       <- samp[[1]]
 Yp.samps    <- samps[,1:np] 

#计算JAGS预测的后验平均值

   <- colMeans(beta.samps)


# 绘制后验预测分布和JAGS预测

 for(j in 1:5)
    # JAGS预测
    y  <- rnorm(20000,mu,)
    plot(density(y),col=2,xlab="Y",main="PPD")

    # 后验预测分布
    lines(density(Yp.samps[,j]))

    # 真值
    abline(v=Yp[j],col=,lwd=2)

代码语言：javascript代码运行次数：0运行复制

 # 95% 置信区间
 +Xp%*% - 1.96*
 +Xp%*% + 1.96*

代码语言：javascript代码运行次数：0运行复制

## [1] 0.9452009

代码语言：javascript代码运行次数：0运行复制

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

代码语言：javascript代码运行次数：0运行复制

## [1] 0.96467

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

本文参与 腾讯云自媒体同步曝光计划，分享自。原始发表：2025-01-14，如有侵权请联系 cloudcommunity@tencent 删除模型数据算法线性回归变量