【R数据处理】GLM（广义线性模型）分析
广义线性模型，是为了克服线性回归模型的缺点出现的，是线性回归模型的推广。首先自变量可以是离散的，也可以是连续的。离散的可以是0-1变量，也可以是多种取值的变量。广义线性模型取消了对残差(因变量)服从正态分布的要求。残差不一定要服从正态分布，可以服从二项、泊松、负二项、正态、伽马、逆高斯等分布，这些分布被统称为指数分布族。（这一段是我在网上找的，想要进一步了解GLM的，请参考R语言实战或者度娘）
在介绍GLM之前，我先说一下为什么我要了解并掌握GLM分析。1）我看到了多篇NC（nature communications）中使用过GLM分析。他们使用GLM要么推断多个自变量对目标变量的解释效应；要么通过算法从很多GLMs中获得最简GLM，然后再根据该GLM预测目标变量的发展趋势；2）看起来这个算法和模型很牛犇。推荐一篇NC供大家在使用该模型时参考“A meta-analysis of global fungal distribution reveals climate-driven patterns.”
模型构建
在前段时间我介绍的随机森林模型的推文中，使用测试数据，我们发现pH是影响物种丰富度（Richness）的主要因素，其它因素对物种的丰富度均没有显著的影响（如下图）。在计算这个随机森林模型的过程中，我们人为的把pH，CN比、P含量、TC（总碳）、Torigin（初始温度）、ECEC（离子交换量）、CP比、NP比和TN（总氮）作为该模型的一个自变量。最终我们发现这些自变量构成的模型对丰富度的解释量为25.45%。现在问题来了，为什么要选择这些自变量而不是那些自变量作为模型中的一个因子？这些自变量的组合是最优组合吗？这个模型是最优最简模型吗？带着这些问题我们来了解广义线性模型。
#load packages
library(tidyverse)
#install.packages("leaps")
library(leaps)
#load data
load("RFdata2.RData")
head(RFdata2)
数据格式如下：
1）计算不同GLMs模型对变量的解释效应
leaps <- regsubsets(Richness~.,data = RFdata2,
nbest=2)
plot(leaps, scale = "adjr2")通过全子集回归分析，我们获得了一批模型及其对应的调整R2（如上图）。这个图的左侧纵坐标为调整R2，横坐标为截距和各个自变量，存在颜色表示包含该自变量，空白表示不包含该自变量。
我们发现当模型仅有一个变量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。该结果表明这两个模型可能都是解释量最高的模型。
为了进一步评估哪个模型是最优模型且同时是最简模型，我们可以看一下每个模型的BIC值，一般来说该值越小则表示模型的拟合度（也就是R2，不是调整R2）越好。
plot(leaps, scale = "bic")我们发现不同GLMs的BIC值排序并不与调整R2一致。结果表明了pH、TN、TC和CN_ratio构成的模型以及pH、P、TC和CP_ratio这两个模型的BIC值最低。查看上一个调整R2的值，它们对应的调整R2分别为0.62和0.62。该结果表明这两个模型都是最简模型。因为它们与最大的拟合度0.66只差0.04，因此，从模型的简单性来说，这两个模型就是最优最简模型。根据自己的科研目的可以选择其中之一。
最终模型如下：
names(RFdata2)
fit <- lm(Richness ~ pH+P+TC+CP_ratio, data = RFdata2)
summary(fit)通过该结果我们发现，该模型显著影响丰富度，且模型中的每个变量都显著影响丰富度，模型的拟合度为0.66，调整拟合度为0.62。
模型的交叉验证
上面我们已经通过算法获得了最优最简模型，那么该模型的稳健性如何呢？下面我们对该模型进行交叉验证。
什么叫交叉验证？
所谓交叉验证指的是将一定比例的样品挑选出来作为训练样本，另一部分样品作为保留样品，先使用训练样品获得回归方程，然后在保留样品上预测。因为保留样品并没有参与模型的构建过程，因此可以用来估测模型的准确性。
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.13053710倍交叉验证的结果表明，我们最终获得的模型对丰富度的实际解释量为0.53；变化性为0.13（这相当于误差）。
然后通过该模型预测因变量的值如下：
#predict valuse
predValue <- predict(fit,RFdata2[,c("pH","P","TC","CP_ratio")],
interval="predict")
predValuefit表示通过该模型预测得到的丰富度值，lwr和upr分别表示下和上边界。
模型中每个变量的重要性
在获得模型后，我们往往还想要知道获得的模型中每一个变量对自变量如何重要，类似于随机森林分析（可以使用随机森林分析预测）也可以通过以下代码预测（参考R语言实战）。代码和结果如下：
#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")
跟我们的随机森林分析的结果对照且相同，GLM模型的结果也表明了pH是影响richness的最主要影响因素。其次是CP比，影响最小的是TC。
希望大家看一下我的群公告，在力所能及的情况下帮一下忙，谢谢。