想发高分文章，不知道有什么冷门且优秀的算法，来看看广义加性模型吧！

2023-07-28 14:19 作者:尔云间 0人读过 | 我要投稿

相信很多小伙伴们都有过类似的苦恼，总感觉自己了解的可以结合生信分析的机器学习算法就那几种，用的人太多了！没有什么吸引力。那么有没有什么冷门点的而且可以结合生信分析的算法呢？

今天，小云就给大家推荐一个，广义加性模型（generalized additive models， GAM）。

GBM是回归家族的一个扩展，且是最强大的统计模型之一，非常自由灵活，可以为任何回归问题建模！传统线性模型简单、直观、便于理解，但是，在现实生活中，变量的作用通常不是线性的，线性假设很可能不能满足实际需求，甚至直接违背实际情况。1985 年 Stone 提出加性模型 (additive models) ，模型中每一个加性项使用单个光滑函数来估计，在每一加性项中可以解释因变量如何随自变量变化而变化，解决了模型中自变量数目较多时，模型的估计方差会加大的问题。1990 年，Hastie 和 Tibshirani 扩展了加性模型的应用范围，提出了GAM。

传统线性回归方程的表达式可以定义为：y=a+bx1+cx2+…+zxn+C，而GAM可以定义为g(y)=W0+W1(x1)+W2(x2)+…Wn(xn)+C，其中wi是非参数光滑函数，g(y)是线性预测值，正是由于非参数光滑函数可以选择多种形式如核函数、局部回归光滑函数等，使得模型非常灵活，揭示出自变量的非线性效应。

概念太多不好理解，接下来小云带大家用一个拟合泊松回归的广义加性模型案例分析来学习GAM的实现原理。

首先我们介绍一下用于分析的数据，具体如下：

本次分析节选了马里兰州河流生物资源调查的部分数据，生物学目的是探索可能影响鱼类物种丰度的环境因素，并对物种丰度变化的原因作出解释。每行为各种鱼类，各列依次为水域中的个体数量(fish)、水域流域面积(acre)、水域溶解氧含量(do2)、水域最大深度(depth)、水域硝酸盐浓度(no3)、水域硫酸盐浓度(so4)、水域温度(temp)。

介绍完数据，就开始我们的整体了，上代码！

1、导入相关R包及数据

#使用 mgcv 包拟合广义加性模型

library(mgcv)

#读取鱼类物种丰度和水体环境数据

dat<- read.delim('fish_data.txt', sep = '\t', row.names = 1)

2、GAM模型构建

#通过 family='poisson' 指定泊松加性模型

#需要在 gam() 函数中指定自变量的局部平滑器类型，如 s() 的样条平滑、lo() 的 LOESS 平滑等，此处以样条平滑 s() 为例，参数 k 可用于指定平滑程度，值越小约平滑，越高扭动越高

gam_poisson <- gam(fish~s(acre)+s(do2)+s(depth)+s(no3)+s(so4)+s(temp), data = dat, family = 'poisson')

summary(gam_poisson)

3、平滑回归曲线参数选择

#平滑回归曲线图，默认显示 95% 置信区间，通过 select 参数选择第 n 个自变量展示

par(mfrow = c(2, 3))

plot(gam_poisson, select = 1, pch = 20, shade = TRUE, residuals = TRUE)

plot(gam_poisson, select = 2, pch = 20, shade = TRUE, residuals = TRUE)

plot(gam_poisson, select = 3, pch = 20, shade = TRUE, residuals = TRUE)

plot(gam_poisson, select = 4, pch = 20, shade = TRUE, residuals = TRUE)

plot(gam_poisson, select = 5, pch = 20, shade = TRUE, residuals = TRUE)

plot(gam_poisson, select = 6, pch = 20, shade = TRUE, residuals = TRUE)

#更改平滑参数，将导致不同的生物学意义解释

#例如更改氧含量的平滑参数，在样条平滑 s() 中调试参数 k

gam_poisson1<-gam(fish~s(acre)+s(do2,k=3)+s(depth)+s(no3)+s(so4)+ s(temp), data = dat, family = 'poisson')

gam_poisson2<-gam(fish~s(acre)+s(do2,k=5)+s(depth)+s(no3)+s(so4)+ s(temp), data = dat, family = 'poisson')

par(mfrow = c(1, 3))

plot(gam_poisson, select = 2, pch = 20, shade = TRUE, residuals = TRUE)

plot(gam_poisson1, select = 2, pch = 20, shade = TRUE, residuals = TRUE)

plot(gam_poisson2, select = 2, pch = 20, shade = TRUE, residuals = TRUE)

4、确定模型及可视化

#考虑了偏大离差问题的泊松加性模型

gam_quasipoisson<-gam(fish~s(acre)+s(do2)+s(depth)+s(no3)+s(so4)+ s(temp), data = dat, family = 'quasipoisson')

summary(gam_quasipoisson)

plot(gam_quasipoisson)

结果展示

如下图所示，图1中在水域流域面积（acre）增加的梯度上，纵轴描述为fish丰度的对数随水域流域面积的增加而增加；但从曲线特征来看似乎是分段式的，15000英亩时存在一个拐点，可能水域面积此时突破了种群密度的限制。类似地理解，fish丰度的对数随水域溶解氧含量（do2）刚开始时升高，氧含量有助于鱼类生存；但随后在大范围区间内波动，表明后续水中溶氧量增加不再是种群数量的限制因素。其他的指标也是可以照这个思路来解读。当然这只是曲线图的趋势，我们还可以更改曲线的平滑程度，从而产生其他的结果。例如图二中，我们更改了do2的平滑参数k，第一个图自动确认平滑曲线得到的，后两个个图分别为k=3和k=5的回归曲线。但在实际的生物学意义解读中，只允许一种正确的解释，因此解释不同参数的平滑回归输出时必须谨慎。