R语言mgcv包详解：广义可加模型（GAM）从入门到实践

R语言mgcv包详解：广义可加模型（GAM）从入门到实践

广义可加模型（GAM）最早由 Hastie 等人提出，是广义线性模型（GLM）的一种非参数扩展形式，与其他统计学模型相比有很多独特的优点，如能够灵活处理 解释变量与响应变量之间的非线性关系；没有规定预测变量的形式，连续型变量和离散型变量均适用；有多种连接函数可选，因此能够对多种数据分布类型（如正态分布、 泊松分布、二项分布等）的响应变量进行处理，应用范围更广。GAM 是最常用的一种描述变量间非线性关系的模型，常用于分析生物与环境之间的关系， 其一般表达式为：

式中，𝑔()为连接函数，𝑓𝑖 为解释变量的平滑函数，𝑌为响应变量，𝑥𝑖 为解释变量， 𝛼截距项，n 为解释变量的数量，ε为残差。

R语言中常用于实现GAM模型的包为mgcv。其主要参数如下：

• formula: GAM的公式，指定了响应变量和解释变量之间的关系。可以使用s()函数来指定光滑项。
• family: 指定模型的分布族。例如，gaussian表示高斯分布，binomial表示二项分布，poisson表示泊松分布等。
• data: 包含数据的data.frame对象。如果未指定，gam将使用当前的环境中的数据。
• weights: 用于加权回归的权重向量。
• subset: 用于指定用于模型拟合的数据子集。
• na.action: 指定如何处理缺失值。默认情况下，gam会删除包含缺失值的观测。
• offset: 偏移向量，用于在模型中添加一个线性项，但不对其进行惩罚。
• method: 指定用于选择光滑参数的方法。例如，GCV.Cp表示使用广义交叉验证。
• optimizer: 指定优化算法。"outer"表示使用外部循环优化光滑参数和线性参数，"newton"表示使用牛顿方法。
• control: 一个包含控制拟合过程参数的列表。例如，可以设置收敛标准和迭代次数。
• scale: 用于对响应变量进行缩放的比例因子。默认为0，表示自动选择。
• select: 逻辑值，指示是否自动选择光滑参数。如果为TRUE，gam将尝试找到最佳的光滑参数。
• knots: 用于指定样条结点的位置。
• sp: 光滑参数的初始值。
• min.sp: 光滑参数的最小值。
• H: 用于指定先验光滑度的矩阵。
• gamma: 用于控制光滑样条惩罚的强度。
• fit: 逻辑值，指示是否立即拟合模型。如果为FALSE，将返回一个未拟合的模型对象。
• paraPen: 用于指定线性项的先验惩罚。
• G: 用于指定分组变量，以便为不同的组拟合不同的光滑参数。
• in.out: 用于指定内部和外部光滑参数的分组。
• drop.unused.levels: 逻辑值，指示是否从模型中删除未使用的因子水平。
• drop.intercept: 逻辑值，指示是否从模型中删除截距项。
• nei: 用于指定邻居矩阵，用于空间光滑。
• discrete: 逻辑值，指示是否将某些变量视为离散变量。

接下来，通过一系列代码示例详细介绍mgcv包的使用方法：

模拟数据（部分）

library(mgcv)
set.seed(2)
dat <- gamSim(1,n=400,dist="normal",scale=2)

拟合GAM模型

b <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3),data=dat)
summary(b)
plot(b,pages=1)
gam.check(b)

自定义光滑参数

G <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3),data=dat,sp=c(0.001,0.01,0.1,1))
plot(G)

使用REML方法

b <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3),data=dat,method="REML")
plot(b)

添加交互项

b <- gam(y ~ te(x0,x1) + s(x2) + s(x3),data=dat)
AIC(b)

平滑ANOVA分解

b <- gam(y ~ s(x0) + s(x1) + s(x2) + ti(x1,x2),data=dat)
AIC(b)

大数据集处理

set.seed(0)
dat <- gamSim(1,n=20000,scale=15)
b6 <- gam(y ~ s(x0,bs="cr",k=100) + s(x1,bs="cr",k=100) + s(x2,bs="cr",k=100) + s(x3,bs="cr",k=100),data=dat,knots=list(x0=seq(0,1,length=100),x1=seq(0,1,length=100),x2=seq(0,1,length=100),x3=seq(0,1,length=100)))
vis.gam(b6,view=c("x0","x1"),plot.type="contour")

b7 <- gam(y ~ s(x0,bs="cr",k=100) + s(x1,bs="cr",k=100) + s(x2,bs="cr",k=100) + s(x3,bs="cr",k=100),data=dat,xt=c(100,100,100,100),seed=123)
vis.gam(b7,view=c("x0","x1"),plot.type="contour")

这些代码片段展示了如何使用mgcv包中的gam函数来拟合各种类型的GAM模型，包括2D平滑、大数据集的处理、用户定义的结点以及不同的可视化技术。通过这些例子，可以学习如何使用mgcv来处理和分析数据，以及如何诊断和调整GAM模型。