含缺失协变量的混合效应模型的多元单边检验

齐培艳，仵旭靓，段西发

（太原科技大学 应用科学学院，山西 太原，030024）

0 引言

人类免疫缺陷病毒（Human Immunodeficiency Virus，HIV）病毒动力学研究是近年来艾滋病研究的热门，非线性混合效应（NLME）模型［1-2］、半参数非线性混合效应（SNLME）模型［3］被广泛应用于HIV病毒动力学、药代动力学分析和消长研究。这些模型通常基于生成数据的底层机制，其参数具有一定意义的物理解释，如抗HIV治疗期间病毒衰减率或随机效应的方差等参数应该是非负的，要满足这些条件就需要对参数进行自然限制或约束。因此，应该在包含限制的前提下对参数进行单边或序约束检验，在假设检验中加入自然参数限制或约束会产生比相应的无约束检验更有效的结果［4］。

约束统计推断一直以来都受到极大的关注，Silvapulle和Sen［5］对其研究进展进行了综述。然而，大多数文献关注的是方差分量和均值的检验［6-7］，对回归模型的单边检验的研究较为有限。Brice等［8］对潜变量模型（LVMs）进行了修正的 Wald检验；Verbeke和 Molenberghs［9］讨论了线性混合效应（LME）模型的约束推理；Qeadan等［10］证明了一类线性混合模型的广义似然比检验（LRT）和Wald检验在检验方差分量时的等价性；Zhou等［11］从HIV病毒动力学模型出发，提出了具有截尾响应的NLME模型的基于近似似然的多元单边检验，仿真结果表明所提出的检验比相应的双边检验或无约束检验更有效；Wang［12］对NLME模型中的多元单边检验发展了似然比检验（LRT）和Wald检验，并证明了它们的有效性要强于传统的双边检验。本文主要考虑SNLME模型中固定参数的单边假设检验。

在实践中，数据往往会在应答或协变量中出现缺失，使用完全情况（CC）方法丢弃所有的不完全观测数据，可能会使结果出现偏差［13］。为了减少这种误差，Wang等［14］和Zhou等［15］分别讨论了基于多重插补和经验似然的具有缺失协变量的NLME模型的多元单边检验。本文针对具有缺失协变量的半参数非线性混合效应模型中固定参数的多元单边检验问题，提出了基于多重插补的似然比检验（LRT）和Wald检验，并通过实例分析对HIV病毒动力学研究中双相指数病毒衰变模型的病毒衰减率是否严格为正进行了检验。这里假设缺失数据是随机缺失的，即缺失可能与观察到的数据有关，而不是与缺失的数据有关［16］。

本文第1节对SNLME模型及多元单边检验进行了介绍；第2节讨论了具有缺失协变量的SNLME模型的多重插补方法；第3节，基于多重插补的多个估计提出了两种检验：似然比检验（LRT）和Wald检验；第4节通过HIV实例分析对提出的方法进行了有效性检验；第5节对本文进行了讨论和总结。

1 SNLME模型及多元单边检验

1.1 SNLME模型

其中μp和ξiq分别是固定系数和随机系数的未知向量。基于hi(t)的假设，可以将ξiq看作独立同分布零均值随机向量的实现。考虑基于百分位结点的自然三次样条基，通过Akaike信息准则（AIC）或贝叶斯信息准则（BIC）来确定结点数（即最佳模型）。

用wp(t)和hiq(t)代替w(t)和hi(t)，则SNLME模型可以用下面的参数NLME模型来近似：

1.2 多元单边检验

2 具有缺失协变量的SNLME模型的多重插补

在HIV病毒动力学研究中，CD4细胞计数（协变量）往往不能完全观察到，这就会导致协变量存在测量误差和缺失，因此需要对协变量过程进行建模。本文考虑以下LME模型来描述协变量过程：

2.1 贝叶斯框架下的多重插补

设θ=(α，ψ)≡(α，β，δ，R，A，D)是模型（3）和（5）未知参数的集合，用a和b表示协变量模型（5）中的随机效应ai和响应模型（3）中的bi的集合。设f(·|·)和π(·)分别为条件密度函数和先验密度函数，并设(X|Y)表示给定Y的X的条件分布，则未观测到的“真”协变量值=Uijα+Vijai完全由平均参数α和随机效应ai决定，的多重插补可 (α，ai)从的多重插补中获得，可由以下给定观测数据的(α，a)的后验预测分布生成：

