INLA and inlabru model comparison

[Brian]

I'm encountering a discrepancy between INLA and inlabru when fitting spatial models with beta-binomial family. For binomial models with formulas like y ~ 0 + b0 + f(s, model = spde) and y ~ 0 + b0 + f(s, model = spde) + f(rr, model = "iid"), both INLA and inlabru produce nearly identical range estimates for the spatial field (differences only in decimal places). 

However, when I switch to family = "betabinomial" with identical data, mesh, and priors, inlabru consistently gives range estimates that are approximately half of what INLA produces (e.g., INLA: ~125 km vs inlabru: ~64 km). Is this a known issue with different numerical approximations being used for the beta-binomial likelihood between the two approaches, and which should be considered more reliable?

betabinomial <- inla(
  formula =  y ~ 0 + b0 + f(s, model = spde),
  family = "betabinomial",
  Ntrials = numtrials,
  data = inla.stack.data(stk.full),
  control.fixed = list(mean = 0, prec = 1),
  control.predictor = list(compute = TRUE, link = 1, A = inla.stack.A(stk.full)),
  control.compute = list(config = TRUE, return.marginals.predictor = TRUE, cpo = T, waic = T, dic = T),
  control.inla = list(tolerance = 1e-7),
  num.threads = 7,
  verbose = F
)

betabinomial_inlabru <- bru(
  n ~ 0 +
    intercept(
      1,
      model = "linear",
      mean.linear = 0,
      prec.linear = 1
    ) +
    field(as.matrix(coords.clusters), model = spde),
  Ntrials = dm2$N,
  data = dm2,
  family = "betabinomial"
)

I have ensured the `coord.clusters` are in the required units (KM) and CRS.

head( coords.clusters)
       km_x     km_y
      <num>    <num>
1: 337.4873 10011.49
2: 335.7025 10010.76
3: 334.3113 10022.43
4: 340.8206 10009.99
5: 343.4788 10007.99
6: 343.8337 10018.47

The spde was created as;

mesh <- fm_mesh_2d(
  boundary = kenya_shapefile2,
  loc = coords.clusters,
  max.edge = c(20, 40),
  cutoff = 10
)

spde <- inla.spde2.pcmatern(
  mesh = mesh,
  alpha = 2,
  prior.range = c(r0, .01),  
  prior.sigma = c(1, .01)  
)

[Finn Lindgren]

Although I can’t see all the details of the data (e.g. I would double check that you’re actually using the same values for Ntrials!), I do see one discrepancy in how you’re doing the estimation; for inla(), you’re setting the optimisation tolerance, but not for inlabru. Please supply the same control.inla argument as options=list(control.inla=…) when calling bru(), to rule out that causing the difference. Since the model predictor is linear, and you’re setting the same priors for both models, inlabru will simply call inla(), and doesn’t do anything special with the likelihood. The only difference should be for nonlinear models, where inlabru might set different tolerances in the subsequent iterations. Also, the actual estimation summaries would be useful, as well as the actual value of r0.

Finn

On 5 Aug 2025, at 09:52, Brian <brian...@gmail.com> wrote:

I'm encountering a discrepancy between INLA and inlabru when fitting spatial models with beta-binomial family. For binomial models with formulas like y ~ 0 + b0 + f(s, model = spde) and y ~ 0 + b0 + f(s, model = spde) + f(rr, model = "iid"), both INLA and inlabru produce nearly identical range estimates for the spatial field (differences only in decimal places). 

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/1db4697e-bf79-4210-b6e0-8fb83859997an%40googlegroups.com.

[Brian]

Thank you Finn for the response. 

I have supplied the  control.inla argument in the bru() function. I have noted that the range has increased slightly to arounf 88Km. The number of trials for inla() was supplied to inla.stack() as shown in the code below, but calls from data frame named 'dm2'

stk.e1 <- inla.stack(
  tag = 'point',
  data = list(y = dm2$n, numtrials = dm2$N),
  A = list(A, 1, 1),
  effects = list(
    s = indexs,
    rr = 1:length(dm2$improved_water),
    b0 = rep(1, nrow(dm2))
  )
)

