fixed effects with spde sub-models
20 views
Skip to first unread message
Courtney Taylor
May 26, 2026, 10:24:17 AM (9 days ago) May 26
to R-inla discussion group
I am wondering if there is a difference between the following model structures in inla.
Model 1: Explicitly defines the fixed effects. spde and fixed components share the same weights so constr = T in the definition of the mesh inside inla.spde2.pcmatern()
components <- ~ Intercept(1) +
fixed(con, model = 'fixed') +
lambda(geometry, model = matern) +
conspecific(geometry, weights = con, model = matern)
Model 2: Considers the fixed effects as a part of the spde
components <- ~ -1 +
lambda(geometry, model = matern) +
conspecific(geometry, weights = con, model = matern)
I have fitted both of these models to my data and I get different results (big differences in the range and stdev hyperparameters) but very similar WAIC scores. Is there a more “correct” way to specify the model?
Model 1: Explicitly defines the fixed effects. spde and fixed components share the same weights so constr = T in the definition of the mesh inside inla.spde2.pcmatern()
components <- ~ Intercept(1) +
fixed(con, model = 'fixed') +
lambda(geometry, model = matern) +
conspecific(geometry, weights = con, model = matern)
Model 2: Considers the fixed effects as a part of the spde
components <- ~ -1 +
lambda(geometry, model = matern) +
conspecific(geometry, weights = con, model = matern)
I have fitted both of these models to my data and I get different results (big differences in the range and stdev hyperparameters) but very similar WAIC scores. Is there a more “correct” way to specify the model?
Finn Lindgren
May 26, 2026, 10:53:08 AM (8 days ago) May 26
to Courtney Taylor, R-inla discussion group
Hi,
So if you used constr=TRUE, then your model 1 can be written as
it's not clear what you mean by "spde and fixed components share the same weights so constr=T".
constr=TRUE for inla.spde2.pcmatern() means "the field should integrate to zero".
constr=TRUE for inla.spde2.pcmatern() means "the field should integrate to zero".
So if you used constr=TRUE, then your model 1 can be written as
eta(s) = Intercept * 1 + fixed_con * con(s) + lambda(s) + conspecific(s) * con(s)
with constraints
\int lambda(s) ds = 0
with constraints
\int lambda(s) ds = 0
\int conspecific(s) sd = 0
and
Intercept and fixed_con are latent coefficients.
Model 2, assuming no constr (i.e. constr=FALSE, the default):
eta(s) = lambda(s) + conspecific(s) * con(s)
without constraints.
In Model 2, lambda and conspecific are forced to include large scale phenomena that in model 1 are handled by the Intercept and fixed_con.
Therefore, one should expect the range and sigma values in Model 2 to be larger than in Model 1, as they handle different aspects.
Both models are valid, but the interpretation of range and sigma are of course very different.
If we let lambda1 = Intercept + lambda in Model 1, and lambda2 = lambda from model 2, we would normally expect
and
Intercept and fixed_con are latent coefficients.
Model 2, assuming no constr (i.e. constr=FALSE, the default):
eta(s) = lambda(s) + conspecific(s) * con(s)
without constraints.
In Model 2, lambda and conspecific are forced to include large scale phenomena that in model 1 are handled by the Intercept and fixed_con.
Therefore, one should expect the range and sigma values in Model 2 to be larger than in Model 1, as they handle different aspects.
Both models are valid, but the interpretation of range and sigma are of course very different.
If we let lambda1 = Intercept + lambda in Model 1, and lambda2 = lambda from model 2, we would normally expect
lambda1 and lambda2 to behave similarly, but the priors are different; lambda1 has a separate prior for its spatial average (and a constrained version of the Matern model for the spatial variation),
whereas the prior for the spatial average of lambda2 is implicitly defined via the Matern model.
Which one is "best"/"most useful" is extremely problem dependent.
In fact, the third option, Model 3, where one takes Model 1 and removes the constraints, is often a practical alternative, where again
lambda3 = Intercept + lambda (with no explicit constraints) typically behaves like lambda1, but where the Intercept and lambda component are mutually non-identifiable.
The main reason this Model 3 can be useful is that it's not always clear over which region one should apply the integrate-to-zero constraint;
the whole mesh, only the domain of interest, or only the observed subregion or set of observation points.
Model 3 removes that arbitrary choice, at the expense of component non-identifiability.
Which one is "best"/"most useful" is extremely problem dependent.
In fact, the third option, Model 3, where one takes Model 1 and removes the constraints, is often a practical alternative, where again
lambda3 = Intercept + lambda (with no explicit constraints) typically behaves like lambda1, but where the Intercept and lambda component are mutually non-identifiable.
The main reason this Model 3 can be useful is that it's not always clear over which region one should apply the integrate-to-zero constraint;
the whole mesh, only the domain of interest, or only the observed subregion or set of observation points.
Model 3 removes that arbitrary choice, at the expense of component non-identifiability.
Model 2 also removes that arbitrary choice, but at the expense of potential issues with the range & sigma parameters.
In summary, "there's no free lunch".
Finn
Finn
--
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/b2ee1ada-47d8-457c-9973-acfa251e1e91n%40googlegroups.com.
Courtney Taylor
May 26, 2026, 11:20:15 AM (8 days ago) May 26
to R-inla discussion group
Thanks very much for this helpful response. I have a follow up question.
When comparing the results of model 1 versus model 2 i also noticed that the spatial fields were more precisely estimated in model 2. That is, the credible intervals were much narrower for the range and sigma hyperparameters. Is there an obvious reason that this would be the case?
Finn Lindgren
May 26, 2026, 11:41:12 AM (8 days ago) May 26
to Courtney Taylor, R-inla discussion group
I would say that
" the spatial fields were more precisely estimated"
does _not_ mean the same thing as
"the credible intervals were much narrower for the range and sigma hyperparameters"
The first statement should refer to the credible intervals for "lambda1" and "lambda2" values (in the notation from my reply), and those are likely to be similar in the two models, as they are linked to the observations in the same way in both models, and at least in principle can represent the same spatial functions. (whereas "lambda" in model 1 and "lambda" in model 2 are fundamentally different, conceptually)
The hyperparameters on the other hand _are not the same in the two models_, as the constraints in model 1 _fundamentally_ changes their interpretation,
so it isn't really meaningful to compare them in that way. They have the same _names_, but they are not the same _conceptually_.
Finn
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/8476db40-d03c-463b-9a11-e0d05a193c34n%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages