[R-sig-ME] Error from glmmTMB().
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Thu Mar 19 03:52:57 CET 2020
Hijacking the stream back slightly in the direction of glmmTMB (not
intentional, just lazy): I made a small change that appears to allow the
original model to run correctly. If you want to try it out, and have
development tools (compilers etc.) and the remotes package installed,
remotes::install_github("glmmTMB/glmmTMB/glmmTMB using cloglog_fix")
should install the patched version (I think).
This will presumably be merged with the master branch sometime soon but
I want to add some tests etc.
On 2020-03-18 10:15 p.m., Rolf Turner wrote:
>
> 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?
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list