一般上式积分是高维的，计算十分复杂，因此，不能直接根据观测数据计算(α，a)的后验预测分布。作为替代方法，马尔可夫链蒙特卡罗方法（MCMC）方法可用于从表达式（6）中取样，且不需要计算积分，同时可以利用Gibbs抽样和Metropolis-Hasting（M-H）算法来实现MCMC方法。

为了避免计算表达式（6）中的积分，从(α，a)的后验预测分布中生成所需的(α，a)的估算，这里使用数据增强算法引入潜在变量(b，ψ)。具体来说，从初始值开始，通过迭代以下两步来创建 (α，a)的估算：对于t=1，2，…，

2.2 算法的实现

对于上一小节提出的贝叶斯框架下的多重插补方法，需要指定参数θ的先验分布。首先假设参数 α，β，δ，R，A，D 是独立的，即 π(θ)≡π(α，β，δ，R，A，D)=π(α)π(β)π(δ)π(R)π(A)π(D)，其中平均参数α和β服从正态先验分布，而精度参数δ，R，A和D服从伽马（Gamma）或威沙特（Wishart）先验分布，如下所示：

其中，超参数 γ，λ，p，q，r，Σα，Σβ，ΣR，ΣA，ΣD已知，Gamma 分布 G(γ，λ)和 Wishart分布 W(ΣR，p)的均值分别为γλ和pΣR。在实践中，先验信息可以从以前的研究或参考文献中获得。对于某些具有可靠先验信息的参数，可以使用方差较小的强先验；对于其他没有足够先验信息的参数，则可以采用非信息性先验分布（差异较大）。

在确定了观测数据的模型和模型中未知参数的先验分布后，可以使用数据增强算法和Gibbs抽样与M-H算法模拟多个独立的条件分布(α，a)。具体来说，使用Gibbs抽样，从条件分布f(α|D，a，b，ψ)，f(ai|yi，zi，bi，θ)，f(bi|yi，zi，ai，θ)，(i=1，…，n)和 f(ψ|D，α，a，b)中迭代样本。

用 β-表示参数 θ中除去 β 外其他的参数，δ-2，A-1，R-1，D-1同理，则精度参数 δ-2，A-1，R-1，D-1的全条件分布为：

3 具有缺失协变量的SNLME模型的多元单边假设检验

如第2节所讨论的，将SNLME模型通过3次样条基函数近似的参数NLME模型生成适当的多重插补之后，我们获得了m个完整的数据集，接下来将这m个结果进行组合便可产生总体推断。对于假设检验问题这里提供了两种方法来组合多个完整数据分析结果：1.可以组合多个完整数据检验统计量，即LRT统计量，然后提出总体检验；2.可以组合完整数据的参数估计，然后基于组合估计提出单个检验统计量，即Wald型检验统计量。

4 HIV实例分析

本节采用SNLME模型来建立HIV病毒动力学模型，经验多项式LME模型来建立协变量过程。将所提出的基于多重插补的检验TLR和TW应用于真实数据的多元单边检验问题中，并与朴素CC法进行比较。

这些数据来自包含45名HIV感染患者的抗HIV治疗研究，他们在研究开始时接受了一种有效的抗逆转录病毒方案的治疗，在开始治疗后第0，2，7，10，14，28天和第6，8，12周测量病毒载量（RNA），在整个研究中也记录了CD4细胞数目，每位患者的观察次数从3到9不等，一共观察了296次。其中，缺失数据出现在随时间变化的协变量CD4中，缺失率大概10%左右。图1显示了45名病人治疗期间病毒载量的轨迹，可以看到，不同患者之间的轨迹不同，但病毒载量的总体趋势都随着时间的推移在下降，这表明抗HIV的治疗可能是有效的。为了避免估计的不稳定性（过大或过小估计），本文标准化了CD4计数并重新定标原始时间（以天为单位），以使新的时间刻度在0和1之间。

图1 抗HIV治疗期间病人的病毒载量轨迹Fig.1 Viral load trajectory of patients during anti-HIV treatment

考虑以下HIV病毒动力学模型来描述病毒载量：

