1
$\begingroup$

I've got a dataset that has temperature (21c or 29c), inoculation (mock (m),single inoculations (c or r), or coinoculation (rc)), and age group (y or o). I am trying to model the interactive effects of temperature, age group, and inoculation on the area under the disease progression stairs (audps) using the glmmTMB package in R as follows:

library(glmmTMB)
glm_model3 <- glmmTMB(audps ~ temp * age_gp * inoc+(1|replicate),
    family = t_family, data = audps_data2)
glm_model3

I had gotten the model to run previously, but then added in the final replicate's data (with more zeroes) today, it produced the following error:

Warning in fitTMB(TMBStruc) : Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
Warning in fitTMB(TMBStruc) : Model convergence problem; function evaluation limit reached without convergence (9). See vignette('troubleshooting')

histogram of audps

I've tried to change the optimizer and fiddled with using a zero-inflated family, but no luck in getting the Hessian issue resolved/getting the model to converge. I've also attempted to run the model using a Gamma family, but it tells me that non-positive values aren't allowed (which I have a hunch is because of the zeroes). I have attached a picture of the histogram of audps, for reference. Any help would be much appreciated!

$\endgroup$
1
  • 1
    $\begingroup$ Could you please substitute an image that includes only the histogram of audps (rather than a screenshot of your entire screen)? $\endgroup$
    – Ben Bolker
    Commented Jul 10 at 15:00

1 Answer 1

2
$\begingroup$

Hard to say, exactly. If your response variable is positive with exact zeros then indeed using a zero-inflated Gamma or a Tweedie distribution as the response would make more sense and might work better.

If you specify zero-inflation explicitly, family = "Gamma" shouldn't complain about zero values in your response, e.g.

glmmTMB(..., family = "Gamma", ziformula = ~1)

(or you could use a more complicated zero-inflation model if you wanted)

family = "Tweedie" should also handle zero-inflated positive responses; it is less computationally efficient and makes different assumptions about how the zero-response component is related to the rest of the responses.


Your solution in comments (fixing the df for the Student-t to 887, i.e. fixing the psi parameter to log(887)) seems odd for a few reasons: (1) it is a strangely specific value [how was it chosen?] and (2) if you set the df this large, you are for all practical purposes fitting a Gaussian.

One possible reason for a non-positive-definite fit would be if glmmTMB were really trying to make your t-distribution into a Gaussian, in which case the df parameter would be going to infinity and the log-likelihood surface would be getting flat ( == non-positive definite Hessian). You could check family_params(glm_model3) to evaluate this possibility ... if it's very large (e.g. > 30) you might consider switching to Gaussian. (I wasn't able to replicate that with a simple example, although I did get lots of "false convergence" warnings ...)

library(glmmTMB)
r <- numeric(20)
for (i in 1:20) {
   set.seed(100+i)
   dd <- data.frame(x = rnorm(10000))
   dd$y <- simulate_new(~ 1 + x,
                  family = gaussian,
                  newdata = dd,
                  newparams = list(beta = c(1, 1), betadisp = 0))[[1]]
   m <- glmmTMB(y ~ 1 + x, family = t_family, data = dd)
   r[i] <- family_params(m)
}
summary(r)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
      90      807  1517740  5074272  5671397 46921288 
$\endgroup$
1
  • $\begingroup$ Ok, so I tried the zero-inflation specification and using the tweedie family and both violated the model assumptions. Somebody in my lab suggested that I try the following edits to fix the shape parameter, and it seems to be up and running! I used a regular glm to get the degrees of freedom and added 1 for my random effect. Does this seem like a reasonable solution? glm_model3<-glmmTMB(audps ~ temp * age_gp * inoc+(1|replicate), family = t_family(),start=list(psi=log(887)),map=list(psi=factor(NA)),data = audps_data2) $\endgroup$
    – user26322319
    Commented Jul 10 at 15:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.