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

10

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!

Meagan
la source

Réponses:

20

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.

Gavin Simpson
la source
Merci Gavin pour les exemples! En travaillant à travers ceux-ci, je l'ai compris. Vous aviez raison - l'erreur était dans le bloc de données newdata. Une fois que j'ai supprimé les guillemets autour de «0» pour la variable «dum» par variable, j'ai pu prédire sans aucune erreur. Erreur de débutant, mais j'avais eu du mal avec ça toute la journée et je pensais que c'était un problème avec le facteur VesselID étant un bon fonctionnement. Merci beaucoup!
Meagan
Comment spécifier plus d'un effet aléatoire à exclure exclude? J'ai essayé d'utiliser c()mais cela ne semble pas fonctionner.
Stefano
L'utilisation d'un vecteur de termes pour exclure fonctionne pour moi: exclude = c("s(x0)", "s(x2)")dites du modèle suivant à b<-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat)partir d' ?predict.gamexemples. Vous devez spécifier les chaînes dans le vecteur transmis excludeavec la notation utilisée par summary()lors de l'affichage des informations sur chaque terme lisse
Gavin Simpson