高维混合效应模型的统计推断：一种新的方法

Statistical Inference of High-Dimensional Mixed Effects Models: A New Approach

摘要: 线性混合效应模型广泛应用于分析聚类或重复测量数据。本文提出了一种拟似然结合弹性网的方法来估计高维线性混合模型中的未知参数，包括固定效应及随机效应的方差分量部分。在此基础上，也提出了相关的统计推断。所提出的方法适用于一般设置，其中随机效应的维度和簇可能很大。关于固定效应，我们提供的方法不依赖于方差分量的结构信息，即对方差分量中所涉及到的复杂未知参数使用代理矩阵进行简化。并且对所提出的方法在各种模拟设置中分别进行了固定效应的误差，假设检验的性能以及方差分量的估计误差的评估，均表现出较优的结果。

Abstract: The linear mixed-effects model is widely used for analyzing clustered or repeated measurements data. This paper proposes a pseudo-likelihood method combined with the elastic net to estimate unknown parameters in high-dimensional linear mixed models, including the variance components of both fixed and random effects. Furthermore, relevant statistical inferences are also presented. The proposed method is applicable to general settings where the dimension and clusters of random effects can be substantial. Regarding fixed effects, our approach does not rely on structural infor-mation about the variance components. Instead, it simplifies the complex unknown parameters in-volved in the variance components using surrogate matrices. The performance of the proposed method is evaluated for fixed effects errors, hypothesis testing, and variance component estimation errors in various simulation settings, all of which demonstrate superior results.

1. 引言

集群数据通常出现在许多领域，如生物学、遗传学和经济学。线性混合效应模型提供了一个灵活的工具来分析这样的集群数据，其中包括重复测量数据，纵向数据和多级数据等 [1] [2] 。线性混合效应模型包含固定和随机效应。在许多基因和经济研究中，协变量的维数可以很大，可能远远大于样本的大小。许多学者提出了各种各样的统计模型和方法，来研究和分析高维数据。然而，大多数方法都局限于处理独立观察的数据，如线性模型和广义线性模型。对高维线性混合效应模型进行统计推断仍然是一个具有挑战性的问题。在这项工作中，我们考虑高维混合效应模型中未知参数的估计和推断。

当维度固定时，很多方法被提出用来估计固定效应和方差项，如Gumedze和Dunne在2011年应用极大似然估计和限制性极大似然估计于混合效应模型的参数估计 [3] ，但遗憾的是，它们不适用于高维情况。Zhang和Tong在2007年提出了随机效应变系数模型的固定效应和方差参数的矩估计 [4] 。Peng和Lu (2012)考虑了固定维线性混合效应模型的矩估计 [5] 。他们都偏向于固定维设置。Ahmn，Zhang和Lu (2012)提出了另一种基于矩的方法 [6] ，该方法用于在固定维设置中估计和选择随机效应的方差分量。当各个集群大小相同时，此方法很适用。

对于固定维设置中方差分量的推断，可以采用似然比、得分和Wald检验等 [7] [8] [9] [10] (Stram & Lee, 1994; Lin 1997; Verbeke & Molenberghs, 2003; Demidenko, 2004)。但是这些方法均基于极大似然估计以及限制性极大似然估计，故也存在缺陷。

在高维中，Fan和Li (2012) [11] 研究了高维线性混合效应模型中，当聚类大小平衡时，固定效应和随机效应的选择无问题，选择变量的一致性要求关于固定效应和随机效应的最小信号强度条件。Bradic、Claeskens和Gueuning (2020) [12] 考虑在具有固定聚类大小、固定数量的随机效应和亚高斯设计的高维线性混合效应模型中测试固定效应的单个系数。但是他们的理论分析都要求随机效应的协方差矩阵具有正定性。并且在高维情况下，对方差分量的估计仍是未知的。基于此，Li等人 [13] 在2022年提出了一种针对于高维混合效应模型的拟似然参数估计方法，该方法基于Lasso并加入了拟似然来处理随机效应的影响。

