【原】【R數(shù)據(jù)處理】GLM（廣義線性模型）分析

科白君 2021-10-26

展開全文

“ 知其然也要知其所以然。” --free傻孩子

"R實戰(zhàn)"專題·第15篇

編輯 | free傻孩子

4306字 |10分鐘閱讀

本期推送內(nèi)容

在數(shù)據(jù)分析的過程中，很多分析方法和模型往往要求目標變量（數(shù)據(jù)）服從某些假設(shè)如正態(tài)分布、方差齊次等。一般來說，如果數(shù)據(jù)不能服從這些假設(shè)，那么采用對應(yīng)的方法或模型獲得的結(jié)果往往不可信。例如，我們經(jīng)常使用的經(jīng)典模型，即形如y = kx +b（在R中形如 lm(y ~ x, data)）的一般線性模型就要求數(shù)據(jù)（目標變量）必須滿足正態(tài)分布和殘差的方差齊次。然而，在實際科研工作中，很多數(shù)據(jù)往往不能滿足以上條件。這種情況就要求我們尋找一種沒有以上假設(shè)的方法來替代存在假設(shè)的模型如：一般線性模型。這種方法之一就是本節(jié)我想給大家推薦的廣義線性模型（GLM）。

廣義線性模型，是為了克服線性回歸模型的缺點出現(xiàn)的，是線性回歸模型的推廣。首先自變量可以是離散的，也可以是連續(xù)的。離散的可以是0-1變量，也可以是多種取值的變量。廣義線性模型取消了對殘差(因變量)服從正態(tài)分布的要求。殘差不一定要服從正態(tài)分布，可以服從二項、泊松、負二項、正態(tài)、伽馬、逆高斯等分布，這些分布被統(tǒng)稱為指數(shù)分布族。（這一段是我在網(wǎng)上找的，想要進一步了解GLM的，請參考R語言實戰(zhàn)或者度娘）

在介紹GLM之前，我先說一下為什么我要了解并掌握GLM分析。1）我看到了多篇NC（nature communications）中使用過GLM分析。他們使用GLM要么推斷多個自變量對目標變量的解釋效應(yīng)；要么通過算法從很多GLMs中獲得最簡GLM，然后再根據(jù)該GLM預(yù)測目標變量的發(fā)展趨勢；2）看起來這個算法和模型很牛犇。推薦一篇NC供大家在使用該模型時參考“A meta-analysis of global fungal distribution reveals climate-driven patterns.”

模型構(gòu)建

在前段時間我介紹的隨機森林模型的推文中，使用測試數(shù)據(jù)，我們發(fā)現(xiàn)pH是影響物種豐富度（Richness）的主要因素，其它因素對物種的豐富度均沒有顯著的影響（如下圖）。在計算這個隨機森林模型的過程中，我們?nèi)藶榈陌裵H，CN比、P含量、TC（總碳）、Torigin（初始溫度）、ECEC（離子交換量）、CP比、NP比和TN（總氮）作為該模型的一個自變量。最終我們發(fā)現(xiàn)這些自變量構(gòu)成的模型對豐富度的解釋量為25.45%。現(xiàn)在問題來了，為什么要選擇這些自變量而不是那些自變量作為模型中的一個因子？這些自變量的組合是最優(yōu)組合嗎？這個模型是最優(yōu)最簡模型嗎？帶著這些問題我們來了解廣義線性模型。

#load packages
library(tidyverse)
#install.packages("leaps")
library(leaps)

#load data
load("RFdata2.RData")
head(RFdata2)

數(shù)據(jù)格式如下：

1）計算不同GLMs模型對變量的解釋效應(yīng)

leaps <- regsubsets(Richness~.,data = RFdata2,
nbest=2)
plot(leaps, scale = "adjr2")

通過全子集回歸分析，我們獲得了一批模型及其對應(yīng)的調(diào)整R2（如上圖）。這個圖的左側(cè)縱坐標為調(diào)整R2，橫坐標為截距和各個自變量，存在顏色表示包含該自變量，空白表示不包含該自變量。

我們發(fā)現(xiàn)當模型僅有一個變量Torigin時（最下方），GLM模型的調(diào)整R2為0.26，而當模型包含Torigin、pH、P、TC、CN_ratio、CP_ratio和NP_ratio時模型的調(diào)整R2最大為0.66；相同的當模型包含Torigin、pH、P、TN、CN_ratio、CP_ratio和NP_ratio時模型的調(diào)整R2也是最大值0.66。該結(jié)果表明這兩個模型可能都是解釋量最高的模型。

