非正态纵向高维数据下基于伪估计量的快速变量选择方法

本发明涉及非正态数据处理方法领域，具体是一种非正态纵向高维数据下基于伪估计量的快速变量选择方法。

背景技术：

1、大数据分析是数据科学重要研究方向，通过对海量数据进行分析，可充分挖掘出高价值的数据，能够为决策、结果预测等提供数据支持，特别是在临床医学、流行病学、社会学等领域具有重要意义。

2、临床医学、流行病学、社会学等领域中纵向观测数据的响应变量往往是非正态分布的，数据分析师们常使用广义线性混合模型对其进行拟合分析。广义线性混合模型由固定效应和随机效应两部分组成，综合了线性模型、线性混合模型和广义线性模型的特性，可以分析满足指数分布族的数据，同时有效处理纵向数据中的簇内异质性，具有广泛的应用。

3、而随着数据收集与存储技术的发展，广义线性混合模型的固定和随机效应中常出现高维预测因子。高维数据存在大量的冗余信息，因此需要有效且快速的变量选择方法，对数据进行一定程度的降维。高维随机效应的存在，会导致基于广义线性混合模型的变量选择法的目标函数涉及复杂的高维积分。由于存在这样的计算难度，传统的算法通常只考虑低维的随机效应，对高维固定效应进行变量选择。然而，改变一组效应的结构可能会对另一组效应的变量选择结果产生影响，因此，固定和随机效应中重要协变量的同时选择就变得尤为重要。

4、在正态数据的广义线性混合模型下，固定效应和随机效应中重要协变量的同时选择已经得到了深入研究。但是在非正态数据的广义线性混合模型下，由于似然函数中存在对随机效应的积分，随着纵向数据中的随机效应维度的增加，现有的基于蒙特卡洛、h似然、惩罚拟似然等正则化方法，在计算效率和变量选择的精确度方面存在很大挑战。

技术实现思路

1、本发明提供了一种非正态纵向高维数据下基于伪估计量的快速变量选择方法，以解决现有技术正侧化方法用于非正态数据的广义线性混合模型时，对于固定效应和随机效应中变量选择精度差、计算效率低的问题。

2、为了达到上述目的，本发明所采用的技术方案为：

3、非正态纵向高维数据下基于伪估计量的快速变量选择方法，包括以下步骤：

4、步骤1、获取样本数据，所述样本数据中包含n个观测个体，每个观测个体有mi个观测数据，每个观测数据在广义线性混合模型中的响应变量均为非正态分布，并记第i个观测个体的第j个观测数据的非正态分布响应变量为yij；

5、步骤2、设第i个观测个体的第j个观测数据在广义线性混合模型的固定效应中的高维协变量为xij，xij为p维向量；设第i个观测个体的第j个观测数据在广义线性混合模型的随机效应中的高维协变量为zij，zij为q维向量；建立广义线性混合模型中非正态分布响应变量yij的期望、方差及固定效应、随机效应中高维协变量的关系模型，如公式(1)、(2)所示：

6、

7、var(yij|bi)＝φv(μij)       (2)，

8、公式(1)、(2)中：

9、μij表示非正态分布响应变量yij的数学期望，记作μij＝e(yij|bi)；

10、var(yij|bi)表示非正态分布响应变量yij的方差；

11、g()和v()是分别作用于数学期望和方差的已知的连接函数；

12、表示固定效应中高维协变量xij的转置矩阵；

13、表示随机效应中高维协变量zij的转置矩阵；

14、φ为待估的离散度参数；

15、β为固定效应中的待估参数；

16、bi为随机效应中服从一个零均值的多元正态分布的随机项，bi的协方差矩阵d为随机效应中的待估参数，且d为q×q的对称矩阵；

17、步骤3、将第i个观测个体所有观测数据的非正态分布响应变量，按照次序排列成向量的形式yi，则有

18、引入伪估计量yi作为向量yi经线性变换后的近似，如公式(3)所示：

19、yi＝vi-1(yi-μi)+xiβ+zibi   (3)；

20、引入加权矩阵w如公式(4)所示：

21、

