[R] Predicting random effects in GAM using mgcv

Meagan Dunphy-Daly m.dunphydaly at gmail.com
Sat Jan 3 22:47:11 CET 2015


Hi all, 

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 (which should allow me to predict with the vessel effects at their predicted values or zero), where "dum" is a vector of 1.

The model runs, but I am having problems predicting.

grid.bin.GOM_A_Swordfish <- data.frame("Distance"=seq(min(GOM$Distance),max(GOM$Distance),length = 100),#range of distances over which to predict
                                 "SetYear" = '2006',
				 "SetMonth" = '6',
                                 "TimePeriod" = 'A',
                                 "SST" = mean(GOM$SST),
                                 "VesselID" = 'MEDIHO', #I picked one of the levels of the random effect, Vessel MEDIHO
                                 "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. Can you provide any advice on how to predict this model without the VesselID term (but still include it in fitting)?

Thanks!
Meagan


More information about the R-help mailing list