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



rstan中实现
统一先验分布
如果模型没有明确指定先验分布,默认情况下,Stan将在参数的合适范围内发出一个统一的先验分布。注意这个先验可能是不合适的,但是只要数据创建了一个合适的后验值就可以了。
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预测的覆盖率仍然可以。