在R的mgcv::gam中验证模型公式、计算AUC及变量重要性

数据集结构 #

structure(list(pa = structure(c(2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 
1L, 2L, 1L, 1L, 2L, 1L, 1L), .Label = c("0", "1"), class = "factor"), 
    lon = c(26.953632, 26.914444, 26.854655, 26.377477, 26.653273, 
    26.739085, 26.732233, 26.67895, 26.6691, 26.925116, 26.771316, 
    26.952233, 26.934466, 26.9493, 26.948333), lat = c(37.65571, 
    37.658056, 37.548262, 37.714353, 37.670897, 37.652183, 37.664717, 
    37.672083, 37.6934, 37.63755, 37.41155, 37.65095, 37.661533, 
    37.65825, 37.652166), distance = c(2664.205501, 2188.408657, 
    1309.509802, 2931.223857, 443.7116677, 83.4248179, 1162.349952, 
    1025.302461, 1447.284772, 156.3081952, 1718.49796, 2120.230705, 
    2940.015299, 2859.658249, 2179.706853), N = c(2L, 3L, 3L, 
    4L, 1L, 3L, 3L, 4L, 8L, 7L, 2L, 0L, 10L, 0L, 0L), nh4 = c(0.0911071189102672, 
    0.0912837530530634, 0.0887604283967188, 0.0809833919295647, 
    0.0806452852518153, 0.0873989977309376, 0.0854938036251452, 
    0.0837840217003991, 0.113291559368372, 0.139553981108798, 
    0.136305334431029, 0.149872598116116, 0.14975582563108, 0.149872598116116, 
    0.149872598116116), ppn = c(3.13649814951996, 3.38222779366539, 
    2.5790228332411, 1.68392748415672, 2.80087243875361, 3.2346900728285, 
    3.17393288172866, 2.63412894585215, 3.14572940860351, 4.80038520203728, 
    5.83457531216185, 5.10820325640801, 5.14342739916075, 5.10820325640801, 
    5.10820325640801)), row.names = c(1L, 2L, 3L, 5L, 6L, 7L, 
8L, 9L, 10L, 11L, 13L, 16L, 17L, 18L, 19L), class = "data.frame")

拟合的GAM模型 #

mod <- mgcv::gam(data=db, family=binomial(link="logit"), method="REML",
           cbind(pa, N) ~ s(lon) + s(lat) + ti(lon, lat, distance, bs = "re") +
             s(nh4)  + s(ppn, k = 10) )

注：pa是二项式存在/不存在变量，N为pa=1时的个体数量

AUC计算代码及报错 #

library(mgcv)   # library for GAM
library(ggplot2)  # for beautiful plots
library(cdata)    # data wrangling
library(sigr)     # AUC calculation

data <- dplyr::select(db, pa, lon, lat, distance, nh4, ppn, N, season) 
randn=runif(nrow(data))
train_idx=randn<=0.8
train=data[train_idx,]
test=data[!train_idx,]

performance=function(y,pred){
  confmat_test=table(truth=y,predict=pred>0.5)
  acc=sum(diag(confmat_test))/sum(confmat_test)
  precision=confmat_test[2,2]/sum(confmat_test[,2])
  recall=confmat_test[2,2]/sum(confmat_test[2,])
  auc=calcAUC(pred,y)
  c(acc,precision,recall,auc)
}
# Posterior probability
train$pred=predict(gam_model,newdata = train,type = "response")
test$pred=predict(gam_model,newdata=test,type="response")

# model performance evaluated using training data
perf_train=performance(train$pa_dd,train$pred)
perf_test=performance(test$pa_dd,test$pred)
perf_mat=rbind(perf_train,perf_test)
colnames(perf_mat)=c("accuracy","precision","recall","AUC")
round(perf_mat,4)

技术疑问 #

1. 上述模型公式是否正确？
2. 如何正确计算该GAM模型的AUC？
3. 如何计算模型中各变量的重要性？

1. 模型公式的正确性分析 #

模型存在三处关键问题：