而在本文中，我们开发了一种新的基于拟似然及弹性网的方法来推断高维线性混合效应模型中的未知参数。我们的方法适用于随机效应的数量可以是大的，集群大小可以是固定的或增长，平衡或不平衡的设置。该方法易于实现，每一步的优化都是解析的或凸的。基于真实协方差矩阵的代理，我们开发了一个惩罚准似然方法的固定效应估计。进一步开发了一个去偏估计的假设检验和建设的固定效应的置信区间。并且基于拟似然方法估计了方差分量。

本文的剩余部分组织如下：在第2节，我们介绍了混合效应模型且进行了相应的符号说明。第3节我们提出了基于弹性网的相关参数的拟似然估计方法。第4节我们展示了仿真实验结果，并将其进行了对比，最后，在第5节进行了总结和讨论。

2. 模型及符号

我们使用聚类数据的设置来呈现线性混合效应模型。对于重复测量数据，重复测量形成聚类。令 $i=1,\cdots ,n$ 表示类群标签。对于第i个类，响应向量为 $y_i\in R^{m_i}$ ，固定效应的设计矩阵为 $X^i\in R^{m_i\times p}$ ，随机效应的设计矩阵为 $Z_i\in R^{m_i\times q}$ 。从而线性混合效应模型 [14] 形式如下所示：

$y_i=X^i{\beta }^{\ast }+Z^i{\gamma }_i+{\epsilon }_i,i=1,\cdots ,n$ , (1)

其中， ${\beta }^{\ast }\in R^p$ 是固定效应的系数向量， ${\gamma }_i\in R^q$ 是第i个类的随机效应的系数向量，而 ${\epsilon }_i$ 是第i个类的噪声向量。对于 $i=1,\cdots ,n$ ，我们假设 ${\gamma }_i$ 和 ${\epsilon }_i$ 独立同分布，且均值为0，方差分别为 $\psi \in R^{q\times q}$ 和 ${\sigma }_e^2{I}_{m_i}$ 。定义 $N={\displaystyle {\sum }_{i=1}^n{m}_i}$ 表示总体样本量。我们将 ${\gamma }_i$ 和 ${\epsilon }_i$ 统称为随机部分。

本文，我们使用i表示第i个聚类，k表示每个聚类中的第k个观测值。 $y,\gamma ,\epsilon$ 和X分别由 $y_i,{\gamma }_i,{\epsilon }_i$ 和 $X^i$ 得到。令 $Z\in R^{N\times nq}$ 是对角块矩阵，第i个块为 $Z^i$ 。令 ${\Sigma }_{\theta }^i=Z^i\psi Z^i+{\sigma }_e^2{I}_{m_i}$ ，且 ${\Sigma }_{\theta }\in R^{N\times N}$ 是第i个块为 ${\Sigma }_{\theta }^i$ 的对角块矩阵。有 ${\Sigma }_z^i={\left(Z^i\right)}^{\text{T}}{Z}^i/m_i$ 和 ${\Sigma }_{z,x}^i={\left(Z^i\right)}^{\text{T}}{X}^i/m_i,i=1,\cdots ,n$ 。对于随机变量 $u\in R$ ，定义它的亚

高斯范数为 ${\Vert u\Vert }_{{\psi }_2}={\sup }_{1\le l}^{-1/2}{E}^{1/l}\left[{\left|u\right|}^l\right]$ 。我们称 ${\Vert u\Vert }_{{\psi }_2,Z}={\sup }_{1\le l}^{-1/2}{E}^{1/l}\left[{\left|u\right|}^l|Z\right]$ 为u的条件亚高斯范数。对于随机向量 $U\in R^{n_0}$ ，定义其亚高斯范数为 ${\Vert U\Vert }_{{\psi }_2}={\sup }_{{\Vert v\Vert }_2=1,v\in R^{n_0}}{\Vert \langle U,v\rangle \Vert }_{{\psi }_2}$ ，而相应的条件亚高斯范数为 ${\Vert U\Vert }_{{\psi }_2,Z}={\sup }_{{\Vert v\Vert }_2=1,v\in R^{n_0}}{\Vert \langle U,v\rangle \Vert }_{{\psi }_2,Z}$ 。

