具有“gp”更平滑的 GAM：在新位置进行预测

Question

我正在使用以下地理附加模型

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。

Answer 1

predict仍然要求模型中使用的所有变量都在 中呈现newdata，但您可以将一些任意值（例如 s）传递0给您没有的协变量，然后使用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"进行 a并加上截距（可能使用）；rowSumsoffset
• 如果type = "terms"、 和"terms"or"exclude"未指定，则按原样返回结果；
• 如果type = "terms"和 您指定了"terms"and/or "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