文章快速检索     高级检索
  北京化工大学学报(自然科学版)  2019, Vol. 46 Issue (3): 123-128   DOI: 10.13543/j.bhxbzr.2019.03.019
0

引用本文  

耿俏, 李志强, 陈少东. 海量数据下线性混合效应模型的估计算法[J]. 北京化工大学学报(自然科学版), 2019, 46(3): 123-128. DOI: 10.13543/j.bhxbzr.2019.03.019.
GENG Qiao, LI ZhiQiang, CHEN ShaoDong. Estimation algorithm for a linear mixed effect model with massive data[J]. Journal of Beijing University of Chemical Technology (Natural Science), 2019, 46(3): 123-128. DOI: 10.13543/j.bhxbzr.2019.03.019.

第一作者

耿俏, 女, 1993年生, 硕士生.

通信联系人

李志强, E-mail:li-zhiqiang2000@163.com

文章历史

收稿日期:2018-09-15
海量数据下线性混合效应模型的估计算法
耿俏 , 李志强 , 陈少东     
北京化工大学 理学院, 北京 100029
摘要:基于以往文献提出线性混合效应模型参数的三步估计方法,避免了繁杂的极大似然估计迭代步骤。同时为进一步解决海量数据下计算估计量时存在的存储瓶颈及计算时间过长问题,在海量纵向数据的两种不同数据格式下,分别基于三步估计方法利用分治算法计算模型参数的估计量。数值模拟和实证分析结果表明,本文所提出的三步估计方法和估计量的分治算法可以减轻计算负担,减少占用内存,解决内存不足的问题,并提高计算速度。
关键词海量数据    纵向数据    线性混合效应模型    三步估计方法    分治算法    
Estimation algorithm for a linear mixed effect model with massive data
GENG Qiao , LI ZhiQiang , CHEN ShaoDong     
Faculty of Science, Beijing University of Chemical Technology, Beijing 100029, China
Abstract: Based on the previous literature, a three-step method for the estimation of the parameters of a linear mixed effect model has been proposed, which avoids the complicated iterative steps of maximum likelihood estimation. At the same time, in order to further solve the storage bottleneck and calculation time when calculating the estimator with massive data, the estimator of the model parameters has been calculated using the divide-and-conquer algorithm based on the three-step estimation method for two different data formats of massive vertical data. The results of numerical simulation and empirical analysis show that the three-step estimation method and the estimator divide-and-conquer algorithm proposed in this paper can reduce the computational burden, reduce the memory consumption, solve the problem of insufficient memory, and improve the calculation speed.
Key words: massive data    longitudinal data    linear mixed effect model    three-step estimation method    divide-and-conquer algorithm    
引言