令 $A\in R^{n_0\times n_0}$ 是对称矩阵， $A\succcurlyeq 0$ 表示矩阵A非负定， $A\succ 0$ 表示正定， ${\Lambda }_{\max }\left(A\right)$ 和 ${\Lambda }_{\min }\left(A\right)$ 分别表示矩阵A的最大和最小特征值，令 ${\Vert A\Vert }_2$ 定义 ${\Lambda }_{\max }\left(A\right)$ ， ${\Vert A\Vert }_1={\max }_j{\displaystyle {\sum }_{i=1}^{n_0}\left|a_{i,j}\right|}$ ， ${\Vert A\Vert }_F=Tr\left(A^{\text{T}}A\right)$ ，其中 $Tr\left(A\right)$ 表示A的迹。令 $c,c_0,c_1,\cdots ,C,C_0,C_1,\cdots$ 表示在不同情况下可能变化的一些通用的正的常数。

3. 参数估计

引理1 (Lasso的收敛速度)：假设响应 $y_i$ 是由模型(1)生成的，X的每一行都是在以Z为条件下由协方差矩阵 ${\Sigma }_{x|z}$ 独立生成的。则对任意固定的 $j\in \left\{1,\cdots ,p\right\}$ ，有

$E\left[{\left|\frac{1}{N}{X}_{.,j}^{\text{T}}\left(Z\gamma +\epsilon \right)\right|}^2|Z\right]=\frac{{\left({\Sigma }_{x|z}\right)}_{j,j}{\sigma }_e^2}{N}+\frac{{\left({\Sigma }_{x|z}\right)}_{j,j}{\displaystyle {\sum }_{i=1}^n{m}_{i}Tr}\left(\Psi {\Sigma }_z^i\right)}{N^2}+\frac{{\displaystyle {\sum }_{i=1}^n{m}_i^2}{\Vert {\Psi }^{1/2}E\left[{\Sigma }_{z,x}^i|Z\right]\Vert }_2^2}{N^2}.$

考虑一般的弹性网作用于模型(1)，我们有：

${\beta }^{\left(lm\right)}=\arg {\min }_{b\in R^p}\left\{\frac{1}{2N}{\Vert y-Xb\Vert }_2^2+{\lambda }_2^{\left(lm\right)}{\Vert b\Vert }_2+{\lambda }_1^{\left(lm\right)}{\Vert b\Vert }_1\right\}.$ (2)

由引理1，当q或者 $m_i$ 增大且X和Z相关时，Lasso的收敛速度对于聚类数据不是最优，从而，弹性网的收敛速度也非最优。因此，我们需要考虑引入新的估计方法。

基于此，我们考虑引入一种新的拟似然方法，将随机部分的方差记为： ${\Sigma }_{{\theta }^{\ast }}^i=Z^i\psi Z^i+{\sigma }_e^2{I}_{m_i}$ ，该式包含较难估计的未知参数。因此，考虑 ${\Sigma }_{{\theta }^{\ast }}^i$ 的代理矩阵：

${\Sigma }_a^i=a{Z}^i{Z}^i+I_{m_i},$

其中， $a>0$ 是某预定常数。记 ${\Sigma }_a\in R^{N\times N}$ 为第i个块为 ${\Sigma }_a^i$ 的对角块矩阵。

引理2：若矩阵Ψ为正定矩阵，则 $\forall a>0$ ，

$\min \left\{\frac{1}{{\sigma }_e^2},\frac{a}{{\Lambda }_{\max }\left(\Psi \right)}\right\}{\Sigma }_a^{-1}\preccurlyeq {\Sigma }_{{\theta }^{\ast }}^{-1}\preccurlyeq \max \left\{\frac{1}{{\sigma }_e^2},\frac{a}{{\Lambda }_{\min }\left(\Psi \right)}\right\}{\Sigma }_a^{-1},$

因此，当Ψ的特征值有界且为正数时， ${\Sigma }_a^{-1}$ 和 ${\Sigma }_{{\theta }^{\ast }}^{-1}$ 有相同的收敛速度。

