相信很多小伙伴們都有過類似的苦惱，總感覺自己了解的可以結(jié)合生信分析的機器學習算法就那幾種，用的人太多了！沒有什么吸引力。那么有沒有什么冷門點的而且可以結(jié)合生信分析的算法呢？

今天，小云就給大家推薦一個，廣義加性模型（generalized additive models， GAM）。

GBM是回歸家族的一個擴展，且是最強大的統(tǒng)計模型之一，非常自由靈活，可以為任何回歸問題建模！傳統(tǒng)線性模型簡單、直觀、便于理解，但是，在現(xiàn)實生活中，變量的作用通常不是線性的，線性假設很可能不能滿足實際需求，甚至直接違背實際情況。1985 年 Stone 提出加性模型 (additive models) ，模型中每一個加性項使用單個光滑函數(shù)來估計，在每一加性項中可以解釋因變量如何隨自變量變化而變化，解決了模型中自變量數(shù)目較多時 ，模型的估計方差會加大的問題。1990 年，Hastie 和 Tibshirani 擴展了加性模型的應用范圍 ，提出了GAM。

傳統(tǒng)線性回歸方程的表達式可以定義為：y=a+bx1+cx2+…+zxn+C，而GAM可以定義為g(y)=W0+W1(x1)+W2(x2)+…Wn(xn)+C，其中wi是非參數(shù)光滑函數(shù)，g(y)是線性預測值，正是由于非參數(shù)光滑函數(shù)可以選擇多種形式如核函數(shù)、局部回歸光滑函數(shù)等，使得模型非常靈活，揭示出自變量的非線性效應。

概念太多不好理解，接下來小云帶大家用一個擬合泊松回歸的廣義加性模型案例分析來學習GAM的實現(xiàn)原理。

首先我們介紹一下用于分析的數(shù)據(jù)，具體如下：

本次分析節(jié)選了馬里蘭州河流生物資源調(diào)查的部分數(shù)據(jù)，生物學目的是探索可能影響魚類物種豐度的環(huán)境因素，并對物種豐度變化的原因作出解釋。每行為各種魚類，各列依次為水域中的個體數(shù)量(fish)、水域流域面積(acre)、水域溶解氧含量(do2)、水域最大深度(depth)、水域硝酸鹽濃度(no3)、水域硫酸鹽濃度(so4)、水域溫度(temp)。

介紹完數(shù)據(jù)，就開始我們的整體了，上代碼！

1、導入相關R包及數(shù)據(jù)

#使用 mgcv 包擬合廣義加性模型

library(mgcv)

#讀取魚類物種豐度和水體環(huán)境數(shù)據(jù)

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

2、GAM模型構(gòu)建

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

#需要在 gam() 函數(shù)中指定自變量的局部平滑器類型，如 s() 的樣條平滑、lo() 的 LOESS 平滑等，此處以樣條平滑 s() 為例，參數(shù) 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、平滑回歸曲線參數(shù)選擇

#平滑回歸曲線圖，默認顯示 95% 置信區(qū)間，通過 select 參數(shù)選擇第 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)

?

#更改平滑參數(shù)，將導致不同的生物學意義解釋

#例如更改氧含量的平滑參數(shù)，在樣條平滑 s() 中調(diào)試參數(shù) 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)

結(jié)果展示

如下圖所示，圖1中在水域流域面積（acre）增加的梯度上，縱軸描述為fish豐度的對數(shù)隨水域流域面積的增加而增加；但從曲線特征來看似乎是分段式的，15000英畝時存在一個拐點，可能水域面積此時突破了種群密度的限制。類似地理解，fish豐度的對數(shù)隨水域溶解氧含量（do2）剛開始時升高，氧含量有助于魚類生存；但隨后在大范圍區(qū)間內(nèi)波動，表明后續(xù)水中溶氧量增加不再是種群數(shù)量的限制因素。其他的指標也是可以照這個思路來解讀。當然這只是曲線圖的趨勢，我們還可以更改曲線的平滑程度，從而產(chǎn)生其他的結(jié)果。例如圖二中，我們更改了do2的平滑參數(shù)k，第一個圖自動確認平滑曲線得到的，后兩個個圖分別為k=3和k=5的回歸曲線。但在實際的生物學意義解讀中，只允許一種正確的解釋，因此解釋不同參數(shù)的平滑回歸輸出時必須謹慎。

整個分析的流程到這里就結(jié)束啦！現(xiàn)在是不是感覺茅塞頓開！原來高分文章離大家是如此之近，心動的小伙伴一定要自己動手試試哦！

（文章結(jié)尾推薦一下小云新開發(fā)的零代碼云生信分析工具平臺包含超多零代碼小工具，上傳數(shù)據(jù)一鍵出圖，感興趣的小伙伴歡迎來參觀喲，網(wǎng)址：http://www.biocloudservice.com/home.html）