【原】【R数据处理】GLM（广义线性模型）分析

科白君 2021-10-26

展开全文

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

"R实战"专题·第15篇

编辑 | free傻孩子

4306字 |10分钟阅读

本期推送内容

在数据分析的过程中，很多分析方法和模型往往要求目标变量（数据）服从某些假设如正态分布、方差齐次等。一般来说，如果数据不能服从这些假设，那么采用对应的方法或模型获得的结果往往不可信。例如，我们经常使用的经典模型，即形如y = kx +b（在R中形如 lm(y ~ x, data)）的一般线性模型就要求数据（目标变量）必须满足正态分布和残差的方差齐次。然而，在实际科研工作中，很多数据往往不能满足以上条件。这种情况就要求我们寻找一种没有以上假设的方法来替代存在假设的模型如：一般线性模型。这种方法之一就是本节我想给大家推荐的广义线性模型（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.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语言实战）。代码和结果如下：

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