在前段時間我介紹的隨機森林模型的推文中，使用測試數(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%?，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模型對變量的解釋效應(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)健性如何呢,？下面我們對該模型進行交叉驗證。

什么叫交叉驗證,？

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

k重交叉驗證,，指的是將樣品分為k個子集,，輪流將k-1個子樣品作為訓(xùn)練集，另外一個子集作為保留集,，最終獲得平均預(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")

跟我們的隨機森林分析的結(jié)果對照且相同，GLM模型的結(jié)果也表明了pH是影響richness的最主要影響因素,。其次是CP比,，影響最小的是TC。

希望大家看一下我的群公告,，在力所能及的情況下幫一下忙,，謝謝。