论文部分内容阅读
摘要：本文主要讨论了广义线性混合模型（即GLMM）模型的两种不同的参数估计方法，基于惩罚拟似然参数估计方法和MCMC贝叶斯参数估计方法。本文考虑借助开源软件R对癫痫病数据集epil进行统计推断，并与其他模型和参数估计方法进行比较。
关键词：广义线性混合模型;方差分量;参数估计;glmmPQL;MCMCglmm
Abstract： This paper mainly discusses two different parameter estimation methods of GLMM model， one is the penalty quasi likelihood parameter estimation method and the other is the MCMC parameter estimation method based on Bayesian， and this paper considers the use of open source software R to estimate the parameters of a specific data set （epil） Compared with other models and parameter estimation methods， it is concluded that the application of GLMM is very extensive and necessary.
Key Words：GLMM;variance component;parameter estimation;glmmPQL;MCMCglmm
在对实际问题建模时，当我们感兴趣的变量是一个正态分布随机变量时，通常考虑使用线性模型（LM），要求响应变量满足正态性和独立性，而实际数据会违背相互独立的假设，这些数据之间往往存在某种相关结构，由此发展了线性混合模型（LMM），在LMM中，放松了独立性要求，考虑响应变量之间的相关性。在后续的研究中发现除了满足正态分布的随机变量之外，还存在很多非正态分布，如泊松分布，二项分布等。對于这类分布，考虑对其均值进行建模，发展成广义线性模型（GLM），相比LM，GLM放松了要求响应变量的正态性，但是仍然需要满足独立性，但实际中数据之间的相关性不可忽略。忽略相关性将会导致偏差或者得出错误的结论。为了体现出数据之间的相关性，在广义线性模型的连接函数上加上了一个随机变量，也就是本文所研究的广义线性混合模型（GLMM）[3]。理论上，广义线性混合模型是广义线性模型和线性模型的结合，同时它也是线性混合模型的拓展。广义线性混合模型适合用来模拟存在相关性的区组效应的数据或纵向数据，此时不再要求响应变量满足正态性和独立性，模型设定更加灵活，更能反映出响应变量之间的相关性，这是广义线性混合模型的一个显著优点。
对于GLMM的参数估计，以往很多学者提出各种方法。针对GLM的参数估计首先构造似然函数，对随机效应进行处理，无论是对条件似然函数的积分处理还是对边际似然函数的极大似然估计[2]的过程中都存在很大的困难，尤其是在对随机效应的积分处理往往伴随各种高维积分。后来的研究中，有学者提出拟似然（quasi—likelihood）函数[1]，详细参考McCullagh[3]，拟似然的提出这类模型提供了推断的基础，后来的学者提出了基于惩罚你似然函数的参数估计方法。对于多数非正态模型，使用Monte-Carlo积分方法[2]。除此之外，还可以在先验分布[4]中使用重要性抽样或者Gibbs抽样避免数值积分，而对于极大似然估计方法最常用的是EM算法和MCEM算法。
参数估计方法比较丰富，但目前针对广义线性模型的参数估计方法主要是基于拟似然或者惩罚拟似然角度出发的极大似然估计，而从贝叶斯后验角度出发进行的参数估计目前还并不多，一般进行贝叶斯参数估计借助于WinBUGS或OpenBUGS等软件进行，遇到较大数据集的情况时运行速度较慢，本文借助于R软件中的MCMCglmm包中的MCMCglmm函数[5]可以很好的避免运行速度过慢的缺陷，同时可以达到一个很好的参数估计的效果。本文主要针对R中自带的数据集epil进行惩罚拟似然参数估计和马尔可夫蒙特卡洛参数估计两种估计方法的结果比较。
1 模型介绍和分析
广义线性混合模型（GLMM）是广义线性模型（GLM）的拓展，GLM需要包含三个部分： （1）响应变量的分布 ;（2）线性预测器 ;（3）连接函数 。模型可表示为：
其中， 为Y的分布函数， 为Y的期望， 是单调可导的连接函数， 为线性预测器。而对广义线性模型中参数 进行参数估计时比较常用的一种方法是迭代加权最小二乘（IWLS）算法，具体算法的细节可以参考McCullagh 和 Nelder（1978）的工作。
而广义线性混合模型（GLMM）与广义线性模型相似，需要以下四个组成部分，除了GLM的三部分还需要添加一个随机效应： ，其中 为 的方差;更具体的说： ，其中 ，且同时有 。令 ，则可以得到样本的似然函数：
其中 表示给定随机项 时Y的条件概率密度函数， 为b的概率密度函数，感兴趣的是模型中参数 的估计，即估计出固定效应 和方差分量 。
2 数据集介绍和描述性分析
本文所考虑使用的数据集为R中MASS包中的癫痫病数据集epil，由Thall和Vail（1990）年给出的59个癫痫病患者的2周内的发病次数，患者被随机分配到两个不同的实验组，一个组进行新药治疗（Trt=1），一个组进行对照的安慰剂治疗（Trt=0），具体介绍可在R中进行查看。
对数据集epil进行描述性分析，新药治疗组的样本相关性比安慰剂组的样本相关性更大，同时从药物治疗组和安慰剂对照组不同时期发病次数箱线图（图1）无法断定癫痫犯病次数随时期的变化产生而产生明显的变化，但是可以明显地看出新药治疗组的发病次数比安慰剂对照组的发病次数水平更低。 3 参数估计
接下来考虑两种不同的参数估计方法来对广义线性混合模型（GLMM）的参数进行估计，基于惩罚拟似然方法（PQL）的极大似然估计方法，具体的细节可参考干晓蓉（2007）的论文工作[6];基于贝叶斯后验分布出发的参数估计方法，具体的细节可以参考Jarrod Hadfield的工作。
3.1 PQL参数估计法
拟似然方法并不要求响应变量是某一个具体已知的分布，只需知道响应变量的均值和方差，在大样本情形下可近似为正态分布。然而为了减少方差估计的无偏性，采用添加惩罚项的方法来提高估计的精度。假设已知响应变量 的条件分布的均值和方差为：
，
连接函数 的逆函数为 ，即： ，其中 ，有 ， ，有 ，其中I为对角阵，对角线元素为单位阵， 。得到拟似然函数：
借助于拉普拉斯近似，进一步可以得到惩罚拟似然函数为：
针对上式进行求导可以得到方程并进一步得到参数的极大似然估計量。接下来考虑在癫痫病数据集上进行实例应用。我们考虑在epil数据集上应用GLMM模型通过PQL方法来估计参数，首先考虑建立模型：
利用R软件中的glmmPQL进行参数估计，考虑两种情况，将period的四个阶段其为数值型变量和因子型变量。
3.2 贝叶斯参数估计
贝叶斯统计推断是从后验分布出发，通常对未知参数有一定的先验信息或无信息先验，当给定先验先验信息，可求后验分布，进而完成贝叶斯参数计。
参考Jarrod Hadfield的工作，考虑第i个数据的概率为： ，假设 服从泊松分布，则有：
其中 为泊松分布规范参数。解释变量之间的线性有下列线性关系：
其中X为与固定效应有关的设计矩阵，Z是与随机效应有关的设计矩阵，e为残差，假设 来自多元正态分布，满足：
利用MCMC抽样方法我们可以得到：
按照 更新参数 ， ，其中C为固定模型的系数矩阵，满足：
其中 ，B为固定效应的先验方差， 和 为从多元正态分布中随机抽取的样本：
且 ， 即为从下列分布中抽样的样本： 。
可借助Jarrod Hadfield[5]编写的MCMCglmm包中的MCMCglmm函数实现广义线性混合模型的参数估计，由于对癫痫病数据集缺乏先验信息，故不设置先验信息，或者设置一个无信息先验信息。利用MCMCglmm函数对数据集进行拟合，考虑每10个间隔进行一次抽样，共抽取13000个样本，除去预烧期的3000个抽样样本，进行贝叶斯参数估计，同样考虑period的两种情形，模型的DIC值分别为1153.947和1158.283，两种情形下随机效应的参数估计结果为：0.252和0.249.
根据表1和表2中的glmmPQL参数估计结果得出结论：新药的回归系数是显著的，年龄并不显著，基础发病次数是显著的，但基础发病次数与新药的交叉作用不显著，同时从整体上看阶段是一个显著的回归项，具体到四个阶段的每一阶段，只有第四阶段的作用是显著的，其它的回归项的显著性和估计结果基本相似。同时两种不同情形下的随机效应估计为0.443和0.444，两者相差不大。这与前面的数据描述性分析是相互对照的。
表1和表2中同样给出了贝叶斯框架下的MCMCglmm参数估计结果，当将period视为数值型变量时，此时随机效应的估计值为 ：0.252，可信区间为（0.129，0.402）。当将period视为因子型变量，此时随机效应的估计值为 ：0.250，可信区间为（0.131，0.409）。
看出无论是将period看成数值型变量还是因子型变量时，各个时期都不显著，这与PQL方法下的结果是不同的。同时在MCMCglmm包中还可以给出随机效应的迭代轨迹图像以及随机效应的概率密度函数图像，可以通过迭代轨迹图判断是否收敛以及基于后验分布的统计推断是否合理。
4 结论
文章主要考虑了两种广义线性混合模型的参数估计方法，PQL方法下的广义线性模型方差分量比MCMC方法下的要略微偏大一些，但两种方法下大部分的回归结果是相似的，这也再次证明了这两种方法的参数估计结果是比较准确的，同时第二种方法下的贝叶斯参数估计方法的效率较高，运行时间比利用WinBUGS或OpenBUGS等软件的效率要高很多，这也是这个方法的一个显著的优势。同时本文也将一般的参数估计方法上拓宽了，将贝叶斯的参数估计方法应用过来，且参数估计的结果表明二者估计的结果是相似且准确的。
参考文献
[1] 曹红艳，刘桂芬，曾平，等.预测性伪似然法和贝叶斯法广义线性混合模型估计[J].中国药物与临床，2008，8（011）：861-863.
[2] 陈永成，STEPHEN W.RAUDENBUSH.用蒙特卡罗方法求取广义线性混合模型之最大似然估计：应用于求取乳癌死亡率之小区域估计[J].应用概率统计，2006，22（1）：69-80.
[3] McCulloch CE，Searle SR.Generalized，linear and mixed models.New York：Wiley，2001.
[4] Gelman A.Prior distribution for variance parameters in hierar-
chical models. Bayesian Analysis， 2006， 1（3）：515-534.
[5] Hadfield J D . MCMC Methods for Multi-Response Generalized Linear Mixed Models： The MCMCglmm R Package[J]. Journal of Statistical Software， 2010， 33（02）：1-22.
[6] 干晓蓉.广义线性混合模型[J].昆明理工大学学报（理工版），2007，32（4）：107-113.
关键词：广义线性混合模型;方差分量;参数估计;glmmPQL;MCMCglmm
Abstract： This paper mainly discusses two different parameter estimation methods of GLMM model， one is the penalty quasi likelihood parameter estimation method and the other is the MCMC parameter estimation method based on Bayesian， and this paper considers the use of open source software R to estimate the parameters of a specific data set （epil） Compared with other models and parameter estimation methods， it is concluded that the application of GLMM is very extensive and necessary.
Key Words：GLMM;variance component;parameter estimation;glmmPQL;MCMCglmm
在对实际问题建模时，当我们感兴趣的变量是一个正态分布随机变量时，通常考虑使用线性模型（LM），要求响应变量满足正态性和独立性，而实际数据会违背相互独立的假设，这些数据之间往往存在某种相关结构，由此发展了线性混合模型（LMM），在LMM中，放松了独立性要求，考虑响应变量之间的相关性。在后续的研究中发现除了满足正态分布的随机变量之外，还存在很多非正态分布，如泊松分布，二项分布等。對于这类分布，考虑对其均值进行建模，发展成广义线性模型（GLM），相比LM，GLM放松了要求响应变量的正态性，但是仍然需要满足独立性，但实际中数据之间的相关性不可忽略。忽略相关性将会导致偏差或者得出错误的结论。为了体现出数据之间的相关性，在广义线性模型的连接函数上加上了一个随机变量，也就是本文所研究的广义线性混合模型（GLMM）[3]。理论上，广义线性混合模型是广义线性模型和线性模型的结合，同时它也是线性混合模型的拓展。广义线性混合模型适合用来模拟存在相关性的区组效应的数据或纵向数据，此时不再要求响应变量满足正态性和独立性，模型设定更加灵活，更能反映出响应变量之间的相关性，这是广义线性混合模型的一个显著优点。
对于GLMM的参数估计，以往很多学者提出各种方法。针对GLM的参数估计首先构造似然函数，对随机效应进行处理，无论是对条件似然函数的积分处理还是对边际似然函数的极大似然估计[2]的过程中都存在很大的困难，尤其是在对随机效应的积分处理往往伴随各种高维积分。后来的研究中，有学者提出拟似然（quasi—likelihood）函数[1]，详细参考McCullagh[3]，拟似然的提出这类模型提供了推断的基础，后来的学者提出了基于惩罚你似然函数的参数估计方法。对于多数非正态模型，使用Monte-Carlo积分方法[2]。除此之外，还可以在先验分布[4]中使用重要性抽样或者Gibbs抽样避免数值积分，而对于极大似然估计方法最常用的是EM算法和MCEM算法。
参数估计方法比较丰富，但目前针对广义线性模型的参数估计方法主要是基于拟似然或者惩罚拟似然角度出发的极大似然估计，而从贝叶斯后验角度出发进行的参数估计目前还并不多，一般进行贝叶斯参数估计借助于WinBUGS或OpenBUGS等软件进行，遇到较大数据集的情况时运行速度较慢，本文借助于R软件中的MCMCglmm包中的MCMCglmm函数[5]可以很好的避免运行速度过慢的缺陷，同时可以达到一个很好的参数估计的效果。本文主要针对R中自带的数据集epil进行惩罚拟似然参数估计和马尔可夫蒙特卡洛参数估计两种估计方法的结果比较。
1 模型介绍和分析
广义线性混合模型（GLMM）是广义线性模型（GLM）的拓展，GLM需要包含三个部分： （1）响应变量的分布 ;（2）线性预测器 ;（3）连接函数 。模型可表示为：
其中， 为Y的分布函数， 为Y的期望， 是单调可导的连接函数， 为线性预测器。而对广义线性模型中参数 进行参数估计时比较常用的一种方法是迭代加权最小二乘（IWLS）算法，具体算法的细节可以参考McCullagh 和 Nelder（1978）的工作。
而广义线性混合模型（GLMM）与广义线性模型相似，需要以下四个组成部分，除了GLM的三部分还需要添加一个随机效应： ，其中 为 的方差;更具体的说： ，其中 ，且同时有 。令 ，则可以得到样本的似然函数：
其中 表示给定随机项 时Y的条件概率密度函数， 为b的概率密度函数，感兴趣的是模型中参数 的估计，即估计出固定效应 和方差分量 。
2 数据集介绍和描述性分析
本文所考虑使用的数据集为R中MASS包中的癫痫病数据集epil，由Thall和Vail（1990）年给出的59个癫痫病患者的2周内的发病次数，患者被随机分配到两个不同的实验组，一个组进行新药治疗（Trt=1），一个组进行对照的安慰剂治疗（Trt=0），具体介绍可在R中进行查看。
对数据集epil进行描述性分析，新药治疗组的样本相关性比安慰剂组的样本相关性更大，同时从药物治疗组和安慰剂对照组不同时期发病次数箱线图（图1）无法断定癫痫犯病次数随时期的变化产生而产生明显的变化，但是可以明显地看出新药治疗组的发病次数比安慰剂对照组的发病次数水平更低。 3 参数估计
接下来考虑两种不同的参数估计方法来对广义线性混合模型（GLMM）的参数进行估计，基于惩罚拟似然方法（PQL）的极大似然估计方法，具体的细节可参考干晓蓉（2007）的论文工作[6];基于贝叶斯后验分布出发的参数估计方法，具体的细节可以参考Jarrod Hadfield的工作。
3.1 PQL参数估计法
拟似然方法并不要求响应变量是某一个具体已知的分布，只需知道响应变量的均值和方差，在大样本情形下可近似为正态分布。然而为了减少方差估计的无偏性，采用添加惩罚项的方法来提高估计的精度。假设已知响应变量 的条件分布的均值和方差为：
，
连接函数 的逆函数为 ，即： ，其中 ，有 ， ，有 ，其中I为对角阵，对角线元素为单位阵， 。得到拟似然函数：
借助于拉普拉斯近似，进一步可以得到惩罚拟似然函数为：
针对上式进行求导可以得到方程并进一步得到参数的极大似然估計量。接下来考虑在癫痫病数据集上进行实例应用。我们考虑在epil数据集上应用GLMM模型通过PQL方法来估计参数，首先考虑建立模型：
利用R软件中的glmmPQL进行参数估计，考虑两种情况，将period的四个阶段其为数值型变量和因子型变量。
3.2 贝叶斯参数估计
贝叶斯统计推断是从后验分布出发，通常对未知参数有一定的先验信息或无信息先验，当给定先验先验信息，可求后验分布，进而完成贝叶斯参数计。
参考Jarrod Hadfield的工作，考虑第i个数据的概率为： ，假设 服从泊松分布，则有：
其中 为泊松分布规范参数。解释变量之间的线性有下列线性关系：
其中X为与固定效应有关的设计矩阵，Z是与随机效应有关的设计矩阵，e为残差，假设 来自多元正态分布，满足：
利用MCMC抽样方法我们可以得到：
按照 更新参数 ， ，其中C为固定模型的系数矩阵，满足：
其中 ，B为固定效应的先验方差， 和 为从多元正态分布中随机抽取的样本：
且 ， 即为从下列分布中抽样的样本： 。
可借助Jarrod Hadfield[5]编写的MCMCglmm包中的MCMCglmm函数实现广义线性混合模型的参数估计，由于对癫痫病数据集缺乏先验信息，故不设置先验信息，或者设置一个无信息先验信息。利用MCMCglmm函数对数据集进行拟合，考虑每10个间隔进行一次抽样，共抽取13000个样本，除去预烧期的3000个抽样样本，进行贝叶斯参数估计，同样考虑period的两种情形，模型的DIC值分别为1153.947和1158.283，两种情形下随机效应的参数估计结果为：0.252和0.249.
根据表1和表2中的glmmPQL参数估计结果得出结论：新药的回归系数是显著的，年龄并不显著，基础发病次数是显著的，但基础发病次数与新药的交叉作用不显著，同时从整体上看阶段是一个显著的回归项，具体到四个阶段的每一阶段，只有第四阶段的作用是显著的，其它的回归项的显著性和估计结果基本相似。同时两种不同情形下的随机效应估计为0.443和0.444，两者相差不大。这与前面的数据描述性分析是相互对照的。
表1和表2中同样给出了贝叶斯框架下的MCMCglmm参数估计结果，当将period视为数值型变量时，此时随机效应的估计值为 ：0.252，可信区间为（0.129，0.402）。当将period视为因子型变量，此时随机效应的估计值为 ：0.250，可信区间为（0.131，0.409）。
看出无论是将period看成数值型变量还是因子型变量时，各个时期都不显著，这与PQL方法下的结果是不同的。同时在MCMCglmm包中还可以给出随机效应的迭代轨迹图像以及随机效应的概率密度函数图像，可以通过迭代轨迹图判断是否收敛以及基于后验分布的统计推断是否合理。
4 结论
文章主要考虑了两种广义线性混合模型的参数估计方法，PQL方法下的广义线性模型方差分量比MCMC方法下的要略微偏大一些，但两种方法下大部分的回归结果是相似的，这也再次证明了这两种方法的参数估计结果是比较准确的，同时第二种方法下的贝叶斯参数估计方法的效率较高，运行时间比利用WinBUGS或OpenBUGS等软件的效率要高很多，这也是这个方法的一个显著的优势。同时本文也将一般的参数估计方法上拓宽了，将贝叶斯的参数估计方法应用过来，且参数估计的结果表明二者估计的结果是相似且准确的。
参考文献
[1] 曹红艳，刘桂芬，曾平，等.预测性伪似然法和贝叶斯法广义线性混合模型估计[J].中国药物与临床，2008，8（011）：861-863.
[2] 陈永成，STEPHEN W.RAUDENBUSH.用蒙特卡罗方法求取广义线性混合模型之最大似然估计：应用于求取乳癌死亡率之小区域估计[J].应用概率统计，2006，22（1）：69-80.
[3] McCulloch CE，Searle SR.Generalized，linear and mixed models.New York：Wiley，2001.
[4] Gelman A.Prior distribution for variance parameters in hierar-
chical models. Bayesian Analysis， 2006， 1（3）：515-534.
[5] Hadfield J D . MCMC Methods for Multi-Response Generalized Linear Mixed Models： The MCMCglmm R Package[J]. Journal of Statistical Software， 2010， 33（02）：1-22.
[6] 干晓蓉.广义线性混合模型[J].昆明理工大学学报（理工版），2007，32（4）：107-113.