一类线性混合模型的广义估计方程方法应用

余俊飞,赵慧秀

(南京理工大学 理学院,江苏 南京 210094)

0 引言

含有两个方差分量的线性混合模型是常用的纵向数据统计分析模型，它在线性模型中引入随机效应，建模了试验个体的不均匀性，并且刻画了同一个体内部观测值之间的相关性.该模型在生物学、医学、经济学等领域具有广泛的应用，因此对于该类模型参数估计的研究具有重要意义.

线性混合模型又称为随机效应模型，已有文献中有许多关于模型参数估计的研究.Liard and Ware[1]提出两阶段随机效应模型，讨论了结合经验贝叶斯和极大似然估计的重复测量模型拟合方法.Dempster and Selwyn[2]阐述了含有两个方差分量线性混合模型的统计计算方法，并比较了基于限制极大似然估计的EM算法和牛顿迭代法的效率.王松桂[3]概述了线性混合模型在参数估计、假设检验等方面的重要成果，并且提出了模型参数的一种谱分解估计.吴密霞等[4]在一组简单条件下证明了线性混合模型固定效应和方差分量可以同时达到最优估计.Gumedze and Dunne[5]综述了线性混合模型不同组成部分的参数估计和推断，重点讨论了方差参数的估计以及随机效应的推理过程.以上方法多是基于数据分布的假设，对于纵向数据的统计分析，Liang and Zeger[6-7]提出的广义估计方程具有广泛的应用.该方法不假设响应变量的分布，通过给定的工作相关阵替代真实的个体内部相关矩阵，其中工作相关阵的选择不影响回归参数及其方差估计的相合性[8].

本文结合广义估计方程处理一类线性混合模型的参数估计问题，该类模型中含有两个方差分量，其中随机效应的方差用来刻画个体内部的相关性.在此基础上可以计算出方差分量的矩估计，并应用广义估计方程估计模型参数.文中以一组大鼠繁殖数据进行了实例分析.

1 基本模型和参数估计

1.1 线性混合模型

考虑纵向数据下的线性混和模型，假设有m个被试个体,xij和yij分别是第i(i=1,2,…，m)个个体在第j(j=1,2,…，ni)次观测的协变量和响应变量.Yi是ni×1的个体响应向量，Xi为已知的ni×p协变量矩阵，Zi是已知的ni×q维设计矩阵，含有两个方差分量的线性混合模型表示为

Yi=Xiβ+Zibi+εi，

(1)

E(Yi)=Xiβ，

(2)

则模型中个体的响应变量服从多元正态分布.这样的模型在纵向数据分析中经常遇到，通过最大化对数似然函数可以求出模型的极大似然估计.

1.2 广义估计方程

var(yij)=φv(μij)，

其中：φ是散布参数；v是已知的方差函数.

基于以上假设，Weddueburn[9]提出了如下的拟似然估计方程

(3)

其中：

对于个体内部具有相关性的纵向数据，Vi的结构是非对角矩阵.Liang and Zeger在此基础上提出了广义估计方程方法.

定义一个对角矩阵

(4)

(5)

方程(5)不同工作相关阵下相关系数的估计方法相同.根据α估计的不同方式，可以迭代求出相应的参数估计值.

则

(6)

1.3 模型参数估计

含有两个方差分量线性混合模型的边际均值和协方差矩阵计算公式为

E(Yi)=E[E(Yi|bi)]，

Vi=cov[E(Yi|bi)]+E[cov(Yi|bi)],

(7)

由此可以得到式(2)的相应结果.

首先对随机截距模型进行分析，即模型中只含有单个随机效应，其协方差矩阵为

Jni是元素全为1的ni×ni矩阵.从Vi的结构可以看出该模型刻画了个体内部观测值之间的等相关性，这与广义估计方程中假设的可交换相关矩阵是等价的.在广义估计方程框架下可以采用矩估计和拟加权最小二乘等方法估计等相关矩阵下的相关系数，并通过对应关系求出模型的方差分量.令

N=n1+n2+…+nm，rij=

(Yij-μij)/v(μij)1/2，

则α的矩估计为

(8)

根据最小二乘中估计方差的思想，可以采用如下公式估计散布参数φ.

(9)

根据方程(5)(8)(9)进行迭代,可以求出回归参数β、相关系数α和散布参数φ.显然β是对应模型(1)中固定效应参数的估计，而方差分量估计的对应结果为

(10)

对含有q个随机效应参数的随机系数模型进行分析，这时模型协方差矩阵的形式为

令

(11)

(12)

然后通过方程(3)、(11)、(12)进行迭代，即可求出模型未知参数估计值.

2 实例分析

本文研究的大鼠繁殖数据可在文献[2]中获得，主要目的是评估试验剂量对大鼠繁殖能力的影响.试验中30只母鼠被随机分配到3个不同的处理组，分别是对照组、高剂量组、低剂量组，每组均有10只母鼠.母鼠生产后记录每只幼鼠的体重，最后得到27窝幼鼠数据，并以幼鼠的体重作为评估试验效应的指标.

每只母鼠的情况具有一定的差异，一般认为来自同一窝的幼鼠之间存在等相关性，运用含有两个方差分量的线性混合模型来分析该数据是合适的.模型的响应变量是幼鼠的体重(weight)，相应的协变量为所在窝中母鼠是否在高剂量组(high)、母鼠是否在低剂量组(low)、幼鼠性别(sex)、所在窝中幼鼠只数(high).其中幼鼠平均体重和所在窝中幼鼠只数的关系如图1所示.

图1 幼鼠平均体重与幼鼠只数的散点图

基于以上数据特征建立如下线性混合模型

(13)

应用广义估计方程方法得到模型参数，结果如表1所列.为评估该方法的估计效果，表中包含了极大似然估计的结果用于对比分析.

从表1结果来看，两种方法的参数估计结果大体上相同，结果表明药物剂量、幼鼠性别、每窝幼鼠只数对幼鼠的体重有明显影响，这些情况与数据的直观分析结果是一致的.表1结果显示广义估计方程方法回归参数估计的标准差更小，该方法不假设响应变量的联合分布，因此能够得到稳健的估计结果.

随机模拟能够评估方法的有效性，以表1中应用极大似然估计得到的结果为真值，即

表1 实例数据参数估计结果

应用统计软件依据公式(13)生成模拟数据，首先随机生成一组数据，应用两种方法估计相关参数.而后生成100组数据，计算各参数估计的平均值，其结果如表2所列.

表2的模拟数据分析结果表明两种方法固定效应的参数估计与真值差异不大.广义估计方程固定效应估计的标准差较小，对于方差分量，在已知分布情形下，极大似然估计更接近真值.

表2 模拟数据参数估计结果

3 结语

文中运用随机截距模型拟合了大鼠繁殖数据，采用广义估计方程方法得到的参数估计取得了较好的效果，这说明提出的方法是可行的.实际应用中，随机系数模型常用来建模时间相关的纵向数据，该模型同样可以应用广义估计方程方法进行分析.