PSYCH统计实验室
前言
Stan 语言是一种概率程序语言，一个完整的 Stan 程序一般由若干个模块构成，包括
数据模块（data）：定义构建贝叶斯模型的观察变量类型、取值范围等。
参数模块（parameters）：定义贝叶斯模型参数的类型、取值范围。
参数转换模块（transformed parameters）（可选）：通过表达式与赋值语句将参数模块的参数转换为新的参数变量。
模型模块（model）：定义模型的核心部分，包括参数的先验分布和数据的似然函数。
预测值模块（generated quantities）（可选）：块通过变量定义与赋值语句定义预测变量。
与 BUGS 语言不同，Stan 程序是基于 C++语言的。因此，各变量使用前必须定义，各模块需按照顺序排列。
接下来我们会用一个简单的示例，来教大家如何进行模型构建以及相应的stan代码的编写。
数据介绍
假设我们有这样一批数据，分别采集了30只小鼠分别在第8、15、22、29、36天的体重，这时，我们想要知道小鼠日龄和体重的关系，此时我们以日龄为自变量、体重为因变量，去构建线性-正态层次模型。
构建模型
明确研究问题
首先，我们明确一下我们的研究问题，在这个案例中，我们的问题是“小鼠年龄和体重的关系如何”。我们的数据是30只小鼠，采集了每一只在不同的日龄（8、15、22、29、36 天）下的体重。这是一种纵向数据，我们对同一个体（小鼠）的体重随时间的变化进行了多次测量。此时，我们关注的核心问题一共是两个：
1.小鼠的平均体重随年龄增长的速度（变化率）
2.不同小鼠的个体差异（即初始体重和体重增长率的差异）
选择模型框架：分层线性模型
由于我们想建模的是每只小鼠的体重随时间的变化，我们可以使用线性模型来描述体重（因变量）和日龄（自变量）之间的关系。这里适合用 分层线性模型，因为：
我们有多只小鼠，每只小鼠都可以看作一个层级。
每只小鼠的初始体重（截距）和体重增长率（斜率）可能不一样。
构建似然函数
为了描述每只⼩⿏在每个⽇龄的体重，我们假设：
•体重 yij 是正态分布的，因为体重数据通常符合正态分布。
•对每只⼩⿏ i，体重与年龄的关系可以⽤线性关系表⽰： αi + βi(xij − xˉ)
yij ：第 i 只⼩⿏在第 j 天的体重。
αi ：第 i 只⼩⿏的截距（初始体重）。
βi ：第 i 只⼩⿏的斜率（体重随⽇龄的变化率）。
xij − xˉ：⽇龄的中⼼化处理（⽤⽇龄减去平均⽇龄），这样可以减少模型的计算复杂度并让估计更稳定。
于是，似然函数可以写成： αi ∼ N(μα,σα),
βi∼N(μβ,σβ)
• σy 表⽰测量体重时的观测误差的标准差。
构建层次先验（分层结构）
由于我们要考虑每只⼩⿏的个体差异，我们假设每只⼩⿏的截距αi和斜率βi都来⾃各⾃的总体分布，即：
•每只⼩⿏的初始体重 αi 是从⼀个总体均值为 μα ，标准差为 σα 的正态分布中抽取的。
•每只⼩⿏的体重变化率 βi 是从⼀个总体均值为 μβ ，标准差为 σβ 的正态分布中抽取的。
因此，层次先验可以写成： αi ∼ N(μα , σα ),
βi ∼ N(μβ , σβ)
给定⾼层次先验
在⻉叶斯模型中，每个参数都需要⼀个先验分布。在这个模型⾥，除了 αi 、βi 有⾃⼰的层次先验，我们还需要给：μα 、μβ（总体均值） 、 σα 、σβ （总体⽅差）以及观测误差σy分配先验。
常⻅的选择是：
•对均值参数 μα 、 μβ 使⽤正态分布先验： μα ∼ N(0, 100),μβ ∼ N(0, 100)
•对⽅差参数 σy 、 σα 、 σβ 使⽤逆伽⻢分布或半正态分布作为先验：σy, σα, σβ ∼ IG(0.001, 0.001)
总结完整模型
将所有这些步骤组合起来，我们得到了这个分层线性模型的完整形式：
根据上⾯的模型结构，我们可以逐步把它转换成 Stan 代码。我们会在 Stan 中定义数据块、参数块、模型块等。
Step 1: 数据块
在 Stan 的 data 块中，我们定义需要输⼊的数据，包括样本数量、⼩⿏数量、⽇龄、体重数据等。
Step 2: 参数块
在 parameters 块中，我们定义需要估计的参数，包括每只⼩⿏的截距和斜率、总体均值和⽅差等。
Step 3: 参数转换模块
transformed parameters ⽤于定义由原始参数转换⽽来的新参数。它在 Stan 的计算流程中会在 parameters 块之后执⾏，也就是说它可以使⽤ parameters 中的参数值。在这⾥，我们将之前 parameters 块中定义的⽅差参数（ sigmasq_y 、 sigmasq_alpha 和sigmasq_beta ）转换为标准差，以便在模型的其他部分更直观地使⽤这些参数。
Step 4: 模型块
在 model 块中，我们定义似然函数以及参数的先验分布。
Step 5: 预测值模块
generated quantities 块⽤于定义从模型参数计算得到的衍⽣量。这些衍⽣量不会参与模型的采样过程，但会基于采样结果计算并存储下来，便于后续分析。在模型中， alpha0 是⼀个计算得到的衍⽣参数，⽤来表⽰基准体重。
数据以 R 列表( list) 形式存储,其中 x 为⼩⿏的⽇龄,xbar 为⽇龄均值,N 为⼩⿏数⽬,T 为观察次数,y 为⼩⿏体重。(前面的数字是行数)
后验分布抽样
通过 library 语句载⼊ rStan 软件包，并通过 Stan 函数执⾏ HMC 抽样，结果存储在fit_rats，该变量存储了模型构建以及后验样本等信息，后续的收敛诊断以及统计推断的数据均来⾃于该变量。该变量存储初始值、Stan 函数的参数设置、Stan 模型、以及后验分布抽样样本等信息。
收敛诊断
通过 traceplot 函数可以对估计参数绘制踪迹图, 以对后验样本是否收敛进⾏诊断,图 3 显⽰了 4 个参数的踪迹图,并⽤灰⾊表⽰预热状态的抽样。
由图可⻅,在抽样 250 次后,4 条⻢尔科夫链很好的重合在⼀起,并保持平稳、⽔平。说明模型已经很好的收敛,更多的收敛诊断⽅法可以使⽤ Stan 提供的 bayesplot 软件包进⾏收敛诊断。
后验分布推断
使⽤ print 函数展⽰后验样本的统计推断结果，可通过指定 pars 参数指定需要展⽰的变量参数，本例展⽰ α0、μβ、μα、σy 4 个参数的后验分布推断结果。
结果分别显⽰了参数估计值、标准误、标准差、95% 可信区间、中位数以及等效样本量和收敛诊断指标 R^ 估计值。结果表明，⼤⿏基线体重为105.13g，在实验期间每天体重平均增加6.2g。对于每⼀个变量，n_ eff 是有效样本量的估计值，各参数 R^ 均为 1，说明⻢尔科夫链已经收敛于⽬标分布。
PSYCH统计实验室
参考文献
刘晋, 汪秀琴, 李天萍, 徐⽂华, & 陈峰. (2019). ⻉叶斯统计分析的新⼯具—Stan. 中国卫⽣统计, 36(3),462-465.
PSYCH统计实验室
通知公告
网络分析课程目前开放视频课啦
单次课200元/讲（学生），250元/讲（非学生）
共有四讲内容：
①横断面网络分析简介与基础
②网络分析与因子分析
③交叉滞后网络分析
④时间序列网络分析
购买后开放视频权限14天，可多次申请。
并赠送所有课程相关资料（无PPT）
如果想申请购买，请联系M18812507626
更多资讯
关注我们
文稿：Travel
排版：Little Star
责编：Wink
审核：摘星
本文由“Psych统计自习室”课题组原创，欢迎转发至朋友圈。如需转载请联系后台，征得作者同意后方可转载。
