混合效应模型在葡萄糖耐糖量测试中的应用研究

Application of Mixed Effect Model in Glucose Tolerance Test

摘要: 本文是对一组来自标准葡萄糖耐糖量实验的平衡纵向数据构建线性混合效应模型、lasso回归和广义线性混合模型，利用AIC准则以及lasso的变量选择选取最优的模型，并对所得模型进行模拟预测，然后计算模型的均方误差并进行比较，选出均方误差相对较小的模型。发现线性混合模型在这组数据中具有较好的预测效果。本文所有计算均用R软件完成。

Abstract: This article is about a group of glucose sugar levels in experimental equilibrium longitudinal resistance from standard data to construct linear mixed effects models, lasso regression and generalized linear mixed models, using the AIC criterion as well as the lasso variable selection to select the optimal model and forecast the model simulated, then the mean square error of the model is calculated and compared, and the model with relatively small mean square error is selected. It is found that the linear mixed model has a good predictive effect in this set of data. All calculations in this paper are completed by R software.

1. 引言

在标准葡萄糖耐糖量实验中 [1]，得到一组平衡纵向数据。我们的研究目的是：评估一下，在这个实验中对照组和肥胖组是否有显著差异。初步的探索表明构建合适的混合效应模型是分析此类数据的关键。

针对我们的数据，我们首先基于探索性数据分析，给出了模型的基本结构 [2]，并根据AIC，BIC准则选择合适的协方差结构和固定效应 [3] [4]，选出最佳模型。我们还用lasso方法进行变量选择，看这两种方法选择的变量是否一致。参数模型的缺点是对分布假设高度敏感，如果分布假设不能满足，那么结论会相差很多，甚至得到完全错误的结论，所以我们进一步用广义线性混合效应模型来拟合这组数据。最后，对三个模型进行400次模拟预测，比较一下哪种方法可以更好地拟合该数据。

2. 模型理论简介

2.1. 线性混合效应模型

线性混合效应模型是Laird和Ware于1982年提出的，它的一般形式为：

$\{\begin{array}{l}{Y}_i=X_i\beta +Z_i{b}_i+{\epsilon }_i\\ b_i~N_q\left(0,D\right)\\ {\epsilon }_i~N_{n_i}\left(0,R_i\right)\\ b_1,\cdots ,b_m,{\epsilon }_1,\cdots ,{\epsilon }_m独立\end{array}$ (2.1)

这里的 $i=1,2,\cdots ,m$， $X_i$ 和 $Z_i$ 是已知的设计矩阵， $\beta$ 是 $p\times 1$ 维固定效应， $b_i$ 是个体随机效应向量。 ${\epsilon }_i$ 为随机误差向量。一般而言， $Z_i$ 的列是 $X_i$ 的子集。而 $Y_i$ 是受试者响应向量，它服从以下这样一个多元正态分布：

$\{\begin{array}{l}E\left[Y_i\right]=X_i\beta \\ Cov\left[Y_i\right]=R_i+Z_{i}G{Z}_i^{\text{T}}={\Sigma }_i\\ Y_i~N_{n_i}\left(X_i\beta ,{\Sigma }_i\right)\end{array}$ (2.2)

向量Y代表完全观察到的数据集合，即 $Y={\left(Y_1,\cdots ,Y_m\right)}^{\text{T}}$。因为Y是一个多元正态随机向量的线性组合，故Y的分布如下：

$Y~N_n\left(X\beta ,\Sigma \right)$

上式中， $N={\displaystyle {\sum }_{i=1}^m{n}_i}$， $X={\left(X_1,X_2,\cdots ,X_m\right)}^{\text{T}}$， $\Sigma =diag\left({\Sigma }_1,{\Sigma }_2,\cdots ,{\Sigma }_m\right)$， $\Sigma$ 的这种形式假定了不同个体之间的所有观察结果都是独立的。

2.2. lasso方法

lasso (Least absolute shrinkage and selection operator)是由Tibshirant于1996年提出的 [5]，该方法可以压缩模型的系数，并将一些系数置为零。该方法的主要特点是可以同时进行参数估计和变量选择。

假定数据为 $\left(x^i,y_i\right),i=1,\cdots ,N$，其中， $x^i={\left(x_{i1},\cdots ,x_{ip}\right)}^{\text{T}}$ 是预测变量， $y_i$ 是响应。在通常的回归中，我们假设观测值是独立的，或者给定 $x_{ij}$ 时， $y_i$ 是有条件独立的。我们假设 $x_{ij}$ 是标准化的，所以 ${\displaystyle {\sum }_i{x}_{ij}}/N=0$， ${\displaystyle {\sum }_i{x}_{ij}^2}/N=1$。

令 $\hat{\beta }={\left({\hat{\beta }}_1,\cdots ,{\hat{\beta }}_p\right)}^{\text{T}}$，lasso估计 $\left(\hat{\alpha },\hat{\beta }\right)$ 定义为：

$\left(\hat{\alpha },\hat{\beta }\right)=\arg \min \left\{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(y_i-\alpha -{\displaystyle \underset{j}{\sum }{\beta }_j{x}_{ij}}\right)}^2}\right\}$，约束为 ${\displaystyle \underset{j}{\sum }\left|{\beta }_j\right|\le t}$

