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

[Rolf Turner]

Wow.  This is impressive, and fascinating, albeit somewhat mystifying.

I tried your code and yes indeed it runs without complaint; no errors, 
no warnings.

I conclude from this that mixed_model() is capable of fitting models 
using the beta binomial distribution, a fact of which I was previously 
unaware.  (I previously thought that only the glmmTMB package provided 
this capability.)

I gather, after reading TFM a bit more, that mixed_model() can 
accommodate custom-built family functions, such as the beta.binomial()
function that you used to produce a fit to my example.

However writing such custom-built family functions is a bit beyond my
very limited capabilities.

I am currently struggling to find suitable models to fit to a number of
rather recalcitrant data sets.  The model presented in my question to
r-sig-mixed-models was just one of a fairly large number with which I 
would like to experiment.  The data in question appear to exhibit 
over-dispersion w.r.t. the binomial distribution.  The random effects 
models that I am fitting seem to produce problematic results; I want to 
try using the beta binomial distribution in addition to (or perhaps 
instead of) random effects modelling.

With glmmTMB() I can (with *some* data sets) fit models like

     glmmTMB(cbind(Dead, Alive) ~ (Trt + 0)/Dose + (Dose | Rep),
             data = X, family = betabinomial(link = "cloglog"),
                dispformula = ~Dose)

I can also set dispformula equal to ~poly(Dose,2) or ~splines::ns(Dose,2).

Can custom-built family functions be created to do this sort of thing 
using mixed_model()?

If so, what would the appropriate value of n_phis be?  I am *guessing* 
that for dispformula = ~Dose, n_phis would be set equal to 2, and for
dispformula = ~poly(Dose,2) it would be 3.  What about for
dispformula = ~splines::ns(Dose,2) ?

Also (sorry for going on!) with glmmTMB() there is the possibility of
omitting random effects entirely and just trying to accommodate 
over-dispersion by means of the beta binomial distribution:

glmmTMB(cbind(Dead, Alive) ~ (Trt + 0)/Dose,data = X,
        family = betabinomial(link = "cloglog"),
                dispformula = ~Dose)

Omitting random effects does not seem to be possible with mixed_model().
Is that correct?  Or have I simply not been using the right syntax?

Thanks for any words of wisdom that you can spare me.

cheers,

Rolf

On 19/03/20 10:03 am, D. Rizopoulos wrote:

> 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?