INLA and inlabru providing different results

[Filhomn]

I’m trying to replicate my inla code with inlabru for a spde model. However, I’m not being successful.

• INLA code

group = as.matrix(as.numeric(as.factor(data[,6])))

A_point <- inla.spde.make.A(mesh, loc = coords, group = group)

spde1 <- inla.spde2.matern(mesh, alpha = 2)

s.index <- inla.spde.make.index("spatial.field", n.spde = spde1$n.spde, n.group = 9)

stack <- inla.stack(data = list(y = data$tot), A = list(A_point, 1), effects = list(c(s.index,list(Intercept = rep(1,spde1$n.spde))), list(alag=data$alag)), tag="obs")

data.inla <- inla(y ~ -1 + Intercept + alag +f(spatial.field, model = spde1, group= spatial.field.group, control.group=list(model="iid")) , family='Gaussian', data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack)), control.fixed = list(prec.intercept = 0.0001), control.compute = list(cpo = T, waic=T), verbose=FALSE)

• Inlabru code

data4= SpatialPointsDataFrame( coords = coords, data = data, proj4string = CRS("+proj=longlat +datum=WGS84") )

data4$Ano = as.numeric(as.factor(data$Ano)) comp_inlabru <- as.formula(paste0('tot ~ -1 + se(1) + alag + s(main = coordinates, model = spde1, group=Ano, ngroup = 9, control.group=list(model="iid"))'))

inlabru <- bru(comp_inlabru, data = data4, family = "Gaussian"

I think the problem is probably the intercept. Any ideas of how I could replicate this code using inlabru?

Thank you very much,

[finnlindgren]

Hi, we need more information about what the problem you're experiencing is, e.g. if there are error messages, or if it's jut the result that has some surprising properties, etc.

[Filhomn]

Hi Finn, Thank you very much for the quick reply. I’m new to inlabru and I’m trying to translate this INLA code to inlabru, which is a spatial-temporal model.

group = as.matrix(as.numeric(as.factor(data[,6])))

A_point <- inla.spde.make.A(mesh, loc = coords, group = group)

spde1 <- inla.spde2.matern(mesh, alpha = 2)

s.index <- inla.spde.make.index("spatial.field", n.spde = spde1$n.spde, n.group = 9)

**stack <- inla.stack(data = list(y = data$tot), A = list(A_point, 1), effects = list(c(s.index,list(Intercept = rep(1,spde1$n.spde))), list(alag=data$alag)), tag="obs")

data.inla <- inla(y ~ -1 + Intercept + alag +f(spatial.field, model = spde1, group= spatial.field.group, control.group=list(model="iid")), family='Gaussian', data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack)), control.fixed = list(prec.intercept = 0.0001), control.compute = list(cpo = T, waic=T), verbose=FALSE)**

• Here is how far I went…. But it just not there yet…

data4= SpatialPointsDataFrame( coords = coords, data = data, proj4string = CRS("+proj=longlat +datum=WGS84") )

data4$Ano = as.numeric(as.factor(data$Ano)) comp_inlabru <- as.formula(paste0('tot ~ Intercept(1) + alag + s(main = coordinates, model = spde1, group=Ano, ngroup = 9, control.group=list(model="iid"))'))

inlabru <- bru(comp_inlabru, data = data4, family = "Gaussian")

Thank you very much!! Kind Regards,

[finnlindgren]

I'd just write the component specification as

comp_inlabru <- tot ~ Intercept(1) + alag + s(main = coordinates, model = spde1, group=Ano, ngroup = 9, control.group=list(model="iid"))

but I think your code should have the same effect as that.

However, that still doesn't explain what the problem you're experiencing is. Are you asking which parts of the "**"-marked inla code you need to convert?

[finnlindgren]

Most of the inla code isn't needed. From what I can see, all you need to do is this:

comp_inlabru <- tot ~ Intercept(1, prec.linear = 0.0001) + alag + s(main = coordinates,
  model = spde1, group=Ano, ngroup = 9, control.group=list(model="iid"))
inlabru <- bru(comp_inlabru, data = data4, family = "Gaussian",
   options = list(control.compute = list(cpo = TRUE, waic=TRUE)))

inlabru handles all the A-matrices and inla.stack calls internally.

[Filhomn]

I have created a toy example to show the problem. The intercept from the models are different. And this difference is very expressive in my real life dataset (not the example below). Am I calculating the intercepts differently?

library(INLA) library(inlabru)

generate dataset

data("shrimp") coords = shrimp$hauls@coords

shrimp$hauls$group= c(rep('a',43), rep('b',43))

mesh <- inla.mesh.2d(loc=coords,max.n.strict = c(128000,128000), max.edge = c(0.05,0.5),max.n = c(48000,16000), min.angle = c(30,21), cutoff = 0.01)

INLA CODE

group = as.matrix(as.numeric(as.factor(shrimp$hauls$group)))

A_point <- inla.spde.make.A(mesh, loc = shrimp$hauls@coords, group = group)

spde1 <- inla.spde2.matern(mesh, alpha = 2)

s.index <- inla.spde.make.index("spatial.field", n.spde = spde1$n.spde, n.group = 2)

stack <- inla.stack(data = list(y = shrimp$hauls$catch), A = list(A_point,1), effects = list(c(s.index,list(Intercept = rep(1,spde1$n.spde))), list(depth=shrimp$hauls$depth)), tag="obs")

data.inla <- inla(y ~ -1 + Intercept + depth +f(spatial.field, model = spde1,group= spatial.field.group, control.group=list(model="iid")) , family='Gaussian', data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack)), control.fixed = list(prec.intercept = 0.0001, prec=0.0001), control.compute = list(cpo = T, waic=T), verbose=FALSE) summary(data.inla)