其中， $t\ge 0$，是调整参数，对于所有的t， $\alpha$ 的所有解 $\hat{\alpha }=\bar{y}$，不失一般性，假设 $\bar{y}=0$，因此可以省略 $\alpha$。

2.3. 广义线性混合效应模型

广义线性混合模型(GLMM)的模型结合了广义线性模型和线性混合模型的优点，进一步扩展了广义线性混合模型，包括随机变量，以说明纵向或聚类数据的变化和相关性。另一方面，GLMM也是对线性混合模型的扩展，允许非正态分布的响应数据。使用GLMM的一个优点是响应数据不必是正态分布的数据。GLMM利用数据的真实本质，而不是依赖于数据的正态分布或转换假设。

广义线性混合效应模型由三部分组成：

(1) 随机组成：

1) 响应； $Y_{ij}|b_i~\exp \left\{\frac{y\theta -b\left(\theta \right)}{a\left(\varphi \right)}+c\left(y,\theta \right)\right\}$

2) 随机效应： $b_i~N_q\left(0,G\right)$

(2) 系统组成： ${\eta }_{ij}=X_{ij}^{\text{T}}\beta +Z_{ij}^{\text{T}}{b}_i$

(3) 链接函数： $g\left\{E\left[y_{ij}|b_i\right]\right\}=X_{ij}^{\text{T}}\beta +Z_{ij}^{\text{T}}{b}_i\left(i=1,\cdots ,m;j=1,\cdots ,n_i\right)$

其中， $y_{ij}$ 是第i个人在第j个时间点的响应，而 $x_{ij}$ 和 $z_{ij}$ 分别是与组合和随机变量相关的协变量向量。 $\beta$ 和 $b_i$ 分别是固定效应和随机效应。

3. 数据说明

本文数据来源于《Randomization Analysis of the Completely Randomized Design Extended to Growth and Response Curves》这篇文献。在一项有关高血糖和相对高胰岛素血症的研究中，对科罗拉多大学医学中心的儿科临床研究病房的13名对照组和20名肥胖组患者进行了标准的葡萄糖耐糖量试验。下表记录了在他们口服标准剂量葡萄糖后，分别在0，0.5，1，1.5，2，3，4，5小时取其血液样本，测其血浆无机磷值得到的，每个病人都测量了8次 [6]。其中数据包括了编号(subject)、组别(treatment)、测量时间(time)、及每个时刻的血浆无机磷测量值(plasma inorganic phosphate measurement)。表1给出了部分病人的数据。

4. 数据分析

4.1. 线性混合效应模型

本文首先运用线性混合效应模型来拟合数据，由于线性混合效应模型假定数据是服从正态分布的，所以分析数据时先进行正态性检验，原始测量值数据直方图和QQ图如图1所示。

Table 1. Data of some patients

表1. 部分病人数据

注：treatment为0表示对照组，1表示肥胖组。

Figure 1. Data histogram and QQ chart of plasma inorganic phosphorus

图1. 血浆无机磷测量值数据直方图和QQ图

为了研究不同组别对血浆无机磷测量值的结果是否有所不同，文章中画出了0，0.5，1，1.5，2，3，4，5小时这8个时间点对应的响应变量(measurement)均值折线图如图2所示 [7]。

从图象中可以看出来测量值在2小时前后变化很大，不同组别有显著的差异，并且测量值对于时间的变化规律不是线性的，在模型中可能需要加入时间的高次项。基于以上分析，我们建立如下线性混合效应模型：

$y_{ij}={\beta }_0+{\beta }_1\cdot {\text{time}}_{ij}+{\beta }_2\cdot {\text{time}}_{ij}^2+{\beta }_3\cdot {\text{time}}_{ij}^3+{\beta }_4\cdot \text{treatment}+{\beta }_5\cdot \text{treatment}\cdot \text{time}+{\epsilon }_{ij}$

