drop1 с объектом GLMMtmb возвращает ошибку: «область не является подмножеством меток терминов» - PullRequest
Новая
       20

drop1 с объектом GLMMtmb возвращает ошибку: «область не является подмножеством меток терминов»

0 голосов
/ 29 октября 2019

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

Из прочтения других источников может показаться, что drop1(model,.~.) даст мне то, что мне нужно, но я получу ошибку: 

"Error in 'drop1.default(model, . ~ .)' : 
  scope is not a subset of term labels"

Я сделаю фиктивный набор данных видовЗначения богатства на графиках в пределах двух уровней обработки в течение трех лет, чтобы продемонстрировать. 

set.seed(1234)
spRich = c(rpois(10,5),rpois(10,8),rpois(10,12),
    rpois(10,12),rpois(10,14),rpois(10,15))
Treatment = c(rep("A",times=30),rep("B",times=30))
Year = rep(rep(c("One","Two","Three"),each=10),times=2)
Site = as.factor(rep(c(1,2,3,4,5,6,7,8,9,10),times=6))
df = data.frame(Site = Site, Plot = paste0(Site,Treatment),
    Year = Year, Treatment = Treatment, Richness = spRich)

Простой график показывает, что, по-видимому, существует взаимодействие между годами и обработками 

ggplot(df,aes(x=Year,y=Richness,color=Treatment)) + theme_classic() + 
    geom_boxplot()

А затем я могу выполнить glmm, используя обобщенное распределение Пуассона 

library(glmmTMB)
library(car)
options(contrasts = c('contr.sum','contr.poly'))
model = glmmTMB(Richness ~ Year + Treatment + Year:Treatment + (1|Site/Plot),
    data=df, family = "genpois")
summary(model)
Anova(model,type=3,test="Chisq")

Вывод из таблицы Anova выглядит следующим образом: 

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: Richness
                  Chisq Df Pr(>Chisq)    
(Intercept)    1652.699  1  < 2.2e-16 ***
Year             20.951  2  2.822e-05 ***
Treatment        62.375  1  2.838e-15 ***
Year:Treatment   13.569  2   0.001131 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ах, я чувствую подтверждение того, что был эффект взаимодействия! Но это статистика теста Вальда, и Вартон (2016) указывает, что регрессия Пуассона со статистикой Вальда по данным о богатстве видов в GLM может иметь коэффициент ошибок типа I ~ 0,20, что будет означать, что да, это кажется значительным, но достаточно близковызывать беспокойство.

Теперь, если я пытаюсь получить статистику отношения правдоподобия, я делаю: 

drop1(model,.~., test="Chisq")

Я получаю сообщение об ошибке.

Я ожидаю, что получу таблицу, аналогичную вызову Anova, как в случае с lm, переданным через drop1. Это смешная вещь, чтобы ожидать? Есть ли способ провести тест для каждого эффекта в отдельности, чтобы я мог объединить результаты в одну таблицу? Если не для проверки отношения правдоподобия, то существует ли процедура повторной выборки или начальной загрузки, которая дает аналогичные результаты для основных эффектов при использовании glmmTMB?

...