GAM с gp smoother: как получить параметры вариограммы?

Я использую следующую геоаддитивную модель

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)). Как получить из результатов параметры вариограммы (самородок, порог)? В выводе модели ничего не нашел.


r regression gam mgcv geostatistics
person user3036416    schedule 02.08.2018    source источник
comment
@李哲源 да, ты прав. Некоторые детали и ссылки приведены в Kammann and Wand (2003).   -  person user3036416    schedule 03.08.2018

Ответы (1)


arrow_upward
1
arrow_downward

В подходе сглаживания корреляционная матрица указывается, поэтому вы оцениваете только параметр дисперсии, то есть порог. Например, вы установили 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 намного выше.

person Zheyuan Li    schedule 09.08.2018