R bam 适用于非常大数据集的广义加性模型

说明

将广义加性模型 (GAM) 拟合到非常大的数据集，术语“GAM”包括任何二次惩罚 GLM(也可以使用 family.mgcv 中列出的扩展系列)。模型项的平滑程度作为拟合的一部分进行估计。在使用中，该函数与 gam 非常相似，只是数值方法是为包含数万个以上数据的数据集而设计的(参见 Wood、Goude 和 Shaw，2015)。 bam 的优点是内存占用比 gam 低得多，但对于大型数据集来说，它也可以更快。 bam 还可以在 parallel 包设置的集群上进行计算。

discrete==TRUE 方法提供了另一种拟合方法(Wood 等人，2017 年；Li 和 Wood，2019 年)。在这种情况下，使用基于协变量值离散化和 C 代码级并行化(由 nthreads 参数而不是 cluster 参数控制)的方法。这扩展了实用的数据集和模型大小。

用法

bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
    na.action=na.omit, offset=NULL,method="fREML",control=list(),
    select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL,
    paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE,
    cluster=NULL,nthreads=1,gc.level=0,use.chol=FALSE,samfrac=1,
    coef=NULL,drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,
    in.out=NULL,...)

参数

formula	GAM 公式(请参阅 formula.gam 和 gam.models )。这与 GLM 的公式完全相同，只是可以将平滑项 s 和 te 添加到右侧以指定线性预测器依赖于预测器的平滑函数(或这些函数的线性函数)。
family	这是一个族对象，指定在拟合等中使用的分布和链接。有关更多详细信息，请参阅glm 和family。也可以使用 family.mgcv 中列出的扩展系列。
data	包含模型响应变量和公式所需的协变量的 DataFrame 或列表。默认情况下，变量取自 environment(formula) ：通常是调用 gam 的环境。
weights	数据对对数似然贡献的先验权重。请注意，例如，权重 2 相当于对完全相同的观察进行两次。如果您想重新加权每个数据的贡献而不改变对数似然的总体大小，那么您应该对权重进行归一化(例如 weights <- weights/mean(weights) )。
subset	一个可选向量，指定要在拟合过程中使用的观测子集。
na.action	一个函数，指示当数据包含“NA”时应该发生什么。默认值由 ‘options’ 的“na.action”设置设置，如果未设置，则为“na.fail”。 “factory-fresh” 默认为“na.omit”。
offset	可用于提供模型偏移以用于拟合。请注意，在预测时，此偏移量将始终被完全忽略，这与 formula 中包含的偏移量不同(这用于符合 lm 和 glm 的行为)。
method	平滑参数估计方法。 "GCV.Cp" 对于未知尺度参数使用 GCV，对于已知尺度使用 Mallows' Cp/UBRE/AIC。 "GACV.Cp" 等效，但使用 GACV 代替 GCV。 "REML" 用于 REML 估计，包括未知尺度，"P-REML" 用于 REML 估计，但使用 Pearson 尺度估计。 "ML" 和 "P-ML" 类似，但使用最大似然代替 REML。默认"fREML" 使用快速 REML 计算。
control	用于替换 gam.control 返回的默认值的拟合控制参数列表。任何未提供的控制参数都保留其默认值。
select	是否应该将选择惩罚添加到平滑效果中，以便原则上可以将它们从模型中惩罚掉？请参阅gamma 以增加惩罚。平滑的副作用不再具有固定效应分量(从贝叶斯角度来看是不正确的先验)，允许对具有相同固定效应结构的模型进行 REML 比较。
scale	如果这是正数，则将其视为已知的尺度参数。负值表示比例参数未知。 0 表示泊松和二项式的尺度参数为 1，否则未知。请注意，对于泊松和二项式情况，(RE)ML 方法只能使用尺度参数 1。
gamma	增加到 1 以上以强制更平滑的配合。 gamma 用于乘以 GCV/UBRE/AIC 分数中的有效自由度(因此 log(n)/2 类似于 BIC)。 n/gamma可以被视为有效样本量，这使得它可以在RE/ML平滑参数估计中发挥类似的作用。
knots	这是一个可选列表，包含用户指定的用于基础构造的结值。对于大多数底座，用户只需提供要使用的结，该结必须与提供的 k 值匹配(请注意，结的数量并不总是只是 k )。请参阅 tprs 了解 "tp"/"ts" 情况下发生的情况。不同的项可以使用不同数量的结，除非它们共享协变量。
sp	此处可以提供平滑参数向量。必须按照平滑项在模型公式中出现的顺序提供平滑参数。负元素表示应该估计参数，因此固定参数和估计参数的混合是可能的。如果平滑共享平滑参数，则 length(sp) 必须对应于底层平滑参数的数量。请注意，discrete=TRUE 可能会导致张量积平滑中的变量重新排序以提高效率，并且 sp 必须按重新排序的顺序提供。
min.sp	可以为平滑参数提供下限。请注意，如果使用此选项，则返回对象中的平滑参数 full.sp 需要添加到此处提供的内容中，以获得实际乘以惩罚的平滑参数。 length(min.sp) 应始终与惩罚总数相同(因此，如果平滑共享平滑参数，它可能比 sp 长)。
paraPen	可选列表，指定应用于参数模型项的任何惩罚。 gam.models 解释更多。
chunk.size	模型矩阵是按这种大小的块创建的，而不是整个形成的。如果 chunk.size < 4*p 则重置为 4*p，其中 p 是系数的数量。
rho	AR1 误差模型可用于 Gaussian-identity 链路模型的残差(基于数据帧顺序)。这是 AR1 相关参数。如果非零，则在 std.rsd 中返回标准化残差(在正确模型下近似不相关)。当 discrete=TRUE 时也可与其他模型一起使用，在这种情况下，AR 模型应用于工作残差并对应于 GEE 近似。
AR.start	与数据长度相同的逻辑变量，TRUE 首先观察 AR1 相关性的独立部分。 DataFrame 中的第一次观察不需要这个。如果NULL则AR1相关性没有中断。
discrete	使用method="fREML"，出于存储和效率原因可以离散化协变量。如果 discrete 是 TRUE ，每个平滑项的数字或数字向量，则会发生离散化。如果提供数字，它们会给出离散化箱的数量。参数项使用指定的最大数量。
cluster	bam 可以使用 parallel 包中的 parLapply 并行计算计算主导的 QR 分解，如果它提供了执行此操作的集群(这里的集群可以是单台机器的一些核心)。查看详细信息和示例代码。
nthreads	用于非集群计算的线程数(例如，组合来自集群节点的结果)。如果 NA 设置为 max(1,length(cluster)) 。查看具体信息。
gc.level	为了减少内存占用，经常调用垃圾Collector会有所帮助，但这需要大量时间。将其设置为零意味着垃圾收集仅在 R 决定应该进行时发生。设置为 2 会频繁进行垃圾回收。 1 介于两者之间。问题不像以前那么严重了。
use.chol	默认情况下，bam 使用非常稳定的 QR 更新方法来获取模型矩阵的 QR 分解。对于条件良好的模型，另一种方法是累积模型矩阵的叉积，然后在最后找到其 Choleski 分解。这在计算上更加有效。
samfrac	对于非常大的样本量广义加性模型，可以通过首先将模型拟合到数据的随机样本，然后使用结果提供起始值来减少模型拟合所需的迭代次数。这种初始拟合是在不严格的收敛公差下运行的，因此通常成本非常低。 samfrac 是要使用的采样分数。 0.1 通常是合理的。
coef	模型系数的初始值
drop.unused.levels	默认情况下，未使用的级别会在拟合之前从因子中删除。对于某些涉及因子变量的平滑，您可能需要将其关闭。仅当您知道自己在做什么时才这样做。
G	如果不是 NULL 则这应该是先前使用 fit=FALSE 调用 bam 返回的对象。导致除 sp , chunk.size , gamma , nthreads , cluster , rho , gc.level , samfrac , use.chol , method 和 scale 之外的所有其他参数被忽略(如果 >0)。
fit	如果 FALSE 则模型设置为拟合但不估计，并返回一个对象，适合作为 G 参数传递给 bam 。
drop.intercept	设置为 TRUE 以强制模型在参数模型部分中实际上没有常量，即使存在因子变量也是如此。
in.out	如果提供，则这是初始值的两个项目列表。 sp 是初始平滑参数估计，scale 是初始尺度参数估计(如果家庭没有，则设置为 1)。
...	传递的进一步论点，例如到gam.fit(例如mustart)。

