General linear mixed model in R which will fit quasi family I am trying to run a GLMM with a quasibinomial family (my data is 0 inflated and I have a negative min x value), but am receiving this error message as quasi families cannot be used in glmer:
Error in lme4::glFormula(formula = slope_prop ~ temp * gly * day + (1 |  : 
  "quasi" families cannot be used in glmer

Is there any way of running a GLMM which will use quasi families? Or am I doing the wrong thing by trying to use a quasibonomial family in the first place?

 A: You could instead use a zero-inflated binomial mixed model. This can be fitted using the GLMMadaptive package and by defining a user-specified family object; see here for more information.
For an example, check this:
# install the development version from GitHub
devtools::install_github("drizopoulos/GLMMadaptive")
# Load the package
library("GLMMadaptive")

# Simulate some zero-inflated Binomial data
set.seed(1234)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 15 # maximum follow-up time

# we constuct a data frame with the design: 
# everyone has a baseline measurment, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
                 time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
                 sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
# design matrices for the fixed and random effects zero-part
X_zi <- model.matrix(~ sex, data = DF)
Z_zi <- model.matrix(~ 1, data = DF)

betas <- c(-2.13, -0.25, 0.24, -0.05) # fixed effects coefficients
gammas <- c(-1.5, 0.5) # fixed-effects zero-part
D11 <- 0.48 # variance of random intercepts
D22 <- 0.1 # variance of random slopes
D_zi <- 0.4 # variance of random intercepts zero-part

# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)), rnorm(n, sd = sqrt(D_zi)))
# linear predictor
eta_y <- drop(X %*% betas + rowSums(Z * b[DF$id, 1:2]))
# we simulate binomial longitudinal data with 10 trials
DF$y <- rbinom(n * K, 10, plogis(eta_y))
# linear predictor zero-part
eta.zi <- drop(X_zi %*% gammas + rowSums(Z_zi * b[DF$id, 3]))
# set the extra zeros
DF$y[as.logical(rbinom(n * K, 1, plogis(eta.zi)))] <- 0


# Define the zero-inflated family object
zi.binomial <- function () {
    stats <- make.link("logit")
    log_dens <- function (y, eta, mu_fun, phis, eta_zi) {
        # the log density function
        N <- y[, 1] + y[, 2]
        y <- y[, 1]
        # Binomial part
        mu <- mu_fun(eta)
        out <- as.matrix(dbinom(y, N, mu, log = TRUE))
        # ZI part
        ind_y0 <- y == 0
        ind_y1 <- y > 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
    }
    structure(list(family = "zero inflated binomial", link = stats$name, linkfun = stats$linkfun,
                   linkinv = stats$linkinv, log_dens = log_dens), class = "family")
}


# Fit the model
fm <- mixed_model(cbind(y, 10 - y) ~ sex * time, random = ~ 1 | id, data = DF,
                  family = zi.binomial(), zi_fixed = ~ sex, 
                  initial_values = list(betas = binomial()))

summary(fm)

