updates

.Rproj.user/shared/notebooks/paths

@@ -8,6 +8,7 @@ C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/R/mixed_model.R="8C72A7F
 C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/README.md="CC30862A"
 C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/_pkgdown.yml="A5519C8D"
 C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/man/GLMMadaptive.Rd="D77150EC"
+C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/man/extra_fams.Rd="2542856A"
 C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/man/marginal_coefs.Rd="C137EA93"
 C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/man/methods.Rd="BE90FBEA"
 C:/Users/Dimitris/Documents/PackagesGitHub/GLMMadaptive/man/mixed_model.Rd="F4B70FF9"

NAMESPACE

@@ -3,7 +3,7 @@ export(mixed_model, marginal_coefs, coef, fixef, ranef, confint, anova, fitted,
        simulate, model.matrix, model.frame, terms, formula, nobs, hurdle.poisson, 
        hurdle.negative.binomial, hurdle.lognormal, beta.fam, hurdle.beta.fam, students.t,
        scoring_rules, cr_setup, cr_marg_probs, VIF, unit.lindley, cooks.distance,
-       beta.binomial, Gamma.fam, censored.normal)
+       beta.binomial, Gamma.fam, censored.normal, zi.binomial)
 
 importFrom("stats", "AIC", "BIC", "asOneSidedFormula", "cov2cor",
              "dbinom", "dnorm", "dpois", "fitted", "glm.fit", "logLik",

NEWS.md

@@ -1,3 +1,8 @@
+# GLMMadaptive 0.8.5
+
+## General
+* Added the `zi.binomial()` family object for fitting zero-inflated Binomial mixed-effects models.
+
 # GLMMadaptive 0.8.0
 
 ## General

R/Fit_Funs.R

@@ -763,6 +763,69 @@ hurdle.negative.binomial <- function () {
               class = "family")
 }
 