接下来，提出具体的拟似然估计方法。使用 ${\Sigma }_a$ 代替 ${\Sigma }_{{\theta }^{\ast }}$ ，并且对响应向量y和固定效应矩阵X进行标准化处理，得到：

$\left(X_a,y_a\right)=\left({\Sigma }_a^{-1/2}X,{\Sigma }_a^{-1/2}y\right).$

3.1. 固定效应的估计

首先基于标准化后的数据 $\left(X_a,y_a\right)$ 对固定效应进行估计，对固定的a，定义：

$\hat{\beta }={\text{argmmin}}_{\beta \in R^p}\left\{\frac{1}{2Tr\left({\Sigma }_a^{-1}\right)}{\Vert y_a-X_a\beta \Vert }_2^2+{\lambda }_2{\Vert \beta \Vert }_2+{\lambda }_1{\Vert \beta \Vert }_1\right\},$ (3)

其中， ${\lambda }_l,l=1,2$ 是调优参数， $Tr\left({\Sigma }_a^{-1}\right)$ 为有效样本量，可以看作是Lasso和岭回归的线性组合。根据Zou等人(2017年) [15] ，我们定义 $\lambda ={\lambda }_1+{\lambda }_2$ ， $\alpha =\frac{{\lambda }_2}{{\lambda }_1}$ ，从而(3)式变为，

$\hat{\beta }={\text{argmmin}}_{\beta \in R^p}\left\{\frac{1}{2Tr\left({\Sigma }_a^{-1}\right)}{\Vert y_a-X_a\beta \Vert }_2^2+\lambda \left(\alpha {\Vert \beta \Vert }_2+\left(1-\alpha \right){\Vert \beta \Vert }_1\right)\right\}.$ (4)

进一步，为了对 ${\beta }^{\ast }$ 进行统计推断，定义如下的无偏估计，

${\hat{\beta }}_j^{\left(db\right)}={\hat{\beta }}_j+\frac{{\hat{\omega }}_j^{\text{T}}\left(y_a-X_a\hat{\beta }\right)}{{\hat{\omega }}_j^{\text{T}}{\left(X_a\right)}_{.,j}},$ (5)

定义 ${\hat{\omega }}_j={\left(X_a\right)}_{.,j}-{\left(X_a\right)}_{.,-j}{\hat{k}}_j$ ，其中，

${\hat{k}}_j=\arg {\min }_{{\hat{k}}_j\in R^{p-1}}\left\{\frac{1}{2Tr\left({\Sigma }_a^{-1}\right)}{\Vert {\left(X_a\right)}_{.,j}-{\left(X_a\right)}_{.,-j}{\hat{k}}_j\Vert }_2^2+{\lambda }_{2j}{\Vert k_j\Vert }_2+{\lambda }_{1j}{\Vert k_j\Vert }_1\right\}$ . (6)

从而， ${\beta }_j^{\ast }$ 双边置信区间为：

${\hat{\beta }}_j^{\left(db\right)}\pm z_{\alpha /2}\sqrt{{\hat{V}}_j},$

这里， $z_{\tau }$ 为标准正态分布的第 $\tau$ 分位数，而 ${\hat{V}}_j$ 为 ${\hat{\beta }}_j^{\left(db\right)}$ 方差的一个估计值，有，

${\hat{V}}_j=\frac{{\displaystyle {\sum }_{i=1}^n{\left[{\left({\hat{\omega }}_j^i\right)}^{\text{T}}\left(y_a^i-X_a^i\hat{\beta }\right)\right]}^2}}{{\left({\hat{\omega }}_j^{\text{T}}{\left(X_a\right)}_{.,j}\right)}^2}.$ (7)

3.2. 方差分量的估计

接下来，我们对模型中的方差分量进行估计。考虑到Ψ是一个对称矩阵，故有如下分解：

$\Psi ={\Psi }_{{\eta }^{\ast }}={\displaystyle {\sum }_{j=1}^d{\eta }_j^{\ast }{G}_j},$ (8)

其中， $G_1,\cdots ,G_d$ 是线性无关的对称基矩阵，且满足：

