GAM与"gp“平滑器:如何检索变差函数参数？

Question

我正在使用下面的geoadditive模型

library(gamair)
library(mgcv)

data(mack)    
mack$log.net.area <- log(mack$net.area)

gm2 <- gam(egg.count ~ s(lon,lat,bs="gp",k=100,m=c(2,10,1)) +
                       s(I(b.depth^.5)) +
                       s(c.dist) +
                       s(temp.20m) +
                       offset(log.net.area),
                       data = mack, family = tw, method = "REML")

这里我使用了一个范围= 10，幂=1 (m=c(2,10,1))的指数协方差函数。如何从结果中检索变差函数参数(块金、底座)？我在模型输出中找不到任何东西。

Answer 1

在平滑法中，指定了相关矩阵，因此您只需估计方差参数，即窗台。例如，您将m = c(2, 10, 1)设置为s(, bs = 'gp')，从而给出一个带有范围参数phi = 10的指数相关矩阵。请注意，除球面相关外，phi与range不同。对于许多相关模型，实际范围是phi的函数。

方差/窗台参数与惩罚回归中的平滑参数密切相关，可以通过将scale参数除以平滑参数得到：

with(gm2, scale / sp["s(lon,lat)"])
#s(lon,lat) 
#  26.20877 

是这样的吗？不是的。这里有一个陷阱：smoothing parameters returned in $sp are not real ones，我们需要以下代码：

gm2_sill <- with(gm2, scale / sp["s(lon,lat)"] * smooth[[1]]$S.scale) 
#s(lon,lat) 
#    7.7772 

然后我们复制您指定的range参数：

gm2_phi <- 10

块必须为零，因为光滑函数是连续的。使用geoR软件包中的lines.variomodel函数，可以可视化s(lon,lat)模拟的潜在高斯空间随机场的半方差函数。

library(geoR)
lines.variomodel(cov.model = "exponential", cov.pars = c(gm2_sill, gm2_phi),
                 nugget = 0, max.dist = 60)
abline(h = gm2_sill, lty = 2)

​

​

然而，对这个变异函数要持怀疑态度。mgcv不是一个解释地统计学的简单环境。低阶平滑器的使用表明，上面的方差参数用于新参数空间中的参数，而不是原始参数空间中的参数。例如，mack数据集的空间场中有630个唯一的空间位置，因此相关矩阵应该是630x630，并且完整的随机效果应该是长度为- 630的向量。但是，通过在s(, bs = 'gp')中设置k = 100，截断的特征分解和随后的低秩近似将随机效应减少到长度为-100。方差参数实际上是针对此向量的，而不是原始向量。这可能解释了为什么基准面和实际范围与数据和预测的s(lon,lat)不一致。

## unique locations
loc <- unique(mack[, c("lon", "lat")])

max(dist(loc))
#[1] 15.98

数据集中两个空间位置之间的最大距离为15.98，但离变差函数的实际范围似乎在40到60之间，这太大了。

## predict `s(lon, lat)`, using the method I told you in your last question
## https://stackoverflow.com/q/51634953/4891738
sp <- predict(gm2,
              data.frame(loc, b.depth = 0, c.dist = 0, temp.20m = 0,
                         log.net.area = 0),
              type = "terms", terms = "s(lon,lat)")

c(var(sp))
#[1] 1.587126

预测的s(lon,lat)只有1.587的方差，但7.77时的基准值要高得多。