Prédire avec des effets aléatoires dans mgcv gam

Je m'intéresse à la modélisation des prises totales de poisson en utilisant gam en mgcv pour modéliser des effets aléatoires simples pour des navires individuels (qui effectuent des déplacements répétés au fil du temps dans la pêche). J'ai 98 sujets, j'ai donc pensé utiliser gam au lieu de gamm pour modéliser les effets aléatoires. Mon modèle est:

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

J'ai codé l'effet aléatoire avec bs = "re" et by = dum (j'ai lu que cela me permettrait de prédire avec les effets du vaisseau à leurs valeurs prédites ou zéro). "dum" est un vecteur de 1.

Le modèle fonctionne, mais j'ai du mal à prévoir. J'ai choisi l'un des navires pour les prévisions (Vessel21) et les valeurs moyennes pour tout le reste, sauf le prédicteur d'intérêt pour les prévisions (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)

L'erreur que je reçois est:

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

Je pense que cela est appelé parce que VesselID est un facteur, mais je l'utilise de manière fluide pour les effets aléatoires.

J'ai pu prédire avec succès l'utilisation de gam sans les effets aléatoires simples (bs = "re").

Pouvez-vous fournir des conseils sur la façon de prédire ce modèle sans le terme VesselID (mais toujours l'inclure dans l'ajustement)?

Je vous remercie!

De la version 1.8.8 de mgcv predict.gam a gagné un excludeargument qui permet la mise à zéro des termes dans le modèle, y compris les effets aléatoires, lors de la prédiction, sans l'astuce factice qui a été suggérée précédemment.

• predict.gamet predict.bammaintenant accepter un 'exclude'argument permettant de mettre à zéro les termes (par exemple les effets aléatoires) pour la prédiction. Pour plus d'efficacité, les termes lisses qui ne sont ni dans termsni dans excludene sont plus évalués et sont plutôt définis sur zéro ou non renvoyés. Tu vois ?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

Une approche plus ancienne

Simon Wood a utilisé l'exemple simple suivant pour vérifier que cela fonctionne:

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

Ce qui fonctionne pour moi. Également:

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

fonctionne également.

Donc, je vérifierais que les données que vous fournissez newdatasont ce que vous pensez que le problème peut ne pas être VesselID- l'erreur vient de la fonction qui aurait été appelée par les predict()appels dans les exemples ci-dessus, et Rail est un facteur dans la premier exemple.