在 mgcv gam 中使用随机效应进行预测

机器算法验证 预言 随机效应模型 广义加法模型 毫克CV
2022-03-17 00:48:34

我有兴趣使用 mgcv 中的 gam 对总渔获量进行建模，以模拟单个船只的简单随机效应（随着时间的推移在渔业中重复旅行）。我有 98 个科目，所以我想我会使用 gam 而不是 gamm 来模拟随机效应。我的模型是：

modelGOM <- gam(TotalFish ~ factor(SetYear) + factor(SetMonth) + factor(TimePeriod) +     
s(SST) + s(VesselID, bs = "re", by = dum) + s(Distance, by = TimePeriod) + 
offset(log(HooksSet)), data = GOM, family = tw(), method = "REML")

我已经用 bs = "re" 和 by = dum 对随机效应进行了编码（我读到这将允许我用它们的预测值或零来预测血管效应）。“dum”是 1 的向量。

模型运行，但我在预测时遇到问题。我选择了其中一艘船进行预测（Vessel21），并为除了预测感兴趣的预测器（距离）之外的所有其他东西选择平均值。

data.frame("Distance"=seq(min(GOM$Distance),max(GOM$Distance),length = 100),
                             "SetYear" = '2006',
                             "SetMonth" = '6',
                             "TimePeriod" = 'A',
                             "SST" = mean(GOM$SST),
                             "VesselID" = 'Vessel21', 
                             "dum" = '0', #to predict without vessel effect
                             "HooksSet" = mean(GOM$HooksSet))

pred_GOM_A_Swordfish <- predict(modelGOM, grid.bin.GOM_A_Swordfish, type = "response", 
se = T)

我得到的错误是：

Error in Predict.matrix.tprs.smooth(object, dk$data) : 
    NA/NaN/Inf in foreign function call (arg 1)
    In addition: Warning message:
    In Ops.factor(xx, object$shift[i]) : - not meaningful for factors

我认为这是因为 VesselID 是一个因素而被调用，但我正在使用它来平滑随机效果。

我已经能够成功地预测使用 gam 而没有简单的随机效应（bs =“re”）。

您能否就如何在没有 VesselID 术语的情况下预测该模型提供任何建议（但仍将其包含在拟合中）？

谢谢！

1个回答

从mgcv predict.gam的 1.8.8 版开始，获得了一个exclude参数，该参数允许在预测时将模型中的项（包括随机效应）归零，而无需使用之前建议的虚拟技巧。

  • predict.gam现在predict.bam接受一个'exclude'参数，允许将术语（例如随机效应）归零以进行预测。为了提高效率，不再评估不在terms或不在的平滑项，而是将其设置为零或不返回。exclude见?predict.gam。
library("mgcv")
require("nlme")
dum <- rep(1,18)
b1 <- gam(travel ~ s(Rail, bs="re", by=dum), data=Rail, method="REML")
b2 <- gam(travel ~ s(Rail, bs="re"), data=Rail, method="REML")

head(predict(b1, newdata = cbind(Rail, dum = dum)))    # ranefs on
head(predict(b1, newdata = cbind(Rail, dum = 0)))      # ranefs off
head(predict(b2, newdata = Rail, exclude = "s(Rail)")) # ranefs off, no dummy

> head(predict(b1, newdata = cbind(Rail, dum = dum)))    # ranefs on
       1        2        3        4        5        6 
54.10852 54.10852 54.10852 31.96909 31.96909 31.96909  
> head(predict(b1, newdata = cbind(Rail, dum = 0)))      # ranefs off
   1    2    3    4    5    6 
66.5 66.5 66.5 66.5 66.5 66.5
> head(predict(b2, newdata = Rail, exclude = "s(Rail)")) # ranefs off, no dummy
   1    2    3    4    5    6 
66.5 66.5 66.5 66.5 66.5 66.5

较旧的方法

Simon Wood 使用以下简单示例来检查它是否有效：

library("mgcv")
require("nlme")
dum <- rep(1,18)
b <- gam(travel ~ s(Rail, bs="re", by=dum), data=Rail, method="REML")
predict(b, newdata=data.frame(Rail="1", dum=0)) ## r.e. "turned off"
predict(b, newdata=data.frame(Rail="1", dum=1)) ## prediction with r.e

这对我有用。同样地：

dum <- rep(1, NROW(na.omit(Orthodont)))
m <- gam(distance ~ s(age, bs = "re", by = dum) + Sex, data = Orthodont)
predict(m, data.frame(age = 8, Sex = "Female", dum = 1))
predict(m, data.frame(age = 8, Sex = "Female", dum = 0))

也有效。

因此，我会检查您提供的数据是否是newdata您认为的数据，因为问题可能不存在VesselID- 错误来自predict()上面示例中的调用将调用的函数，并且 Rail是导致第一个例子。

其它你可能感兴趣的问题
上一篇如何通过反向传播训练 SVM？ 下一篇是否有标准化的偏度和峰度等价物？