Predicting with random effects in mgcv gam

Question

I am interested in modeling total fish catch using gam in mgcv to model simple random effects for individual vessels (that make repeated trips over time in the fishery). I have 98 subjects, so I thought I would use gam instead of gamm to model the random effects. My model is:

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

I have coded the random effect with bs = "re" and by = dum (I read that this would allow me to predict with the vessel effects at their predicted values or zero). "dum" is a vector of 1.

The model runs, but I am having problems predicting. I picked one of the vessels for the predictions (Vessel21) and average values for everything else except the predictor of interest for predictions (Distance).

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)

The error that I'm getting is:

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

I think this is being called because VesselID is a factor, but I'm using it a smooth for the random effects.

I have been able to successfully predict using gam without the simple random effects (bs = "re").

Can you provide any advice on how to predict this model without the VesselID term (but still include it in fitting)?

Thank you!

Gavin Simpson · Accepted Answer · 2016-02-17T15:40:53.993

From version 1.8.8 of mgcv predict.gam has gained an exclude argument which allows for the zeroing out of terms in the model, including random effects, when predicting, without the dummy trick that was suggested previously.

predict.gam and predict.bam now accept an 'exclude' argument allowing terms (e.g. random effects) to be zeroed for prediction. For efficiency, smooth terms not in terms or in exclude are no longer evaluated, and are instead set to zero or not returned. See ?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

Older approach

Simon Wood has used the following simple example to check this is working:

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

Which works for me. Likewise:

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

also works.

So I would check the data you are supplying in newdata is what you think it is as the problem may not be with VesselID — the error is coming from the function that would have been called by the predict() calls in the examples above, and Rail is a factor in the first example.

Thank you, Gavin, for the examples! In working through those, I figured it out. You were correct - the error was in the newdata data frame. Once I removed the quotation marks around '0' for the "dum" by variable, I was able to predict without any errors. Rookie mistake, but I had been struggling with it all day and thought it was a problem with the VesselID factor being a smooth. Thank you so much! — Meagan, Jan 04 '15 at 02:07
How can one specify more than one random effect to exclude with `exclude`? I tried using `c()` but it doesn't seem to work. — Stefano, Jan 01 '18 at 21:17
Using a vector of terms to exclude works for me: `exclude = c("s(x0)", "s(x2)")` say from the following model `b — Gavin Simpson, Jan 01 '18 at 21:54

Predicting with random effects in mgcv gam

1 Answers1

Older approach

Linked