纵向监测连续非随机缺失数据变系数模型及其应用*

山西医科大学公共卫生学院卫生统计教研室(030001) 徐丽红 刘志永 刘桂芬 罗天娥

纵向研究是指同一观察对象在不同时间点连续监测的一种方法。通常情况下，观测时间多是相对固定的(如隔月检查、每周测量并记录等);但在医学监测过程中，常会由于休假、病人有事等实际问题，在原计划规定的时间点未能获取到监测值。若忽视时间点不统一，或收集资料的时间是在一段时间内获取的具体问题，采用常规的重复测量资料的方差分析，删除有缺失数据的观察单位，有可能引致估计偏差加大。变系数模型就是将其均值参数和方差成分看作是随观测时间变化的函数，利用模式混合模型的基本原理，通过贝叶斯惩罚样条构建变系数模型，从而获得更精确的参数估计。该方法不但允许均值参数随时间变化，而且允许模型中方差成分随缺失时间变化，实际应用中具有较大的灵活性。

变系数模型的基本原理

假定监测资料{Yi(t)}和{Xi(t)}分别表示第i个体在t时刻的反应变量与协变量值，且Y为连续型随机变量，则有:

{Yi(t)|Xi(t)}∈N(μi(t)，Vi(s，t))，

式中 μi(t)=Xi(t)β，Vi(s，t)=cov(Yi(s)，Yi(t)|Xi(s)，Xi(t))，(s≤t)。令 θ为模型参数，则{Yi(t)|Xi(t)}∈F(θ;Xi(t))，设Ui为研究对象规定的退出时间，此时可将反应变量的联合密度函数分解为:

f(Yi(t)，Xi(t)，Ui)=f(Yi(t)|Xi(t)，Ui)f(Ui|Xi(t))(1)模型中{Yi(t)|Xi(t)，Ui}∈F{θ(U);Xi(t)}，若反应变量为连续型随机变量，则分布函数记作:

{Yi|Xi，Zi，bi，Ui}∈N［Xiβ(u)+Zibi，Ri{φ(u)}］(2)

{bi|Ui，δi}∈N(0，Gi{φ(u)})，式中Xi:ni×p维固定效应矩阵，Zi:ni×q维随机效应矩阵。可见变系数模型中不仅位置参数随时间而变化，且离散参数也可以是一个时间依赖参数，即βi和φi随u的变化而变化。

1．参数估计

令随访时间分布为f(u|x，π)，Θ为参数集合，则第i个体的似然函数记作

2．贝叶斯惩罚样条

惩罚样条函数是非参数法的一种，它在参数设定方面具有很大灵活性。贝叶斯法可利用样本信息结合先验分布得到后验，利用后验分布进行推断。贝叶斯惩罚样条结合了以上两方面优点。θ(·)的三次样条函数记作:

式中 φ =(η0，η1，γ1，γ2，…，γk)T，γk为平滑参数，而 vk为节点，它是U的第k/(k+1)百分位数。为保证曲线的灵活性应尽量选择较多节点，但过多的节点数也可能会造成过度拟合，Ruppert认为采用5-20个节点效果较好〔1〕。

3．敏感性分析

假定包含缺失信息的观测值与未包含缺失信息的观测值，建模时都具有相同的参数估计结果〔2〕。为了验证在适当延长观测时间后是否仍能得到相似的结论，实际问题分析中有必要进行敏感性分析，因为该假定无法采用已观测到的数据进行验证。若在模型后增加敏感项ω(u)(tij－u)+，令模型中无法证明的部分为 ω(u)，ω(u)=a{max(U) －u}/range(U)，则称 a为敏感性参数。

高血压社区规范化管理的变系数模型实现

1．随机抽取全国高血压社区规范化管理项目山西分中心负责的太原市四个社区卫生服务中心1 275名高血压患者的监测资料。根据危险因素量化分层标准，将社区规范化管理的高血压患者分为1级、2级和3级，分别实施不同的管理程序。1级和2级高血压患者分别在基线调查后每隔三个月随访一次，一年内随访四次;3级高血压患者隔月随访一次，一年内随访六次。本文选取3级高血压患者301名作为研究对象，共有264名3级高血压患者完成了六次随访，失访率为12.3%。由于本次调查采取患者主动前来接受调查的方法，许多研究对象未能按照研究方案规定的两个月一次接受调查，表现为不规则的观测时间和连续的缺失时间。分别以收缩压或舒张压(mmHg)作为反应变量，年龄(岁)、性别、高血压病程(年)和社区中心作为自变量，构建不同参数随观测时间变化的变系数模型。

