Support for INLA

[florianhartig]

via Email:

In my previous work, I enjoy using the DHARMa package to test spatial and temporal residual autocorrelation and need to run similar tests on the current sp model. However, the model runs using INLA which is an unsupported package with DHARMa.

I found this example to help with code development : https://aosmith.rbind.io/2017/12/21/using-dharma-for-residual-checks-of-unsupported-models/#just-the-code-please but without success.

I wonder if you have some codes that could help us or if you would be keen to support the development. The idea is to run a basic INLA model without spatial and/or temporal components and perform tests to inform the need of these components into the modelling framework based on the detection of spatial and/or temporal autocorrelation. I am also interested to use them to ensure absences of such autocorrelation into the residuals from the final model. The current model uses a binomial probability distribution with some covariates. There is also a second model using beta distribution that works very well and it’s a bit more advanced in the development. I am happy to share model outputs from the second model to help the development.

[florianhartig]

Hello,

I'm not a frequent INLA user, but from a superficial look, I don't see an easy way to couple INLA to DHARMa without major coding efforts in the DHARMa wrappers.

However, there is a function to create posterior predictive simulations https://rdrr.io/github/inbo/INLA/man/posterior.sample.html, so it should be possible to create those and then read them into DHARMa via createDHARMA().

I have made an attempt to do so (see below), but the result doesn't seem to be correct. In particular I encountered the following problems:

1) Predictions ("Predictor") seems to be after the link function, but are not strictly positive 2) I didn't find a way to simulate the response directly, thus the rpois 3) Not quite sure how the RE (if introduced) are handled, i.e. if you simulate starting from hyper parameters, or if you condition on the latents 4) Most importantly, however, I didn't manage to get a good residual pattern, which I would have expected in this example (data is created as assumed in the model).

However, this may all be due to my imperfect understanding of INLA. I would suggest you have to contact the INLA developers (in particular @hrue) to check how to create correct posterior predictive simulations, then it should work with DHARMa.

library(INLA)
library(DHARMa)

dat = createData(n = 500, randomEffectVariance = 0)
str(dat) # 100 obs

fit = inla(observedResponse ~ 1 + Environment1 , family = "poisson", data = dat, 
           control.family=list(link='log'),
           control.predictor=list(link=1, compute=TRUE),
           control.compute=list(config = TRUE))

summary(fit)

sim <- inla.posterior.sample(50, m)
sim = inla.posterior.sample.eval("Predictor", sim)
poi <- function(x){
  x = ifelse(x<0, 0, x)
  rpois(1,x)
} 
sim = apply(sim, c(1,2), poi)

res = createDHARMa(sim, observedResponse = dat$observedResponse)
plot(res)

[JulieVercelloni]

Thanks a lot for your answer @florianhartig. We will check how to create the correct posterior predictive simulations and keep you posted if we have positive results.

@hrue, if possible, your guidance on the best way to create the posterior predictive simulations will be very appreciated here.

Cheers

[florianhartig]

If you do so, make sure to understand how you can control on what latents you want to condition on. See these comments in the DHARMa Vignette for Bayesian https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMaForBayesians.html#conditional-vs.-unconditional-simulations-in-hierarchical-models

[tmeeha]

If anyone else is interested in this, here is some commented code, modified from Florian, above, that uses DHARMa to evaluate the fit of a simple inla ZIP model.

Make data

library(INLA)
library(DHARMa)
n <- 100
p <- 0.75
lam <- 2
y <- rpois(n, lambda = lam)
y[rbinom(n, size=1, prob=p) == 1] = 0
dat1 <- data.frame(y=y)

Make model For the family used in this example, the right side of the formula relates to the count component of the model. Covariates of p cannot be specified with this family. There are ways to do it with inla, but they are more complicated.

mod1 <- inla(y ~ 1, 
             family = "zeroinflatedpoisson1", 
             data = dat1,
             control.predictor=list(compute=TRUE),
             control.compute=list(config=TRUE))
summary(mod1)

Get posterior samples Get 50 posterior samples for p and lambda. Note that lambda needs to be exponentiated before simulation.

ss1 <- 50
samps1 <- inla.posterior.sample(ss1, mod1)
get_p <- function(...){return(theta[1])}
sim_p <- inla.posterior.sample.eval(get_p, samps1)
sim_lam <- exp(inla.posterior.sample.eval("Predictor", samps1))

Simulate counts Use p and lambda samples to simulate counts

out_sim <- sim_lam
for(i in 1:nrow(sim_lam)){
  for(j in 1:ncol(sim_lam)){
    p <- sim_p[,j]
    p <- 1 - p
    y <- rpois(1, sim_lam[i,j])
    y <- y * rbinom(n=1, size=1, prob=p)
    out_sim[i, j] <- y
  }
}

Run diagnostics using DHARMa

res <- createDHARMa(simulatedResponse=out_sim, 
                    observedResponse=dat1$y,
                    integerResponse = T)
plot(res)