GAM с gp более плавным: предсказывайте в новых местах

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

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")

Как я могу использовать его для прогнозирования значения egg.count в новых местоположениях (lon/lat), где у меня нет ковариатных данных, как в kriging?

Например, скажем, я хочу предсказать egg.count в этих новых местах.

    lon lat
1  -3.00  44
4  -2.75  44
7  -2.50  44
10 -2.25  44
13 -2.00  44
16 -1.75  44

но здесь я не знаю значений ковариат (b.depth, c.dist, temp.20m, log.net.area).


r regression gam mgcv kriging
person user3036416    schedule 01.08.2018    source источник

Ответы (1)


arrow_upward
4
arrow_downward

predict по-прежнему требует, чтобы все переменные, используемые в вашей модели, были представлены в newdata, но вы можете передать некоторые произвольные значения, например 0s, в те ковариаты, которых у вас нет, а затем использовать type = "terms" и terms = name_of_the_wanted_smooth_term для продолжения. Использовать

sapply(gm2$smooth, "[[", "label")
#[1] "s(lon,lat)"        "s(I(b.depth^0.5))" "s(c.dist)"        
#[4] "s(temp.20m)"

чтобы проверить, какие гладкие члены есть в вашей модели.

## new spatial locations to predict
newdat <- read.table(text = "lon lat
                             1  -3.00  44
                             4  -2.75  44
                             7  -2.50  44
                             10 -2.25  44
                             13 -2.00  44
                             16 -1.75  44")

## "garbage" values, just to pass the variable names checking in `predict.gam`
newdat[c("b.depth", "c.dist", "temp.20m", "log.net.area")] <- 0

## prediction on the link scale
pred_link <- predict(gm2, newdata = newdat, type = "terms", terms = "s(lon,lat)")
#   s(lon,lat)
#1  -1.9881967
#4  -1.9137971
#7  -1.6365945
#10 -1.1247837
#13 -0.7910023
#16 -0.7234683
#attr(,"constant")
#(Intercept) 
#   2.553535 

## simplify to vector
pred_link <- attr(pred_link, "constant") + rowSums(pred_link)
#[1] 0.5653381 0.6397377 0.9169403 1.4287511 1.7625325 1.8300665

## prediction on the response scale
pred_response <- gm2$family$linkinv(pred_link)
#[1] 1.760043 1.895983 2.501625 4.173484 5.827176 6.234301

Обычно я не использую predict.gam, если хочу сделать прогноз для определенного гладкого термина. Логика predict.gam заключается в том, чтобы сначала сделать прогноз для всех терминов, то есть так же, как вы делаете type = "terms". потом

  • если type = "link", сделать rowSums для всех прогнозов по срокам плюс перехват (возможно, с offset);
  • если type = "terms" и "terms" или "exclude" не указаны, вернуть результат как есть;
  • если type = "terms" и вы указали "terms" и/или "exclude", выполняется некоторая постобработка, чтобы удалить термины, которые вам не нужны, и дать вам только те, которые вам нужны.

Таким образом, predict.gam всегда будет выполнять вычисления для всех терминов, даже если вам нужен только один термин.

Зная неэффективность этого, я сделаю следующее:

sm <- gm2$smooth[[1]]  ## extract smooth construction info for `s(lon,lat)`
Xp <- PredictMat(sm, newdat)  ## predictor matrix
b <- gm2$coefficients[with(sm, first.para:last.para)]  ## coefficients for this term
pred_link <- c(Xp %*% b) + gm2$coef[[1]]  ## this term + intercept
#[1] 0.5653381 0.6397377 0.9169403 1.4287511 1.7625325 1.8300665
pred_response <- gm2$family$linkinv(pred_link)
#[1] 1.760043 1.895983 2.501625 4.173484 5.827176 6.234301

Видите ли, мы получаем тот же результат.

Не будет ли результат зависеть от значения, присвоенного ковариатам (здесь 0)?

Некоторое предсказание мусора будет сделано для этих значений мусора, но predict.gam отбрасывает их в конце.

Спасибо, вы правы. Я не совсем уверен, что понимаю, почему тогда есть возможность добавлять значения ковариатов в новых местах.

Сопровождение кода, насколько я понимаю, очень сложно для такого большого пакета, как mgcv. Код необходимо значительно изменить, если вы хотите, чтобы он соответствовал потребностям каждого пользователя. Очевидно, что описанная здесь логика predict.gam будет неэффективной, когда люди вроде вас просто хотят, чтобы она предсказывала определенное сглаживание. И теоретически, если это так, проверка имен переменных в newdata может игнорировать те термины, которые не нужны пользователям. Но это требует значительного изменения predict.gam и потенциально может привести к множеству ошибок из-за изменений кода. Кроме того, вы должны отправить журнал изменений в CRAN, и CRAN может просто не обрадоваться такому радикальному изменению.

Саймон однажды поделился своими чувствами: Многие говорят мне, что я должен писать mgcv так или так, но я просто не могу. Да, посочувствуйте такому автору/сопровождающему пакета, как он.

Спасибо за ответ об обновлении. Однако я не понимаю, почему прогнозы не зависят от значений ковариат в новых местах.

Это будет зависеть от того, предоставите ли вы значения ковариатов для b.depth, c.dist, temp.20m, log.net.area. Но так как у вас нет их в новых местах, прогноз состоит в том, чтобы предположить, что эти эффекты будут 0.

Хорошо, спасибо, теперь я вижу! Итак, правильно ли было бы сказать, что в отсутствие значений ковариат в новых местах я только предсказываю реакцию на основе пространственной автокорреляции остатков?

Вы только предсказываете пространственное поле/гладкое. В подходе GAM пространственное поле моделируется как часть среднего, а не как дисперсия-ковариация (как в кригинге), поэтому я думаю, что ваше использование «остатков» здесь неверно.

Да, ты прав. Просто чтобы понять, что делает этот код: правильно ли будет сказать, что я предсказываю, как отклик меняется в пространстве, но не его фактические значения в новых местах (поскольку для этого мне понадобятся значения ковариат в этих местах)?

Правильный. Вы можете попробовать predict.gam с terms = "s(lon,lat)" или без него, чтобы помочь вам переварить вывод. Посмотрите, как она меняется, когда вы меняете значения мусора, передаваемые другим ковариатам.

## a possible set of garbage values for covariates
newdat[c("b.depth", "c.dist", "temp.20m", "log.net.area")] <- 0

predict(gm2, newdat, type = "terms")
#   s(lon,lat) s(I(b.depth^0.5)) s(c.dist) s(temp.20m)
#1  -1.9881967          -1.05514 0.4739174   -1.466549
#4  -1.9137971          -1.05514 0.4739174   -1.466549
#7  -1.6365945          -1.05514 0.4739174   -1.466549
#10 -1.1247837          -1.05514 0.4739174   -1.466549
#13 -0.7910023          -1.05514 0.4739174   -1.466549
#16 -0.7234683          -1.05514 0.4739174   -1.466549
#attr(,"constant")
#(Intercept) 
#   2.553535 

predict(gm2, newdat, type = "terms", terms = "s(lon,lat)")
#   s(lon,lat)
#1  -1.9881967
#4  -1.9137971
#7  -1.6365945
#10 -1.1247837
#13 -0.7910023
#16 -0.7234683
#attr(,"constant")
#(Intercept) 
#   2.553535 

## another possible set of garbage values for covariates
newdat[c("b.depth", "c.dist", "temp.20m", "log.net.area")] <- 1
#   s(lon,lat) s(I(b.depth^0.5))  s(c.dist) s(temp.20m)
#1  -1.9881967        -0.9858522 -0.3749018   -1.269878
#4  -1.9137971        -0.9858522 -0.3749018   -1.269878
#7  -1.6365945        -0.9858522 -0.3749018   -1.269878
#10 -1.1247837        -0.9858522 -0.3749018   -1.269878
#13 -0.7910023        -0.9858522 -0.3749018   -1.269878
#16 -0.7234683        -0.9858522 -0.3749018   -1.269878
#attr(,"constant")
#(Intercept) 
#   2.553535 

predict(gm2, newdat, type = "terms", terms = "s(lon,lat)")
#   s(lon,lat)
#1  -1.9881967
#4  -1.9137971
#7  -1.6365945
#10 -1.1247837
#13 -0.7910023
#16 -0.7234683
#attr(,"constant")
#(Intercept) 
#   2.553535
person Zheyuan Li    schedule 01.08.2018