纵向数据是通过对每个个体在不同的时间点上进行观测得到的数据集,广泛应用在经济学、医学等领域,具有组内数据相关、组间数据独立的特点。线性混合效应模型[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维的固定效应向量,biq×1维随机效应向量,XijZij分别表示与之相关的协变量。假定随机效应bi~N(0, Σ),模型误差εij~N(0, σ2),二者均满足独立同分布条件,且biεij相互独立。

1.2 三步估计方法

为了避免迭代计算,提出三步估计方法分别对线性混合效应模型的方差分量和固定效应进行估计,利用与文献[9]类似的证明步骤可以得到估计的相合性和渐近正态性,因此本文只考虑估计算法。

首先利用最小二乘法对系数进行初步估计。令$\mathit{\boldsymbol{Y}}={{\left( \mathit{\boldsymbol{Y}}_{1}^{\text{T}}, \cdots , \mathit{\boldsymbol{Y}}_{n}^{\text{T}} \right)}^{\text{T}}}$${{\mathit{\boldsymbol{Y}}}_{i}}={{\left( {{\mathit{\boldsymbol{Y}}}_{i1}}, \cdots , {{\mathit{\boldsymbol{Y}}}_{i{{n}_{i}}}} \right)}^{\text{T}}}$, Xε类似定义,$\mathit{\boldsymbol{b}}={{\left( \mathit{\boldsymbol{b}}_{1}^{\text{T}}, \cdots , \mathit{\boldsymbol{b}}_{n}^{\text{T}} \right)}^{\text{T}}}$${{\mathit{\boldsymbol{Z}}}_{i}}={{\left( {{\mathit{\boldsymbol{Z}}}_{i1}}, \cdots , {{\mathit{\boldsymbol{Z}}}_{i{{n}_{i}}}} \right)}^{\text{T}}}$Z =diag(Z1, …, Zn), 则模型(1)可重新表示为

$ \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,进一步有${{{\hat{U}}}_{ij}}={{Y}_{ij}}-\mathit{\boldsymbol{X}}_{ij}^{\text{T}}{{\widehat{\mathit{\boldsymbol{ }}\!\!\beta\!\!\text{ }}}_{0}}$。记Ui=(Ui1, …, Uini)T,则${{{\mathit{\boldsymbol{\hat{U}}}}}_{i}}={{\left( {{{\hat{U}}}_{i1}}, \cdots , {{{\hat{U}}}_{i{{n}_{i}}}} \right)}^{\text{T}}}$,进一步有Ui= Zi bi+ εi。从而利用最小二乘法进行估计可得随机效应bi的估计值${{\widetilde{\mathit{\boldsymbol{b}}}}_{i}}={{\left( \mathit{\boldsymbol{Z}}_{i}^{\text{T}}{{\mathit{\boldsymbol{Z}}}_{i}} \right)}^{-1}}\mathit{\boldsymbol{Z}}_{i}^{\text{T}}{{\mathit{\boldsymbol{U}}}_{i}}$,其中Ui是未知的,故用Ui的估计值${{{\mathit{\boldsymbol{\hat{U}}}}}_{i}}$代替Ui可得${{\widehat{\mathit{\boldsymbol{b}}}}_{i}}={{\left( \mathit{\boldsymbol{Z}}_{i}^{\text{T}}{{\mathit{\boldsymbol{Z}}}_{i}} \right)}^{-1}}\mathit{\boldsymbol{Z}}_{i}^{\text{T}}{{\widehat{\mathit{\boldsymbol{U}}}}_{i}}$。进一步记${{{\mathit{\boldsymbol{\tilde{U}}}}}_{i}}={{\mathit{\boldsymbol{Z}}}_{i}}{{\widetilde{\mathit{\boldsymbol{b}}}}_{i}}={{\mathit{\boldsymbol{Z}}}_{i}}{{\left( \mathit{\boldsymbol{Z}}_{i}^{\text{T}}{{\mathit{\boldsymbol{Z}}}_{i}} \right)}^{-1}}\mathit{\boldsymbol{Z}}_{i}^{\text{T}}{{\mathit{\boldsymbol{U}}}_{i}}\equiv {{\mathit{\boldsymbol{H}}}_{i}}{{\mathit{\boldsymbol{U}}}_{i}}$,则模型${{{\mathit{\boldsymbol{\tilde{U}}}}}_{i}}\text{=}{{\mathit{\boldsymbol{Z}}}_{i}}{{\mathit{\boldsymbol{b}}}_{i}}+{{\mathit{\boldsymbol{ }}\!\!\varepsilon\!\!\text{ }}_{i}}$的残差平方和为${{S}_{i}}={{\left( {{\mathit{\boldsymbol{U}}}_{i}}-{{{\mathit{\boldsymbol{\tilde{U}}}}}_{i}} \right)}^{\text{T}}}\left( {{\mathit{\boldsymbol{U}}}_{i}}-{{{\mathit{\boldsymbol{\tilde{U}}}}}_{i}} \right)=\mathit{\boldsymbol{U}}_{i}^{\text{T}}{{\mathit{\boldsymbol{U}}}_{i}}-\mathit{\boldsymbol{U}}_{i}^{\text{T}}{{\mathit{\boldsymbol{H}}}_{i}}{{\mathit{\boldsymbol{U}}}_{i}}$,从而有

$ {{\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]的估计步骤,进一步考虑系数的估计值${\mathit{\boldsymbol{\widehat \beta }}_0}$ 的影响,可以类似构造出方差的估计

$ \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)

式中$N=\sum\limits_{i=1}^{n}{{{n}_{i}}}$qbi的维数,pβ的维数。下一步估计Σ。根据Ui的定义可知

$ {\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)中最后两项的阶为${O_p}\left( {{n^{\frac{1}{2}}}} \right)$,因此相对于其他项可将其忽略不计,由此可得

$ \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)
2 海量数据下线性混合效应模型的分治算法

考虑海量数据的两种情形:①个体数据量较大,但组内数据量较小;②个体数据量较小,但组内数据量较大。此时会遇到存储瓶颈及计算低效的问题,1.2节中的估计方法不再适用。因此在此两种情形下分别考虑利用分治算法对前面提出的三步估计方法进行调整。

2.1 个体数据量较大,组内数据量较小情形

$\mathit{\boldsymbol{D}} = \left\{ {\left( {{\mathit{\boldsymbol{X}}_{ij}}, {\mathit{\boldsymbol{Z}}_{ij}}, {\mathit{\boldsymbol{Y}}_{ij}}} \right)} \right\}_{i = 1j = 1}^{{n_i}}$为全部数据集。根据分治算法,将数据集D划分为K个子集D1, D2, …, DK,每个子集的样本数m=n/K,每个子集包含的数据为Xkij, Zkij, Ykij,其中1≤jnki,1≤im,1≤kK且每个子集包含的数据不能超过单机处理能力。

首先计算初步估计。