stk.p1 <- inla.stack(
  tag = "pred",
  data = list(y = NA, numtrials = NA),
  A = list(Ap, 1, 1),
  effects = list(
    s = indexs,
    rr = (length(pred$improved_water) + 1):(length(pred$improved_water) + nrow(pred)),
    b0 = rep(1, nrow(pred))
  )
)

> any(inla.stack.data(stk.e1)$numtrials != dm2$N)
[1] FALSE

r0 is set to be 5% of kenya in north-south direction

> r0 <- as.vector(0.05 * (st_bbox(kenya_shapefile2)["ymax"] - st_bbox(kenya_shapefile2)["ymin"]))
> r0
[1] 54.11497

betabinomial_inlabru <- bru(
  n ~ 0 +
    intercept(
      1,
      model = "linear",
      mean.linear = 0,
      prec.linear = 1
    ) +
    field(as.matrix(coords.clusters), model = spde),

  options = list(control.inla = list(tolerance = 1e-7)),

  Ntrials = dm2$N,
  data = dm2,
  family = "betabinomial"
)

Here is a summary from inla()

Time used:
    Pre = 0.942, Running = 206, Post = 3.06, Total = 210
Fixed effects:
     mean   sd 0.025quant 0.5quant 0.975quant   mode kld
b0 -4.696 0.26     -5.207   -4.696     -4.185 -4.696   0

Random effects:
  Name   Model
    s SPDE2 model

Model hyperparameters:
                                                    mean    sd 0.025quant 0.5quant 0.975quant    mode
overdispersion for the betabinomial observations   0.103 0.002      0.099    0.103      0.106   0.103
Range for s                                      124.667 0.892    122.840  124.691    126.349 124.808
Stdev for s                                        1.278 0.025      1.225    1.279      1.322   1.287

Deviance Information Criterion (DIC) ...............: 3138.02
Deviance Information Criterion (DIC, saturated) ....: 1684.69
Effective number of parameters .....................: 103.01

Watanabe-Akaike information criterion (WAIC) ...: 3128.23
Effective number of parameters .................: 87.60

Marginal log-Likelihood:  -1662.89
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Here is a summary from bru()

inlabru version: 2.13.0
INLA version: 25.03.13
Components:
intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
field: main = spde(as.matrix(coords.clusters)), group = exchangeable(1L), replicate = iid(1L), NULL
Observation models:
  Family: 'betabinomial'
    Tag: <No tag>
    Data class: 'tbl_df', 'tbl', 'data.frame'
    Response class: 'integer'
    Predictor: n ~ .
    Additive/Linear: TRUE/TRUE
    Used components: effects[intercept, field], latent[]
Time used:
    Pre = 1.11, Running = 180, Post = 0.727, Total = 182
Fixed effects:
            mean    sd 0.025quant 0.5quant 0.975quant   mode kld
intercept -4.934 0.203     -5.329   -4.935     -4.532 -4.935   0

Random effects:
  Name   Model
    field SPDE2 model

Model hyperparameters:
                                                   mean    sd 0.025quant 0.5quant 0.975quant  mode
overdispersion for the betabinomial observations  0.132 0.023      0.104    0.127       0.19  0.11
Range for field                                  88.757 9.699     68.437   89.131     105.83 93.24
Stdev for field                                   1.070 0.160      0.729    1.082       1.33  1.18

Marginal log-Likelihood:  -1648.46
 is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

[Finn Lindgren]

Check the results after running Inla.rerun() and bru_rerun(), respectively. Might be the optimisation is having difficulties.

Finn

On 5 Aug 2025, at 12:12, Brian <brian...@gmail.com> wrote:

Thank you Finn for the response. 

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/3f5023cc-436b-4f2b-b851-e9a3fc6da5c0n%40googlegroups.com.

[Brian]

Hi Finn,  

Thank you for the suggestion! I have now run inla.rerun() and bru_rerun() on both models, and I can confirm that the results are now identical between INLA and inlabru. Both approaches now give a range of ~64 km (95% CI: 48-85), resolving the previous discrepancy where INLA gave ~125 km. However, when rerunning the models, I encountered the following warnings:

