From 48ebd1798a784a359aff9d6b3be8bd05329fdc50 Mon Sep 17 00:00:00 2001 From: Dimitris Rizopoulos Date: Mon, 7 Feb 2022 15:42:03 +0100 Subject: [PATCH] updates --- .Rproj.user/shared/notebooks/paths | 1 + NAMESPACE | 2 +- NEWS.md | 5 +++ R/Fit_Funs.R | 63 ++++++++++++++++++++++++++++++ R/methods.R | 13 ++++++ R/mixed_model.R | 5 ++- README.md | 2 + man/extra_fams.Rd | 3 ++ 8 files changed, 92 insertions(+), 2 deletions(-) diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths index e057c89..0d6b9cb 100644 --- a/.Rproj.user/shared/notebooks/paths +++ b/.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" diff --git a/NAMESPACE b/NAMESPACE index 293e05e..f3a174c 100644 --- a/NAMESPACE +++ b/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", diff --git a/NEWS.md b/NEWS.md index 7d21b23..0501695 100644 --- a/NEWS.md +++ b/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 diff --git a/R/Fit_Funs.R b/R/Fit_Funs.R index eec03a8..ef6337c 100644 --- a/R/Fit_Funs.R +++ b/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) { diff --git a/R/methods.R b/R/methods.R index 1ebcc4d..b8d0668 100644 --- a/R/methods.R +++ b/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) diff --git a/R/mixed_model.R b/R/mixed_model.R index fa91c06..e6fd7a8 100644 --- a/R/mixed_model.R +++ b/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 { diff --git a/README.md b/README.md index 74ce61e..45529bf 100644 --- a/README.md +++ b/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. diff --git a/man/extra_fams.Rd b/man/extra_fams.Rd index 453175b..7085f7b 100644 --- a/man/extra_fams.Rd +++ b/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()