+zi.binomial <- function () {
+    stats <- make.link(link = "logit")
+    log_dens <- function (y, eta, mu_fun, phis, eta_zi) {
+        # the log density function
+        # Binomial part
+        mu <- mu_fun(eta)
+        y <- as.matrix(y)
+        N <- if (ncol(y) == 2L) y[, 1L] + y[, 2L] else rep(1L, nrow(y))
+        out <- as.matrix(dbinom(y[, 1L], N, mu, TRUE))
+        # ZI part
+        ind_y0 <- y[, 1L] == 0
+        ind_y1 <- y[, 1L] > 0
+        pis <- as.matrix(plogis(eta_zi))
+        # combined
+        out[ind_y0, ] <- log(pis[ind_y0, ] + (1 - pis[ind_y0, ]) * exp(out[ind_y0, ]))
+        out[ind_y1, ] <- log(1 - pis[ind_y1, ]) + out[ind_y1, ]
+        attr(out, "mu_y") <- mu
+        out
+    }
+    score_eta_fun <- function (y, mu, phis, eta_zi) {
+        # Binomial part
+        mu <- as.matrix(mu)
+        y <- as.matrix(y)
+        N <- if (ncol(y) == 2L) y[, 1L] + y[, 2L] else rep(1L, nrow(y))
+        out <- y[, 1L] * (1 - mu) - (N - y[, 1L]) * mu
+        # ZI part
+        ind_y0 <- y[, 1L] == 0
+        eta_zi <- as.matrix(eta_zi)
+        pis <- plogis(eta_zi[ind_y0, ])
+        mu0 <- mu[ind_y0, ]
+        N0 <- N[ind_y0]
+        pis1 <- 1 - pis
+        den <- pis + pis1 * (1 - mu0)^N0
+        out[ind_y0, ] <- - (N0 * mu0 * pis1 * (1 - mu0)^N0) / den
+        out
+    }
+    score_eta_zi_fun <- function (y, mu, phis, eta_zi) {
+        y <- as.matrix(y)
+        N <- if (ncol(y) == 2L) y[, 1L] + y[, 2L] else rep(1L, nrow(y))
+        ind_y0 <- y[, 1L] == 0
+        ind_y1 <- y[, 1L] > 0
+        pis <- as.matrix(plogis(eta_zi))
+        mu <- as.matrix(mu)
+        # Binomial part
+        out <- mu
+        out[ind_y1, ] <- - pis[ind_y1, ]
+        # ZI part
+        mu0 <- mu[ind_y0, ]
+        N0 <- N[ind_y0]
+        pis1 <- 1 - pis[ind_y0, ]
+        FF <- (1 - mu0)^N0
+        den <- pis[ind_y0, ] + pis1 * FF
+        out[ind_y0, ] <- pis[ind_y0, ] * pis1 * (1 - FF) / den
+        out
+    }
+    structure(list(family = "zero-inflated binomial", link = stats$name, 
+                   linkfun = stats$linkfun, linkinv = stats$linkinv, log_dens = log_dens,
+                   score_eta_fun = score_eta_fun,
+                   score_eta_zi_fun = score_eta_zi_fun,
+                   score_phis_fun = NULL),
+              class = "family")
+}
+
 hurdle.lognormal <- function () {
     stats <- make.link("identity")
     log_dens <- function (y, eta, mu_fun, phis, eta_zi) {

R/methods.R

@@ -1325,6 +1325,19 @@ simulate.MixMod <- function (object, nsim = 1, seed = NULL,
                 rbinom(n = n, size = .N, prob = mu)
             }
             environment(sim_fun) <- env
+        } else if (object$family$family == "zero-inflated binomial") {
+            N <- if ((NCOL(y <- model.response(object$model_frames$mfX))) == 2) 
+                y[, 1] + y[, 2] else 1
+            .N <- N
+            env <- new.env(parent = .GlobalEnv)
+            assign(".N", N, envir = env)
+            sim_fun <- function (n, mu, phis, eta_zi) {
+                out <- rbinom(n = n, size = .N, prob = mu)
+                extra_zeros <- as.logical(rbinom(n, 1, plogis(eta_zi)))
+                out[extra_zeros] <- 0
+                out
+            }
+            environment(sim_fun) <- env
         } else if (object$family$family == "poisson") {
             sim_fun <- function (n, mu, phis, eta_zi) {
                 rpois(n = n, lambda = mu)

R/mixed_model.R

@@ -19,7 +19,8 @@ mixed_model <- function (fixed, random, data, family, weights = NULL,
              "  use 'family = GLMMadaptive::negative.binomial()'.")
     }
     if (family$family %in% c("zero-inflated poisson", "zero-inflated negative binomial",
-                             "hurdle poisson", "hurdle negative binomial", "hurdle beta") && 
+                             "hurdle poisson", "hurdle negative binomial", "hurdle beta", 
+                             "zero-inflated binomial") && 
         is.null(zi_fixed)) {
         stop("you have defined a family with an extra zero-part;\nat least argument ",
              "'zi_fixed' needs to be defined, and potentially also argument 'zi_random'.")
@@ -169,6 +170,8 @@ mixed_model <- function (fixed, random, data, family, weights = NULL,
                 glm.fit(X, y, family = Gamma(), offset = offset)$coefficients
             else if (family$family == "censored normal")
                 glm.fit(X, y, family = gaussian(), offset = offset)$coefficients
+            else if (family$family == "zero-inflated binomial")
+                glm.fit(X, y, family = binomial(), offset = offset)$coefficients
             else
                 glm.fit(X, y, family = family, offset = offset)$coefficients
         } else {

README.md

@@ -39,6 +39,8 @@ object.
 - Hurdle Poisson and negative binomial models using the `hurdle.poisson()` and 
 `hurdle.negative.binomial()` family objects.
 
+- Zero-inflated binomial models using the `zi.binomial()` family objects.
+
 - Two-part/hurdle mixed models for semi-continuous normal data using the 
 `hurdle.lognormal()` family object.
 

man/extra_fams.Rd

@@ -3,6 +3,7 @@
 \alias{beta.fam}
 \alias{beta.binomial}
 \alias{zi.poisson}
+\alias{zi.binomial}
 \alias{zi.negative.binomial}
 \alias{hurdle.poisson}
 \alias{hurdle.negative.binomial}
@@ -13,6 +14,7 @@
 \alias{censored.normal}
 
 
+
 \title{
 Family functions for Student's-t, Beta, Zero-Inflated and Hurdle Poisson and Negative Binomial,
 Hurdle Log-Normal, Hurdle Beta, Gamma, and Censored Normal Mixed Models}
@@ -28,6 +30,7 @@ mixed-effects model, using
 students.t(df = stop("'df' must be specified"), link = "identity")
 beta.fam()
 zi.poisson()
+zi.binomial()
 zi.negative.binomial()
 hurdle.poisson()
 hurdle.negative.binomial()