适用于大数据集的广义可加模型①②，本文主要内容关键词为：适用于论文,广义论文,模型论文,数据论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
如今，回归问题中的观测值个数经常会达到几万个甚至是几百万个。处理如此庞大的数据，有些问题需要建立全新的模型，但也有些问题，只要计算机可以完成运算，现有的模型依然可行。本文的研究旨在，在使用适中的计算机（modest computer hardware）的情况下，如何利用广义可加模型（Generalized Additive Model，GAM）来处理大数据集，并且在模型拟合过程中能够得到平滑参数的估计值。 本文是在处理一个实际问题时受到的启发。图1显示了自2002年9月1日以来，整个法国每隔半小时的用电量，单位为GW（千兆瓦，十亿瓦特，
瓦）。一直以来，法国能源公司——法国电力公司（Electricit'e de France，EDF）在利用GAMs（选取了一系列变量）对用电量做短期预测上取得了重大成就，特别是基于24小时以前的用电量数据进行的预测。然而，现有的GAM拟合方法不能一次性对全部数据进行计算，而是必须使用48个单独的模型来对每天每隔半小时的数据进行拟合（Pierrot和Goude，2011）。显然，若能只使用一个模型是最好不过的。
尽管对现有的GAM估计方法带来了挑战，但就目前的标准来看，上例中的数据还是相对较小的。比如，已有非常多的统计模型研究空气污染与呼吸死亡率之间的关系。方法之一便是采用GAMs，将死亡率分解成为随时间平滑变动的背景成分和污染因素（光滑或者线性的），并同时假定观察到的死亡数服从泊松分布。Peng和Welty（2004）搜集了美国108个城市近5000天3个年龄段居民每天的污染死亡率数据。Wood（2006）单独对芝加哥的研究表明，臭氧层与温度之间存在非常强的交互作用，如果这种研究能够重复，将会具有重要的意义。然而，这种作用仅对少数日子的数据敏感，因此，应该基于Peng和Welty（2004）搜集的数据集，对美国其他城市臭氧与温度的交互作用进行验证。理想情况下，我们希望能够同时对数据集中120万个观测值进行GAM拟合。显然，这样的拟合已经远远超出了现有GAM估计方法的运算能力，但是基于本文的方法，利用具有适度计算机硬件的电脑就可以进行拟合。对于包含几万个观测值的数据集，通用的GAMs及相关模型（如Wood，2011）的拟合方法都表现出非常强的有效性和稳健性。但对于更大的数据集，这些方法就因需要占用太大的内存而最终导致拟合的失效。这里的困难，简单来说，就是模型矩阵变得太大：记n、p分别为模型矩阵的行数和列数，M是平滑参数的个数，那么GAM拟合方法所需的存储空间通常为
，显然，这处理起来会非常困难。以下我们将给出一种简单的处理方法，该方法可以通过更新模型矩阵的因子分解来避免对GAM模型的整个模型矩阵进行计算。更重要的是，本文给出了在该情形下如何调整平滑参数估计的方法。 本文考虑的一般模型形式为：
其中，
是n个观测值中的一个单变量响应变量，服从指数族分布（或者至少是知道均值方差关系在一定的尺度参数内），i=1，2，…，n，g是已知的平滑单调连接函数，A是n行的模型矩阵，θ是由未知参数组成的列向量，
是已知的线性函数，
是未知的单变量或多变量平滑函数，并且平滑度也未知。对于每一个
，都有一个测度偏离平滑度的函数
（f）。 式（1）最常见的例子就是GAM，此时
是评估函数（Hastie和Tibshirani，1986，1990）。其他的例子还有，变系数模型（varying-coefficient models）、函数性广义线性函数（functional generalied linear models）（例如Marx和Eilers，1999）和结构可加回归模型（structured additive regression models）（例如Fahrmeir等，2004）等。当然，对于不同的模型形式有不同的估计方法。本文主要关注
为中等秩惩罚回归样条函数的情形（例如Parker和Rice，1985；Eilers和Marx，1996）。在这种情形下，采用惩罚迭代再加权最小二乘法进行拟合，采用广义交叉验证方法（GCV）、限制极大似然（RMEL）或类似方法（参考Wood，2011）进行平滑参数选择，可以从整个模型类中获得非常有效的计算方法（例如Wood，2000）。 二、高斯恒等连接情形 （一）基本情形 首先考虑，
服从独立且方差为
的正态分布，g是恒等函数。
是通过线性基（如B-样条基，或者薄板回归样条基）生成的，
是关于基系数的二次方形式。在这种情形下，模型的期望可以改写为： E（y）=Xβ （2） 其中，X为包括A和评估基函数的n×p阶模型矩阵，β包括θ和全部基系数。我们假设p＜n，因为这种方法几乎只有在这种情况下才有意义。对β的估计可以最小化以下方程：
然而，λ的估计要相对复杂一些。方法之一是最小化预测误差的GCV方法，即最小化关于平滑参数的以下方程：
这里的问题是只要我们有了R、f和
，就有了进行拟合所需的所有参数，而此时的X就不再起作用了。因此如果我们可以不通过构建X而得到这些变量的值，就可以不用耗费大量的计算机储存空间而对模型进行估计。当使用Mallows
和RMEL方法时，上述情况仍然成立。本文不再赘余。 事实上，可以仅利用任何时候构建的X的小的子矩阵来计算R、f和
具体的方法可以采用QR分解的迭代更新，或者稳定性稍差的Choleski分解等。 （二）相关误差 本文第四部分对短期用电量的预测问题也需要对自相关残差进行建模，通常采用简单的AR（p）相关结构。首先将高斯恒等连接模型式（2）简化为y=Xβ+e，其中，e的协方差矩阵为
∑，∑是自回归AR（p）的相关矩阵。则
的Choleski因子C是带状的，ε=Ce独立同分布于N（0，
）。综上，如果
，则：
该式是式（2）的变形，因此同样可以利用上述方法来估计β。这里唯一需要变化的就是，如果使用REML来估计ρ，必须对REML的对数值进行修正，以使的可以被C转化，但考虑到C是三角矩阵，因而可以很容易获得其行列式的对数。在计算上，可以使用简单的一维搜索来得到ρ，对于每一个ρ都需要对模型重新进行拟合。注意到，C是带状结构的，并且可以不通过求解
得到，这使得AR（p）在计算机上实现起来更方便：相比对整个矩阵进行昂贵的乘积运算，
仅涉及到对X，y相邻行的加权差分，这显然能节省更多的存储空间。 三、广义可加模型拟合 在广义的情形下，未知函数和他们的惩罚样条都可以表示成上述简单高斯恒等连接的情形。唯一的改变是，模型变成了一个过度参数化的广义线性模型：
然后迭代以下步骤至收敛：
对于中等大小的数据集，对于每个λ通过PIRLS方法迭代至收敛相对比较容易，其中λ可以使用广义的GCV、
、拉普拉斯近似REML等方法估计得到（Wood，2008，2011）。但对于大数据集，使用这些方法就需要大X好几倍的存储空间去有效地计算平滑选择准则的导数。为了避免耗费如此巨大的存储空间，我们可以使用一种更早期的方法来替代，即最早始于Gu（1992）的研究。如此，在PIRLS算法的每一步便能简单地应用GCV，
，REML等方法对线性模型的平滑参数进行选择。Gu（1992，2002）将该方法命名为“绩效导向的迭代”，它与Breslow和Clayton（1993）的惩罚伪似然估计很相似。该方法通常是收敛的，即使不能保证必然收敛，使收敛出现问题的不利条件也会随着n的增大不断减轻。 （一）对于大数据集的绩效导向迭代 对于大数据集的绩效导向迭代，可以通过在PIRLS算法的每一步对WX矩阵使用QR-更新方法实施。这就意味着，对于模型矩阵（子矩阵）的计算在每一步都需要重复，但是这些计算所需的存储空间是O（np）而不是QR-分解方法的
，因此，对那些平滑基来说，运算成本就不再重要。基本算法如下。 1.初始化。
2.迭代。 第1步：令
=D，R是一个p阶0矩阵，f是一个0向量，D=0，r=0。 第2步：重复以下步骤（a）～（f），k=1，…，M。
（f）构造：
第3步：令
。 第4步：若q＞0，通过比较当前异常D与之前的异常
来检验是否收敛。若达到收敛，则步骤停止（若q值已超出了预先设定的范围，也应该停止）。
第6步：令q←q+1。 收敛后最终得到的
和λ即是系数和平滑参数的估计值。进一步的推断常基于贝叶斯近似
其中，
可以通过REML的最优化过程估计得到，或者采用以下方程：
参考Wood（2006）以获得进一步的信息。需要说明的是，如附录B所述，第2步可以简化。 3.初始值的二次抽样。 使用全部的数据集来运行上述PIRLS算法的前几步，显然，在计算上是非常浪费的。这无异于耗费大量的工作去拟合一个本是错误的模型。因此，通常情况下，首先从数据集中随机抽取5%～10%的数据来构成子样本，然后用该子样本估计得到的β和λ作为初始值去对整个数据集进行拟合。在实际操作中，这样处理后一般能节省PIRLS算法拟合整个数据集所需的1～2步的计算时间。 （二）平滑选择步骤的说明 在PIRLS迭代中的每一步：这些准则所需的假设对模型来说也是适用的，因而在将GCV和
应用到模型时不需要特别的理由。
Reissue和Ogden（2009）的研究显示，相比
不太可能有多个局部最小值（或者说与
紧密相关）。在绩效导向的迭代背景下，这意味着REML促进了迭代的收敛，因为它减小了迭代在多个最优值之间的循环。 （三）关于p＜n的假设 采用GAM处理大数据集，计算可行性取决于使用降秩平滑使p远小于n的能力。但一个明显的问题是，假设p的增长速度比n的增长速度慢得多，是否是合理的。显然，如果平滑项是用来作为去除残差自相关性的随机部分，或者在有新数据加入延长了时间轴时用来平滑时间，依然假设p的增长速度远小于n的增长速度便不再合理。否则，所有的理论依据就会变为，平滑基仅需要随着样本的增长而缓慢地增长：例如，考虑有均衡间隔结的三次回归样条情形。三次样条的偏差平方的平均值为
，其中h是结间距，k是结点数（例如de Boor，1978）。由基本回归理论，样条方差的平均值为O（k/n），其中n是数据个数。为了避免偏差与方差之间的互相影响，并且考虑到当n→∞时的次优均方误差，则我们需要选取与偏差平方和方差相同阶的k，比如
。考虑到在任何有限样本集中，我们会选择降低均方误差（相对于纯回归）的惩罚度。因此，考虑惩罚的回归同样可以采用以上比率，即
（当然，其他比率也是可以的）。事实上，Gu和Kim（2002）研究显示，带有惩罚的基本维度应该是一个
的规模，即与1000个数据的数据集相比，100万个数据所需的系数个数仅仅是它的5倍左右。 四、短期用电量预测 EDF成功建立了日用电量预测模型，该模型基于分割数据，即每半小时的数据，对48个子集中的每一个都拟合GAM。尽管以上方法能够用来预测，但是采用单独的48个模型存在以下3点不足： （a）未能有效利用数据，即相邻的每半小时数据之间的相关性未利用。 （b）不好解释，每半小时模型之间的连续性未能体现出来，这不符合实际。 （c）进行预测的模型应必须具有统计上的稳定性，与此同时，模型的预测目的也要求必须采用如GCV的预测误差准则来检验这种模型的平滑估计。但我们都知道，GCV及相关准则会产生过拟合，并且会随着样本量的减小而出现一系列问题（例如参考Reiss和Ogden，2009）。尤其是要对数据的48个子数据建立单独的模型，更会增加这一问题的严重性。以致于每次更新模型估计时，都会降低模型的稳定性，并且还会给预测人员增加大量的模型检验工作。因此，对全部数据一次性拟合1个模型就会显著地减小过拟合现象，并且也会减少模型检验的工作量。 本文要建立模型，最主要的目的就是能够使用1个模型来代替48个单独的模型。以下给出该模型的一种形式。 图1的数据来自EDF的数据库，该数据库中还包括气象数据（如来自法国气象局（Meteo France）的提前1天的摄氏温度的预测值，以及云层覆盖率（以覆盖天空的1/8为单位））、日历和资费等信息。尽管用电量的影响因素非常复杂，以致于只有从统计的角度才能完成预测，EDF坚持预测模型必须能够具有符合实际的解释力。在一定程度上，当数据出现在正常拟合模型的数据集外时，这对实际预测异常事件非常有帮助。 已有的预测经验和探索性数据分析显示，基于24小时之前的用电量，24和48小时之前的预测和实际温度、预测云层覆盖率和整年的日期，能够得到很好的预测效果。一般来说，用电量在温度为15℃左右时最低，低于该温度，随着供暖需求的增加，用电量也会随之提高；高于15℃，由于空调的使用，用电量也会有少量增加。需要注意的是，用于预测的温度是整个法国的平均温度，即以用电量加权得到的温度。其他条件相似时，发现云层覆盖率对用电量也表现出非常重要的直接影响，这可能主要是受建筑的被动式太阳能加热和照明的影响。 另一个重要因素是，EDF通过采用财政刺激降低用电大户需求的策略来控制冬季预计用电需求异常高的日子。EDF服务部门提前一天预测采用这种刺激的“特殊资费日”的晚上7点用电需求降低的量。显然，这一预测值同样可以作为一个预测变量。基于以上考虑，以及变量探索性分析的结果和每半小时的单独的48个模型，并且在对模型进行筛选后，本文提出以下模型：
式（6）中的单项在很大程度上基于EDF研究经验中的预期影响因素，在某些情况下，还需要基于先验基础，比如温度效应的设置。然而，也存在大量的结构性假设，是之前未曾有的经验。一个需要考虑的问题是效应该设为相加还是相乘形式呢？这可以通过估计具有和没有对数连接的模型解决。在估计模型后，假设ρ=0，两个模型的
在3位有效数字上是相等的，也即在预测能力上无差异。考虑到相加形式比相乘形式能够更好地处理相关性，因此本文决定采用相加形式。另外一个问题是如何处理季节波动。本文采用年中的循环效应来处理。当然，也有一种更简洁的方法，即变系数方法（Eilers等，2008），其季节波动是用截断傅里叶级数来处理的：
其中，
（t，I）是随时间和时刻变化的平滑函数，它决定了循环周期T的相位和振幅。我们尝试过用这种形式的模型处理季节波动，但是
和预测误差结果显示，该模型比式（6）的预测效果相对要差一些。 为了说明模型的有效性，本文采用截至2008年8月31日的数据进行模型估计，然后对接下来的年份利用该模型进行下一日用电量的预测，并且利用每天的新数据去更新模型估计。考虑到银行假日比较特殊，并且通常排除在正常的自动预测外，因此本文也进行了相应的处理。考虑到模型的预测特性，采用GCV对平滑参数和AR参数进行估计。注意到，在进行初始模型拟合时，需要设置基使得
有能够运行到2009年9月的域，以确保在整个预测过程中，基都是合适的。这设置起来并不难，但如果我们采用基于Choleski的更新方案则可能会带来数值问题，因为初始的
可能会出现亏秩或者近乎亏秩。 图2是初始拟合数据的残差图。可见，较大的预测误差多出现在周一早上；大的异常值多出现在白天；残差在冬天较大，此时的用电量较高。图2（d）展示的是预测残差最后一年（2008年8月31日-2009年，译者注）的数据。可见预测的残差与初始拟合的残差相差不大。AR参数估计值为0.98，图2也显示了将参数设为0的简化模型的预测残差。当残差之间的相关性被忽略时，会导致预测效果变差。这可由表1中的平均绝对百分比误差（mean absolute percentage error，MAPE）和均方根误差（root-mean-square error，RESE）看出。第22页图3（a）直接给出了利用图1（b）中的全部用电量数据拟合模型得到的用电量。
可见，基于AR模型的预测误差相对更低，并且拟合数据集和预测数据集之间的错配也更小。ρ=0时的较大错配，也几乎完全是由忽略相关性的过渡拟合造成的：带有相关性模型的有效自由度相当于不带相关性模型的83%。这也证明了本文在第二部分提出的方法在应用上的重要性。下页图3（b）和图3（c）分别显示了ρ=0模型残差和ρ=0.98标准残差的自相关函数。显然，AR（1）模型带来了很大的提高，但仍有更大的改进空间。式（6）的预测能力不仅可与Pierrot和Goude（2011）的相媲美，且还有前文所述单一模型拟合的三项实际优势，而且可以简化当新数据生成时模型更新的估计过程。为了证明后一点，我们在3.1GH内核、i3540处理器、3.7GB随机存储器、LINUX运行环境计算机（相当于一台至少价值600美元的个人计算机）上，比较了运算时间。运用mgcv软件包里的gam，对与式（6）等价的48个独立模型进行估计，不加AR残差，需要1个小时左右的时间；而用本文提出的方法对式（6）的初始估计只需要不到半个小时的时间，而且还包括对ρ的计算。随后的更新，本文提出的方法仅仅需要不到2分钟；而利用之前的方法，需要对所有的48个模型完全重新拟合。
至此，又会顺理成章地提出这样的问题：既然模型是对很长的一段时期的拟合，那么需不需要每天更新。本文之所以对每天都进行更新是因为，即使是对很长一段时期数据的抽样，也不能保证该样本里包含了所有的情形，并且预测变量也是高度自相关的。这也就是说，如果有几天异常的情形，那么在这种情况下，最近几天数据的信息就可能会构成用电量信息的一个不可忽略部分，进而就会对这些异常协变量邻近平滑函数的估计产生不可忽略的影响。由此，在预测模型的拟合过程中，通常不会舍弃最近的数据。 本文所提出的模型是对现有GAM运算难题的最直接的解决。虽然式（6）在8GB内存的计算机上，利用Wood（2008）的方法也能运算。但Wood的方法的计算速度比本文所提方法的1/10还低，并且需要通过完全的重新拟合来对模型进行更新。在模型开发过程中考虑得更复杂的模型，需要使用更多的基和更复杂的日期类型分类，因而需要更大的内存，而且未来的目标是对一些地区效应进行建模，如果没有本文提出的这种新方法，这种模型是完全不能建立的。 为了说明可解释力，
的估计结果如图4所示。注意到，在一天的任何一段时间，一个工作日的用电量与前一个工作日用电量是相关的线性关系，从均值回归中更能明显地看出。假期则表现出不同的情况，在低用电量时没有很好的线性特征，但当用电量超过50GW时则会表现出明显的线性特征。假期后的工作日，大多是星期一，则表现出与假期用电模式截然相反的特征，在低用电量时最为突出。 综上，本文所建立的方法提高了EDF预测模型的稳定性和解释能力，通过对全部数据同时进行拟合，使得能对自相关性进行建模，当有新数据时，能更有效地进行更新估计。 为了提高EDF预测模型的运算速度和稳定性，本文展示了如何利用GAMs来处理大数据的问题。本文方法最明显的优势在于，它是在现有方法上的相对直接的扩展，但却能在可建模的数据集规模上和某些情况下的拟合速度上都有实质的提高。第一部分中介绍的空气污染的平滑泊松回归的例子，是实际改善最鲜明的一个例子。在那个例子中，如果同时对全部数据进行拟合，仅模型矩阵就需要7GB的存储空间；而利用本文的方法，只需要1GB的存储空间。就目前所知，还未在其他公开出版物上发现像我们处理这样大的空气污染数据集的模型拟合方法。如果存储空间没问题的话，估计所需的时间（如第四部分所述，在普通的计算机上仅需要不到12分钟的时间）也将会比现有的方法快100倍。
除本文所提到的电力领域，遥感、基因技术、金融和信息技术等领域也都正在产生越来越大的数据集。对于有些问题，处理这些数据集需要全新的分析方法；但也有些问题，已有的方法仍是有效的，解决问题的关键所在是探索可以在计算机上运算的方法。本文所提出的方法就是为了解决这一问题，因此希望能够得到广泛的应用，当然不限于以上所提及的领域。 本文所讨论的方法是用R中的mgcv包（Wood，2009）bam函数运行的。 ①本文翻译自Journal of the Royal Statistical Society:Applied Statistics Series C,Appl.Statist.（2015）,64,Part1,pp.139-155。译文有不恰当之处，还请批评指正。 ②刊出时有删节——编者注。
标签：大数据论文; 预测模型论文; 参数估计论文; 残差分析论文; 数据拟合论文; 空间数据论文; 线性拟合论文; 迭代模型论文; 迭代计算论文; 矩阵分解论文; 时间计算论文;
