Przewidywanie losowych efektów w mgcv gam

10

Interesuje mnie modelowanie całkowitego połowu ryb za pomocą gam w mgcv do modelowania prostych efektów losowych dla poszczególnych statków (które odbywają wielokrotne podróże w czasie na łowisku). Mam 98 przedmiotów, więc pomyślałem, że użyję gam zamiast gamma do modelowania efektów losowych. Mój model to:

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")

Zakodowałem losowy efekt za pomocą bs = "re" i by = dum (czytałem, że to pozwoli mi przewidzieć z efektami naczynia przy ich przewidywanych wartościach lub zero). „dum” jest wektorem 1.

Model działa, ale mam problemy z przewidywaniem. Wybrałem jeden z naczyń do prognoz (Vessel21) i średnie wartości dla wszystkiego innego oprócz predyktora zainteresowania prognozami (Odległość).

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)

Występuje błąd:

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

Myślę, że jest to wywoływane, ponieważ VesselID jest czynnikiem, ale używam go jako płynnego dla losowych efektów.

Udało mi się z powodzeniem przewidzieć za pomocą gam bez prostych efektów losowych (bs = „re”).

Czy możesz udzielić porady, jak przewidzieć ten model bez określenia VesselID (ale nadal uwzględniać go w dopasowaniu)?

Dziękuję Ci!

Meagan
źródło

Odpowiedzi:

20

Od wersji 1.8.8 mgcv predict.gam zyskał excludeargument, który pozwala na zerowanie pojęć w modelu, w tym efekty losowe podczas przewidywania, bez sztuczki sztucznej, która została wcześniej zasugerowana.

  • predict.gami predict.bamteraz zaakceptuj 'exclude'argument umożliwiający wyzerowanie terminów (np. efektów losowych) w celu przewidywania. Dla skuteczności, gładka warunki nie w termslub excludenie są już ocenione, a zamiast tego są ustawione na zero lub nie powrócił. Zobaczyć ?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

Starsze podejście

Simon Wood użył następującego prostego przykładu, aby sprawdzić, czy to działa:

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

Który działa dla mnie. Również:

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))

działa również.

Chciałbym więc sprawdzić, czy dane, które podajesz, są newdatatakie, jak myślisz, ponieważ problem może nie wynikać z VesselID- błąd pochodzi z funkcji, która zostałaby wywołana przez predict()wywołania w powyższych przykładach i Rail jest czynnikiem pierwszy przykład.

Gavin Simpson
źródło
Dziękuję Gavin za przykłady! Opracowując je, doszedłem do tego. Miałeś rację - błąd występował w ramce danych newdata. Po usunięciu cudzysłowu wokół „0” dla zmiennej „dum” według zmiennej mogłem przewidzieć bez żadnych błędów. Błąd nowicjusza, ale zmagałem się z tym przez cały dzień i myślałem, że to problem z płynnością czynnika VesselID. Dziękuję bardzo!
Meagan
Jak określić więcej niż jeden losowy efekt, który należy wykluczyć exclude? Próbowałem użyć, c()ale nie działa.
Stefano
Użycie wektora terminów do wykluczenia działa dla mnie: exclude = c("s(x0)", "s(x2)")powiedzmy z poniższego modelu b<-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat)z ?predict.gamprzykładów. Musisz podać ciągi znaków w wektorze przekazywanym excludeza pomocą notacji używanej summary()podczas wyświetlania informacji o każdym gładkim terminie
Gavin Simpson