ГАМ с "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).
1 ответ
predict все еще требует, чтобы все переменные, используемые в вашей модели, были представлены в newdata, но вы можете передать некоторые произвольные значения, такие как 0 s, для тех ковариат у вас нет, затем используйте 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