${\displaystyle {\sum }_{j=1}^d{c}_j{G}_j}=0,\text{iff}{c}_1=\cdots =c_d=0.$ (9)

用于估计方差分量的方法是高斯最大似然方法。将数据分为两部分， $I_1\cup I_2=\left[n\right]$ ， $I_1\cap I_2=\varnothing$ ，且 $\left|I_1\right|\approx \left|I_2\right|\approx n/2$ 。令 ${\hat{\beta }}^{\left(2\right)}$ 是数据 ${\left\{X^i,Z^i,y_i\right\}}_{i\in I_2}$ 估计得到的初始值，接下来在 $i\in I_1$ 中计算残差 ${\hat{r}}_i=y_i-X^i{\hat{\beta }}^{\left(2\right)}$ ，并且由下式得到 ${\sigma }_e^2$ 的估计，

${\hat{\sigma }}_e^2=\frac{1}{{\displaystyle {\sum }_{i\in I_1}Tr}\left(P_{Z^i}^{\perp }\right)}{\displaystyle {\sum }_{i\in I_1}{\Vert P_{Z^i}^{\perp }{\hat{r}}_i\Vert }_2^2}.$ (10)

现在，我们估计 ${\eta }^{\ast }$ ，

$\hat{\eta }=\arg {\min }_{\eta \in R^d}{\displaystyle {\sum }_{i\in I_1}{\Vert {\left({\Sigma }_a^i\right)}^{-1/2}\left({\hat{r}}_i{\hat{r}}_i^{\text{T}}-Z^i{\Psi }_{\eta }{\left(Z^i\right)}^{\text{T}}-{\hat{\sigma }}_e^2{I}_{n_i}\right){\left({\Sigma }_a^i\right)}^{-1/2}\Vert }_F^2},$ (11)

其中 $K\ge K_2$ 是常数，且 ${\hat{\sigma }}_e^2$ 由(10)式得到。

公式(10)的基本原理是观测值 $P_{Z^i}^{\perp }\left(y_i-X^i{\beta }^{\ast }\right)$ 的协方差矩阵为 ${\sigma }_e^2{P}_{Z^i}^{\perp }$ ，它只涉及目标参数 ${\sigma }_e^2$ 。将 ${\beta }^{\ast }$ 用它的拟似然估计来代替，得到(10)。此估计量仅当 ${\displaystyle {\sum }_{i\in I_1}Tr}\left(P_{Z^i}^{\perp }\right)>0$ 时才成立，即 ${\displaystyle {\sum }_{i\in I_1}{m}_i\max \left\{0,1-q/m_i\right\}}>0$ 。

3.3. 去偏估计量的渐近性质

条件1 (亚高斯随机变量)：随机噪声项 ${\epsilon }_{i,k},i=1,\cdots ,n;k=1,\cdots ,m$ 独立分布于均值为0方差为 $0<{\sigma }_e^2<K_0<\infty$ 。 ${\epsilon }_{i,k}$ 的亚高斯范数以K₀为上界。随机效应 ${\gamma }_i\in R^q,i=1,\cdots ,n$ 独立分布于均值为0，方差为 $\Psi \le K_1{I}_q$ ，其中 $K_1>0$ 是常数。对于 $i=1,\cdots ,n$ ， ${\epsilon }_i$ 和 ${\gamma }_i$ 独立于彼此，且独立于 $\left(X^i,Z^i\right)$ 。 ${\Sigma }_{{\theta }^{\ast }}^{-1/2}\left(Z\gamma +\epsilon \right)$ 的亚高斯范数以K₀为界。

条件2：在给定Z的条件下，X的每一行都是独立的，具有零均值和协方差矩阵 ${\Sigma }_{X|Z}$ 满足 $0<K_{\ast }\le {\Lambda }_{\min }\left({\Sigma }_{X|Z}\right)\le {\Lambda }_{\max }\left({\Sigma }_{X|Z}\right)\le K^{\ast }<\infty$ 。在以Z为条件时， $X_{k,.}^i$ 的条件亚高斯范数以K₀为上界。

