predict
仍然需要模型中使用的所有变量
newdata
,但可以传入一些任意值,如
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
如果我想预测一个特定的平滑项逻辑
预测.gam
先做所有的预测,也就是说,和你做的一样
type=“条款”
. 那么
-
如果
type = "link"
,做一个
rowSums
在所有的术语预测加上截距(可能是
offset
);
-
如果
type=“条款”
,和
"terms"
或
"exclude"
未指定,返回结果;
-
如果
type=“条款”
你已经指定
“条款”
和/或
“排除”
,完成了一些后期处理,以删除不想要的条款,并只给您想要的条款。
所以,
预测.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)?
将在这些垃圾值上进行一些垃圾预测,但是
预测.gam
最后丢弃它们。
谢谢,你说得对。我不完全确定为什么会有在新位置添加协变量值的选项。
就我看来,代码维护对于像
mgcv
. 如果您想让代码适合每个用户的需要,那么需要对其进行显著的更改。很明显
预测.gam
我在这里所描述的逻辑将是低效的,当人们,像你一样,只想预测某个平稳。理论上,如果是这种情况,变量名将签入
新数据
可以忽略那些用户不想要的术语。但是,这需要
预测.gam
,并且可能由于代码更改而引入许多错误。此外,你必须向克兰提交一个变更日志,克兰可能只是不高兴看到这种剧烈的变化。
西蒙曾经分享过他的感受:
有很多人告诉我,我应该写
mgcv公司
不管是这样还是那样,但我不能
. 是的,向他这样的软件包作者/维护者表示同情。
谢谢你的更新回复。然而,我不明白为什么预测不依赖于新位置的协变量值。
它将取决于您是否为
b.depth
,
c.dist
,
temp.20m
,
log.net.area
. 但是由于你没有在新的地点使用它们,所以预测只是假设这些影响是
零
.
好的,谢谢,我现在明白了!所以说在没有新位置的协变量值的情况下,我只是根据残差的空间自相关来预测响应,这是正确的吗?
你只是在预测空间场/平滑。在GAM方法中,空间场被建模为均值的一部分,而不是方差协方差(如kriging),所以我认为这里使用“残差”是不正确的。
是的,你说得对。为了理解这段代码的作用:如果说我预测的是响应在空间上的变化,而不是在新位置的实际值,这是否正确(因为我需要这些位置的协变量值)?
对的。你可以试试
预测.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
