关于数据挖掘:R语言STAN贝叶斯线性回归模型分析气候变化影响北半球海冰范围和可视化检查模型收敛性附代码数据

原文链接：http://tecdat.cn/?p=24334

最近咱们被客户要求撰写对于贝叶斯线性回归的钻研报告，包含一些图形和统计输入。

像任何统计建模一样，贝叶斯建模可能须要为你的钻研问题设计适合的模型，而后开发该模型，使其合乎你的数据假如并运行

1. 理解 Stan

统计模型能够在 R 或其余统计语言的各种包中进行拟合。但有时你在概念上能够设计的完满模型，在限度了你能够应用的散布和复杂性的软件包或程序中很难或不可能实现。这时你可能想转而应用统计编程语言，如 Stan。

Stan 是一种旧式的语言，它提供了一种更全面的学习和实现贝叶斯模型的办法，能够适应简单的数据结构。Stan 开发团队的一个指标是通过清晰的语法、更好的采样器（这里的采样是指从贝叶斯后验散布中抽取样本）以及与许多平台（包含 R、RStudio、ggplot2 和 Shiny）的集成，使贝叶斯建模更易于应用。

在这个入门教程中，咱们将从一个线性模型开始，经验模型建设的迭代过程。在咱们的高级 stan 教程中，咱们将摸索更简单的模型构造。

首先，在建设模型之前，你须要定义你的问题并理解你的数据。摸索它们，绘制它们，计算一些汇总统计。

一旦你对你的数据和你想用统计模型答复的问题有了理解，你就能够开始建设贝叶斯模型的迭代过程。

• 设计你的模型。
• 抉择先验
• 对后验散布进行采样。
• 查看模型收敛（traceplots、rhats）
• 应用后验预测批判性地评估模型并查看它们与您的数据的比拟状况
• 反复…

模仿数据也是很好的做法，以确保你的模型正确，作为测试你的模型的另一种形式。

2. 数据

首先，让咱们找到一个能够拟合简略线性模型的数据集。气候变化对地球最显着的影响之一是北半球每年海冰范畴的缩小。让咱们应用 Stan 的线性模型摸索海冰范畴如何随工夫变动。

通过运行 setwd("your-file-path") 蕴含您本人的文件门路的代码，将您的工作目录设置为您保留数据的文件夹。当初，让咱们加载数据：

# 增加 stringsAsFactors = F 意味着数字变量将不会被
# 作为因子 / 分类变量读入
ece <- red.cv("sv", stinsAsFators = F)

咱们来看一下数据：

咱们能够用这些数据提出什么钻研问题？以下状况如何：

钻研问题： 北半球的海冰范畴是否会随着工夫的推移而缩小？

为了摸索这个问题的答案，首先咱们能够做一个数字。

plot(th ~ yr, data)

图 1. 北半球海冰范畴随工夫的变动。

当初，让咱们应用 lm().

l1 <- lm(exnoh ~ yer, data = sie)
summary(l1)

咱们能够将该模型增加到咱们的绘图中：

ablne(m1, l = 2, ty = 2, w = 3)

图 2. 北半球海冰范畴随工夫的变动（加上线性模型拟合）。

记住线性模型的方程：

y = α + β ∗ x + 误差

在 Stan 你须要指定你想模型。

兴许咱们曾经找到了问题的答案，但本教程的重点是摸索应用编程语言 Stan，所以当初让咱们尝试在 Stan 中编写雷同的模型。

筹备数据

让咱们重命名变量并将年份从 1 索引到 39。对于贝叶斯模型的一个要害是您必须应用信息散布来形容数据中的变动。因而，您心愿确保您的数据合乎这些散布，并且它们将实用于您的模型。在这种状况下，咱们真的想晓得从数据集的开始到数据集完结的海冰是否产生了变动，而不是 1979 年到 2017 年。咱们不须要咱们的模型预计 500 年或 600 年的海冰是什么样的，就在咱们的数据集的持续时间内。因而，咱们将年份数据设置为索引 1 到 30 年。