定理1 (无偏估计量的渐近性质)：假设条件1和条件2成立。令 ${\lambda }_l\wedge {\lambda }_{lj}\ge c_1\sqrt{\log p/Tr\left({\Sigma }_a^{-1}\right)}$ ， $l=1,2$ ， $j=1,2,\cdots ,p-1$ ， $c_1$ 是一个足够大的常数。对于在(5)式中定义的 ${\hat{\beta }}_j^{\left(db\right)}$ ，若 ${\left(s\log p\right)}^2\vee \log n{\max }_i{m}_i\ll Tr\left({\Sigma }_a^{-1}\right){\Lambda }_{\min }\left({\Sigma }_a^{-1/2}{\Sigma }_{{\theta }^{\ast }}{\Sigma }_a^{-1/2}\right)$ 以及 $\left|H_j\right|\log p\ll Tr\left({\Sigma }_a^{-1}\right)$ ，则，

$V_j^{-1/2}\left({\hat{\beta }}_j^{\left(db\right)}-{\beta }_j^{\ast }\right)=R_j+o_p\left(1\right),$ (12)

其中， $R_j\stackrel{D}{\to }N\left(0,1\right)$ ，对

$V_j=\frac{{\hat{\omega }}_j^T{\Sigma }_a^{-1/2}{\Sigma }_{{\theta }^{*}}{\Sigma }_a^{-1/2}{\hat{\omega }}_j}{{\left\{{\hat{\omega }}_j^{\text{T}}{\left(X_a\right)}_{.,j}\right\}}^2}.$

这里， $V_j$ 的大小满足，

$V_j=\frac{{\left({\Sigma }_{X|Z}^{-1}\right)}_{j,j}Tr\left({\Sigma }_a^{-1}{\Sigma }_{{\theta }^{\ast }}{\Sigma }_a^{-1}\right)}{T{r}^2\left({\Sigma }_a^{-1}\right)}\left(1+o_p\left(1\right)\right).$

4. 仿真实验

在本节中，我们进行了仿真模拟来评估所提出的方法的power性能，并将其与相关的方法进行比较。我们考虑和Li等人(2022)相同的模拟设置：令 $N=144$ ， $p=300$ 。 $\left(X,Z\right)$ 的每一行都独立同分布产生于均值为0，方差如 ${\Sigma }_x=I_p$ ， ${\Sigma }_z=I_q$ 以及 ${\left({\Sigma }_{X,Z}\right)}_{k,j}={\rho }^j,1\le j,k\le q$ ， ${\left({\Sigma }_{X,Z}\right)}_{k,.}=0,k>q$ 。也就是说，如果 $j\le p$ ，则 $X_j$ 和Z之间的相关性为非0，如果 $j>q$ ，则相关性为0。随机噪声项独立同分布于 ${\epsilon }_i~N\left(0,0.25I_{m_i}\right)$ 。我们考虑 $q\in \left\{2,8,14\right\}$ ，响应y经由模型(1)生成，其中 $s=5$ 且 ${\beta }_{1:5}=\left(1,0.5,0.2,0.1,0.05\right)$ 和相等的簇大小，即 $m_1=m_2=\cdots =m_n=m$ 。每个设置都使用300次独立的Monte Carlo模拟进行重复。

4.1. 固定效应的统计推断

我们首先检查所提出的估计量的误差。考虑两种随机效应的协方差矩阵形式，一个是正定的Ψ，我们记为“p.d.Ψ”，其中， ${\Psi }_{j,k}={0.56}^{\left|j-k\right|},1\le j,k\le q$ ，以及奇异的Ψ，记为“singularΨ”，其中Ψ是对角矩阵，且 ${\Psi }_{j,j}=0.56,1\le j\le q/2$ ，否则 ${\Psi }_{j,j}=0$ 。对于所提出的方法，首先通过交叉验证选择a，且调优参数 ${\lambda }_i$ 通过 ${\hat{\sigma }}_x\sqrt{2\log p/N}$ 计算，其中， ${\hat{\sigma }}_x$ 由scaled-Lasso和观测值 $\left({\left(X_a\right)}_{.,j},{\left(X_a\right)}_{.,-j}\right)$ 进行计算。

