strengejacke/ggeffects

Problem with ggpredict when fixing a parameter value in glmmTMB

Closed this issue · 2 comments

Hello,

I would like to use ggpredict to make effects plots for a model fitted via glmmTMB. The problem I am running into is that I am fixing one of the parameters using the argument map in glmmTMB and this (I think) creates a mismatch in dimensions that causes problems with ggpredict. Below is a set of examples.

First, let's explore a model that will work

library(glmmTMB)
library(ggeffects)

data("Salamanders")

m1 <- glmmTMB(
  count ~ cover  + (1 | site),
  family = poisson,
  data = Salamanders)
dat1 <- ggpredict(m1, c("cover", "mined"))
plot(dat1)

All is good with the above.

Second, let's look at the same model, but where we fix the standard deviation of the random effect to a set value.

# Set the model, but do not fit it.
m2 <- glmmTMB(
  count ~ cover + (1 | site),
  family = poisson,
  data = Salamanders, doFit = FALSE) 
# Specify the value of the parameter we want to fix.
m2$parameters$theta[1] <- log(1)
# Tell glmmTMB to keep it fix (don't estimate)
m2$mapArg <- list(theta=factor(NA)) 
# Now fit the model
m2.f <- glmmTMB:::fitTMB(m2)

# Try ggpredict on the model
dat2 <- ggpredict(m2.f, "cover")

This will return:

Error in EvalADFunObject(ADFun, theta, order = order, hessiancols = cols,  : 
  Wrong parameter length.

This is due to a mismatch in dimensions. I found it helpful to compare the results of the first and second models.

m1$fit$parfull
m2.f$fit$parfull

The fix parameter value does not show up in the parfull (or par), which makes sense since it's not estimated, but I think may be the source of the problem.

I tried to use the condition argument, to see if that would help.

dat2 <- ggpredict(m2.fm, "cover", condition = c(theta = m2$parameters$theta[1]))

But I get the following error:

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

I note that it works with the package effects, but I would like to use ggeffects.

library(effects)
plot(Effect("cover", m1))
plot(Effect("cover", m2.f))

I'm sorry if I am missing something obvious!

Many thanks!

Thanks to this great discussion (glmmTMB/glmmTMB#535), I found a solution.

m3 <- glmmTMB(
  count ~ cover + (1 | site),
  family = poisson,
  data = Salamanders, start = list(theta = log(1)), map = list(theta = factor(NA)))



dat3 <- ggpredict(m3, "cover")
plot(dat3, add.data=TRUE)

Great you found a working solution!