x <- I(year - 1978)

咱们能够应用新数据从新运行该线性模型。

summary(lm1)

咱们还能够从咱们的简略模型中提取一些要害的汇总统计数据，以便咱们 Stan 稍后能够将它们与模型的输入进行比拟。

coeff[1] # 截距值
coeff[2] # 斜率
sigma(lm1) # 残差 

当初让咱们将其转换为用于输出 Stan 模型的数据框。传递给 Stan 的数据须要是命名对象列表。此处给出的名称须要与模型中应用的变量名称相匹配。

库
–

请确保装置了以下库（这些是本 Stan 教程和下一个教程的库）。rstan 是最重要的，如果您没有 C++ 编译器，则须要一些额定的货色。

3. Stan 程序

咱们将首先用语言编写一个线性模型 Stan。这能够写在你的 R 脚本中，或者独自保留为一个 .stan 文件并调用到 R.

一个 Stan 程序具备三个必须的“块”：

• “数据”块：您能够在其中申明数据类型、它们的维度、任何限度（即 upper = 或 lower =，用作查看 Stan）及其名称。
• “参数”块：您能够在此处指明要建模的参数和名称。对于线性回归，咱们心愿对回归线四周的误差的截距、任何斜率和标准偏差进行建模。
• “模型”块：这是蕴含任何抽样语句的中央，包含正在应用的模型。模型块是指明要为参数蕴含的任何先验散布的中央。如果未定义 Stan 先验，则 应用默认先验 uniform(-infinity, +infinity)。您能够在申明参数时应用下限或上限来限度先验（即 lower = 0\> 确保参数为正）。

采样由 ~ 符号示意，并且 Stan 曾经蕴含许多常见的散布作为矢量化函数。

还有四个可选块：

• “性能”
• “ 转化的数据 ”
• “ 转换后的参数
• “ 生成的数量 ”

正文 // 在 Stan 中用 示意。该 write("model code", "file_name") 容许咱们在 R 脚本中编写 Stan 模型并将文件输入到工作目录（或者您能够设置不同的文件门路）。