${\mathit{\boldsymbol{X}}_{ki}} = {\left( {{\mathit{\boldsymbol{X}}_{ki1}}, \cdots , {\mathit{\boldsymbol{X}}_{ki{n_{ki}}}}} \right)^{\rm{T}}}$${\mathit{\boldsymbol{X}}_k} = {\left( {\mathit{\boldsymbol{X}}_{k1}^{\rm{T}}, \cdots , \mathit{\boldsymbol{X}}_{km}^{\rm{T}}} \right)^{\rm{T}}}$, Ykεk类似定义,Zki=(Zki1, …, Zkinki)TZk=diag(Zk1, …, Zkm),bk=(bk1T, …, bkmT)T

则第k个子集对应的数据模型可写作

$ \boldsymbol{Y}_{k}=\boldsymbol{X}_{k} \boldsymbol{\beta}+\boldsymbol{Z}_{k} \boldsymbol{b}_{k}+\boldsymbol{\varepsilon}_{k} $ (10)

对每个子数据集直接作最小二乘估计有

${\mathit{\boldsymbol{\widehat \beta }}_{0k}} = {\left( {\mathit{\boldsymbol{X}}_k^{\rm{T}}{\mathit{\boldsymbol{X}}_k}} \right)^{ - 1}}\mathit{\boldsymbol{X}}_k^{\rm{T}}{\mathit{\boldsymbol{Y}}_k}$,1≤kK,根据分治算法,整合各个子集的结果,保存每组$\mathit{\boldsymbol{X}}_k^{\rm{T}}{\mathit{\boldsymbol{X}}_k}$, ${\mathit{\boldsymbol{\widehat \beta }}_{0k}}$,最后求得初步估计

$ \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 U}}}_{ki}} = {\left( {{{\mathit{\boldsymbol{\hat U}}}_{ki1}}, \cdots , {{\mathit{\boldsymbol{\hat U}}}_{ki{n_i}}}} \right)^{\rm{T}}}$${{\mathit{\boldsymbol{\hat U}}}_k} = {\left( {{{\mathit{\boldsymbol{\hat U}}}_{k1}}, \cdots , {{\mathit{\boldsymbol{\hat U}}}_{km}}} \right)^{\rm{T}}}$与公式(4)类似,可得

$ {\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节可知,式中${\mathit{\boldsymbol{\widehat b}}_{ki}} = {\left( {\mathit{\boldsymbol{Z}}_{ki}^{\rm{T}}{\mathit{\boldsymbol{Z}}_{ki}}} \right)^{ - 1}}\mathit{\boldsymbol{Z}}_{ki}^{\rm{T}}{\mathit{\boldsymbol{\widehat U}}_{ki}}$

公式(11)~(13)中,对应每个数据子集需计算保存的数据为$\left\{ {\mathit{\boldsymbol{\widehat U}}_k^{\rm{T}}{{\mathit{\boldsymbol{\widehat U}}}_k}, \mathit{\boldsymbol{\widehat U}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}{{\mathit{\boldsymbol{\widehat U}}}_k}, \mathit{\boldsymbol{Z}}_{ki}^{\rm{T}}{\mathit{\boldsymbol{Z}}_{ik}}, \mathit{\boldsymbol{Z}}_{ki}^{\rm{T}}{{\mathit{\boldsymbol{\widehat U}}}_{ki}}} \right\}$

最后计算模型系数β的加权估计。根据方差分量的计算公式(12)~(13)可知,第k个子数据模型的协方差矩阵的估计为${\mathit{\boldsymbol{\widehat V}}_k} = {\mathit{\boldsymbol{Z}}_k}\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}\mathit{\boldsymbol{Z}}_k^{\rm{T}} + {{\hat \sigma }^2}{\mathit{\boldsymbol{I}}_{{n_i} \times m}}$。对每个子集用加权最小二乘法进行估计有${\mathit{\boldsymbol{\widehat \beta }}_k} = {\left( {\mathit{\boldsymbol{X}}_k^{\rm{T}}\mathit{\boldsymbol{\widehat V}}_k^{ - 1}{\mathit{\boldsymbol{X}}_k}} \right)^{ - 1}}\mathit{\boldsymbol{X}}_k^{\rm{T}}\mathit{\boldsymbol{\widehat V}}_k^{ - 1}{\mathit{\boldsymbol{Y}}_k}$,由此可知,保留数据$\mathit{\boldsymbol{X}}_k^{\rm{T}}\mathit{\boldsymbol{\widehat V}}_k^{ - 1}{\mathit{\boldsymbol{X}}_k}$${\mathit{\boldsymbol{\widehat \beta }}_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)
2.2 个体数据量较小,组内数据量较大情形

对每个个体运用分治算法。