细节

当 discrete=FALSE 、bam 进行操作时，首先使用数据的代表性子样本设置平滑的基本特征。然后使用 predict.gam 分块构建模型矩阵。对于每个块，来自整个模型矩阵的 QR 分解的因子 R 以及 Q'y 都会更新。以及 y 的平方和。在块处理结束时，进行拟合，而无需形成整个模型矩阵。

在一般情况下，在 PIRLS 的每个步骤中，对加权模型矩阵和加权伪数据使用相同的技巧。在每个阶段(面向性能的迭代)对工作模型执行平滑度选择，以保持较小的内存占用。在基于 GCV 或 Cp/UBRE/AIC 的模型选择的情况下，这很容易证明，而对于 REML/ML，则通过 Q'z 的渐近多元正态性来证明，其中 z 是 IRLS 伪数据。

有关完整方法的详细信息，请参阅 Wood、Goude 和 Shaw (2015)。

请注意，POI 不如 gam 使用的默认嵌套迭代那么稳定，但对于非常大、信息丰富的数据集，这不太重要。

还要注意，如果使用昂贵的基础，则可能将大部分计算时间花费在基础评估上。实际上，这意味着应避免默认的 "tp" 基础：几乎任何其他基础(例如 "cr" 或 "ps" )都可以在一维情况下使用，并且张量积平滑( te )通常要小得多在多维情况下代价高昂。

