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!