[R-sig-ME] Error from glmmTMB().

[D. Rizopoulos]

The model seems to fit with the GLMMadaptive package, i.e.,

library("GLMMadaptive")

fm <- mixed_model(cbind(Dead, Alive) ~ (Trt + 0) / Dose, data = X,
                  random = ~ Dose | Rep, 
                  family = beta.binomial(link = "cloglog"), n_phis = 1)

summary(fm)

where

beta.binomial <- function (link = "logit") {
    .link <- link
    env <- new.env(parent = .GlobalEnv)
    assign(".link", link, envir = env)
    stats <- make.link(link)
    dbbinom <- function (x, size, prob, phi, log = FALSE) {
        A <- phi * prob
        B <- phi * (1 - prob)
        log_numerator <- lbeta(x + A, size - x + B)
        log_denominator <- lbeta(A, B)
        fact <- lchoose(size, x)
        if (log) {
            fact + log_numerator - log_denominator
        } else {
            exp(fact + log_numerator - log_denominator)
        }
    }
    log_dens <- function (y, eta, mu_fun, phis, eta_zi) {
        phi <- exp(phis)
        eta <- as.matrix(eta)
        mu_y <- mu_fun(eta)
        out <- if (NCOL(y) == 2L) {
            dbbinom(y[, 1L], y[, 1L] + y[, 2L], mu_y, phi, TRUE)
        } else {
            dbbinom(y, rep(1L, length(y)), mu_y, phi, TRUE)
        }
        attr(out, "mu_y") <- mu_y
        out
    }
    score_eta_fun <- function (y, mu, phis, eta_zi) {
        phi <- exp(phis)
        mu <- as.matrix(mu)
        if (NCOL(y) == 2L) {
            size <- y[, 1L] + y[, 2L]
            y <- y[, 1L]
        } else {
            size <- rep(1L, length(y))
        }
        phi_mu <- phi * mu
        phi_1mu <- phi * (1 - mu)
        comp1 <- (digamma(y + phi_mu) - digamma(size - y + phi_1mu)) * phi
        comp2 <- (digamma(phi_mu) - digamma(phi_1mu)) * phi
        mu.eta <- switch(.link,
                         "logit" = mu - mu * mu,
                         "cloglog" = - (1 - mu) * log(1 - mu))
        out <- (comp1 - comp2) * mu.eta
        out
    }
    score_phis_fun <- function (y, mu, phis, eta_zi) {
        phi <- exp(phis)
        mu <- as.matrix(mu)
        if (NCOL(y) == 2L) {
            size <- y[, 1L] + y[, 2L]
            y <- y[, 1L]
        } else {
            size <- rep(1L, length(y))
        }
        mu1 <- 1 - mu
        phi_mu <- phi * mu
        phi_1mu <- phi * mu1
        comp1 <- digamma(y + phi_mu) * mu + digamma(size - y + phi_1mu) * mu1 - 
            digamma(size + phi)
        comp2 <- digamma(phi_mu) * mu + digamma(phi_1mu) * mu1 - digamma(phi)
        out <- (comp1 - comp2) * phi
        out
    }
    structure(list(family = "beta binomial", link = stats$name, 
                   linkfun = stats$linkfun, linkinv = stats$linkinv, log_dens = log_dens,
                   score_eta_fun = score_eta_fun,
                   score_phis_fun = score_phis_fun),
              class = "family")
}

Best,
Dimitris


-----Original Message-----
From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> On Behalf Of Rolf Turner
Sent: Sunday, March 15, 2020 4:32 AM
To: r-sig-mixed-models using r-project.org
Subject: [R-sig-ME] Error from glmmTMB().


I am getting an error, that I have no idea what to do about, from glmmTMB():

library(glmmTMB)
fmla <- cbind(Dead, Alive) ~ (Trt + 0)/Dose + (Dose | Rep)
X    <- dget("X.txt")
fit  <- glmmTMB(fmla, data = X, family = betabinomial(link = "cloglog"),
                dispformula = ~1)

> Error in optimHess(par.fixed, obj$fn, obj$gr) : 
>   gradient in optim evaluated to length 1 not 16 In addition: There 
> were 16 warnings (use warnings() to see them)

The warnings are all repetitions of

> 1: In nlminb(start = par, objective = fn, gradient = gr,  ... :
>   NA/NaN function evaluation

The error sounds to me like something is amiss in the code.

Can anyone confirm/deny/suggest what I might do to get this call to
glmmTMB() to run?

Thanks.

cheers,

Rolf Turner

P.S.  The data set is attached.

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276