3

我正在使用以下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)) 的指数协方差函数。如何从结果中检索变异函数参数（块金、窗台）？我在模型输出中找不到任何东西。

4

1 回答 1

1

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

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

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

这是正确的吗？不，这里有一个陷阱：返回的平滑参数$sp不是真实的，我们需要以下内容：

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

我们复制您指定的范围参数：

gm2_phi <- 10

块金必须为零，因为平滑函数是连续的。使用包中的lines.variomodel函数geoR，您可以可视化由 建模的潜在高斯空间随机场的半变异函数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 个唯一空间位置，因此相关矩阵应为 630 x 630，全随机效应应为长度为 630 的向量。但是通过设置截断k = 100的s(, bs = 'gp')特征分解和随后的低秩近似，将随机效应减少到长度为 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 的基数要高得多。

于 2018-08-09T14:26:59.467 回答