Search code examples
rterraglmmtmb

Create a raster with predictions from glmmTMB when a parameter is fixed


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)              

enter image description here

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)

enter image description here 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!


Solution

  • 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)  
    

    Which returns enter image description here

    Thanks @RobertHijmans