22、公式(3)-(4)中：μi为根据非正态分布响应变量yij的数学期望建立的向量，并有

23、xi为根据第i个观测个体所有观测数据在固定效应中的高维协变量，按照先后次序建立的mi×p维矩阵，并有

24、zi为根据第i个观测个体所有观测数据在随机效应中的高维协变量，按照先后次序建立的mi×q维矩阵，并有

25、vi-1为对角矩阵vi的逆矩阵，对角矩阵g-1表示连接函数g()的逆函数；(g-1)′表示逆函数的一阶导数；为第i个观测个体的第j个观测数据在固定效应中的高维协变量xij的转置矩阵；为第i个观测个体的第j个观测数据随机效应中的高维协变量zij的转置矩阵；diag{}表示对角矩阵，其对角元素为{}内的值；

26、步骤4、采用多次迭代的两步正则算法，在每次迭代中更新伪估计量的估计值，基于伪估计量新的估计值，进行正则化运算，多次迭代运算过程如下：

27、(4.1)首先，为固定效应中的待估参数β、随机效应中的待估参数d、待估的离散度参数φ、随机效应中的随机项bi分别设置初始估计值，分别记为这些初始值代入公式(3)，可得到伪估计量yi的初始值；

28、(4.2)然后，在第s步迭代中，s＝1,…,s，将固定效应中的待估参数β、随机效应中的随机项bi以及离散度参数φ在s-1步迭代中得到的估计值代入公式(3)和(4)，得到伪估计量yi的当前估计值yi(s)和加权矩阵w的估计矩阵w(s)，其中第0步迭代的估计值为步骤(4.1)中设置的初始值；

29、(4.3)接着，基于加权矩阵w的估计矩阵w(s)，计算权矩阵p(s)、q(s)，如公式(5)、(6)所示：

30、

31、q(s)＝w(s)-w(s)x(xtw(s)x)-1xtw(s)     (6)；

32、公式(5)、(6)中：(w(s))-1表示加权矩阵w的估计矩阵w(s)的逆矩阵；

33、表示在s-1步迭代中得到的随机效应中的待估参数d估计值；

34、x为样本的所有观测个体在固定效应中的高维协变量组成的n×p维的矩阵，并有x＝(x1t,…,xnt)t；

35、z为样本的所有观测个体在随机效应中的高维协变量组成的n×nq维的分块矩阵，并有z＝diag{z1,…,zn}；

36、

37、(4.4)基于伪估计量的当前估计值yi(s)和权矩阵p(s)、q(s)，建立施加惩罚项的目标函数如公式(7)、(8)所示：

38、

39、

40、公式(7)、(8)中：βj表示待估向量β中第j个元素；

41、b★k表示随机效应中第k个随机项的样本标准差，并有bik表示向量bi中第k个元素；

42、b为所有观测个体在随机效应中的随机项建立的向量，并有b＝(b1t,…,bnt)t；

43、表示第s-1步迭代中得到估计值的逆矩阵；

44、

45、表示目标函数h1中的已知惩罚函数，表示目标函数h2中的已知惩罚函数；λ1为惩罚函数的调优参数，λ2为惩罚函数的调优参数；

46、并对目标函数h1、h2最小化，以得到基于当前伪估计量的近似最大惩罚似然估计，作为此时固定效应中待估参数β的当前估计值随机效应中随机项bi的当前估计值

47、(4.5)得到固定效应中待估参数β的当前估计值随机效应中随机项bi的当前估计值后，根据公式(9)、(10)计算随机效应中待估参数d的当前估计值离散度参数φ的当前估计值公式(9)、(10)如下所示：

48、

49、

50、公式(9)、(10)中，上标t表示转置矩阵；

51、(4.6)重复步骤(4.2)-(4.5)，计算得到的固定效应中待估参数β当前估计值随机效应中随机项bi当前估计值取代各自对应的初始估计值，更新伪估计量和加权矩阵，由此进行多次迭代运算，直至两步正则算法收敛，得到固定效应中待估参数β的近似最大似然估计值、随机效应中随机项bi的近似最大似然估计值、随机效应中随机项的协方差矩阵d的近似最大似然估计值；