采用 SAS9.2整理分析数据集，通过 Win-BUGS14.0建立变系数模型，并通过 R软件中的R2WINBUGS包来调用WinBUGS程序，通过 Gibbs抽样进行MCMC迭代，迭代进行15000次，退火5000次。根据迭代序列图可见所有参数的迭代历史轨迹呈现“毛毛虫”状，概率密度均呈现中间高，两边低的倒钟形分布，且自相关图逐渐趋近于0，说明参数值达到收敛。最终分析结果输出到R软件中。

图1(1)～(7)分别给出了截距、年龄、性别、高血压病程、社区二、社区三和社区四的回归系数随时间的变化情况。由图可见，无论收缩压还是舒张压，年龄参数随观测时间延长呈下降趋势，表明年龄较大的高血压患者血压控制更理想，且舒张压值下降的幅度比收缩压值大;收缩压的性别参数随时间呈上升趋势，而舒张压的性别参数随访期内基本不变;收缩压的高血压病程参数基本保持为0，而舒张压则上升且为正值，表明高血压病程长的患者舒张压控制不理想，而对收缩压的影响则较小;中心2和中心3不论对收缩压还是舒张压，在整个观测时间内基本不变，而中心3则呈现下降趋势，说明中心3的血压控制较好。图1(8)～(10)分别给出随机系数截距、随机系数斜率以及个体间误差随观测时间的变化情况。随机系数截距和斜率随时间延长基本为一个常数，收缩压的个体间误差呈现下降趋势，而舒张压则略呈上升趋势。

图1 固定效应参数与随机效应参数随观测时间(天)变化图

表1 分别以收缩压和舒张压为反应变量的参数估计结果

表1为由变系数模型得到的收缩压和舒张压的各个参数估计值及其可信区间。其中收缩压的参数，性别的可信区间未包括0，表明女性的收缩压低于男性;舒张压的参数，性别和中心3的可信区间未包括0，表明女性舒张压低于男性，社区3的舒张压均值低于社区1。

分别设定敏感性参数a为5和10，比较各个社区分性别的血压估计值变化情况。结果见表2。可以看出三种条件下的估计结果均较为接近。

表2 不同敏感性参数的分析结果

讨 论

随访研究中经常会遇到监测时间不规则问题。完全数据分析中，个体内相关系数的变化一般仅影响方差成分，而对回归参数的影响较小。当数据中包含信息缺失时，个体内相关系数的变化也可能会引致回归系数的变化〔3〕。因此本文建模时也同时建立了方差成分随时间变化的平滑函数，以得到更为精确的分析结果。

由于贝叶斯分析结果对先验选取非常敏感，若错误地选择先验，可能引致错误的参数估计结果。因此，在参数分布未知的情况下，应选择较为保守的无信息先验。Crainiceanu在平滑参数估计中将方差成分的先验分布设定为Gamma(A，B)，并指出当A≤0.01且B≤‖b‖2/2时，先验分布与后验分布之间相互独立〔5〕。当不满足上述条件时，先验分布就会对后验分布有影响。因此本文在先验选取上对位置参数设定了无信息先验，而将方差成分的先验分布设定为 Gamma(10－6，10－6)。

WinBUGS软件是进行贝叶斯分析最常用的工具。该软件通过Gibbs抽样和MCMC方法来进行复杂统计模型推断，其优点是运算速度快，结果准确。但它对数据输入格式要求严格，在数据库操作方面还无法与目前常用的SAS和R软件相提并论，且由于其储存空间有限，当数据量较大时常常会出现超过记忆容量的情况，实际应用中受到较大限制。Sibylle〔5〕认为R软件推出的R2WinBUGS包允许在R环境下运行WinBUGS程序，可解决此问题。本文通过四个社区高血压规范化管理实例，进一步证实并很好地解决了这一问题。

1．Diggle P，Heagerty P，Liang KY，et al．Analysis of Longitudinal Data．New York:Oxford University Press，2002．

2．Daniels M，Hogan，J．Missing data in longitudinal studies:strategies for bayesian modeling and sensitivity analysis．Monographs on Statistics and Applied Probability，2008，101．New York:Chapman ＆ Hall．

3．Li S，Joseph WH．Varying-coefficient models for longitudinal processes with continuous-time informative dropout．Biostatistics，2010，11(1):93-110．

4．Ciprian MC，David R，Wand MP．Bayesian analysis for penalized spline regression using WinBUGS．Johns Hopkins University，Dept．of Biostatistics Working Papers，2007:40-74．

5．Sibylle S，Uwe L，Andrew G．R2WinBUGS:A package for running Win-BUGS from R．Journal of Statistical Software，2005，12(3):1-16．