drop1 с объектом GLMMtmb возвращает ошибку: "область действия не является подмножеством меток терминов"
Я пытался найти способ получить результаты теста отношения правдоподобия, когда использую 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)
Простой график показывает, что, по-видимому, существует взаимодействие лечения x год
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
Ах, я чувствую себя подтвержденным, что был эффект взаимодействия! Но это статистика теста Вальда, и Warton (2016) указывает, что регрессия Пуассона со статистикой Вальда по данным о видовом богатстве в GLM может иметь частоту ошибок типа I ~ 0,20, что означает, что да, это кажется значительным, но достаточно близко. вызывать беспокойство.
Теперь, если я попытаюсь получить статистику отношения правдоподобия, я сделаю:
drop1(model,.~., test="Chisq")
Я получаю сообщение об ошибке.
Я ожидал, что получу таблицу, похожую на вызов Anova, как и в случае с lm прошел через drop1. Разве этого ожидать нелепо? Есть ли способ провести тест для каждого эффекта индивидуально, чтобы я мог объединить результаты в одну таблицу? Если не тесты отношения правдоподобия, то существует ли процедура повторной выборки или начальной загрузки, которая дает аналогичные результаты для основных эффектов при использовании glmmTMB?