如果提供 cluster 作为使用 parallel 包中的 makeCluster (或 makeForkCluster )设置的集群，则使用该集群并行执行模型矩阵的速率限制 QR 分解。请注意，速度通常不会那么快。在多核计算机上，通常最好将集群大小设置为物理核心的数量，该数量通常小于 detectCores 报告的数量。使用超过物理核心数量可能会导致根本没有加速(甚至减慢)。请注意，高度并行的 BLAS 可能会抵消使用核心集群的所有优势。并行计算当然比串行计算需要更多的内存。请参阅示例。

当discrete=TRUE时，协变量数据首先被离散化。离散化是在逐个平滑的基础上进行的，或者在张量积平滑(或任何可以表示为这样的平滑，例如随机效应)的情况下，分别针对每个边平滑进行。然后以离散值评估所需的样条基数，并将其与指示它们与哪个原始观测值相关的索引向量一起存储。拟合是通过在每个迭代工作模型上使用 REML 平滑参数选择的面向性能的迭代/PQL 版本(对于默认方法)。迭代基于 REML 分数的导数，无需计算分数本身，从而允许将昂贵的计算减少为每次迭代一个并行块 Cholesky 分解(加上两个成本相等但易于并行化的基本操作)。与标准 POI/PQL 不同，每一步仅进行工作模型平滑参数更新的一步(而不是迭代到每个工作模型的一组最佳平滑参数)。在每个步骤中，都需要模型矩阵的加权模型矩阵叉积 - 这是根据在离散协变量值处评估的预先计算的基函数有效地计算出来的。使用张量积项进行高效计算意味着张量积内的某些项可以重新排序以获得最大效率。有关完整详细信息，请参阅 Wood 等人 (2017) 以及 Li 和 Wood (2019)。

