AIC() 和 model$aic 在 mgcv::gam() 中给出不同的结果
AIC() and model$aic give different result in mgcv::gam()
我正在使用具有空间平滑的 mgcv::gam() 拟合 GAM,类似于 this example (data here 中的模型)。我发现使用 AIC(model)
会得到与 model$aic
不同的结果。为什么是这样?哪个是正确的?
library('mgcv')
galveston <- read.csv('gbtemp.csv')
galveston <- transform(galveston,
datetime = as.POSIXct(paste(DATE, TIME),
format = '%m/%d/%y %H:%M', tz = "CDT"))
galveston <- transform(galveston,
STATION_ID = factor(STATION_ID),
DoY = as.numeric(format(datetime, format = '%j')),
ToD = as.numeric(format(datetime, format = '%H')) +
(as.numeric(format(datetime, format = '%M')) / 60))
knots <- list(DoY = c(0.5, 366.5))
M <- list(c(1, 0.5), NA)
m <- bam(MEASUREMENT ~
s(ToD, k = 10) +
s(DoY, k = 30, bs = 'cc') +
s(YEAR, k = 30) +
s(LONGITUDE, LATITUDE, k = 100, bs = 'ds', m = c(1, 0.5)) +
ti(DoY, YEAR, bs = c('cc', 'tp'), k = c(15, 15)) +
ti(LONGITUDE, LATITUDE, ToD, d = c(2,1), bs = c('ds','tp'),
m = M, k = c(20, 10)) +
ti(LONGITUDE, LATITUDE, DoY, d = c(2,1), bs = c('ds','cc'),
m = M, k = c(25, 15)) +
ti(LONGITUDE, LATITUDE, YEAR, d = c(2,1), bs = c('ds','tp'),
m = M, k = c(25, 15)),
data = galveston, method = 'fREML', knots = knots,
nthreads = 4, discrete = TRUE)
AIC(m)
[1] 57073.08
m$aic
[1] 57053.21
请注意,我给出的示例使用 bam()
而不是 gam()
,但结果是相同的。
我无法用更简单的模型复制它(来自 here 的示例):
set.seed(2)
dat <- gamSim(1,n=400,dist="normal",scale=2)
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
AIC(b)
[1] 1696.143
b$aic
[1] 1696.143
差异是因为$aic
中存储的作为AIC中复杂度校正的自由度是模型的有效自由度(EDF)。这已被证明过于自由或保守,并且可能导致 AIC 始终根据使用的是边际 AIC 还是条件 AIC 选择更复杂的模型或更简单的模型。
有一些方法可以纠正这种行为,mgcv 实现了 Wood 等人 (2016) 中的一个方法,该方法对自由度进行了校正。这是通过 logLik.gam()
函数完成的,该函数由 AIC.gam()
调用。这也解释了差异,因为 $aic
是没有应用校正的标准 AIC,而 IIRC 是 GAM 对象的一个组成部分,它明显早于 Wood 等人 (2016) 的工作。
至于为什么你不能用简单的例子复制这个,那是因为校正需要使用拟合的组件,只有当用于拟合的 method
是 "REML"
或 "ML"
(包括 "fREML"
用于 bam()
,并且在使用 Extended Fellner Schall 或 BFGS 优化器时也不包括:
> library('mgcv')
Loading required package: nlme
This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.
> set.seed(2)
> dat <- gamSim(1,n=400,dist="normal",scale=2)
Gu & Wahba 4 term additive model
> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat, method = 'REML')
> AIC(b)
[1] 1698.504
> b$aic
[1] 1696.894
method
默认使用 GCV。
参考资料
Wood, S.N., Pya, N., Säfken, B., 2016。一般平滑模型的平滑参数和模型选择。 J. Am。状态。副教授。 111, 1548–1563。 https://doi.org/10.1080/01621459.2016.1180986
我正在使用具有空间平滑的 mgcv::gam() 拟合 GAM,类似于 this example (data here 中的模型)。我发现使用 AIC(model)
会得到与 model$aic
不同的结果。为什么是这样?哪个是正确的?
library('mgcv')
galveston <- read.csv('gbtemp.csv')
galveston <- transform(galveston,
datetime = as.POSIXct(paste(DATE, TIME),
format = '%m/%d/%y %H:%M', tz = "CDT"))
galveston <- transform(galveston,
STATION_ID = factor(STATION_ID),
DoY = as.numeric(format(datetime, format = '%j')),
ToD = as.numeric(format(datetime, format = '%H')) +
(as.numeric(format(datetime, format = '%M')) / 60))
knots <- list(DoY = c(0.5, 366.5))
M <- list(c(1, 0.5), NA)
m <- bam(MEASUREMENT ~
s(ToD, k = 10) +
s(DoY, k = 30, bs = 'cc') +
s(YEAR, k = 30) +
s(LONGITUDE, LATITUDE, k = 100, bs = 'ds', m = c(1, 0.5)) +
ti(DoY, YEAR, bs = c('cc', 'tp'), k = c(15, 15)) +
ti(LONGITUDE, LATITUDE, ToD, d = c(2,1), bs = c('ds','tp'),
m = M, k = c(20, 10)) +
ti(LONGITUDE, LATITUDE, DoY, d = c(2,1), bs = c('ds','cc'),
m = M, k = c(25, 15)) +
ti(LONGITUDE, LATITUDE, YEAR, d = c(2,1), bs = c('ds','tp'),
m = M, k = c(25, 15)),
data = galveston, method = 'fREML', knots = knots,
nthreads = 4, discrete = TRUE)
AIC(m)
[1] 57073.08
m$aic
[1] 57053.21
请注意,我给出的示例使用 bam()
而不是 gam()
,但结果是相同的。
我无法用更简单的模型复制它(来自 here 的示例):
set.seed(2)
dat <- gamSim(1,n=400,dist="normal",scale=2)
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
AIC(b)
[1] 1696.143
b$aic
[1] 1696.143
差异是因为$aic
中存储的作为AIC中复杂度校正的自由度是模型的有效自由度(EDF)。这已被证明过于自由或保守,并且可能导致 AIC 始终根据使用的是边际 AIC 还是条件 AIC 选择更复杂的模型或更简单的模型。
有一些方法可以纠正这种行为,mgcv 实现了 Wood 等人 (2016) 中的一个方法,该方法对自由度进行了校正。这是通过 logLik.gam()
函数完成的,该函数由 AIC.gam()
调用。这也解释了差异,因为 $aic
是没有应用校正的标准 AIC,而 IIRC 是 GAM 对象的一个组成部分,它明显早于 Wood 等人 (2016) 的工作。
至于为什么你不能用简单的例子复制这个,那是因为校正需要使用拟合的组件,只有当用于拟合的 method
是 "REML"
或 "ML"
(包括 "fREML"
用于 bam()
,并且在使用 Extended Fellner Schall 或 BFGS 优化器时也不包括:
> library('mgcv')
Loading required package: nlme
This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.
> set.seed(2)
> dat <- gamSim(1,n=400,dist="normal",scale=2)
Gu & Wahba 4 term additive model
> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat, method = 'REML')
> AIC(b)
[1] 1698.504
> b$aic
[1] 1696.894
method
默认使用 GCV。
参考资料
Wood, S.N., Pya, N., Säfken, B., 2016。一般平滑模型的平滑参数和模型选择。 J. Am。状态。副教授。 111, 1548–1563。 https://doi.org/10.1080/01621459.2016.1180986