基于Weibull回归模型的遗传方差分量模型研究，本文主要内容关键词为：模型论文,方差论文,分量论文,Weibull论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
遗传流行病学家系研究中，由于共享的遗传和环境因素作用，家庭内各成员的疾病表型往往是非独立的，遗传方差分量模型常用于研究遗传和环境因素的作用大小。对数量性状和质量性状，可在广义线性混合模型(GLMM)的框架下用MCMC的方法估计家庭内遗传和环境方差分量[1,2]，对具有删失的生存时间资料，如疾病的发病年龄等，我们将上述方法进一步推广，在生存分析领域内构造生存时间资料的遗传方差分量模型。本文将带有随机效应的Weibull回归模型扩展到遗传方差分量模型，并利用MCMC方法估计固定效应参数和随机效应方差分量，用于分析家系资料中删失性状的遗传和环境因素的作用。
一、模型构造
此为Weibull回归模型，回归系数β和形状参数α是模型中待估参数。
当生存时间数据不满足独立性要求时，生存分析模型中常通过随机效应ω来指定类内相关。随机效应也称为脆弱(frailty)，通常以乘积的形式作用于基准危险率，表示同一类内成员享有共同的基准危险，不同类成员基准危险不同。脆弱项的分布研究较多的为gamma分布或对数正态分布，当假定ω为对数正态分布时，令T=ln（ω），上述Weibull回归模型可以方便地扩展到家系生存时间数据的遗传方差分量模型[3,4]。
可以看到，基于Weibull回归的遗传方差分量模型可看作单个脆弱模型的推广，其基准危险函数指定为Weibull分布，属参数生存模型。
和数量性状和质量性状的遗传方差分量模型一样，模型假设：①Hardy-Weinberg平衡；②没有异位显性；③没有基因环境交互作用；④随机婚配；⑤随机效应之间没有协方差。
二、参数估计方法
对数正态分布的脆弱模型边际似然很难明确地表达，往往通过惩罚偏似然函数得到近似边际似然。但目前基于Weibull回归的脆弱模型在常规统计软件如SAS，SPSS中难以实现，Bayes统计中马尔可夫链蒙特卡罗(markov chain monte carlo，MCMC)模拟技术的发展为许多复杂模型的参数估计提供了可能。在Bayes统计中，对未知参数的推断基于后验分布，许多后验统计量一般都可以归纳为对后验分布的积分计算。为避免复杂的高维积分的困难，可以随机模拟服从该后验分布的一个独立样本，用该样本来描述各种统计量，这就是蒙特卡罗(monte carlo，MC)积分；而该独立样本可以通过展开一个马氏链(markov chain，MC)来获得，即MCMC方法。为构造具有特定平稳分布的链，可采用Gibbs抽样的方法，依次从各参数的全条件分布中连续抽样直至收敛，而Gibbs抽样的方法可在WinBUGS软件中实现。本研究将Gibbs抽样方法应用到有多个随机效应的遗传方差分量模型，进行Weibull分布形状参数、回归系数和各方差分量的估计。对基于Weibull回归的遗传方差分量模型，可在WinBUGS软件中编写程序指定模型，生存时间服从Weibull分布，删失的生存时间指定为截断的Weibull分布，对数尺度参数为固定效应和随机效应之和。由于各随机效应服从均数向量为0的多元正态分布，在程序中难于指定，故需要对随机效应重新参数化，对核心家庭中的父亲、母亲和子女，分别指定不同的随机效应设计矩阵。由于重新参数化后的随机矩阵为不满秩状态，导致随机效应方差分量无明确的解析解，因此在WinBUGS软件中采用一种循环估计的方法。随机效应重新参数化和循环估计的方法可参见数量性状和质量性状的遗传方差分量模型研究[1,2]，不再赘述。参数先验分布的指定均采用无信息无验，回归系数指定为正态分布N(0，
)；随机效应的精度
指定为Pareto(1，0.01)分布，相当于方差在间隔[0，100]上的均匀分布；Weibull分布的形状参数指定为Gamma(1，)分布。
三、模拟研究
为了解基于Weibull回归的遗传方差分量模型参数估计的准确性，本文进行模拟研究，模拟核心家系资料，包括父亲、母亲和两个子女。模拟数据集在SAS软件中产生，协变量
来自二项分布(1，0.3)，
来自均数为0，方差为1.2的正态分布；随机效应项来自均数向量为0，方差协方差满足模型(3)的多元正态分布。生存时间资料来自Weibull分布，用均匀分布指定删失的比例为10%。指定回归系数
和方差分量
的真值分别为(2.0，0.5)和(0.2，1.5，0.8)，威布尔分布的形状参数α为4.3。模拟家系数为1000和250时，分别拟合20个家系数据集。
在软件WinBUGS中，采用循环估计方法，设定前后两次方差分量估计值相差容许值小于0.1。当循环过程中方差分量出现负估计时，将其设为
，继续循环。家系数为1000时每次运算共迭代10000次，burn-in=5000，用于“退火”，后5000次用于模拟抽样，thin=5。家系数为250时每次运算共迭代20000次，burn-in=10000，后10000次用于模拟抽样。最后根据每个数据集的参数估计值计算参数估计值的均数和标准误。方法评价从两方面：点估计和区间估计。点估计通过计算参数估计值的均数（标准误）及均数距真值的偏差，区间估计通过计算估计值的可信区间对真值的覆盖率来评价。
表1显示了家系数为1000和250时20个模拟数据集的参数估计结果。可以看到，Weibull分布的形状参数、固定效应和随机效应的方基分量，20个模拟数据集的参数估计值的均数都接近真值，绝对偏差均小于0.1，和家系数为1000时相比，家系数为250时参数估计值均数的标准误较大；用中位数和均数估计参数统计量没有明显差别，除Weibull分布形状参数外，95%可信区间对参数真值的覆盖率均大于80%，随机效应方差分量的估计准确率和固定效应相比没有明显不同。
表1 家系数为1000和250时20个数据集的WinBUGS模拟结果
*a：20次参数估计均数的均数（均数的标准误）和均数的最小最大值；b：20次参数估计中位数的中位数和中位数的最小最大值。下同。
四、讨论
生存分析领域中，研究生存时间和相关因素之间关系常用的回归模型有参数回归模型和半参数回归模型（如Cox比例危险模型），Weibull回归模型属参数回归模型的一种。Weibull分布用来拟合生存时间非常灵活，其危险率曲线可以单调上升（形状参数α＞1时）、单调下降（当α＜1时）或保持不变（当α=1时）等各种情形；另一方面，威布尔分布的生存函数、危险率函数和概率密度函数的形式相对简单；此外，Weibull回归模型是唯一具有比例危险和加速时间表达式的参数模型。由于这些性质，使得Weibull回归模型在生存分析领域中的应用极为广泛[5]。本研究在Weibull回归模型的框架下，将随机效应项合并到固定效应项中，用于分析家系资料中的生存时间资料，在估计协变量向量对生存时间影响的同时，估计随机效应的方差分量。不失一般性地，该思路可推广到临床或流行病学研究中重复测量或纵向研究中生存时间资料的分析。
近年来，非独立生存时间数据的统计方法已受到统计学家越来越多的关注。边际模型和脆弱模型是近年来研究较多的多变量生存模型。通常边际模型将类内生存时间的依赖性看做多余参数而不予指定，比较调整依赖性后的参数标准误和独立假设下的参数标准误。如果关联程度本身是研究的主要兴趣，则更常用到脆弱模型。脆弱模型通过引入一个共享的随机效应来描述类内生存时间之间的依赖性，假定每类中有一个随机效应（脆弱），以乘积的形式作用于基准危险率。脆弱的分布可以有多种形式，研究较多的有gamma分布和对数正态分布，本研究采用对数正态分布的形式，将单个随机效应的生存模型推广到多个随机效应，根据家系成员间的亲缘关系构造遗传方差分量模型，用于研究遗传因素和环境因素对删失的生存时间资料的影响。
Weibull回归模型基础上的遗传方差分量模型也可看作数量性状和质量性状遗传方差分量模型的推广，只要方差分量的指定符合生物学解释。因此在参数估计过程中，尤其是随机效应方差分量的估计中，本研究和以往研究[1,2]存在共同的问题。首先是模型的重新参数化问题。与通常的非独立数据不同，家系结构数据有需要特殊指定的相关关系，每个随机效应均服从多元正态分布，因此在数据分析时，需要将模型的随机效应重新参数化，将服从多元正态分布的随机效应之和分解为若干服从单变量正同分布的各随机效应之和，并且保证新的随机效应项之和满足随机效应的方差协方差结构。重新参数化的方法很多，本研究和质量性状遗传方差分量模型一样，为了得到更一致的参数估计，重新参数化的模型允许各方差分量的最终估计值为负，从而更符合生物学上的解释。从模拟结果中也可看出，当家系数为250时，对第一套参数真值的估计中，20个数据集中有2个数据集遗传方差分量
的最终估计值为负。另一个问题是采用循环的估计方法，和质量性状遗传方差分量模型类似，在Weibull回归模型中，生存时间的误差由Weibull分布指定，因此随机效应项中缺少随机误差效应项，随机效应方差中也没有随机误差方差，这样重新参数化后就存在部分“多余”方差分量，导致随机效应方差分量没有明确的解析解。在WinBUGS软件中采用循环估计的方法找到方差分量的近似解，但循环估计方法容易出现由于多参数间负相关导致参数值无法最终收敛的问题，但本次模拟实验中尚未碰到此类问题。