从表1可以看出，我们所提出的估计量具有较小的估计误差，且随着m和q的增加，误差变化不大，说明我们提出的方法具有较好的稳定性。同时，随着 ${\beta }_j^{\ast }$ 的值递减，误差也呈现递减的规律。

Table 1. Standard error for p.d.Ψ and singularΨ when β j ∗ ∈ { 1 , 0.5 , 0.2 , 0 }

表1. ${\beta }_j^{\ast }\in \left\{1,0.5,0.2,0\right\}$ 时在正定Ψ以及奇异Ψ下的标准误

接下来，我们检查基于 ${\hat{\beta }}_j^{\left(db\right)}$ 的假设检验的power性能。考虑两个随机效应的协方差矩阵，第一个是“p.d.Ψ”，另一个为“singularΨ”，其余参数设置和之前相同。

在表2中，我们展示了本文提出的方法和Li等人(2023)提出的方法的I类误差和power。我们的运行时间和Li等人的时间相差不大。在理想情况下，原假设H₀下的拒绝率应该接近5%，并且对{1, 0.5, 0.2}的拒绝率应该大于5%。从中可以看到我们的方法和Li等人的方法在控制第I类误差方面都是有效的。然而，在保证第I类错误的情况下，我们的方法相对于Li等人的power性能较好。

4.2. 方差分量的统计推断

真实的固定效应和数据生成步骤与4.1节相同。我们使用全部数据来估计 ${\sigma }_e^2$ 和 ${\eta }^{\ast }$ 。令 ${\sigma }_e^2=0.25$ ，我们考虑 $d=2$ 的对角矩阵Ψ，基矩阵的设置如下所示：

$G_1=\left(\begin{array}{cc}{I}_{q/2}& 0\\ 0& 0\end{array}\right),G_2=\left(\begin{array}{cc}0& 0\\ 0& I_{q/2}\end{array}\right),$

对于正定矩阵Ψ， ${\eta }^{\ast }={\left(0.56,0.56\right)}^{\text{T}}$ ，而奇异矩阵Ψ，有 ${\eta }^{\ast }={\left(0.56,0\right)}^{\text{T}}$ 。表3显示了 ${\sigma }_e^2\left(\text{mae}.{\sigma }_e^2\right)$ 、 ${\eta }_1\left(\text{mae}.{\eta }_1\right)$ 和 ${\eta }_2\left(\text{mae}.{\eta }_2\right)$ 的平均绝对误差。我们可以看到，估计的方差分量的平均绝对误差保持在较小的水平，且在Ψ是正定的情况下，估计的平均绝度误差小于奇异情况下的。

Table 2. The rejection rate for testing H0: β j ∗ = 0 at 95% level for β j * ∈ { 1 , 0.5 , 0.2 , 0 } with positive definite (p.d.) and singular when ρ = 0

表2. 测试H₀： ${\beta }_j^{\ast }=0$ 在95%水平上的power， ${\beta }_j^{*}\in \left\{1,0.5,0.2,0\right\}$ ，且Ψ具有正定矩阵(p.d)以及 $\rho =0$ 时的奇异矩阵(singular)两种形式

Table 3. Estimation of the variance components with the proposed method for positive definite and singular when ρ = 0

表3. 当 $\rho =0$ 时，利用所提出的方法对方差分量进行了正定和奇异估计

5. 总结与讨论

本文我们基于高维线性混合模型的框架，考虑了未知参数的估计和推断问题。并且，提出了一种新的结合拟似然和弹性网的方法，该方法考虑使用拟似然消去未知参数的影响，并且通过弹性网进行变量选择及降维。在建模重复测量和纵向数据，特别是当集群大小是大的或异构的时候，具有普遍适用性。我们提出的估计过程计算效率高，不需要很强的对于随机效应和误差的分布假设。仿真实验表明，我们的方法在假设检验方面的性能优于之前的方法，并且，估计量的误差也处于合理的范围内。