Как генерировать основные эффекты для обобщенных линейных смешанных моделей glmmTMB - PullRequest
Новая
       84

Как генерировать основные эффекты для обобщенных линейных смешанных моделей glmmTMB

0 голосов
/ 09 января 2019

Я провел серию обобщенных линейных смешанных моделей с использованием glmmTMB, чтобы оценить различия в поведении животных на разных объектах. Насыщенная модель включает полный набор основных эффектов, один случайный эффект (случайное пересечение), один вложенный эффект и два непрерывных предиктора (см. Код ниже). Я хочу оценить значимость основных эффектов в модели (три основных эффекта, представляющих два основных эффекта и их эффект взаимодействия: поведение, объект, поведение: объект) и использую anova () в базовом пакете для сравнения моделей друг с другом но продолжайте получать значения хи-квадрат и степени свободы, которые, кажется, не имеют смысла. Ни одна из моделей не выдает ошибок при запуске самостоятельно, но сравнения с использованием anova () не имеют смысла. Я не понимаю, почему или что я должен делать, чтобы обойти это.

Учитывая, что модели были запущены с использованием пакета glmmTMB в R, я не могу использовать функцию Anova () из пакета car, поскольку объекты модели glmmTMB не поддерживаются, и при этом генерируются сообщения об ошибках. Я также пытался использовать пакет glmmADMB, но обнаружил, что он намного медленнее (некоторые модели буквально отнимают часы), и сравнение моделей было таким же бессмысленным.

Мой код: 

rm()
detach()
library(glmmTMB)
dataset<-read.csv("Dur18Bletters.csv", header=T)
attach(dataset)
names(dataset)

mod1<-glmmTMB(Count~Facility+Beh+Facility:Beh+(1|Sex)+(1|ID/Day)+Age+Attendance,family=nbinom2)
mod2<-glmmTMB(Count~Facility+Beh+(1|Sex)+(1|ID/Day)+Age+Attendance,family=nbinom2)
mod3<-glmmTMB(Count~Facility+Facility:Beh+(1|Sex)+(1|ID/Day)+Age+Attendance,family=nbinom2)
mod4<-glmmTMB(Count~Beh+Facility:Beh+(1|Sex)+(1|ID/Day)+Age+Attendance,family=nbinom2)

anova(mod1,mod2,test="LRT")
anova(mod1,mod3,test="LRT")
anova(mod1,mod4,test="LRT")

Результаты трех тестов на отклонения: 

Data: NULL
Models:
mod2: Count ~ Facility + Beh + (1 | Sex) + (1 | ID/Day) + Age + People, zi=~0, disp=~1
mod1: Count ~ Facility + Beh + Facility:Beh + (1 | Sex) + (1 | ID/Day) + , zi=~0, disp=~1
mod1:     Age + People, zi=~0, disp=~1
     Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
mod2 26 19672 19823 -9810.0    19620                             
mod1 60 19447 19795 -9663.6    19327 292.69     34  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(mod1,mod3,test="LRT")
Data: NULL
Models:
mod1: Count ~ Facility + Beh + Facility:Beh + (1 | Sex) + (1 | ID/Day) + , zi=~0, disp=~1
mod1:     Age + People, zi=~0, disp=~1
mod3: Count ~ Facility + Facility:Beh + (1 | Sex) + (1 | ID/Day) + , zi=~0, disp=~1
mod3:     Age + People, zi=~0, disp=~1
     Df   AIC   BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
mod1 60 19447 19795 -9663.6    19327                        
mod3 60 19447 19795 -9663.6    19327     0      0          1
> anova(mod1,mod4,test="LRT")
Data: NULL
Models:
mod1: Count ~ Facility + Beh + Facility:Beh + (1 | Sex) + (1 | ID/Day) + , zi=~0, disp=~1
mod1:     Age + People, zi=~0, disp=~1
mod4: Count ~ Beh + Facility:Beh + (1 | Sex) + (1 | ID/Day) + Age + , zi=~0, disp=~1
mod4:     People, zi=~0, disp=~1
     Df   AIC   BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
mod1 60 19447 19795 -9663.6    19327                            
mod4 60 19447 19795 -9663.6    19327     0      0  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Когда я проверяю сравнение первой модели (т.е. anova (mod1, mod2)), результат кажется разумным и заслуживающим доверия. Таким образом, я бы пришел к выводу, что основной эффект, который был удален в mod2, внес существенный вклад в модель, о чем свидетельствует существенное различие между двумя моделями. Невозможно, чтобы значение хи-квадрат 0 генерировало значение р <0,0001. </p>

Любая помощь, которую вы можете оказать, будет принята с благодарностью! Я не особенно привержен этому способу оценки основных эффектов (это единственный метод, который я мог найти при просмотре форумов), поэтому, если у кого-то есть альтернативный пакет / метод, который я мог бы попробовать, это было бы здорово! Я в полной растерянности.

...