Сообщение об ошибке при использовании метлы для получения коэффициентов из моделей нулевой инфляции glmmTMB

Используя данные Owls и пакет glmmTMB, я хочу визуально сравнить коэффициенты регрессии из разных моделей с нулевым завышением, которые отличаются используемым семейством (ZIPOISS, ZINB1, ZINB2) и без/без смещения (logBroodSize).

Однако моя первая проблема - получить коэффициенты. Функция tidy из пакета broom должна предоставить вам коэффициенты для их последующего построения с помощью ggplot, но при попытке получить их я получаю следующую ошибку:

modList= list(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)  
coefs= ldply(modList,tidy,effect="fixed",conf.int=TRUE,
                  .id="model") %>%
    mutate(term=abbfun(term)) %>%
    select(model,term,estimate,conf.low,conf.high) %>%
    filter(!term %in% c("Intercept","Intercept.1","NCalls","zi_NCalls"))

Error in as.data.frame.default(x) :
cannot coerce class ""glmmTMB"" to a data.frame
Also: Warning message:
In tidy.default(X[[i]], ...) :
No method for tidying an S3 object of class glmmTMB , using as.data.frame

Есть идеи, что может быть не так? Мне уже сказали, что причиной может быть отсутствие правильной версии broom, однако я установил ее правильную версию... Далее приведен код для воспроизводимого примера:

    # Packages and dataset
    library(glmmTMB)
    library(broom) # devtools::install_github("bbolker/broom")
    library(plyr)
    library(dplyr)
    data(Owls,package="glmmTMB")
    Owls = plyr::rename(Owls, c(SiblingNegotiation="NCalls"))
    Owls = transform(Owls,
             ArrivalTime=scale(ArrivalTime, center=TRUE, scale=FALSE),
             obs=factor(seq(nrow(Owls))))

    # Models
    zipoiss<-glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
    offset(logBroodSize)+(1|Nest),
    data=Owls,
    ziformula = ~ 1,
    family="poisson")

    zinb2<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
    offset(logBroodSize)+(1|Nest),
    data=Owls,
    ziformula = ~ 1,
    family="nbinom2")

    zinb1 <- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
    offset(logBroodSize)+(1|Nest),
    data=Owls,
    ziformula = ~ 1,
    family="nbinom1")

    zinb1_bs<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+ 
    BroodSize+(1|Nest),
    data=Owls,
    ziformula = ~ 1,
    family="nbinom1")

    zinb2_bs<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
    BroodSize+(1|Nest),
    data=Owls,
    ziformula = ~ 1,family="nbinom2")

    # Get coefficients ("coefs" does not work yet...)
    modList = list(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
    coefs = ldply(modList,tidy,effect="fixed",conf.int=TRUE,
          .id="model") %>%
    mutate(term=abbfun(term)) %>%
    select(model,term,estimate,conf.low,conf.high) %>%
    filter(!term %in% c("Intercept","Intercept.1","NCalls","zi_NCalls"))

r broom coefficients
person Sergio Ramos    schedule 05.03.2018    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Теперь это работает (на сегодняшний день) с новым/разрабатываемым пакетом broom.mixed, например

devtools::install_github("bbolker/broom.mixed")

Надеюсь, это будет в CRAN в ближайшее время, но для меня это лишь средний приоритет, и я не хочу выпускать его в CRAN, пока оно не будет готово хотя бы на 90%. Пулл-реквесты приветствуются!

Последний шаг немного изменился (во-первых, у меня, кажется, нет abbfun):

modList = lme4:::namedList(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
coefs = ldply(modList,tidy,effect="fixed",conf.int=TRUE,
              .id="model") %>%
  select(model,term,component,estimate,conf.low,conf.high) %>%
  filter(!term %in% c("(Intercept)","NCalls"))

Предостережение: разработка этих инструментов для glmmTMB моделей является довольно новой/экспериментальной; вам следует (1) тщательно проверить результаты и (2) сообщить о проблеме если что-то не работает, как ожидалось.

person Ben Bolker    schedule 25.03.2018