时间戳

首页 > 解决方案 > 带有“gp”的 GAM 更平滑：在新位置进行预测

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

问题描述

我正在使用以下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")

如何使用它来预测我没有协变量数据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）。

标签： rregressiongammgcvkriging

解决方案

predict仍然需要将模型中使用的所有变量都显示在 中newdata，但是您可以将一些任意值（例如0s）传递给那些您没有的协变量，然后使用type = "terms"andterms = 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", and "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

推荐阅读