• 响应变量格式错误：cbind(pa, N)不符合二项式模型的计数格式（要求是cbind(成功数, 失败数)）。你的数据是存在/缺失+丰度结构，若仅做存在性预测，直接用pa作为响应变量即可；若要同时考虑丰度，需改用零膨胀类GAM（如ziplss族）。
• 张量交互项错误：ti(lon, lat, distance, bs = "re")中，bs="re"是随机效应平滑，仅适用于分组变量，而lon、lat、distance均为连续变量，应改用te(lon, lat, distance)（张量乘积平滑）来构建连续变量的交互效应。
• 变量冗余：同时加入s(lon)、s(lat)和包含这两个变量的交互项，会导致效应重叠，建议移除单独的主效应平滑项，或通过分层张量项（ti）明确区分主效应与交互效应。

修正后的模型示例（仅做存在性预测）：

db$pa_num <- as.numeric(as.character(db$pa))
mod <- mgcv::gam(data=db, family=binomial(link="logit"), method="REML",
           pa_num ~ te(lon, lat, distance) + s(nh4) + s(ppn, k=10))

2. 正确计算AUC的方法 #

首先修正代码中的变量名错误（模型名、响应变量名），再确保响应变量为0/1数值型，最后计算AUC：

修正后的完整代码 #

# 加载依赖包
library(mgcv)
library(sigr)
library(dplyr)

# 预处理数据
db$pa_num <- as.numeric(as.character(db$pa))
data <- dplyr::select(db, pa_num, lon, lat, distance, nh4, ppn) 

# 划分训练/测试集
set.seed(123)
train_idx <- sample(nrow(data), floor(0.8*nrow(data)))
train <- data[train_idx,]
test <- data[-train_idx,]

# 拟合模型（用修正后的公式）
mod <- mgcv::gam(data=train, family=binomial(link="logit"), method="REML",
           pa_num ~ te(lon, lat, distance) + s(nh4) + s(ppn, k=10))

# 预测概率
train$pred <- predict(mod, newdata=train, type="response")
test$pred <- predict(mod, newdata=test, type="response")

# 性能计算函数（处理边界情况）
performance <- function(y, pred){
  y <- as.numeric(y)
  # 固定类别顺序避免混淆矩阵错误
  confmat <- table(truth=factor(y, levels=c(0,1)),
                   predict=factor(as.numeric(pred>0.5), levels=c(0,1)))
  acc <- sum(diag(confmat))/sum(confmat)
  precision <- ifelse(sum(confmat[,2])==0, NA, confmat[2,2]/sum(confmat[,2]))
  recall <- ifelse(sum(confmat[2,])==0, NA, confmat[2,2]/sum(confmat[2,]))
  auc <- sigr::calcAUC(pred, y)
  c(accuracy=acc, precision=precision, recall=recall, AUC=auc)
}

# 计算并输出结果
perf_train <- performance(train$pa_num, train$pred)
perf_test <- performance(test$pa_num, test$pred)
round(rbind(Train=perf_train, Test=perf_test), 4)

备选方案（用pROC包） #

library(pROC)
# 训练集AUC
auc(roc(train$pa_num, train$pred))
# 测试集AUC
auc(roc(test$pa_num, test$pred))

3. 计算模型变量重要性 #

GAM中常用三种变量重要性计算方式：

方法1：基于平滑项的偏差解释占比 #

通过summary.gam输出的%dev explained，判断每个平滑项解释的偏差占总解释偏差的比例，占比越高变量越重要：

sum_mod <- summary(mod)
imp_df <- data.frame(
  Variable = rownames(sum_mod$s.table),
  Deviance_Explained = sum_mod$s.table[,"%dev"],
  P_value = sum_mod$s.table[,"p-value"]
)
# 按重要性排序
imp_df[order(-imp_df$Deviance_Explained),]

方法2：置换重要性 #

通过置换某一变量的取值，观察模型性能的下降幅度，下降越多变量越重要：

library(permimp)
# 计算置换重要性
perm_imp <- permimp(mod, method="permute", nperm=100, progress="text")
# 可视化结果
plot(perm_imp)

方法3：基于有效自由度（edf） #

平滑项的有效自由度越高，说明变量的非线性效应越强，可作为重要性的辅助参考（需结合偏差解释占比）：

data.frame(
  Variable = rownames(sum_mod$s.table),
  EDF = sum_mod$s.table[,"edf"]
)

内容的提问来源于stack exchange，提问作者Mauri21