在 mgcv gam 中使用隨機效應進行預測

我有興趣使用 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 術語的情況下預測該模型提供任何建議（但仍將其包含在擬合中）？

謝謝！

從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是第一個例子。

引用自：https://stats.stackexchange.com/questions/131106