使用nthreads自变量控制discrete=TRUE并行计算时。对于此方法，不使用集群计算，并且不需要 parallel 包。请注意，并行化的实际速度取决于安装的 BLAS 和您的硬件。使用(R 默认)参考 BLAS 使用多个线程可以产生显著的差异，但是使用单线程调整的 BLAS(例如 openblas)，效果不太明显(因为缓存使用通常针对一个线程进行优化，然后是次优的)给几个人的)。然而，无论您使用多少线程，调整后的 BLAS 通常比使用参考 BLAS 快得多。如果您安装了多线程 BLAS，那么您应该将 nthreads 保留为 1，因为从多个线程调用多线程 BLAS 通常会减慢速度：唯一的例外是您可能选择形成离散矩阵叉积(拟合例程中的主要成本)以多线程方式，但使用单线程代码进行其他计算：这可以通过例如nthreads=c(2,1) ，它将使用 2 个线程用于离散内积，1 个线程用于大多数调用 BLAS 的代码。这并不是说多线程性能经常令人失望的根本原因是大多数计算机的内存带宽受到严重限制，而不是翻转率受到限制。很难足够快地将数据传输到一个核心，更不用说尝试同时将数据传输到多个核心了。

discrete=TRUE 通常会产生与没有离散化的方法相同的结果，因为协变量通常只采用适度数量的离散值，因此离散化过程中根本不涉及任何近似。即使涉及一些近似，差异通常也很小，因为算法会尽可能地进行边际离散。例如，张量积平滑的每个边际都是单独离散化的，而不是离散化到协变量值的网格上(对于等效的各向同性平滑，我们必须离散化到网格上)。边际方法允许非常精细的尺度离散化，因此近似误差非常低。请注意，当使用平滑 id 机制链接平滑参数时，离散方法无法强制链接的碱基相同，因此与非离散方法的一些差异将会很明显。

也可以使用 family.mgcv 中给出的扩展系列。这些额外参数是通过最大化惩罚似然来估计的，而不是像 gam 中的限制边际似然。因此估计值可能与 gam 返回的值略有不同。估计是通过牛顿迭代来找到额外参数(例如负二项式的 theta 参数或缩放 t 的自由度和尺度)来完成的，在拟合过程的每次迭代中给定模型系数时最大化对数似然。

值

类 "gam" 的对象，如 gamObject 中所述。

警告

如果使用默认的 "tp" 基础，则例程可能会比最佳运行速度慢。

对于相同的数据集，此例程的稳定性不如‘gam’。

使用 discrete=TRUE 时，可以有效计算 te 项，但 t2 则不然。

例子

library(mgcv)
## See help("mgcv-parallel") for using bam in parallel

## Sample sizes are small for fast run times.

set.seed(3)
dat <- gamSim(1,n=25000,dist="normal",scale=20)
bs <- "cr";k <- 12
b <- bam(y ~ s(x0,bs=bs)+s(x1,bs=bs)+s(x2,bs=bs,k=k)+
           s(x3,bs=bs),data=dat)
summary(b)
plot(b,pages=1,rug=FALSE)  ## plot smooths, but not rug
plot(b,pages=1,rug=FALSE,seWithMean=TRUE) ## `with intercept' CIs

 
ba <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+
            s(x3,bs=bs,k=k),data=dat,method="GCV.Cp") ## use GCV
summary(ba)

## A Poisson example...

k <- 15
dat <- gamSim(1,n=21000,dist="poisson",scale=.1)

system.time(b1 <- bam(y ~ s(x0,bs=bs)+s(x1,bs=bs)+s(x2,bs=bs,k=k),
            data=dat,family=poisson()))
b1

## Similar using faster discrete method...

 
system.time(b2 <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+
            s(x3,bs=bs,k=k),data=dat,family=poisson(),discrete=TRUE))
b2

作者

Simon N. Wood simon.wood@r-project.org

参考

Wood, S.N., Goude, Y. & Shaw S. (2015) Generalized additive models for large datasets. Journal of the Royal Statistical Society, Series C 64(1): 139-155. doi:10.1111/rssc.12068

Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association. 112(519):1199-1210 doi:10.1080/01621459.2016.1195744

Li, Z & S.N. Wood (2020) Faster model matrix crossproducts for large generalized linear models with discretized covariates. Statistics and Computing. 30:19-25 doi:10.1007/s11222-019-09864-2

也可以看看

mgcv.parallel , mgcv-package , gamObject , gam.models , smooth.terms , linear.functional.terms , s , te predict.gam、plot.gam、summary.gam、gam.side、gam.selection、gam.control gam.check、linear.functional.terms negbin 、magic、vis.gam