*** WARNING *** GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=297
*** WARNING *** GMRFLib_2order_approx: reset counter for 34 NAN/INF values in logl

I have three follow-up questions:

1. What do these warnings mean and should I be concerned about them for my results?

2. Should I routinely use inla.rerun()/bru_rerun() after fitting models, as a best practice?

betabinomial_rerun <- inla.rerun(betabinomial)
 *** WARNING *** GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=297

> betabinomial_inlabru_rerun <- bru_rerun(betabinomial_inlabru)
 *** WARNING *** GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=645
 *** WARNING *** GMRFLib_2order_approx: reset counter for 34 NAN/INF values in logl

> summary(betabinomial_rerun)
Time used:
    Pre = 0.929, Running = 12.8, Post = 1.47, Total = 15.2

Fixed effects:
     mean    sd 0.025quant 0.5quant 0.975quant   mode kld

b0 -5.591 0.289     -6.172   -5.586     -5.034 -5.585   0

Random effects:
  Name   Model
    s SPDE2 model

Model hyperparameters:
                                                   mean    sd 0.025quant 0.5quant 0.975quant

overdispersion for the betabinomial observations  0.099 0.011       0.08    0.099      0.121
Range for s                                      64.121 9.435      47.84   63.344     84.908
Stdev for s                                       1.872 0.162       1.58    1.864      2.214
                                                   mode
overdispersion for the betabinomial observations  0.098
Range for s                                      61.652
Stdev for s                                       1.846

Deviance Information Criterion (DIC) ...............: 3120.84
Deviance Information Criterion (DIC, saturated) ....: 1680.37
Effective number of parameters .....................: 251.78

Watanabe-Akaike information criterion (WAIC) ...: 3103.86
Effective number of parameters .................: 194.00

Marginal log-Likelihood:  -1620.93

CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

> summary(betabinomial_inlabru_rerun)

inlabru version: 2.13.0
INLA version: 25.03.13
Components:
intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
field: main = spde(as.matrix(coords.clusters)), group = exchangeable(1L), replicate = iid(1L), NULL
Observation models:
  Family: 'betabinomial'
    Tag: <No tag>
    Data class: 'tbl_df', 'tbl', 'data.frame'
    Response class: 'integer'
    Predictor: n ~ .
    Additive/Linear: TRUE/TRUE
    Used components: effects[intercept, field], latent[]
Time used:

    Pre = 0.905, Running = 8.37, Post = 0.26, Total = 9.54

Fixed effects:
            mean    sd 0.025quant 0.5quant 0.975quant   mode kld

intercept -5.591 0.289     -6.171   -5.586     -5.034 -5.585   0

Random effects:
  Name   Model
    field SPDE2 model

Model hyperparameters:
                                                   mean    sd 0.025quant 0.5quant 0.975quant

overdispersion for the betabinomial observations  0.099 0.010       0.08    0.099      0.121
Range for field                                  64.120 9.356      47.97   63.350     84.730
Stdev for field                                   1.872 0.162       1.58    1.864      2.213
                                                   mode
overdispersion for the betabinomial observations  0.098
Range for field                                  61.665
Stdev for field                                   1.846

Marginal log-Likelihood:  -1620.95

[Finn Lindgren]

Did you also try setting the same tolerance option for bru to see if that would make it give the same result as a inla() call? And vice versa, if inla() would match bru() if you _dont_ set the tolerance? Its clearly an optimiser issue, so the tolerance argument is likely to matter…

Finn

On 5 Aug 2025, at 16:55, Brian <brian...@gmail.com> wrote:

Hi Finn,  

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/79b19276-6519-4463-8fca-918314c326b7n%40googlegroups.com.

[Brian]

Yes. If set the same tolerance option for bru(), the range is still different if I don't use bru_rerun(), However, when I remove the tolerance option in both, the range is the same even without using bru_rerun() and inla.rerun(). I agree that the tolerance argument matters in this case.