write("// 简略线性回归的模型

数据

 int < lower = 1 > N; // 样本大小
vector[N] x;// 预测
vecor[N] y;// 后果


参数 

 real alha; // 截距
 real beta; // 斜率 (回归系数)
 real < loer = 0 > sima; // 误差 SD



模型 

 y ~ nrmal(alpha + x * beta , siga);



产生的数量 


// 后验预测散布 "。"md1.stan"

首先，咱们应该查看咱们的 Stan 模型以确保咱们编写了一个文件。

当初让咱们保留该文件门路。

 "model1.stan"

在这里，咱们隐式地将 uniform(-infinity, +infinity) 先验用于咱们的参数。

4. 运行 Stan 模型

Stan 程序 C++ 在被应用之前被恪守。这意味着在 R 能够应用模型之前须要运行 C++ 代码。为此，您必须 C++ 装置编译器。编译后，您能够在每个会话中屡次应用模型，但在开始新 R 会话时必须从新编译。有许多 C++ 编译器，而且它们在不同零碎中通常是不同的。如果您的模型一堆谬误，请不要放心。只有模型能够与该 stan() 函数一起应用，它就能够正确编译。如果咱们想应用以前编写的 .stan 文件，咱们在 file 函数中应用 参数 stan_model()。

咱们通过应用 stan() 函数拟合咱们的模型，并为它提供模型、数据，并批示预热的迭代次数（这些迭代稍后不会用于后验散布，因为它们只是模型“预热””），总迭代次数，咱们要运行的链数，咱们要应用的内核数（Stan 为并行化而设置），它示意同时运行的链数（即，如果您的计算机有四个内核），您能够在每个链上运行一个链，同时创立四个链）和细化，这是咱们想要存储咱们的预热后迭代的频率。“thin = 1”将保留每次迭代，“thin = 2”将保留每一秒，依此类推……

Stan 如果 warmup = 未指定参数，则主动应用一半的迭代作为预热。

stan(fie =moel1, data = data, wrmup = 500, ier = 1000, chais = 4, cres = 2, thn = 1)

拜访 fit 对象的内容 

后果 stan() 保留为 stanfit 对象（S4 类）。

咱们能够通过执行对象的名称来获取参数估计和采样器诊断的汇总统计信息：

fit

模型输入展现了什么？你怎么晓得你的模型曾经收敛了？您能看到批示您的 C++ 编译器已运行的文本吗？

从这个输入中，咱们能够通过查看 Rhat 每个参数的值来疾速评估模型收敛性。当这些值等于或靠近 1 时，链曾经收敛。还有许多其余诊断办法，但这对 Stan 来说很重要。

咱们还能够通过从模型对象中提取参数来查看参数的残缺后验。有很多办法能够查看后验。

poteir <- exrat(fit)

extract() 将每个参数的后验预计放入一个列表中。

让咱们与咱们之前应用“lm”的预计进行比拟：

plot(y ~ x)

图 3. 北半球海冰范畴随工夫的变动（比拟 Stan 线性模型拟合和个别 lm 拟合）。

后果与 lm 输入雷同。这是因为咱们应用了一个简略的模型，并且在咱们的参数上搁置了非信息先验。

将回归线预计中的可变性可视化的一种办法是绘制来自后验的多个预计。

plot(y ~ x, pch = 20)

图 4. 北半球海冰范畴随工夫的变动（Stan 线性模型拟合）。

5. 扭转先验

咱们再试一次，但当初对海冰和工夫之间的关系采纳更有信息量的预设。咱们将应用小标准差的正态先验。如果咱们应用标准差十分大的正态先验（比方 1000，或者 10000），它们的作用与平均先验十分类似。

 参数 

 real alha; // 截距
 real bta; // 斜率 (回归系数)
 real < lowr = 0 > sima; // 误差 SD



模型

 beta ~ nomal(1, 0.1);
 y ~ normal(apha + x * beta , siga);

咱们将拟合该模型并将其与应用平均先验的均值预计进行比拟。

 stan(stn_oel)



plot(y ~ x)

图 5. 北半球海冰范畴随工夫的变动（Stan 线性模型拟合）。

后验预测产生了什么变动？模型是否更好地拟合数据？为什么模型拟合产生了变动？通过制作十分窄的先验散布，咱们的模型扭转了什么？

点击题目查阅往期内容

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

左右滑动查看更多

01

02

03

04

尝试本人将先验更改为一些不同的数字，看看会产生什么，这是贝叶斯建模中的一个常见问题，如果您的先验散布十分窄，但不合乎您对系统或数据分布的了解，您能够运行无奈有意义地解释数据变动的模型。然而，这并不是说您不应该抉择一些信息丰盛的先验，您的确心愿应用先前的剖析和对您的钻研零碎的了解来告知您的模型先验和设计。您只须要认真思考您做出的每个建模决策！

6. 收敛诊断

在持续之前，咱们应该再次查看模型参数的 Rhat 值、无效样本大小 (n_eff) 和跟踪图，以确保模型已收敛且牢靠。

n_f 是无效样本大小的粗略度量。通常只须要放心这个数字小于迭代次数的 1/100 或 1/1000。

对于跟踪图，咱们能够间接从后验查看它们：

plot(alpha, type = "l")
plot(beta, type = "l")
plot(sigma, type = "l")

图 6. alpha 的迹线图，截距。

对于更简略的模型，收敛通常不是问题，除非您的代码中有谬误，或者运行采样器的迭代次数太少。

收敛性差

尝试仅运行 50 次迭代的模型并查看跟踪图。

fid <- stan(iter = 50)

这在预热后也有一些“转换”，表明模型指定谬误，或者采样器未能齐全采样后验。

plot(alpha, type = "l")
plot(beta, type = "l")
plot(sigma, type = "l")

图 7. alpha 的谬误轨迹图，截距。

参数汇总

咱们也能够间接通过后验失去参数的汇总。让咱们还绘制非贝叶斯线性模型值，以确保咱们的模型运行正确

par(mfrow = c(1,3))

plot(dnsty(alpha)

图 8.Stan 模型拟合的密度图散布 与个别 lm 拟合的估计值比拟。

从后验咱们能够间接计算任何参数超过或低于某个感兴趣值的概率。

beta > 0 的概率：

sum(beta>0)/lenth(beta)
# 0

beta > 0.2 的概率：

sum(bea>0.2)/legth(beta)
# 0

诊断图 rstan

尽管咱们能够间接应用后验，但 rstan 内置了许多有用的性能。

plot(fit)

图 9.Stan 模型不同链的跟踪图。

咱们还能够查看后验密度和直方图。

dens(it)
hist(ft)

图 10.Stan 模型中截距、斜率和残差方差的后验密度图和直方图。

咱们能够生成批示均匀参数估计值和咱们可能感兴趣的任何置信区间的图。请留神，beta 和 sigma 参数的 95% 置信区间 十分小，因而您只能看到点。依据您本人数据的差别，当您进行本人的剖析时，您可能会看到更小或更大的置信区间。

plot(fit)

图 11.Stan 模型的参数估计。

后验预测查看

对于预测和作为模型诊断的另一种模式，Stan 能够应用随机数生成器在每次迭代中为每个数据点生成预测值。通过这种形式，咱们能够生成预测，这些预测也代表了咱们模型和数据生成过程中的不确定性。可用于获取咱们想要的对于后验的任何其余信息，或对新数据进行预测。

 参数

 real alpha; // 截距
 real beta; // 斜率 (回归系数)




 y ~ nomal(x * eta + alpa, sgma);



产生的数量 


 for (n in 1:N) {yp[n] = normlrg(x[n] * bta + apha, sima)。}

请留神，GQ（生成量）块不反对矢量化，因而咱们必须将其放入循环中。然而因为它被编译为 C++，循环实际上十分快，并且 Stan 每次迭代只评估一次 GQ 块，因而它不会为您的采样减少太多工夫。通常，数据生成函数将是您在模型块中应用的散布，但带有 _rng 后缀。

 stan(modl2GQ, data , ier = 1000, hans = 4, cres = 2, tin = 1)

y_rep 从后验中提取 值。

解决 y_rep 值有很多抉择。

每一行都是模型的一次迭代（繁多后验预计）。

咱们能够制作一些更丑陋的图。这个包是 ggplot2。

在 200 次后验抽样中，比拟 y 的密度和 y 的密度。

poy(y, yrep[1:200,])

图 12. 比拟随机后验抽取的估计值。

在这里，咱们看到数据（深蓝色）与咱们的后验预测十分吻合。

咱们还能够应用它来比拟汇总统计的估计值。

pp(y = y, yep = yrep, tat = "mean")

图 13. 比拟汇总统计的估计值。

咱们能够更改传递给 stat 函数的函数，甚至能够本人编写！

咱们能够考察每个数据点的均匀后验预测与每个数据点的察看值（默认线为 1:1）

scttrg(y = y, yrp = yrep)

图 14. 每个数据点的均匀后验预测与每个数据点的观测值。

所以当初您曾经学习了如何运行线性模型 Stan 并查看模型收敛性。

如有任何问题，请分割咱们！

点击文末 “浏览原文”

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

本文选自《R 语言 STAN 贝叶斯线性回归模型剖析气候变化影响北半球海冰范畴和可视化查看模型收敛性》。

点击题目查阅往期内容

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