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