Prever com efeitos aleatórios em mgcv gam

10

Estou interessado em modelar a captura total de peixes usando gam em mgcv para modelar efeitos aleatórios simples para navios individuais (que fazem viagens repetidas ao longo do tempo na pesca). Eu tenho 98 sujeitos, então pensei em usar gam em vez de gamm para modelar os efeitos aleatórios. Meu modelo é:

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

Eu codifiquei o efeito aleatório com bs = "re" e by = dum (eu li que isso me permitiria prever com os efeitos dos vasos em seus valores previstos ou zero). "dum" é um vetor de 1.

O modelo é executado, mas estou tendo problemas para prever. Escolhi um dos navios para as previsões (navio 21) e os valores médios para todo o resto, exceto o preditor de interesse pelas previsões (distância).

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)

O erro que estou recebendo é:

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

Acho que isso está sendo chamado porque o VesselID é um fator, mas estou usando-o de maneira suave para os efeitos aleatórios.

Consegui prever com sucesso o uso de gam sem os simples efeitos aleatórios (bs = "re").

Você pode fornecer algum conselho sobre como prever esse modelo sem o termo VesselID (mas ainda incluí-lo no ajuste)?

Obrigado!

Meagan
fonte

Respostas:

20

A partir da versão 1.8.8 do mgcv predict.gam ganhou um excludeargumento que permite zerar os termos no modelo, incluindo efeitos aleatórios, ao prever, sem o truque fictício sugerido anteriormente.

  • predict.game predict.bamagora aceite um 'exclude'argumento que permita que os termos (por exemplo, efeitos aleatórios) sejam zerados para a previsão. Para eficiência, termos suaves que não estão dentro termsou excludenão são mais avaliados e, em vez disso, são definidos como zero ou não retornados. Veja ?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

Abordagem mais antiga

Simon Wood usou o seguinte exemplo simples para verificar se está funcionando:

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

O que funciona para mim. Da mesma forma:

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

também funciona.

Portanto, eu verificaria se os dados que você está fornecendo newdatasão o que você acha que é, pois o problema pode não estar ocorrendo VesselID- o erro está vindo da função que seria chamada pelas predict()chamadas nos exemplos acima e Rail é um fator no primeiro exemplo.

Gavin Simpson
fonte
Obrigado, Gavin, pelos exemplos! Ao trabalhar com isso, eu descobri. Você estava correto - o erro estava no quadro de dados de novos dados. Depois que removi as aspas em torno de '0' para o "dum" por variável, fui capaz de prever sem erros. Erro de novato, mas eu estava lutando com isso o dia todo e pensei que era um problema com o fator VesselID sendo suave. Muito obrigado!
Meagan
Como se pode especificar mais de um efeito aleatório para excluir exclude? Eu tentei usar, c()mas parece não funcionar.
Stefano
Usar um vetor de termos para excluir obras para mim: exclude = c("s(x0)", "s(x2)")diga a partir do seguinte modelo a b<-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat)partir de ?predict.gamexemplos. Você precisa especificar as cordas do vetor passado para excludecom a notação usada por summary()ao exibir as informações sobre cada termo suave
Gavin Simpson