R语言和Stan,JAGS：用rstan,rjag建立贝叶斯多元线性回归预测数据

By tecdat4月 7, 2021大数据部落, 数理统计, 机器学习, 计算机科学与技术JAGS, rjag, Rstan, stan, 回归, 多元, 线性回归, 贝叶斯, 预测

本文将介绍如何在R中做贝叶斯回归分析，R中有不少包可以用来做贝叶斯回归分析，比如最早的（同时也是参考文献和例子最多的）R2WinBUGS包。

由Kaizong Ye，Sherry Deng撰写

这个包会调用WinBUGS软件来拟合模型，后来的JAGS软件也使用与之类似的算法来做贝叶斯分析。然而JAGS的自由度更大，扩展性也更好。

近来，STAN和它对应的R包rstan一起进入了人们的视线。STAN使用的算法与WinBUGS和JAGS不同，它改用了一种更强大的算法使它能完成WinBUGS无法胜任的任务。同时Stan在计算上也更为快捷，能节约时间。

可下载资源

完整程序、数据和文档（word）

例子

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

JAGS，全称是Just another Gibbs sampler，是基于BUGS语言开发的利用MCMC来进行贝叶斯建模的软件包。它没有提供建模所用的GUI以及MCMC抽样的后处理，这些要在其它的程序软件上来处理，比如说利用R包（rjags）来调用JAGS并后处理MCMC的输出。JAGS相对于WinBUGS/OpenBUGS的主要优点在于平台的独立性，可以应用于各种操作系统，而WinBUGS/OpenBUGS只能应用于windows系统；JAGS也可以在64-bit平台上以64-bit应用来进行编译。

JAGS和R的交互非常好，下面我们使用rjags包来实现R对JAGS的调用。运行一个JAGS模型是指在参数的后验分布中生成抽样，需要这样5个步骤：

定义模型
初始化
编译
适应
监测
后续的MCMC收敛诊断、模型评价等等工作是要由R来完成的。当然，在使用rjags之前，要保证JAGS已经安装在你的电脑上。（JAGS下载：http://sourceforge.net/projects/mcmc-jags/files/）

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

加载并标准化选举数据

1

2

3

4

5

6

7

8

9

10

# 加载数据

load("elec.RData")

Y <- Y[!is.na(Y+rowSums(X))]

X <- X[!is.na(Y+rowSums(X)),]

n <- length(Y)

p <- ncol(X)

视频

R语言中RStan贝叶斯层次模型分析示例

1	## [1] 3111

1

p

1

## [1] 15

1

2

3

4

5

6

X <- scale(X)

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

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

table(test)

1

2

3

## test

## FALSE TRUE

## 100 3011

1

2

3

4

5

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

Xo <- X[!test,]

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

Xp <- X[test,]

选举数据的探索性分析

1	boxplot(X, las = 3

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

最受欢迎的见解

1.使用R语言进行METROPLIS-IN-GIBBS采样和MCMC运行

2.R语言中的Stan概率编程MCMC采样的贝叶斯模型

3.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

4.R语言BUGS JAGS贝叶斯分析马尔科夫链蒙特卡洛方法（MCMC）采样

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

6.R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

7.R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数

8.R语言使用Metropolis- Hasting抽样算法进行逻辑回归

9.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

rstan中实现

统一先验分布

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

R语言中的Stan概率编程MCMC采样的贝叶斯模型

阅读文章 ➜

1

2

3

4

5

6

7

8

9

10

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

1

2

3

rstan_options(auto_write = TRUE)

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

1	print(fit)

1	hist(fit, pars = pars)

随时关注您喜欢的主题

1

dens(fit)

1	traceplot(fit)

rjags中实现

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

1

2

3

4

5

6

7

8

9

10

11

12

13

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中编译模型

1

2

3

# 注意：Yp不发送给JAGS

jags.model(model,

data = list(Yo=Yo,no=no,np=np,p=p,Xo=Xo,Xp=Xp))

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

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

1

2

3

4

5

6

7

8

9

10

11

#提取每个参数的样本

samps <- samp[[1]]

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

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

beta.mn <- colMeans(beta.samps)

1

2

3

4

5

6

7

8

9

10

11

12

# 绘制后验预测分布和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)

1

2

3

# 95% 置信区间

alpha.mn+Xp%*%beta.mn - 1.96*sigma.mn

alpha.mn+Xp%*%beta.mn + 1.96*sigma.mn

1	## [1] 0.9452009

1

2

3

# PPD 95% 置信区间

apply(Yp.samps,2,quantile,0.025)

apply(Yp.samps,2,quantile,0.975)

1	## [1] 0.9634673

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

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

非常感谢您阅读本文，如需帮助请联系我们！

相关文章