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

科白君 2021-10-26

展開全文

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

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

編輯 | free傻孩子

4306字 |10分鐘閱讀

本期推送內容

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

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

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

模型構建

在前段時間我介紹的隨機森林模型的推文中，使用測試數(shù)據(jù)，我們發(fā)現(xiàn)pH是影響物種豐富度（Richness）的主要因素，其它因素對物種的豐富度均沒有顯著的影響（如下圖）。在計算這個隨機森林模型的過程中，我們人為的把pH，CN比、P含量、TC（總碳）、Torigin（初始溫度）、ECEC（離子交換量）、CP比、NP比和TN（總氮）作為該模型的一個自變量。最終我們發(fā)現(xiàn)這些自變量構成的模型對豐富度的解釋量為25.45%?，F(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模型對變量的解釋效應

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

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

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

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

plot(leaps, scale = "bic")

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

最終模型如下：

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

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

模型的交叉驗證

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

什么叫交叉驗證？

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

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

代碼如下：

#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倍交叉驗證的結果表明，我們最終獲得的模型對豐富度的實際解釋量為0.53；變化性為0.13（這相當于誤差）。

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

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

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

模型中每個變量的重要性

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

#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")