I'm trying to create a raster showing the predictions of a model fitted with glmmTMB
. I can generally do it with no problem using the suggestion from here: https://stackoverflow.com/a/72736029
However, I am running into problems when I fix one of the parameter values, specifically the standard deviation of the random effect.
Here is an example (that is completely meaningless but explains the problem). Let's start with a simple model that works.
library(terra)
library(glmmTMB)
logo <- rast(system.file("ex/logo.tif", package="terra"))
names(logo) <- c("red", "green", "blue")
logo$blue <- as.factor(logo$blue)
logo.df <- as.data.frame(logo)
m1 <- glmmTMB(red ~ green + (1|blue),
family=poisson, data = logo.df)
# There is a warning, and clearly this is a bad model
# But it's irrelevant for this example
#predict to a raster
r1 <- predict(logo, m1, const=data.frame(blue=NA))
plot(r1)
The above works great and makes a plot that makes sense.
Now let me use the same model, but let's fix the standard deviation of the random effect to a specific value.
m2 <- glmmTMB(red ~ green + (1|blue),
family=poisson, data = logo.df, doFit=FALSE)
m2$parameters$theta[1] <- log(1e3) # Set large value
m2$mapArg <- list(theta=factor(NA)) # Tell TMB to keep it fix (don't estimate)
m2.f <- glmmTMB:::fitTMB(m2)
# There is a warning, and clearly this is a bad model
# But it's irrelevant for this example
r2 <- predict(logo, m2.f, const=data.frame(blue=NA))
plot(r2)
This creates a plot that is just one colour because the predictions are all 0s. It's a bit strange since there is an estimated relationship with green.
Not sure how to tackle the problem.
Thanks!
As indicated by @RobertHijmans this does not have to do with the raster, since the problem occurs when using predict()
on values. I have found this very helpful discussion: https://github.com/glmmTMB/glmmTMB/issues/535 which indicates that we need to use the new map
argument instead of maprArg
when fixing a value. So here is now a version that works.
library(glmmTMB)
library(terra)
logo <- rast(system.file("ex/logo.tif", package="terra"))
names(logo) <- c("red", "green", "blue")
logo$blue <- as.factor(logo$blue)
logo.df <- as.data.frame(logo)
m4 <- glmmTMB(red ~ green + (1|blue),
family=poisson, data = logo.df,
map= list(theta=factor(NA)) ,
start = list(theta=log(1e3)))
summary(m4)
predict(m4, data.frame(red=2,green=100,blue=NA))
r4 <- predict(logo$green, m4, const=data.frame(blue=NA))
plot(r4)
Thanks @RobertHijmans