这是我的第一篇文章，我对 GAM 还很陌生；道歉。

我运行了一系列 16 个广义加性负二项式模型（gam、family=nb、mgcv 包），其复杂性不断增加，将 12 种鱼类的丰度（独立）建模为几个环境指标（及其组合）的函数：

这是一个大型数据集；几个模型示例：

m4(s): ccount ~ s(sal, k=6, bs="tp")

m12(ost): ccount ~ s(do, k=6, bs="tp") + s(sal, k=6, bs="tp") + s(temp, k=6, bs="tp")

我进行了一系列模型检查，然后绘制了每个模型的估计自由度 (edf)、偏差解释和 AIC（按照从最小到最复杂的近似顺序），以帮助每个物种的模型选择。

法国电力公司

偏差解释

AIC

似乎解释的偏差只是 AIC 的倒数，并且 AIC 似乎不会因模型 EDF 的变化而受到惩罚。我绘制了偏差解释与 AIC，并惊讶于它们几乎完全相关。这似乎违反直觉。

AIC 与 Dev.Expl

问题/评论 1 - 我发现了三种不同的方式来调用 AIC：

1) model$aic
2) AIC(model)
3) extractAIC(model)

(3) 仅针对参数模型（即非 GAM）进行了记录。目前尚不清楚（1）或（2）是否更适合基于 mgcv pkg 的 gam 对象？两者都给出了几乎相同的结果（见下文），但似乎使用了不同的 df（edf 与参考 df）。我假设（1）是最好的；但是，尚不清楚使用了什么 df，并且有 logLik 与 (2) 一起用于 gam 对象的文档。目前还不清楚这是否已经内置在（1）中？https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/logLik.gam.html

问题 2 - 无论模型复杂性如何，期望解释的偏差和 AIC 具有如此强的相关性是否合理？我不这么认为。dev.expl 的改进是否足够大，以至于 edf 的变化无关紧要？模型中的 edf 是否已经对偏差解释进行了惩罚？

谢谢你的帮助。

这是一个来自 mtcars 的简单可重现示例，显示了相同的模式：

library(mgcv)

# 3 GAM models
b1 <-gam(mpg~s(hp), data=mtcars)
b2 <-gam(mpg~s(wt), data=mtcars)
b3 <-gam(mpg~s(hp)+s(wt), data=mtcars)

# AIC1
a1.1 <- b1$aic
a1.2 <- b2$aic
a1.3 <- b3$aic
a1 <- c(a1.1,a1.2,a1.3)

# AIC2
a2.1 <- AIC(b1)
a2.2 <- AIC(b2)
a2.3 <- AIC(b3)
a2 <- c(a2.1, a2.2, a2.3)

# dev.explained
d1 <- summary(b1)$dev.expl
d2 <- summary(b2)$dev.expl
d3 <- summary(b3)$dev.expl
d <- c(d1,d2,d3)

par(mfrow=c(2,2))
plot(a1~d, type=c("b"), xlab="deviance explained", ylab="model$aic")
plot(a1~a2, col="black", type="b", lwd=1, cex=1, xlab = "AIC(model)", ylab="model$aic")

MTCARS GAM 示例：AIC vs Dev.Expl & model$aic vs AIC(model)

根据 Jon（下）的有用评论，我进一步探索了如何使用 mtcars 的 2-factor gam（b3，上）在 mgcv 中计算 gam 项。我能够计算出与模型输出相匹配的偏差解释和 AIC（尽管 df 的计算似乎双重计算了截距以获得与 gam 输出相同的 df；奇怪）。

MTCARS GAM 平滑

我无法进行的唯一计算是 Deviance。提供的模型表明偏差=2logL（模型|数据）-常数）。I，假设与 AIC 相同的常数 2*p，不能产生与 gam 模型输出相同的结果。这只能证明，虽然掌握了AIC和dev.expl，但我仍然不明白Deviance是如何计算的，以及它与AIC的关系。

> # GAM based on mtcars

> library(mgcv)
> b <-gam(mpg~s(hp)+s(wt), data=mtcars)
> summary(b)

Family: gaussian 
Link function: identity 

Formula:
mpg ~ s(hp) + s(wt)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  20.0906     0.3723   53.97   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df      F  p-value    
s(hp) 2.409  3.016  6.305  0.00218 ** 
s(wt) 2.085  2.523 15.277 1.19e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.878   Deviance explained = 89.6%
GCV = 5.3542  Scale est. = 4.4349    n = 32
> plot(b)



> # LOG LIKELIHOOD OF GAM MODEL (AND DF)

> LL <- logLik(b)
> LL
'log Lik.' -66.22412 (df=6.494385)



> # MODEL DF = no. of parameters (p) = sum of smooth (edf) + parametric (nsdf) terms

>     # note: intercept gets counted twice (in b$edf and b$nsdf) to attain same df as logLik(b) and aic = AIC(b)
>     # note: summary(b)$pTerms.pv(does not include intercept as a parameter term; gives zero result for unspecified param param)
> p <-  sum(b$edf)+sum(b$nsdf) 
> p
[1] 6.494385



> # AIC = penalty - twice the log likelihood (2p-LL)

> 2*p - 2*LL
'log Lik.' 145.437 (df=6.494385)
> AIC(b)
[1] 145.437
> b$aic
[1] 145.437



> #DEVIANCE EXPLAINED = 1- (residual dev / total (null) deviance)

> 1-(b$deviance/b$null.deviance)  
[1] 0.895608
> summary(b)$dev.expl
[1] 0.895608



> #DEVIANCE = the unexplained error in the model (residual deviance)

> 2*LL - 2*p    # DOES NOT MATCH b$deviance!
'log Lik.' -145.437 (df=6.494385)
> b$deviance
[1] 117.5503
>

任何进一步的想法将不胜感激。

GAM (mgcv)：AIC 与偏差解释

AIC 与 Dev.Expl