软件X：Walker包——时变系数贝叶斯广义线性模型

=ϵ∼i，∼×=软件X 18（2022）101016原始软件出版物具有时变系数的有效贝叶斯广义线性模型：R中的Walker包朱尼·赫尔斯克芬兰于韦斯屈莱大学数学与统计系，FI-40014 Jyväskyläar t i cl e i nf o文章历史记录：收到2020年2020年8月10日收到修订版，2022年保留字：贝叶斯推理时变回归R马尔可夫链蒙特卡罗a b st ra ctR包步行器将标准贝叶斯一般线性模型扩展到解释变量的影响可以随时间变化的情况。例如，这可以模拟干预措施的影响，如税收政策的变化，随着时间的推移逐渐增加其影响支持贝叶斯推理的马尔可夫链蒙特卡罗算法基于Stan软件提供的汉密尔顿蒙特卡罗，使用模型的状态空间表示对回归系数进行边缘化，以实现有效的低维采样。©2022作者（S）。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本1.0.3-1用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2020_212Code Ocean compute capsule none法律代码许可证GPL 3使用git的代码版本控制系统软件代码语言，工具和服务使用R，Stan，C++。编译要求，运行环境依赖R 3.4.0及以上版本，C++14，R packages bayesplot，BH，coda，dupiter，ggplot2，Hmisc、KFAS、loo、方法、RcppParallel、rlang、rstan、rstantools、StanHeaders、Rcpp、RcppArmadillo、RcppEigen如果可用，链接到开发人员文档/手册https://cran.r-project.org/web/packages/walker/walker.pdf技术支持电子邮件jouni. iki.fi1. 动机和意义假设感兴趣的时间序列y1，. . .，yT，其线性地依赖于某个其它预测器时间序列X1，. . .，XT，其中XT′（x1 t，. . .，x pt）“是时间点t处的预测变量，并将线性高斯时间序列回归模型定义为y t= Xt′βt+ βt，t = 1，. . . ，T，（1）其中，N（0，σ2）和β是p个未知回归系数的向量（第一个通常是截距项）。假设y t和某个预测变量x之间的关系在t = 1，. . . ，T，改变测量变量之间关系的混杂因素。在某些情况下，允许回归系数随时间变化可以缓解这些问题。我们对感兴趣的现象的知识也可能已经导致怀疑的时变关系（参见，例如，[1]中的第1章）。基本时间序列回归模型（1）可以扩展为允许未知回归系数βt随时间变化。这可以通过各种方式来完成，例如，通过基于核平滑[2]构建半参数模型，或者通过使用动态贝叶斯网络[3]或状态空间模型[4]来参数化。我们遵循状态空间建模方法，定义了一个具有随机游走的高斯时间序列回归模型系数作为yt=X′βt+βt，感兴趣的时间段方便的线性关系近似可能只适用于分段，如果基本的不βt+1=βt+ηt，（二）关系是非线性的，或者当存在一些不可测量的电子邮件地址：jouni. jyu.fi。其中ηtN（0，D），D是含有对角元素σ2，i1，.的p p对角矩阵。. . ，p，并定义先验分布对于第一时间点β1，作为N（μ β，σ2）。中的底部方程https://doi.org/10.1016/j.softx.2022.1010161β12352-7110/©2022作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx朱尼·赫尔斯克软件X 18（2022）1010162==∼|||联系我们|=|||=（2）定义了回归系数的随机游走过程，D定义了系数的可变程度（D 0时，模型崩溃为基本回归模型）。我们的目标是对未知回归系数β1，. . . ，βT和标准差σ（σ1，σ2，η，. . . ，σp，η）。虽然原则上我们可以使用通用马尔可夫链蒙特卡罗（MCMC）软件（如Stan[5]或BUGS[6]）来估计这些，但由于我们感兴趣的模型的性质，这些标准实现可能在计算上效率低下，并且容易出现与底层MCMC算法收敛相关的严重问题。例如，在（块）吉布斯采样方法中，我们通过从p（β σ，y）和p（σ β，y）采样来瞄准联合后验p（β，σy）。但由于不同时间点的系数β之间以及与相关标准差参数之间的强自相关性，该MCMC方案可能导致缓慢混合。此外，要采样的参数的总数随着数据点T的数量而增加。虽然由Stan提出的Hamilton MonteCarlo算法通常比Gibbs型算法更有效地探索高维后验，但我们仍然遇到类似的问题。R [7]包遍历器使用的另一种解决方案基于模型（2）可以被写为线性高斯状态空间模型（参见，例如，[8]中的第3.6.1节）。这允许我们在MCMC采样期间通过使用卡尔曼滤波器边缘化回归系数β，从而快速准确地推断边缘后验p（σy）。然后，通过模拟给定标准差边际后验的回归系数，得到相应的联合后验p（σ，βy）p（β σ，y）p（σy）。该采样可以例如通过模拟平滑算法[9]来执行。回归系数的边际化不能直接推广到广义线性模型，如泊松回归，因为边际对数似然是棘手的。然而，可以使用该指数状态空间模型的高斯近似[10]，然后可以使用重要性采样类型校正[11]对来自近似后验的所得样本进行加权，再次导致渐近精确推断。当将回归系数建模为简单的随机游走时，这些系数的后验估计值可能具有较大的短期变化，这在实践中可能不现实。一种使估计更加平滑的方法是从随机游走系数切换到积分的二阶随机游走系数，定义为βt+1=βt+νt，νt+ 1=νt+αt，n=N（0，D）。这是一个局部线性趋势模型[12]（也称为积分随机游走），限制是β水平上没有噪声。对于该模型，我们可以应用与随机游走系数情况相同的估计技术。更复杂的β模式也是可能的。例如，我们可以定义βt1N（βt，γtσ n），其中γ1，. . . ，γT是已知的单调递减的值序列，导致βt随时间逐渐收敛到常数的情况。2. 软件描述Walker的稳定版本可以在CRAN，1上获得，而当前的开发版本可以从Github安装。2MCMC采样由rstan（Stan的R接口）处理，1 https://cran.r-project.org/package=walker网站。2 https://github.com/helske/walker网站。而模型定义和结果分析都是用R语言完成的，这就导致了快速和灵活的建模。为了定义模型，walker包使用类似的Wilkinson–Rogers model formulation syntax [例如，通过键入fit<- walker（y ~ 0 + x +rw1（~z， beta =c（0，1）， sigma =c（2，1））， sigma_y =c（2，1））我们定义一个模型，该模型具有用于预测器x的时不变系数、用于z的一阶随机游走系数和时变截距项。β（包括截距）和σ的先验定义为长度为2的向量，分别定义正态分布的平均值和标准差以及Gamma分布的形状和Functionwalker基于公式和先前的定义创建模型，然后从rstan包调用采样函数。然后可以使用函数as.data.frame将所得后验样本转换为数据帧格式，以便于可视化和进一步分析。除 了 用 于 高 斯 情 况 的 walker 和 用 于 泊 松 和 负 二 项 式 的walker_glm之外，模型，该包包含用于结果可视化的附加功能（例如，plot_coef和pp_check）和样本外预测（函数predict）。此外，由于建模函数返回MCMC采样中使用的完整stanfit，因此可以使用几个Stan相关R软件包（如ShinyStan [14]）提供的许多通用诊断和图形工具来分析该对象。3. 说明性实例作为说明性示例，让我们考虑由随机游走（即时变截距）和两个预测器生成的长度为n100的观测值y。首先，我们模拟系数、预测变量和观测值：set.seed（1）n-100β 1-累积和（c（0.5，rnorm（n- 1，0，sd=0.05）β 2-累积和（c（-1，rnorm（n-1，0， sd=0.15） x1- rnorm（n，平均值= 2）x2-cos（1：n）rw<- cumsum（rnorm（n，0，0.5））信号- rw+ beta1 * x1 +beta2 * x2 y- rnorm（n，signal，0.5）然后我们可以调用函数walker。如第2节所述，模型被定义为一个简单的公式对象，除了前面的定义外，我们还可以将各种参数传递给rstan的采样方法，例如迭代次数iter和用于MCMC的链的数量（默认值分别为2000和4fit-walker<（y~0+rw1（~x1+x2，beta=c（0，10），sigma=c（2， 4）），sigma_y =c（2，4）， chains =2， seed = 1）walker的输出是walker_fit对象，其本质上是具有来自Stan的采样函数、原始观测值y和协变量矩阵xreg的stanfit的列表。 这允许我们使用stanfit对象的所有后处理函数。朱尼·赫尔斯克软件X 18（2022）1010163Fig. 1. 说 明 性示例的时变回归系数和真实数据生成值（虚线）的后验均值（实线）和95%后验区间（阴影区域）。图图1显示了Walker如何相对较好地恢复真实系数：后验区间包含真实系数，后验均值估计紧密跟随真实值。为了绘制图形，我们首先使用函数as.data.frame来提取我们的系数后验样本作为数据帧，然后使用包dapr[15]和ggplot2[16]分别总结和绘制后验均值和95%后验区间：library（dplr）library（ggplot2）sumr<- as.data.frame（fit，type =“tv”） %>%group_by（variable，time）%>%总结（平均值=平均值（值），lwr =分位数（值，0.025），upr =分位数（值，0.975））sumr$true- c（rw，beta1，beta2）ggplot（sumr，aes（y= mean，x =time，colour=variable））+geom_ribbon（aes（ymin= lwr，ymax=upr，fill=variable），颜色= NA，alpha = 0.2）+ geom_line（aes（线型=“估计值”），lwd = 1）+geom_line（aes（y= true，linetype=“True”），lwd= 1）+scale_linetype_manual（values= c（“solid”，“dashed”））+theme_bw（）+xlab（“Time”）+ylab（“Value”）+theme（legend.position=“bottom”，legend.title= element_blank（））更多的例子可以在软件包vignette和函数文档页面中找到（例如，输入vignette（''walker''）或？walker_glm），包括walker的边缘化方法和Stan的“naive”实现的比较，以及一个可扩展性的例子。4. 影响和结论Walker软件包以计算高效的方式将标准贝叶斯广义线性模型扩展到灵活的时变系数情况，这使得经济学、社会科学和其他领域的研究人员能够放松对响应变量和（某些）预测变量之间稳定、时不变关系的有时不合理的假设。类似的方法以前已经用于最大似然设置，例如在研究广告[17]和对国际储备的需求[18]。有几种方法可以在未来扩展沃克。已经有一些计划用于附加形式的时变系数（例如平稳自回归过程），支持更多先验和响应变量的附加分布负二项式）。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作确认这项工作得到了芬兰科学院研究资助284513，312605，311877和331817的支持。引用[1]莫里森·M回归和状态空间模型中随机游走系数的检验。Springer ScienceBusiness Media; 2012.[2]罗宾逊PM。时变参数的非参数估计。经济结构变化的统计分析和预测。Springer;1989，p. 253比64[3]杨伟华，王晓平，王晓平.动态网络预测模型。在：第八届人工智能不确定性会议论文集。美国加利福尼亚州旧金山：摩根考夫曼出版公司; 1992年，p. 41比8[4]哈维AC，菲利普斯GDA。时变参数回归模型的估计。In：Deistler M，FürstE，Schwödiauer G，editors.游戏、经济动态和时间序列分析。Heidelberg：Physica-VerlagHD; 1982，p. 306-21[5] Stan开发团队Stan C++库2016年，版本2.15.0。网址http://mc-stan.org/。[6]Lunn DJ，Thomas A，Best N，Spiegelhalter D. WinBUGS -贝叶斯建模框架：概念、结构和可扩展性。2000年，http：//dx.doi.org/10.1023/A:1008929526011网站。[7]R核心团队R：统计计算的语言和环境。奥地利维也纳：R Foundation forStatistical Computing; 2020，URLhttps：//www.R-project.org/网站。[8]Durbin J，Koopman SJ.时间序列分析的状态空间方法。第2版纽约：牛津大学出版社，2012.[9]Durbin J，Koopman SJ.一种用于状态空间时间序列分析的简单有效的仿真平滑器。Biometrika2002;89：603-15.[10]Durbin J，Koopman SJ.非高斯状态空间模型的蒙特卡罗极大似然估计。Biometrika1997;84（3）：669-84.[11] Vihola M，Helske J，Franks J.基于近似边际马尔可夫链蒙特卡罗的重要性抽 样 类 型 估 计 。 斯 堪 的 纳 维 亚 河 Stat. 2020;47 （ 4 ） ： 1339-76.http://dx.doi.org/10.1111/sjos.12492网站。[12]哈维AC。预测，结构时间序列模型和卡尔曼滤波器。北京：清华大学出版社.[13]Wilkinson GN，Rogers CE.方差分析析因模型的符号描述。J R Stat Soc Ser CAppl Stat 1973;22（3）：392http://www.jstor.org/stable/2346786[14]Stan开发团队shinystan：贝叶斯模型的交互式视觉和数值诊断与后验分析2017，R包版本2.3.0. 网址http://mc-stan.org/。朱尼·赫尔斯克软件X 18（2022）1010164[15]Wickham H，François R，Henry L，Müller K. dmgr：一种数据操作的语法。2018年，R包版本0.7.6。网址https://CRAN.R-project。org/package= documentr.[16]威克姆·H Ggplot2：用于数据分析的优雅图形。Springer-Verlag New York;2016，URLhttps://ggplot2.tidyverse.org。[17]文卡特斯瓦兰·金努坎通用广告;和结构异质性假说。Can J Agric Econ Rev CanAgroeconomie 1994;42 （ 3 ） ： 381-96. http://dx.doi.org/10.1111/j.1744-7976.1994.tb00032.x网站。[18]Bahmani-Oskooee M，Brown F.卡尔曼滤波法估计国际储备需求。应用经济2004;36（15）：1655网址：//dx.doi.org/10.1080/0003684042000218543网站。