广义可加模型(GAM)作为广义线性模型(GLM)的扩展,摒弃了对变量间预设参数关系的依赖,从而能够敏锐地捕捉数据间的复杂关联。它不仅具备识别线性关系的能力,更能游刃有余地处理数据间的非线性关系,这使得构建更为贴切的非线性模型成为可能。因此,GAM在气候变化、环境科学、生态学等诸多研究领域中得到了广泛的应用。然而,尽管GAM模型在数据处理和分析中展现出显著的优势,但在其结果分析过程中,一个常见的挑战在于缺乏直观且有效的方法来量化各解释变量对模型拟合度贡献率的大小。这一点在共线性或共曲线性的情况下尤为重要,因为它直接关系到我们如何准确评估各解释变量的相对重要性。因此,寻找简单而有效的途径来定量分析各解释变量的贡献率,成为当前GAM模型应用中亟待解决的问题。
图1. 平均共享方差理论示意图. 原理是将共享方差平均分配给相关的解释变量,每个解释变量的贡献率由独自解释(unique)部分和所分配的平均共享方差的和,优势是所有解释变量贡献率总和刚好等于总拟合度(R2).
为了解决这一问题,Lai et al. (2024) 将“平均共享方差”的思想拓展至GAM模型,研发R包 gam.hp 在线发布于CRAN(https://cran.r-project.org/web/packages/gam.hp/index.html),该包主要作用是定量分解由mgcv包给出的GAM模型的校正R2或离差解释率(explained deviance),提供了在共线性或共曲线性情况下定量评估变量拟合度贡献率及相对重要性的新方法。
文中以openair包提供的英国伦敦空气质量数据mydata为例进行演示。
(1)打开R,运行下面命令,下载安装gam.hp软件包,载入gam.hp和相关软件包。
install.packages("gam.hp")
install.packages("openair")
library(gam.hp)
library(openair)
library(mgcv)
(2)载入mydata数据,建立GAM模型,分析变量共曲线性,评估模型效果。
data("mydata")
mydata <- na.omit(mydata)[,c("o3","ws","nox","co")]
mod1<- gam(o3~s(ws)+s(nox)+s(co),data=mydata)
concurvity(mod1)
## para s(ws) s(nox) s(co)
## worst 5.408694e-20 0.04970121 0.7966771 0.7975821
## observed 5.408694e-20 0.04610312 0.7421644 0.6727814
## estimate 5.408694e-20 0.03262776 0.6461008 0.5952891
summary.gam(mod1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## o3 ~ s(ws) + s(nox) + s(co)
##
## Parametric coefficients:
## Estimate Std.Error t value Pr(>|t|)
## (Intercept) 7.12094 0.02585 275.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ws) 6.568 7.520 439.46 <2e-16 ***
## s(nox) 8.947 8.998 1470.64 <2e-16 ***
## s(co) 7.512 8.355 33.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.505 Deviance explained = 50.5%
## GCV = 28.425 Scale est. = 28.409 n = 42524
(3)定量分解summary.gam(mod1)输出的R-sq.(adj)得到各解释变量的独自效应、平均分配值、单独效应、贡献率等,绘图展示单独效应。
hp1<- gam.hp(mod1,type = 'adjR2')
hp1
## $adjusted.R2
## [1] 0.5051936
##
## $hierarchical.partitioning
## Unique Average.share Individual I.perc(%)
## s(ws) 0.0384 0.0029 0.0413 8.17
## s(nox) 0.1524 0.1579 0.3103 61.42
## s(co) 0.0033 0.1503 0.1536 30.40
##
## $variables
## [1] "s(ws) " " s(nox) " " s(co)"
##
## $type
## [1] "hierarchical.partitioning"
##
## attr(,"class")
## [1] "gamhp"
plot(hp1)
(4)定量分解summary.gam(mod1)输出的Deviance explained得到各解释变量的独自效应、平均分配值、单独效应、贡献率等,绘图展示单独效应。
hp2<- gam.hp(mod1,type = 'dev')
hp2
## $Explained.deviance
## [1] 0.5054616
##
## $hierarchical.partitioning
## Unique Average.share Individual I.perc(%)
## s(ws) 0.0385 0.0029 0.0414 8.19
## s(nox) 0.1525 0.1578 0.3103 61.40
## s(co) 0.0034 0.1503 0.1537 30.41
##
## $variables
## [1] "s(ws) " " s(nox) " " s(co)"
##
## $type
## [1] "hierarchical.partitioning"
##
## attr(,"class")
## [1] "gamhp"
plot(hp2)
欢迎大家尝试gam.hp包,也欢迎大家对gam.hp包多提改进意见(请直接联系lai@njfu.edu.cn), 以便使它成为更有用的工具。使用此包请正确引用如下文献:
Jiangshan Lai, Jing Tang, Tingyuan Li, Aiying Zhang, Lingfeng Mao. 2024. Evaluating the relative importance of predictors in Generalized Additive Models using the gam.hp R package. Plant Diversity. https://doi.org/10.1016/j.pld.2024.06.002.
我们也邀请广大研究者使用另外两个包rdacca.hp和glmm.hp分别用于多元数据的典范分析,广义线性混合效应模型(包括常规的多元回归)解释变量相对重要性的评估功能。然而,正如在此文章讨论部分强调的一样,共线性问题在统计学中是一个固有的难题,至今尚未有完美的解决方案。 “平均共享方差”概念旨在为这一问题提供一种简洁且实用的方法而已,虽然它已经在很多领域得到了广泛应用,但并非万能之策。因此,鼓励大家在使用这个方法的时候,既要充分利用它们的优势,也要保持审慎的态度,结合自身的专业知识进行判断。请记住,统计方法只是工具,真正的科学发现应基于对专业知识深入的理解和全面的分析,切勿被统计结果所“绑架”,而是应该让数据和分析服务于我们的研究目的。