采用AIC或BIC准则选择最佳模型。使用自然三次样条基与百分位数结点的线性组合来逼近非参数函数w(t)和hi(t)，令ψ0(t)=ϕ0(t)≡1，取与（2）式中相同的三次样条基函数，其中q≤p。采用AIC/BIC准则由表1可确定模型（2）中p=3，q=1即：

表1 非参数函数w(t)和hi(t)的模型选择Table 1 Model selection for nonparametric functions w(t)and hi(t)

模 型 中 固 定 效 应 为 β=(P1，P2，λ1，λ2，β1，β2，β3)T，随 机 效 应 为 bi=(b1i，b2i，b3i，b4i)T，假 设eiji.i.d～N(0，σ2I)，bi～N(0，D)，随机误差和随机效应相互独立。表示真实的（但不可观察的）CD4计数。众所周知的是，CD4计数的测量通常有很大的误差。因此，这里假设初始病毒衰减率与真实的CD4值有关，而不是与观察到的CD4值有关。

对于协变量CD4，本文采用经验多项式LME模型模拟其轨迹：

在HIV病毒动力学研究中，研究人员经常感兴趣的是治疗的有效性，它可以通过真实衰减率是否严格为正来表明结果。模型中λ1，λ2分别表示生产性感染细胞，长期性感染细胞和潜伏性感染细胞的病毒衰减率，由于第一阶段为快速指数衰变阶段，而第二阶段衰变缓慢［17］，因此，在建立模型时，可以假设 P1＞P2，λ1＞λ2，以使模型具有可辨识性。

对于HIV病毒动力学模型（10），考虑多元单边检验问题（4），假设D是对角矩阵，则可以得到所有参数的集合为 θ=(P1，P2，λ1，λ2，β1，β2，β3，δ2，d11，d22，d33，d44)T，因此，矩阵R有以下形式：

其中矩阵R的秩r=2。

这样，针对病毒衰减率是否严格为正的检验问题可以写成如下形式：

接下来，使用组合LRT检验和总体Wald检验对参数进行检验，并得出结论。

为了解决协变量缺失问题，基于MCMC迭代，去掉初始的5万次迭代之后，再从接下来的15万次迭代中每500个MCMC抽样中选取一个样本，这样就得到了300个用于统计推断的样本。然后基于第四节的介绍，为协变量CD4生成m=5个多重插补数据集，并对这5个数据集进行参数估计。表2给出了基于估计结果得到的两个感兴趣参数的最大似然估计及均值，可以看到，正如之前假设的，患者第一阶段的病毒衰减率远大于第二阶段。

表2 五组多重插补结果的估计值及其平均值Table 2 Estimated and average values of five groups of multiple imputation

按照第四节给出的检验统计量计算可以得到组合LRT检验统计量为TLR=638.7，总体Wald检验统计量为TW=621.1，且两个检验的P值均接近于0。因此，有强有力的证据拒绝H0，从而得出结论：在抗HIV治疗期间，病毒载量显著下降。而基于CC方法的检验统计量为280.61，它的结果也拒绝了H0，但是比本文提出的检验统计量TLR和TW要弱。

5 结论

多重插补方法［18］的主要优点就是可以对输入数据集使用现有的完全数据方法，并且可以对输入数据集执行各种完全数据假设检验。本文基于多重插补方法对具有缺失协变量的半参数非线性混合效应模型的多元单边检验问题进行了讨论，给出了两个检验统计量，似然比检验（LRT）统计量和Wald检验统计量，并推导了其渐近零分布。并且为了验证检验统计量的性能，将给出的统计量应用到真实的HIV治疗研究中，通过计算我们得出了三个检验统计量的值，实例分析表明基于多重插补的似然比检验（LRT）统计量和Wald检验统计量给出的结果相似，而且比朴素的CC方法更提供了更强有力的证据。

在实践中，当缺失数据是信息性的或不可忽略的时候，必须假设缺失数据机制模型，并将缺失数据模型合并到联合模型中。同时，多元单边检验问题可以推广到其他混合效应模型，如广义线性模型、非参数混合效应模型、Logistic回归模型、比例风险模型等。这些问题超出了本文的范围，可以作为下一步的研究。

当然本文使用的SNLME模型不只适用于抗HIV病毒的研究，也可用于其他病毒治疗或药物研究中，如对结核病患者药物治疗的建模，糖尿病患者血浆中亮氨酸动力学的群体模型等。