Figure 2. Mean line graph under different groups

图2. 不同组别下的均值折线图

首先，对上述模型选取不同的协方差结构，基于AIC，BIC准则，我们选择个体内部的协方差结构为AR(1)。这与数据的探索性分析中，与散点图得出的结果一致。然后，我们选取最佳的随机效应，根据AIC，BIC准则，当选择时间(time)为随机效应时模型更优。最后，基于我们选择的个体内部的协方差结构和随机效应的结构可以选择固定效应，固定效应包括：treatment，时间的三次多项式，时间与组别的交互项。根据模型的基本选取准则，接着我们对模型进行参数估计，由于经典最大似然估计方法估计出来的方差是有偏估计，因此我们这里使用限制最大似然估计(REML)方法来获得方差分量的最大似然估计(见表2)。限制最大似然估计的想法是将数据分为两部分，一部分用来估计固定效应，另一部分用来估计协方差参数 [8] [9]。可得模型及其限制极大似然估计为：

$\begin{array}{c}y=3.947508-1.313492\ast \text{time}+0.467245\ast \text{time}2-0.041090\ast \text{time}3\\ +0.667951\ast \text{treatment}-0.114951\ast \text{time}\ast \text{treatment}\end{array}$

Table 2. Restricted maximum likelihood (REML) estimation for linear mixed models

表2. 线性混合模型的限制最大似然(REML)估计

4.2. lasso方法

用R软件里的lars程序包实现lasso变量选择。函数lars()，提供了通过回归变量x和因变量y求解其回归模型。

在对数据用线性混合效应建立模型的基础上，我们用lasso进行模型回归及其变量选择 [10]。如图3表示在进行lasso回归时，自变量被选入的顺序。

可以看到图3中的竖线对应的是lasso中迭代的次数，对应的系数值不为0的自变量即为选入的，竖线的标号与lasso回归中的step相对应。

根据MallowsCp统计量，即Cp值选择Cp值最小的模型及模型对应系数，见表3。

4.3. 广义线性混合效应模型

取链接函数为密度函数，建立如下模型并采用极大似然(ML)估计方法来进行参数估计(见表4)。

${\mu }_{ij}={\beta }_0+{\beta }_1\cdot {\text{time}}_{ij}+{\beta }_2\cdot \text{treatment}+{\beta }_3{\text{time}}_{ij}\cdot \text{treatment}+b_i$

Figure 3. Lasso regression

图3. lasso回归

Table 3. The corresponding coefficient of lasso regression

表3. lasso回归对应系数

Table 4. Maximum likelihood estimation of generalized linear mixed effects model

表4. 广义线性混合效应模型的最大似然估计

通过以上三种方法对该组数据进行拟合，发现每个系数都比较显著，没有剔除的项，三种方法都表明组别(treatment)的系数较大，说明肥胖对血浆的无机磷酸盐含量有较大影响。就线性混合效应模型与lasso方法来看，两者之间的变量选择结果一致，而且参数估计的结果也比较接近。

5. 模型预测与比较

接下来我们用所建立的模型进行模拟预测，首先生成与原数据同分布的数据，然后再用所得模型进行预测，我们用均方误差来评价模型的预测能力。均方误差的定义为：

$\text{MSE}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{\left({\text{observed}}_i-{\text{predicted}}_i\right)}^2}$

其中，N为样本容量。当预测值与真实值完全相同时为0，误差越大，MSE值越大。如图4给出了三个模型的均方误差

从图4可以看出，黑色的水平线代表了线性混合模型的均方误差的平均值，要低于其他两种模型的均方误差的平均值，而且线性混合模型整体的均方误差也低于其他两种的均方误差，所以在该组数据下我们采用线性混合模型来拟合会更好，这是因为这组数据本身接近正态分布，有着很好的性质。在实际应用中，我们要基于数据本身选择更优的模型以便更好的解决问题。

Figure 4. Mean square error chart

图4. 均方误差图

6. 结论

基于对数据的探索性分析，如果数据接近正态分布，我们可以通过线性混合效应模型来对数据进行分析和预测。如果数据不是正态分布，可以对数据进行转换，然后用线性混合效应模型分析数据。在本文中，线性混合效应模型的拟合精度要比lasso回归和广义线性混合模型的拟合精度要高。我们发现肥胖对血浆中的无机磷酸盐含量有显著的影响，进而会导致高血压、糖尿病等多种疾病，所以，对于在儿童中由肥胖引起的血浆中的无机磷酸盐含量变化研究是有必要的。

NOTES

^*通讯作者。