纵向数据是通过对每个个体在不同的时间点上进行观测得到的数据集,广泛应用在经济学、医学等领域,具有组内数据相关、组间数据独立的特点。线性混合效应模型[1-5]能够很好地刻画数据间的组内相关性和组间独立性,因此是分析处理纵向数据的常用模型之一。
随着信息、网络技术的迅速发展,数据增长变化速度惊人,呈现出数据海量化趋势[6]。传统的线性混合模型估计方法在海量数据下遇到一系列的挑战,如存储瓶颈、计算效率等问题[7-8],因此有必要探求新的算法对以往线性混合效应模型的估计方法进行改进。
为了构造线性混合效应模型中系数的最优估计,Wu等[5]提出了基于协方差结构的加权最小二乘法,但该估计方法的效率与纵向数据的协方差结构估计方法密切相关。常用的协方差结构估计方法有极大似然法和限制极大似然法等,然而这些方法需要进行反复的迭代计算及优化步骤,当数据量非常大时会导致很高的计算时间和计算量成本。为了解决此类问题,Sun等[9]在研究随机效应变系数模型时提出了一种方差分量估计的简便方法,该方法不需要进行迭代计算,能够适应海量数据情形。
然而直接利用矩阵形式的加权最小二乘法估计模型系数会遇到存储问题,因此分治算法被提出并在海量数据的统计计算中得到了广泛应用[10-13]。分治算法的核心思想是先把复杂问题分解成几个子问题,找到这几个子问题的解法后,再利用合适的方法把每个子问题的结果组合起来得到整个问题的结果。Chang等[11]提出的分治算法是解决海量数据分析问题的最实用方法;Guha等[12]把最小二乘法与分治算法结合起来解决海量数据下线性回归模型的估计问题:首先将海量数据集分成一系列可管理的数据块,然后对每个数据块独立运行最小二乘方法,最后整合每块数据的结果得到最终的结果。但之前学者们也只是利用分治算法解决海量数据下简单数据模型的估计问题,对更复杂的、应用更广泛的线性混合效应模型在海量数据下的估计算法还并未研究。
本文基于文献[5]的线性混合效应模型系数的加权最小二乘估计方法和文献[9]的方差分量的估计方法以及分治算法,提出海量数据下线性混合效应模型参数的估计算法。首先提出线性混合模型参数的三步估计方法,避免了海量数据下通常所需的基于极大似然估计或限制极大似然估计的迭代计算,再针对海量数据的两种不同情形,基于分治算法分步骤计算估计量。实验表明本文算法不仅可以减少内存占用,解决存储问题,而且可以缩短计算时间,提高运算速度。
1 线性混合效应模型及其估计方法 1.1 线性混合效应模型考虑纵向数据线性混合效应模型
$ {Y_{ij}} = \mathit{\boldsymbol{X}}_{ij}^{\rm{T}}\mathit{\boldsymbol{\beta }} + \mathit{\boldsymbol{Z}}_{ij}^{\rm{T}}{\mathit{\boldsymbol{b}}_i} + {\varepsilon _{ij}}, i = 1, \cdots , n;j = 1, \cdots , {n_i} $ | (1) |
式中,未知参数β是p×1维的固定效应向量,bi是q×1维随机效应向量,Xij和Zij分别表示与之相关的协变量。假定随机效应bi~N(0, Σ),模型误差εij~N(0, σ2),二者均满足独立同分布条件,且bi与εij相互独立。
1.2 三步估计方法为了避免迭代计算,提出三步估计方法分别对线性混合效应模型的方差分量和固定效应进行估计,利用与文献[9]类似的证明步骤可以得到估计的相合性和渐近正态性,因此本文只考虑估计算法。
首先利用最小二乘法对系数进行初步估计。令
$ \mathit{\boldsymbol{Y}} = \mathit{\boldsymbol{X\beta }} + \mathit{\boldsymbol{Zb}} + \mathit{\boldsymbol{\varepsilon }} $ | (2) |
对该模型系数直接求最小二乘估计,可得初步估计
$ {\mathit{\boldsymbol{\widehat \beta }}_0} = {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{Y}} $ | (3) |
其次估计方差分量。利用文献[9]的原理来估计方差分量σ2,Σ。
令Uij= ZijTbi+εij,则有Yij= XijT β +Uij,进一步有
$ {{\mathit{\boldsymbol{\hat S}}}_i} = \mathit{\boldsymbol{\hat U}}_i^{\rm{T}}{{\mathit{\boldsymbol{\hat U}}}_i} - \mathit{\boldsymbol{\hat U}}_i^{\rm{T}}{\mathit{\boldsymbol{H}}_i}{{\mathit{\boldsymbol{\hat U}}}_i} $ | (4) |
根据文献[9]的估计步骤,进一步考虑系数的估计值
$ \begin{array}{l} {{\hat \sigma }^2} = {\left[ {\sum\limits_{i = 1}^n {\left( {{n_i} - q} \right)} - p} \right]^{ - 1}}\sum\limits_{i = 1}^n {{{\hat S}_i}} = (N - qn - \\ p{)^{ - 1}}\sum\limits_{i = 1}^n {{{\hat S}_i}} \end{array} $ | (5) |
式中
$ {\mathit{\boldsymbol{\widetilde b}}_i} = {\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)^{ - 1}}\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{U}}_i} = {\mathit{\boldsymbol{b}}_i} + {\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{{\bf{Z}}_i}} \right)^{ - 1}}\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{\varepsilon }}_i} $ |
进一步有
$ \begin{array}{l} \sum\limits_{i = 1}^n {{{\mathit{\boldsymbol{\widetilde b}}}_i}} \mathit{\boldsymbol{\widetilde b}}_i^{\rm{T}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{b}}_i}} \mathit{\boldsymbol{b}}_i^{\rm{T}} + \sum\limits_{i = 1}^n {{{\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)}^{ - 1}}} \mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{\varepsilon }}_i}\mathit{\boldsymbol{\varepsilon }}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}\\ {\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)^{ - 1}} + \sum\limits_{i = 1}^n {{{\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)}^{ - 1}}} \mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{\varepsilon }}_i}\mathit{\boldsymbol{b}}_i^{\rm{T}} + \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{b}}_i}} \mathit{\boldsymbol{\varepsilon }}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}\\ {\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)^{ - 1}} \end{array} $ | (6) |
通过直接计算各项的一、二阶矩,可证式(6)中最后两项的阶为
$ \begin{array}{l} \frac{1}{n}\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{b}}_i}} \mathit{\boldsymbol{b}}_i^{\rm{T}} \approx \frac{1}{n}\left[ {\sum\limits_{i = 1}^n {{{\mathit{\boldsymbol{\widetilde b}}}_i}} \mathit{\boldsymbol{\widetilde b}}_i^{\rm{r}} - \sum\limits_{i = 1}^n {{{\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)}^{ - 1}}} } \right.\\ \left. {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{\varepsilon }}_i}\mathit{\boldsymbol{\varepsilon }}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}{{\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)}^{ - 1}}} \right] \approx \frac{1}{n}\left[ {\sum\limits_{i = 1}^n {{{\mathit{\boldsymbol{\widetilde b}}}_i}} \mathit{\boldsymbol{\widetilde b}}_i^{\rm{r}} - {\sigma ^2}\sum\limits_{i = 1}^n {} } \right.\\ \left. {{{\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)}^{ - 1}}} \right] \end{array} $ |
从而可得
$ \mathit{\boldsymbol{ \boldsymbol{\hat \varSigma} }} = \frac{1}{n}\sum\limits_{i = 1}^n {{{\mathit{\boldsymbol{\widehat b}}}_i}} \mathit{\boldsymbol{\widehat b}}_i^{\rm{T}} - \frac{1}{n}{\hat \sigma ^2}\sum\limits_{i = 1}^n {{{\left( {\mathit{\boldsymbol{Z}}_i^{\rm{T}}{\mathit{\boldsymbol{Z}}_i}} \right)}^{ - 1}}} $ | (7) |
最后计算β的加权最小二乘估计,记
$ \boldsymbol{V}=\operatorname{diag}\left(\boldsymbol{V}_{1}, \cdots, \boldsymbol{V}_{n}\right) $ |
$ {\mathit{\boldsymbol{V}}_i} = {\mathop{\rm var}} \left( {{\mathit{\boldsymbol{Z}}_i}{\mathit{\boldsymbol{b}}_i} + {\mathit{\boldsymbol{\varepsilon }}_i}} \right) = {\mathit{\boldsymbol{Z}}_i}\mathit{\boldsymbol{ \boldsymbol{\varSigma} Z}}_i^{\rm{T}} + {\mathit{\boldsymbol{\sigma }}^2}{\mathit{\boldsymbol{I}}_{{n_i}}} $ |
由公式(6)、(7)有
$ {\mathit{\boldsymbol{\widehat V}}_i} = {\mathit{\boldsymbol{Z}}_i}\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}\mathit{\boldsymbol{Z}}_i^{\rm{T}} + {\hat \sigma ^2}{\mathit{\boldsymbol{I}}_{{n_i}}}, \quad \widehat {\bf{V}} = {\mathop{\rm diag}\nolimits} \left( {{{\widehat {\bf{V}}}_1}, \cdots , {{\widehat {\bf{V}}}_n}} \right) $ |
根据文献[5]的加权最小二乘估计方法可得
$ \hat{\boldsymbol{\beta}}=\left(\boldsymbol{X}^{\mathit{\boldsymbol{T}}} \hat{\boldsymbol{V}}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\mathit{\boldsymbol{T}}} \hat{\boldsymbol{V}}^{-1} \boldsymbol{Y} $ | (8) |
或利用文献[5]的样本求和公式(9)计算
$ \hat{\boldsymbol{\beta}}=\left(\sum\limits_{i=1}^{n} \boldsymbol{X}_{i}^{\mathit{\boldsymbol{T}}} \hat{\boldsymbol{V}}_{i}^{-1} \boldsymbol{X}_{\boldsymbol{i}}\right)^{-1}\left(\sum\limits_{i=1}^{n} \boldsymbol{X}_{\boldsymbol{i}}^{\mathit{\boldsymbol{T}}} \hat{\boldsymbol{V}}_{i}^{-1} \boldsymbol{Y}_{i}\right) $ | (9) |
考虑海量数据的两种情形:①个体数据量较大,但组内数据量较小;②个体数据量较小,但组内数据量较大。此时会遇到存储瓶颈及计算低效的问题,1.2节中的估计方法不再适用。因此在此两种情形下分别考虑利用分治算法对前面提出的三步估计方法进行调整。
2.1 个体数据量较大,组内数据量较小情形令
首先计算初步估计。
记
则第k个子集对应的数据模型可写作
$ \boldsymbol{Y}_{k}=\boldsymbol{X}_{k} \boldsymbol{\beta}+\boldsymbol{Z}_{k} \boldsymbol{b}_{k}+\boldsymbol{\varepsilon}_{k} $ | (10) |
对每个子数据集直接作最小二乘估计有
$ \begin{array}{l} {\mathit{\boldsymbol{\widehat \beta }}_{0}} = {\left( {\mathit{\boldsymbol{X}}^{\rm{T}}{\mathit{\boldsymbol{X}}}} \right)^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{Y}} = {\left( {\sum\limits_{k = 1}^K {{\bf{X}}_k^{\rm{T}}} {{\bf{X}}_k}} \right)^{ - 1}}\\ \sum\limits_{k = 1}^K {{\bf{X}}_k^{\rm{T}}} {{\bf{X}}_k }{\mathit{\boldsymbol{\widehat \beta }}_{0k}} \end{array} $ | (11) |
其次利用分治算法估计方差分量。记
$ {\mathit{\boldsymbol{\hat S}}_k} = \mathit{\boldsymbol{\hat U}}_k^{\rm{T}}{\mathit{\boldsymbol{\hat U}}_k} - \mathit{\boldsymbol{\hat U}}_k^{\rm{T}}{{\bf{H}}_k}{\widehat {\bf{U}}_k} $ | (12) |
所以有
$ \hat{\sigma}^{2}=(N-q n-p)^{-1} \sum\limits_{k=1}^{K} \hat{S}_{k} $ | (13) |
$ \mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }} = \frac{1}{n}\sum\limits_{k = 1}^K {\sum\limits_{i = 1}^m {{{\mathit{\boldsymbol{\widehat b}}}_{ki}}} } \mathit{\boldsymbol{\widehat b}}_{ki}^{\rm{T}} - \frac{1}{n}{\hat \sigma ^2}\sum\limits_{k = 1}^K {\sum\limits_{i = 1}^m {{{\left( {\mathit{\boldsymbol{Z}}_{ki}^{\rm{T}}{\mathit{\boldsymbol{Z}}_{ki}}} \right)}^{ - 1}}} } $ | (14) |
由1.2节可知,式中
公式(11)~(13)中,对应每个数据子集需计算保存的数据为
最后计算模型系数β的加权估计。根据方差分量的计算公式(12)~(13)可知,第k个子数据模型的协方差矩阵的估计为
$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\widehat \beta }} = {{\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}{{\mathit{\boldsymbol{\widehat V}}}^{ - 1}}\mathit{\boldsymbol{X}}} \right)}^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}{{\mathit{\boldsymbol{\widehat V}}}^{ - 1}}\mathit{\boldsymbol{Y}} = }\\ {{{\left( {\sum\limits_{k = 1}^K {\mathit{\boldsymbol{X}}_k^{\rm{T}}} \mathit{\boldsymbol{\widehat V}}_k^{ - 1}{\mathit{\boldsymbol{X}}_k}} \right)}^{ - 1}}\sum\limits_{k = 1}^K {\mathit{\boldsymbol{X}}_k^{\rm{T}}} \mathit{\boldsymbol{\widehat V}}_k^{ - 1}{\mathit{\boldsymbol{X}}_k}{{\mathit{\boldsymbol{\widehat \beta }}}_k}} \end{array} $ | (15) |
对每个个体运用分治算法。
令Di={(Xij, Zij, Yij)}j=1ni,i=1, …, n, 将数据Di划分为K个子集Di1, Di2, …, DiK,每个子集的样本数为mi=ni/K。每个子集包含的数据为Xikj, Zikj, Yikj,其中1≤j≤mi,1≤k≤K,1≤i≤n,并且每个子集包含的数据不能超过单机处理能力。
首先计算初步估计。记Xik=(Xik1, …, Xikmi)T,Xi=(Xi1T, …, XiKT)T,Yi, Zi, εi类似定义,第i个个体的第k个子集对应的数据模型为
$ \boldsymbol{Y}_{i k}=\boldsymbol{X}_{i k} \boldsymbol{\beta}+\boldsymbol{Z}_{i k} \boldsymbol{b}_{i}+\boldsymbol{\varepsilon}_{i k} $ | (16) |
对每个个体的每个数据集直接作最小二乘估计有
根据分治算法,整合每个个体每个子集的结果,保存每组
$ \begin{array}{l} {\mathit{\boldsymbol{\widehat \beta }}_0} = {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{Y}} = \left[ {\sum\limits_{i = 1}^n {\left( {\sum\limits_{k = 1}^K {\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}} } \right.} } \right.\\ {\mathit{\boldsymbol{X}}_{ik}}){]^{ - 1}}\sum\limits_{i = 1}^n {\left( {\sum\limits_{k = 1}^K {\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}} {\mathit{\boldsymbol{X}}_i}{{\mathit{\boldsymbol{\widehat \beta }}}_{0ik}}} \right)} \end{array} $ | (17) |
其次利用分治算法估计方差分量。记
$ \hat{S}_{i}=\sum\limits_{k=1}^{K} \hat{\boldsymbol{U}}_{i k}^{\mathit{\boldsymbol{T}}} \hat{\boldsymbol{U}}_{i k}-\sum\limits_{k=1}^{K} \hat{\boldsymbol{U}}_{i k}^{\mathit{\boldsymbol{T}}} \boldsymbol{H}_{i k} \hat{\boldsymbol{U}}_{i k} $ | (18) |
所以有
$ \hat{\sigma}^{2}=(N-q n-p)^{-1} \sum\limits_{i=1}^{n} \hat{S}_{i} $ | (19) |
$ \mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }} = \frac{1}{n}\sum\limits_{i = 1}^n {{{\mathit{\boldsymbol{\widehat b}}}_i}} \mathit{\boldsymbol{\widehat b}}_i^{\rm{T}} - \frac{1}{n}{\hat \sigma ^2}\sum\limits_{i = 1}^n {{{\left( {\sum\limits_{k = 1}^K {\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}} {\mathit{\boldsymbol{Z}}_{ik}}} \right)}^{ - 1}}} $ | (20) |
由1.2节可知,式中
根据公式(17)~(18)可知,对应每个个体子集,需计算保存的数据有
最后计算模型系数β的加权估计。根据方差分量的计算公式(18)~(19),第i个个体的第k个子数据模型的协方差矩阵的估计为
$ \hat{\boldsymbol{\beta}}_{i k}=\left(\boldsymbol{X}_{i k}^{\mathit{\boldsymbol{T}}} \hat{\boldsymbol{V}}_{i k}^{-1} \boldsymbol{X}_{i k}\right)^{-1} \boldsymbol{X}_{i k}^{\mathit{\boldsymbol{T}}} \hat{\boldsymbol{V}}_{i k}^{-1} \boldsymbol{Y}_{i k} $ | (21) |
由此可知,对每个个体的每个子集保留数据
$ \begin{array}{l} \mathit{\boldsymbol{\widehat \beta }} = {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}{{\mathit{\boldsymbol{\widehat V}}}^{ - 1}}\mathit{\boldsymbol{X}}} \right)^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}{\mathit{\boldsymbol{\widehat V}}^{ - 1}}\mathit{\boldsymbol{Y}} = \left[ {\sum\limits_{i = 1}^n {\left( {\sum\limits_{k = 1}^k {\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}} \mathit{\boldsymbol{\widehat V}}_{ik}^{ - 1}} \right.} } \right.\\ {\mathit{\boldsymbol{X}}_{ik}}){]^{ - 1}}\sum\limits_{i = 1}^n {\left( {\sum\limits_{k = 1}^K {\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}} \mathit{\boldsymbol{\widehat V}}_{ik}^{ - 1}{\mathit{\boldsymbol{X}}_{ik}}{{\mathit{\boldsymbol{\widehat \beta }}}_{ik}}} \right)} \end{array} $ | (22) |
为了将本文估计算法与极大似然法[5]进行对比,验证本文算法在计算时间上的优势,同时检验算法的可行性,采用了Matlab软件进行数值模拟。
3.1 模拟数据设纵向数据线性混合效应模型如下
$ Y_{i j}=\boldsymbol{X}_{i j}^{\mathit{\boldsymbol{T}}} \boldsymbol{\beta}+\boldsymbol{Z}_{i j}^{\mathit{\boldsymbol{T}}} \boldsymbol{b}_{i}+\varepsilon_{i j}, i=1, \cdots, n ; j=1, \cdots, n_{i} $ | (23) |
分别考虑海量数据下的两种样本情形。
情形(1) 个体数量较大,但组内数据量较小时,考虑两种样本:①n=10000, ni=10,总样本量100000;②n=100000, ni=10,总样本量1000000。
情形(2) 个体有有限个,但组内数据量较大时,样本为n=10, ni=10000。
另外参数部分系数为:β =(1, 2, -3, 1, -2)T;Xij=(Xij1, Xij2, Xij3, Xij4, Xij5)T,其中Xij1~N(0, 1),Xij2~N(0, 2),Xij3~N(0, 1),Xij4~N(0, 2),Xij5~N(0, 1),且Xij1、Xij2、Xij3、Xij4、Xij5相互独立;Zij=(Zij1, Zij2)T,Zij1和Zij2都独立同分布于N(0, 2);bi=(bi1, bi2)T,bi1和bi2都独立同分布于N(0, 1);εij~N(0, σ2),其中σ2=2,且bi与εi相互独立。模拟次数50次。
3.2 情形(1)的模拟结果与分析在情形(1)的两种样本情况下,分别用极大似然法[5]和本文算法估计模型系数和方差分量,并将两种方法所用的计算时间及估计精度进行对比。
3.2.1 估计精度和计算效率利用均方误差(MSE)衡量模拟估计值与真实值之间的偏差,对估计精度进行分析,MSE值越小,估计精度越高。均方误差计算公式如下
$ \gamma_{\rm{{MSE}}}=\frac{1}{m} \sum\limits_{k=1}^{m}\left(x_{{\rm{obs}}, k}-x_{\rm{{real}}}\right)^{2} $ | (24) |
式(24)中,m是模拟次数,xreal为参数真值,xobs, k是第k次模拟得到的参数估计值。
通过模拟得到情形(1)中两种样本量下的计算时间和均方误差的值如表 1所示。
由表 1可知,极大似然法与本文算法估计出的模型参数的MSE值都较小,说明参数的估计效果较好,估计精确度较高。另外,与极大似然法相比,本文算法估计参数可以在几乎不降低估计精确度的同时将运行效率提高4倍,说明该算法可以大幅减少计算时间,提高计算速度。
3.2.2 数据集块数对计算时间的影响为了研究计算时间与子数据集块数之间的关系,以说明利用分治算法进行分块计算的必要性,分别将两种样本量下的数据分成不同的数据集块数,记录估计参数所用时间,结果如表 2、3所示。
由表 2和表 3看出, 当K=1时,直接利用加权最小二乘公式(8)估计模型系数,此时会超出内存而无法计算。样本①、②下,即K分别取10 000和100 000时,利用样本求和公式(9)估计模型系数。K取表中其他值时,按照本文算法进行计算,当K缓慢增加时,所用时间逐渐减少;当K取一个恰当的值,如K分别取样本①、②中2 000和20 000时,计算时间达到最小,并且与样本求和方式相比,运行速度可提高20%以上。这是因为当K增加时,每个数据集块分到的样本数量同时减少,因此计算量和计算时间也随之减少;同时在单机处理能力有限且固定情况下,一定数量级的数据在合适的单机数目下可以达到最佳计算速度。而实际中的数据量更大,通常达百万级甚至千万级,此时本文算法的效果会更好。
综上说明,本文算法可以降低内存开销,减少计算时间,提高计算速度。
3.3 情形(2)的模拟结果与分析情形(2)下,本文算法所用计算时间与数据集块数的关系如表 4所示。K=1时利用样本求和公式(9)计算系数;K取其他值时利用本文算法计算,可以看出时间随数据块数的变化趋势与情形(1)类似,且K=200时t=1.38 s, 计算时间达到最小。同时利用极大似然法计算情形(2)模型参数,所用时间为2 269.71 s,再次证明了本文算法的高效性。
选取2006—2010年全国的285个地级及以上城市的地区生产总值、工业总产值和社会消费品零售总额、固定资产投资与城乡居民储蓄年末余额数据,实证研究工业总产值和社会消费品零售总额对地区生产总值的影响。
以地区生产总值(Y)作为响应变量,工业总产值(X1)和社会消费品零售总额(X2)作为固定效应的协变量;考虑到固定资产投资(Z1)和城乡居民储蓄年末余额(Z2)也会影响地区生产总值,故将这两个因素作为随机效应部分, 建立线性混合效应模型
$ Y_{i j}=\boldsymbol{X}_{i j}^{\mathit{\boldsymbol{T}}} \boldsymbol{\beta}+\boldsymbol{Z}_{i j}^{\mathit{\boldsymbol{T}}} \boldsymbol{b}_{i}+\boldsymbol{\varepsilon}_{i j}, i=1, \cdots, 285 ; j=1, \cdots, 5 $ | (25) |
将数据带入模型(25)可知,此时数据结构符合情形(1),因此利用情形(1)下的估计算法进行参数估计,具体结果见表 5。
由表 5可知,工业总产值和社会消费品零售总额的系数都为正值,说明这两个因素正面影响地区生产总值,且所建模型是有效的。所以要提高地区生产总值,就要努力提高工业生产水平,增加社会消费品零售总额,提高民众收入,改善消费环境,稳定物价水平,解决医疗、卫生、就业等问题。
5 结束语本文针对海量数据下的两种情形,将线性混合效应模型系数和方差分量的估计方法与分治算法相结合,提出了三步估计算法,使得复杂的统计模型也可以应用于海量数据情形。数值模拟表明该算法可以大幅减少计算时间,提高计算效率,并解决海量数据下计算机内存不足的问题。通过实证分析证明了本文模型及估计算法在实际应用中的可行性。
下一步,可以考虑采用分治算法与其他估计方法结合,从而实现其他统计模型在海量数据下的估计与应用。
[1] |
JIANG J. Linear and generalized linear mixed models and their applications[M]. New York: Springer, 2007.
|
[2] |
李光辉.线性混合模型中的参数估计[D].西宁: 青海师范大学, 2010. LI G H. Parameter estimation in linear mixed model[D]. Xining: Qing Hai Normal University, 2010. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10746-1011040283.htm |
[3] |
许王莉. 线性混合效应模型中方差分量的估计[J]. 应用概率统计, 2009, 25(3): 301-308. XU W L. Estimation of variance components in linear mixed effects models[J]. Application Probability Statistics, 2009, 25(3): 301-308. (in Chinese) DOI:10.3969/j.issn.1001-4268.2009.03.009 |
[4] |
孙燕, 柴根象. 纵向数据线性混合效应模型的统计分析[J]. 应用科学学报, 2006, 24(1): 70-73. SUN Y, CHAI G X. Statistical analysis of linear mixed effect model for longitudinal data[J]. Journal of Applied Science, 2006, 24(1): 70-73. (in Chinese) DOI:10.3969/j.issn.0255-8297.2006.01.015 |
[5] |
WU H L, ZHANG J T. Nonparametric regression methods for longitudinal data analysis[M]. Hoboken: John Wiley & Sons, Inc., 2006: 17-25.
|
[6] |
WU X D, ZHU X Q, WU G Q, et al. Data mining with big data[J]. IEEE Transactions on Knowledge & Data Engineering, 2013, 26(1): 97-107. |
[7] |
FAN J Q, HAN F, LIU H. Challenges of big data analysis[J]. National Science Review, 2014, 1(2): 293-314. DOI:10.1093/nsr/nwt032 |
[8] |
ZHOU Z H, CHAWLA N V, WILLIAMS G J, et al. Big data opportunities and challenges:discussions from data analytics perspectives[J]. IEEE Computational Intelligence Magazine, 2014, 9(4): 62-74. DOI:10.1109/MCI.2014.2350953 |
[9] |
SUN Y, ZHANG W Y, TONG H. Estimation of the covariance matrix of random effects in longitudinal studies[J]. Annals of Statistics, 2007, 35(6): 2795-2814. DOI:10.1214/009053607000000523 |
[10] |
LIN N, XI R B. Aggregated estimating equation estimation[J]. Statistics & Its Interface, 2011, 1(1): 73-83. |
[11] |
CHANG X Y, LIN S B, WANG Y. Divide and conquer local average regression[J]. Electronic Journal of Statistics, 2017, 11(1): 1326-1350. DOI:10.1214/17-EJS1265 |
[12] |
GUHA S, HAFEN R, ROUNDS J, et al. Large complex data:divide and recombine (d & r) with rhipe[J]. The ISI's Journal for the Rapid Dissemination of Statistics Research, 2015, 1(1): 53-67. |
[13] |
CHEN X Y, XIE M G. A split-and-conquer approach for analysis of extraordinarily large data[J]. Statistica Sinica, 2014, 24(4): 1655-1684. |