為了進一步評估哪個模型是最優(yōu)模型且同時是最簡模型，我們可以看一下每個模型的BIC值，一般來說該值越小則表示模型的擬合度（也就是R2，不是調(diào)整R2）越好。

plot(leaps, scale = "bic")

我們發(fā)現(xiàn)不同GLMs的BIC值排序并不與調(diào)整R2一致。結(jié)果表明了pH、TN、TC和CN_ratio構(gòu)成的模型以及pH、P、TC和CP_ratio這兩個模型的BIC值最低。查看上一個調(diào)整R2的值，它們對應(yīng)的調(diào)整R2分別為0.62和0.62。該結(jié)果表明這兩個模型都是最簡模型。因為它們與最大的擬合度0.66只差0.04，因此，從模型的簡單性來說，這兩個模型就是最優(yōu)最簡模型。根據(jù)自己的科研目的可以選擇其中之一。

最終模型如下：

names(RFdata2)
fit <- lm(Richness ~ pH+P+TC+CP_ratio, data = RFdata2)
summary(fit)

通過該結(jié)果我們發(fā)現(xiàn)，該模型顯著影響豐富度，且模型中的每個變量都顯著影響豐富度，模型的擬合度為0.66，調(diào)整擬合度為0.62。

模型的交叉驗證

上面我們已經(jīng)通過算法獲得了最優(yōu)最簡模型，那么該模型的穩(wěn)健性如何呢？下面我們對該模型進行交叉驗證。

什么叫交叉驗證？

所謂交叉驗證指的是將一定比例的樣品挑選出來作為訓練樣本，另一部分樣品作為保留樣品，先使用訓練樣品獲得回歸方程，然后在保留樣品上預(yù)測。因為保留樣品并沒有參與模型的構(gòu)建過程，因此可以用來估測模型的準確性。

k重交叉驗證，指的是將樣品分為k個子集，輪流將k-1個子樣品作為訓練集，另外一個子集作為保留集，最終獲得平均預(yù)測值。

代碼如下：

#install.packages("bootstrap")
library(bootstrap)


shrinkage <- function(fit, k = 10){
require(bootstrap)
set.seed(123)


theta.fit <- function(x,y){lsfit(x,y)}
theta.predict <- function(fit,x){cbind(1,x) %*% fit$coef}


x <- fit$model[,2:ncol(fit$model)]
y <- fit$model[,1]


results <- crossval(x, y, theta.fit,theta.predict, ngroup = k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
cat("Original R-square =", r2, "\n")
cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
cat("Change =", r2-r2cv, "\n")
}

shrinkage(fit, k =10)
#Original R-square = 0.6993476
#10 Fold Cross-Validated R-square = 0.527686

#Change = 0.130537

10倍交叉驗證的結(jié)果表明，我們最終獲得的模型對豐富度的實際解釋量為0.53；變化性為0.13（這相當于誤差）。

然后通過該模型預(yù)測因變量的值如下：

#predict valuse
predValue <- predict(fit,RFdata2[,c("pH","P","TC","CP_ratio")],
interval="predict")
predValue

fit表示通過該模型預(yù)測得到的豐富度值，lwr和upr分別表示下和上邊界。

模型中每個變量的重要性

在獲得模型后，我們往往還想要知道獲得的模型中每一個變量對自變量如何重要，類似于隨機森林分析（可以使用隨機森林分析預(yù)測）也可以通過以下代碼預(yù)測（參考R語言實戰(zhàn)）。代碼和結(jié)果如下：

#importance of each variables
relweights <- function(fit,...){
set.seed(123)
options(digits = 3)


R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar,1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda^2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta^2)
rawwgt <- lambdasq %*% beta^2
import <- (rawwgt/rsquare) *100
import <- as.data.frame(import)
rownames(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
dotchart(import$Weights, labels = rownames(import),
xlab = "% of R-Square", pch = 19,
main = "Relative importance of predictor variables",
sub = paste("Total R-Square =",round(rsquare,digits = 2)),
...)
return(import)
}
relweights(fit,col = "blue")