Differences between factor estimates with INLA vs inlabru

[ASeatonSpatial]

I am using inlabru to fit the spatio-temporal model from section 7.1 in the SPDE textbook. This model includes a factor variable with 3 levels and has no intercept. When I fit the model using INLA I get parameter estimates:

> inla_fit$summary.fixed
         mean        sd 0.025quant   0.5quant 0.975quant       mode          kld
wA -1.7945173 0.3231770 -2.4312689 -1.7954229 -1.1536005 -1.7972378 1.737425e-06
wB -0.7939925 0.3231795 -1.4307402 -0.7949011 -0.1530616 -0.7967220 1.739200e-06
wC  0.2050322 0.3231785 -0.4317095  0.2041221  0.8459657  0.2022984 1.743636e-06

When I fit the model using inlabru it is interesting that the factor is not seen as a "fixed effect" but instead viewed as a random effect. There are also only 2 levels. It looks like level A is set to zero by inlabru.

> inlabru_fit$summary.random$fac
  ID     mean          sd 0.025quant 0.5quant 0.975quant     mode kld
1  B 1.000462 0.005475537  0.9897114 1.000462   1.011203 1.000462   0
2  C 1.999512 0.005449700  1.9888128 1.999512   2.010203 1.999512   0

Here is the output of INLA:::summary.inla(inlabru_fit) without the call (it's very long and the result of all the internal inlabru code to convert the inlabru syntax to something INLA understands, I've put the inlabru code below instead)

Time used:
    Pre = 2.56, Running = 49, Post = 1.11, Total = 52.7 
Random effects:
  Name    Model
    fac IID model
   grf SPDE2 model

Model hyperparameters:
                                          mean    sd 0.025quant 0.5quant 0.975quant
Precision for the Gaussian observations 97.316 2.848     91.789   97.291    103.004
Range for grf                            3.026 0.188      2.683    3.016      3.421
Stdev for grf                            1.076 0.058      0.970    1.073      1.197
GroupRho for grf                         0.718 0.019      0.679    0.718      0.756
                                          mode
Precision for the Gaussian observations 97.268
Range for grf                            2.991
Stdev for grf                            1.066
GroupRho for grf                         0.718

Expected number of effective parameters(stdev): 1130.21(0.00)
Number of equivalent replicates : 3.27 

Deviance Information Criterion (DIC) ...............: -5299.44
Deviance Information Criterion (DIC, saturated) ....: 4829.44
Effective number of parameters .....................: 1130.21

Watanabe-Akaike information criterion (WAIC) ...: -5287.93
Effective number of parameters .................: 922.77

Marginal log-Likelihood:  1255.41 
Posterior marginals for the linear predictor and
 the fitted values are computed

It doesn't appear as though an intercept has been included anyway. inlabru_fit$summary.fixed is empty. How can I reconcile these two different estimates for the factors? I should say that the model seemed to fit with no issues both in inlabru and INLA and the predictions from the SPDE component look basically identical across both models so I think I have fitted the same model but it could be I have not.

This is the formula for the INLA model (I've not included the code for the stack or the priors). w is the name of the factor covariate.

# Model formula
formulae <- y ~ 0 + w + f(i, model = spde, group = i.group, 
                          control.group = list(model = 'ar1', hyper = h.spec)) 

# Model fitting
inla_fit <- inla(formulae,  data = inla.stack.data(sdat), 
            control.predictor = list(compute = TRUE,
                                     A = inla.stack.A(sdat)), 
            control.family = list(hyper = list(theta = prec.prior)), 
            control.fixed = list(expand.factor.strategy = 'inla'))

This is what I did in inlabru:

cmp = y ~ 0 + 
  fac(main = w, model = "factor") +
  grf(main = cbind(xcoo, ycoo),
      model = spde,
      group = time,
      control.group = list(model = 'ar1', hyper = h.spec))

inlabru_fit = bru(cmp,
                  data = dat,
                  family = "gaussian",
                  options = list(control.predictor = list(compute = TRUE,
                                                          A = inla.stack.A(sdat)), 
                                 control.family = list(hyper = list(theta = prec.prior)), 
                                 control.fixed = list(expand.factor.strategy = 'inla')))

I am assuming there is some difference in the way inlabru handles factors by default vs whatever INLA does. I couldn't find anything that mentions this in the documentation, could be a gap there.

[ASeatonSpatial]

The reason this happens is because by default inlabru currently sets first level to 0 and treats it as a contrast. Changing this to an "iid" model with fixed precision seems to do the same thing as including w in the linear predictor in INLA:

cmp = y ~ 0 + 
  fac(main = w, 
      model = "iid",
      hyper = list(prec = list(initial = log(1e-6), 
                             fixed = TRUE))) +
  grf(main = cbind(xcoo, ycoo),
      model = spde,
      group = time,
      control.group = list(model = 'ar1', hyper = h.spec))

which gives

> inlabru_fit$summary.random$fac
  ID       mean        sd 0.025quant   0.5quant 0.975quant       mode kld
1  A -1.8032138 0.2921339 -2.3767711 -1.8032220 -1.2301351 -1.8032138   0
2  B -0.8027115 0.2921354 -1.3762718 -0.8027197 -0.2296298 -0.8027115   0
3  C  0.1963100 0.2921343 -0.3772481  0.1963017  0.7693894  0.1963100   0
> inla_fit$summary.fixed
         mean        sd 0.025quant   0.5quant 0.975quant       mode          kld
wA -1.7945173 0.3231770 -2.4312689 -1.7954229 -1.1536005 -1.7972378 1.737425e-06
wB -0.7939925 0.3231795 -1.4307402 -0.7949011 -0.1530616 -0.7967220 1.739200e-06
wC  0.2050322 0.3231785 -0.4317095  0.2041221  0.8459657  0.2022984 1.743636e-06

[finnlindgren]

Fixed in commit 0055058. The factor model is deprecated in favour of factor_contrast, and factor_full can now be used for the non-contrast fixed precision version.

[ASeatonSpatial]

Just in case this is useful to anyone who finds their way here. The way to do this in the latest version is now:

cmp = y ~ 0 +
  fac(main = w,
      model = "factor_full") +
  grf(main = cbind(xcoo, ycoo),
      model = spde,
      group = time,
      control.group = list(model = 'ar1', hyper = h.spec))

This seems to work and gives same results as above. Note, doing something like

cmp = y ~  0 + w

will throw an error if w is a factor variable.