Di={(Xij, Zij, Yij)}j=1nii=1, …, n, 将数据Di划分为K个子集Di1, Di2, …, DiK,每个子集的样本数为mi=ni/K。每个子集包含的数据为Xikj, Zikj, Yikj,其中1≤jmi,1≤kK,1≤in,并且每个子集包含的数据不能超过单机处理能力。

首先计算初步估计。记Xik=(Xik1, …, Xikmi)TXi=(Xi1T, …, XiKT)TYi, 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)

对每个个体的每个数据集直接作最小二乘估计有${\mathit{\boldsymbol{\widehat \beta }}_{{\rm{oik }}}} = {\left( {\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{X}}_{ik}}} \right)^{ - 1}}\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{Y}}_{ik}}$,1≤in, 1≤kK

根据分治算法,整合每个个体每个子集的结果,保存每组$\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{X}}_{ik}}$${\mathit{\boldsymbol{\widehat \beta }}_{{\rm{oik }}}}$,最后求出初步估计

$ \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)

其次利用分治算法估计方差分量。记${{\mathit{\boldsymbol{\hat U}}}_{ik}} = {\left( {{{\hat U}_{ik1}}, \cdots , {{\hat U}_{ik{m_i}}}} \right)^{\rm{T}}}$,与公式(4)类似,可得

$ \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节可知,式中${\mathit{\boldsymbol{\widehat b}}_i} = {\left( {\sum\limits_{k = 1}^K {\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}} {\mathit{\boldsymbol{Z}}_{ik}}} \right)^{ - 1}}\sum\limits_{k = 1}^K {\mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}} {\mathit{\boldsymbol{\widehat U}}_{ik}}$

根据公式(17)~(18)可知,对应每个个体子集,需计算保存的数据有$\left\{ {\mathit{\boldsymbol{\widehat U}}_{ik}^{\rm{T}}{{\mathit{\boldsymbol{\widehat U}}}_{ik}}, \mathit{\boldsymbol{\widehat U}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{H}}_{ik}}{{\mathit{\boldsymbol{\widehat U}}}_{ik}}, \mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{Z}}_{ik}}} \right., \mathit{\boldsymbol{Z}}_{ik}^{\rm{T}}{\mathit{\boldsymbol{\widehat U}}_{ik}}\} $

最后计算模型系数β的加权估计。根据方差分量的计算公式(18)~(19),第i个个体的第k个子数据模型的协方差矩阵的估计为${\mathit{\boldsymbol{\widehat V}}_{ik}} = {{\bf{Z}}_{ik}}\widehat {{\bf{\Sigma }}Z}_{ik}^{\rm{T}} + {{\mathit{\boldsymbol{\hat \sigma }}}^2}{\mathit{\boldsymbol{I}}_{{m_i}}}$,对每个个体的每个子集运用加权最小二乘估计有

$ \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)

由此可知,对每个个体的每个子集保留数据$\mathit{\boldsymbol{X}}_{ik}^{\rm{T}}\mathit{\boldsymbol{\widehat V}}_{ik}^{ - 1}{\mathit{\boldsymbol{X}}_{ik}}$${\mathit{\boldsymbol{\widehat \beta }}_{ik}}$,得到

$ \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)
3 数值模拟

为了将本文估计算法与极大似然法[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)TXij=(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),且Xij1Xij2Xij3Xij4Xij5相互独立;Zij=(Zij1, Zij2)TZij1Zij2都独立同分布于N(0, 2);bi=(bi1, bi2)Tbi1bi2都独立同分布于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所示。

下载CSV 表 1 情形(1)下两种算法的计算时间及估计精度比较 Table 1 Comparison of calculation time and estimated accuracy of two algorithms under case (1)

表 1可知,极大似然法与本文算法估计出的模型参数的MSE值都较小,说明参数的估计效果较好,估计精确度较高。另外,与极大似然法相比,本文算法估计参数可以在几乎不降低估计精确度的同时将运行效率提高4倍,说明该算法可以大幅减少计算时间,提高计算速度。

3.2.2 数据集块数对计算时间的影响

为了研究计算时间与子数据集块数之间的关系,以说明利用分治算法进行分块计算的必要性,分别将两种样本量下的数据分成不同的数据集块数,记录估计参数所用时间,结果如表 23所示。

下载CSV 表 2 样本①计算时间与数据集块数的关系 Table 2 The relationship between the calculation time and the number of data sets in sample ①
下载CSV 表 3 样本②计算时间与数据集块数的关系 Table 3 The relationship between the calculation time and the number of data sets in sample ②

表 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,再次证明了本文算法的高效性。

下载CSV 表 4 情形(2)下计算时间与数据集块数的关系 Table 4 The relationship between the calculation time and the number of data sets in case (2)
4 实证分析

选取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

下载CSV 表 5 β的估计值和标准差 Table 5 Estimated value and standard deviation of β

表 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.