Call c("inla(formula = y ~ -1 + Intercept + depth + f(spatial.field, ", " model = spde1, group = spatial.field.group, control.group = list(model = \"iid\")), ", " family = \"Gaussian\", data = inla.stack.data(stack), verbose = FALSE, ", " control.compute = list(cpo = T, waic = T), control.predictor = list(A = inla.stack.A(stack)), ", " control.fixed = list(prec.intercept = 1e-04, prec = 1e-04))" ) Time used: Pre = 1.39, Running = 17.4, Post = 0.219, Total = 19 Fixed effects: mean sd 0.025quant 0.5quant 0.975quant mode kld Intercept -6.376 18.619 -47.136 -5.638 29.282 -4.652 0 depth 0.067 0.022 0.024 0.067 0.109 0.067 0

Random effects: Name Model spatial.field SPDE2 model

Model hyperparameters: mean sd 0.025quant 0.5quant 0.975quant mode Precision for the Gaussian observations 0.005 0.001 0.003 0.005 0.006 0.005 Theta1 for spatial.field -4.683 0.615 -5.932 -4.664 -3.520 -4.597 Theta2 for spatial.field 1.011 0.646 -0.219 0.995 2.318 0.938

Expected number of effective parameters(stdev): 7.84(2.38) Number of equivalent replicates : 10.97

Watanabe-Akaike information criterion (WAIC) ...: 721.87 Effective number of parameters .................: 9.79

Marginal log-Likelihood: -391.86 CPO and PIT are computed

Posterior marginals for the linear predictor and the fitted values are computed

INLABRU CODE

data = data.frame(catch= shrimp$hauls$catch, depth = shrimp$hauls$depth, group=shrimp$hauls$group) data$group = as.numeric(as.factor(data$group))

data4= SpatialPointsDataFrame( coords = coords, data = data, proj4string = CRS("+proj=longlat +datum=WGS84") )

spde1 <- inla.spde2.matern(mesh, alpha = 2)

comp_inlabru <- catch ~ Intercept(1, prec.linear = 0.0001) + depth + s(main = coordinates, model = spde1, group = group, ngroup =2) inlabru <- bru(comp_inlabru, data = data4, family = "Gaussian", options = list(control.compute = list(cpo = T, waic=T)))

summary(inlabru)

inlabru version: 2.3.1.9000 INLA version: 21.02.23 Components: Intercept: Model types main='linear', group='exchangeable', replicate='iid' depth: Model types main='linear', group='exchangeable', replicate='iid' s: Model types main='spde', group='exchangeable', replicate='iid' Likelihoods: Family: 'Gaussian' Data class: 'SpatialPointsDataFrame' Predictor: catch ~ . Time used: Pre = 1.75, Running = 274, Post = 0.977, Total = 277 Fixed effects: mean sd 0.025quant 0.5quant 0.975quant mode kld Intercept -6.111 23.327 -53.341 -5.911 40.66 -5.419 0 depth 0.068 0.022 0.025 0.068 0.11 0.068 0

Random effects: Name Model s SPDE2 model

Model hyperparameters: mean sd 0.025quant 0.5quant 0.975quant mode Precision for the Gaussian observations 0.005 0.001 0.003 0.005 0.006 0.004 Theta1 for s -4.670 0.661 -6.015 -4.651 -3.416 -4.580 Theta2 for s 0.964 0.755 -0.468 0.942 2.499 0.863 GroupRho for s 0.172 0.612 -0.930 0.279 0.976 0.973

Expected number of effective parameters(stdev): 7.66(2.37) Number of equivalent replicates : 11.23

Deviance Information Criterion (DIC) ...............: 720.53 Deviance Information Criterion (DIC, saturated) ....: 99.11 Effective number of parameters .....................: 9.33

Watanabe-Akaike information criterion (WAIC) ...: 721.98 Effective number of parameters .................: 9.64

Marginal log-Likelihood: -390.57 CPO and PIT are computed

Posterior marginals for the linear predictor and the fitted values are computed

Thank you very much for all the help! Kind Regards,

[finnlindgren]

Please include the summary output so I don't have to run the code myself. Take a look at the reprex package for an easy way to format the example code. https://reprex.tidyverse.org/articles/articles/learn-reprex.html

[finnlindgren]

Note that the INLA default precision for fixed effects is 0.001, not 0.0001. In inlabru, to get the same model for "depth", you may need to add control.fixed=list(prec = 0.0001) to the inlabru options list.

[Filhomn]

I’m sorry for that. I have updated the summary output. Using “control.fixed=list(prec = 0.0001)” worked! Thank you very much for your time and this amazing work with this package! Soon I’m updating the INLA content from all my Bayesian spatial model lectures to inlabru, much easier for the students to learn.

Best Regards,

[finnlindgren]

Good to hear! Please let me know if you run into further issues. I know there may still be a few types of INLA models that are not fully supported by inlabru, or that need special syntax (in particular latent models that have 2n elements in the internal inla storage instead of n elements, and models with special response objects, such as as inla.surv() models; these are next on my list of things to add).

Btw, I notice you forgot to specify control.group in the inlabru s-component definition; that works the same as in INLA. As you can see from the summary output, the default group model in inla (and inlabru) is "exchangeable", not "iid", hence the extra "Rho" parameter output!

[finnlindgren]

You should also make sure to use inla.spde2.pcmatern() wherever possible instead of inla.spde2.matern(), since the PC-priors tend to work much better than the old log-normal priors for the range and variance controlling parameters.

[Filhomn]

Thank you very much! I will try it!

Best Regards,