Search code examples
rcoefficientsbroom

Error message when using broom to get coefficients from glmmTMB zero inflation models


Using the Owls data and the glmmTMB package, I want to visually compare the regression coefficients from different zero-Inflated models that differ in the family used (ZIPOISS, ZINB1, ZINB2) and with/out the offset (logBroodSize).

However my first problem is to get the coefficients. The tidy function from package broom should provide you with the coefficients to plot them later with ggplot, but I get the following error when I try to get them:

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

Any idea of what could be wrong? I was already told that not having a right version of broom could be the reason, however I have had installed the right version of it... Code for a reproducible example is provided next:

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

Solution

  • This now works (as of today) with the new/under-development broom.mixed package, e.g.

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

    Hopefully this will be on CRAN sometime soon, but it's only medium priority for me, and I don't want to release it to CRAN until it's at least 90% baked. Pull requests welcome!

    The last step changed a little bit (for one thing, I don't seem to have 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"))
    

    Caveat: The development of these tools for glmmTMB models is pretty new/experimental; you should (1) sanity-check your results carefully and (2) report an issue if something doesn't work as expected.