52、基于待估参数β、随机项bi、协方差矩阵d的近似最大似然估计值的非零元素所对应的变量，确定固定效应与随机效应中的重要协变量，认为以上参数近似最大似然估计值的零元素所对应的变量与响应变量无关并将其丢弃，由此实现对固定效应和随机效应进行变量选择。

53、进一步的，步骤2中公式(1)、(2)的连接函数g()和v()，分别是根据样本数据响应变量观测值的特征在广义线性混合模型中确定的作用于期望和方差的连接函数。

54、进一步的，步骤4中，固定效应中的待估参数β的初始估计值设为零向量，随机效应中待估参数d的初始估计值设为单位对角矩阵，离散度参数φ的初始估计值设为1，随机效应中,随机项bi的初始估计值设为零向量。

55、进一步的，步骤4中，惩罚函数采用scad惩罚函数。

56、进一步的，步骤4中，惩罚函数采用lasso惩罚函数。

57、进一步的，步骤4中，迭代的两步正则算法的收敛条件为迭代次数达到预设阈值。

58、进一步的，步骤4中，迭代的两步正则算法的收敛条件为最新一次迭代与其相邻的上一次迭代得到的待估参数β、d的近似最大似然估计值的差值的2范数小于预设阈值。

59、与现有技术相比，本发明的有益效果为：

60、本发明提供了一种基于伪估计量的快速变量选择方法，针对非正态响应变量及广义线性混合模型中的高维协变量，将广义线性混合模型中的估计问题与线性混合模型的估计问题关联起来，通过两步正则算法多次最小化估计迭代运算，对固定效应和随机效应进行变量选择。

61、相比限定低维随机效应对高维固定效应的变量选择、或是限定低维固定效应对高维随机效应的变量选择的两种单独选择法，本发明算法可以用来处理固定效应和随机效应中均存在高维协变量的情况。与现有正则化方法相比，由于算法中的两个目标函数分别是关于固定效应中待估参数β和随机效应中随机项bi的函数，两个目标函数的最小化过程完全独立，在固定效应和随机效应中仅估计一组效应，而不需要估计另一组效应，一组效应的变量选择无需限定另一组效应的维度，从而实现对两种效应的变量进行同时选择。这使得本发明方法与其他同时选择两组效应的方法相比，不仅便于计算机实现，计算效率更高，而且能够更具稳定性和计算精度。本发明无论是从数学原理还是融合算法实现都比较简单。

62、非正态响应变量的纵向数据临床医学、流行病学、经济学、社会学等应用领域中广泛存在，例如：具有相同疾病的不同病人的随着时间或疾病发展变化的包括生物标记物浓度、疾病状态等指标、不同城市随时间变化的空气污染物浓度、对不同产品每天记录的消费者投诉次数、学生成绩、股票收益、工资收入等等。而随着现代技术的发展，高维海量数据呈现爆炸式增长。例如，将基因表达数据、临床电子病历数据以及实验室化验数据作为预测因子，分析这些高维预测因子与病人疾病指标的关系，根据快速且有效的数据分析结果，从而辅助医生做出精准的疾病诊疗方案。又如，影响空气污染是因素是多方面的，包括政府政策、经济发展、城市基础设施建设、人口流动、气候气温湿度、地形地貌等等，快速且有效的数据分析结果，可以辅助施策者做出精准决策。

63、本发明基于以上实际应用目的，在用于分析非正态响应变量的纵向数据时，根据算法对参数β、d输出结果，若β的某个元素的值为零，则该位置对应的固定效应的协变量与响应变量无关、可被丢弃；若矩阵d的某个对角线元素为零，则对应的随机效应的协变量与响应变量无关、可被丢弃，从而实现对高维协变量数据进行快速且精准的降维。此外根据参数β的估计值、b预测值的非零元素的正负号及绝对值大小，可以判断协变量对响应变量的影响力。再结合离散度参数φ的估计，将这些参数估计值代入公式(1)-(2)，得到响应变量条件期望和方差的估计值，实现对响应变量的精准预测。