tests/testthat/test-zi.R

In glmmTMB: Generalized Linear Mixed Models using Template Model Builder

stopifnot(require("testthat"),
          require("glmmTMB"),
          require("lme4"))

## simulate something smaller than the Owls data set?


context("ZI models")
data(Owls)

test_that("zi", {
  skip_on_cran()
    ## Fit negative binomial model with "constant" Zero Inflation :
    owls_nb1 <<- glmmTMB(SiblingNegotiation ~ FoodTreatment*SexParent +
                             (1|Nest)+offset(log(BroodSize)),
                         family = nbinom1(),
                         ziformula = ~1, data=Owls)
    owls_nb2 <<- update(owls_nb1,
                        ziformula = ~ FoodTreatment*SexParent + (1|Nest))
    owls_nb3 <<- update(owls_nb1,ziformula=~.)

    expect_equal(fixef(owls_nb2),
  structure(list(cond = structure(c(0.812028613585629, -0.342496105044418,
      -0.0751681324132088, 0.122484981295054),
      .Names = c("(Intercept)", "FoodTreatmentSatiated", "SexParentMale",
                 "FoodTreatmentSatiated:SexParentMale")),
      zi = structure(c(-2.20863281353936, 1.86779027553285, -0.825200772653965,
    0.451435813933449), .Names = c("(Intercept)", "FoodTreatmentSatiated",
     "SexParentMale", "FoodTreatmentSatiated:SexParentMale")),
    disp = structure(1.33089630005212, .Names = "(Intercept)")),
    .Names = c("cond", "zi", "disp"), class = "fixef.glmmTMB"),
  tolerance=1e-5)
    expect_equal(fixef(owls_nb2),fixef(owls_nb3))
})

test_that("zi beta and Gamma", {
    skip_on_cran()
    suppressWarnings(RNGversion("3.5.1"))
    set.seed(101)
    dd <- data.frame(yb=c(rbeta(100,shape1=2,shape2=1),rep(0,10)),
                     yg=c(rgamma(100,shape=1.5,rate=1),rep(0,10)))
    expect_error(glmmTMB(yb~1, data=dd, family=beta_family),
                 "y values must be")
    m1 <- glmmTMB(yb~1, data=dd, family=beta_family, zi=~1)
    expect_equal(unname(plogis(fixef(m1)[["zi"]])),1/11)
    expect_equal(unname(fixef(m1)[["cond"]]), 0.6211636, tolerance=1e-5)
    ## need *both* ziformula and family=ziGamma for gamma-hurdle
    expect_error(glmmTMB(yg~1, data=dd, family=Gamma),
                 "non-positive values not allowed")
    expect_error(glmmTMB(yg~1, zi=~1, data=dd, family=Gamma),
                 "non-positive values not allowed")
    expect_error(glmmTMB(yg~1, data=dd, family=ziGamma),
                 "non-positive values not allowed")
    m2 <- glmmTMB(yg~1, data=dd, family=ziGamma(link="log"), zi=~1)
    expect_equal(unname(plogis(fixef(m2)[["zi"]])),1/11)
    expect_equal(unname(fixef(m2)[["cond"]]), 0.3995267, tolerance=1e-5)
})