inla.spde2.matern

[dan]

Hi

I would like to translate a geospatial model to inla with the appropriate prior distributions. This mostly works as desired, but I do not know how to set my priors properly for the spatial field. I assume a exponential correlation of the form

cor(s1, s2) = sigma^2 * exp(- rho * dist(s1, s2))

and I would like to set priors 1/sigma^2 ~ Gamma(a1, a2) and rho ~ Unif(b1,b2). Is this possible with INLA?
I have read Lindgren and Rue and several books but can not find a solution which comes close. So far I use

inla.spde2.matern(mesh, alpha = 1.5) or

inla.spde2.pcmatern(mesh, alpha = 1.5, prior.range = c(0.3, 0.5), prior.sigma = c(10, 0.01))

For the inla.spde2.matern there are options prior.variance.nomial, prior.range.nominal which sound promising. However I could not find a documentation as to what they do exactly which I could understand.

Any help is appreciated, thanks in advanced.

[Finn Lindgren]

Hi Dan,

the sigma-prior is possible; I'm not sure if the "uniform on a finite interval" prior for rho is possible.

The inla.spde2.matern implements a general family that links a set of theta-parameters to tau and kappa in an spde formulation of the Matern family, (kappa^2-Laplacian)^(alpha/2)(tau x(s)) = W(s) (see the function help text).

The parameter link is user-defined, so allows some flexibility.

Smoothness:

  nu = alpha - dim/2 = 0.5 for alpha=1.5.

Variance:

  sigma^2 = const / kappa^(2*nu) / tau^2

Precision:

  1/sigma^2 = kappa * tau^2 / const

The range controlling parameter rho in your parameterisation is kappa (i.e. rho is proportional to the inverse range), so

  theta1 = log(1/sigma^2) = - log(const) + log(kappa) + 2*log(tau)

  theta2 = log(rho) = log(kappa)

Solving for log(tau) and log(kappa),

  log(tau) = -log(const) / 2 + theta1 / 2 + theta2

  log(kappa) = theta2,

so you need

  B.tau = cbind(-log(const) / 2, 1/2, 1)

  B.kappa = cbind(0, 0, 1)

in order to get theta1 and theta2 to be the log() of 1/sigma^2 and rho.

The constant "const" can be derived from the expression for the variance in Lindgren et al (2011, JRSSB) or Lindgren and Rue 2015 (https://www.jstatsoft.org/article/view/v063i19).

You can then set the prior for theta1 to a loggamma(a1,a2) distribution to get your desired prior for 1/sigma^2, via

  hyper=list(theta1=list(...))

For theta2, you'd be looking for a prior with density  exp(theta2)/(b2-b1) on the interval (log(b1), log(b2)),

which isn't directly available. I see there's an "expression" prior model that might be able to specify a general density expression, but the finite interval support is bound to cause problems.

So in summary, yes, it's partly possible, but in particular the rho-prior is a bit problematic, and from a modelling perspective it's probably better to use the pc-prior option via inla.spde2.pcmatern instead.

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 on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/9762bfc7-adfb-45ee-ba10-731eba9ec82fn%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com

[dan]

Hi Finn

Thank you very much for this quick and thorough response, that helps me a lot. Just to clarify details in you calculation:

  theta1 = log(1/sigma^2) = - log(const) + log(kappa) + 2*log(tau)

  theta2 = log(rho) = log(kappa)

Solving for log(tau) and log(kappa),

  log(tau) = -log(const) / 2 + theta1 / 2 + theta2

  log(kappa) = theta2,

so you need

  B.tau = cbind(-log(const) / 2, 1/2, 1)

  B.kappa = cbind(0, 0, 1)

I got:

Solving for log(tau) and log(kappa),

  log(tau) = + log(const) / 2 + theta1 / 2 - theta2 /2

  log(kappa) = theta2,

so you need

  B.tau = cbind(+ log(const) / 2, 1/2,  - 1/2)

  B.kappa = cbind(0, 0, 1)

Do you agree?

Assuming that I would like a Gamma prior also for the kappa (alias rho), the call to inla.spde2.matern is then

INLA::inla.spde2.matern(

  mesh = mesh,

  alpha = 1.5,

  B.tau = cbind(log(const) / 2, .5, -.5),

  B.kappa = cbind(0, 0, 1),

  hyper = list(theta1 = list(param = c(a1, a2), prior = "loggamma"),

               theta2 = list(param = c(b1, b2), prior = "loggamma"))

)?

I just ask because there is no hyper argument specified for the inla.spde2.matern and I have learned not to trust the ... argument too much in R.

I have to wrap my head around how to get a suitable value for this const and will consider the resources that you have suggested, thanks again. 

Dan

[dan]

Sorry, I hit send to early. For the second part you probably mean:

spde <- INLA::inla.spde2.matern(

  mesh = mesh,

  alpha = 1.5,

  B.tau = cbind(log(const) / 2, .5, -.5),

  B.kappa = cbind(0, 0, 1)

)

formula <- y ~ 0 + f(spatial.field, model = spde, 

      hyper = list(theta1 = list(param = c(a1, a2), prior = "loggamma"),

                           theta2 = list(param = c(b1, b2), prior = "loggamma")))

?

[Finn Lindgren]

Hi,

yes, you're right, I missed one of the 1/2 factors, as well as the sign on the log(const).

And yes, the hyper argument is for the f() specification, as you wrote in your other reply.

(Or you can alter the spde object by replacing the contents of spde$f$hyper.default instead, after the inla.spde2.matern call and before setting up the inla model.)

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/90680bf2-17fc-444c-8603-ea2ea5b89a0fn%40googlegroups.com.

[dan]

Great! Thanks for the clarifications.

[dan]

Hello again

Unfortunately, the said solution does not work. I put together a reproducible example borrowed from https://becarioprecario.bitbucket.io/spde-gitbook/ch-intro.html :

rm(list = ls())

library(INLA)

pl.dom <- cbind(c(0, 1, 1, 0.7, 0), c(0, 0, 0.7, 1, 1))

mesh5 <- inla.mesh.2d(loc.domain = pl.dom, max.e = c(0.092, 0.2))

coords <- as.matrix(SPDEtoy[, 1:2])

A5 <- inla.spde.make.A(mesh5, loc = coords)

spde5 <- inla.spde2.matern(

  mesh = mesh5,

  alpha = 1.5,

  B.tau = matrix(c(-0.5 * log(4 * pi), 0.5,-0.5), 1, 3),

  B.kappa =  matrix(c(0, 0, 1), 1, 3)

)

stk5 <- inla.stack(

  data = list(resp = SPDEtoy$y),

  A = list(A5, 1),

  effects = list(i = 1:spde5$n.spde,

                 beta0 = rep(1, nrow(SPDEtoy))),

  tag = 'est'

)

# working default ---------------------------------------------------------

res5 <- inla(

  resp ~ 0 + beta0 + f(i, model = spde5),

  data = inla.stack.data(stk5),

  control.predictor = list(A = inla.stack.A(stk5))

)

# try1 to reset prior -----------------------------------------------------

res5 <- inla(

  resp ~ 0 + beta0 + f(

    i,

    model = spde5,

    hyper = list(

      theta1 = list(param = c(2.01, 1.01), prior = "loggamma"),

      theta2 = list(param = c(2.01, 1.01), prior = "loggamma")

    )

  ),

  data = inla.stack.data(stk5),

  control.predictor = list(A = inla.stack.A(stk5))

)

# try2 to reset prior -----------------------------------------------------

spde5$f$hyper.default <- list(

  theta1 = list(param = c(2.01, 1.01), prior = "loggamma"),

  theta2 = list(param = c(2.01, 1.01), prior = "loggamma")

)

res5 <- inla(

  resp ~ 0 + beta0 + f(i, model = spde5),

  data = inla.stack.data(stk5),

  control.predictor = list(A = inla.stack.A(stk5))

)

I get the following error:

Error in inla.set.hyper(model = model, section = "latent", hyper = hyper,  : 

  Setting hyperparameter `theta2'. Key `prior' = 'none' is 'read-only', cannot change to `loggamma'.

I could not find the inla.set.hyper in the documentation for the users, so I can not tell what is wrong, but telling from the 'read-only' hint, is it forbidden to change priors for theta1 and theta2? Is there a way to work around this, or have I misunderstood some of the hints above?

Best, Dan

[Finn Lindgren]

Hi,

is the error from both the "try1" and "try2" approaches, or only from one of them?

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/fddf8e04-40b7-4a06-9640-376e33538c62n%40googlegroups.com.

[dan]

Hi

Try1: "Error in inla.set.hyper(model = model, section = "latent", hyper = hyper,  : 

  Setting hyperparameter `theta2'. Key `prior' = 'none' is 'read-only', cannot change to `loggamma'."

Try2: "Error in inla.set.hyper(model, section, hyper = hyper.default, debug = debug) : 

  Setting hyperparameter `theta2'. Key `prior' = 'none' is 'read-only', cannot change to `loggamma'."

slight difference in the first lines.

Dan

[Finn Lindgren]

Hi Dan,

I checked this with Haavard, and it turns out that the internals for the spde2 models can only use special multivariate priors; the default log-normal prior and the joint PC-prior for range and sigma are what's implemented.

Either one would need to add the gamma-product as a bivariate prior, or implement the model separately, via inla.rgeneric.

Since for ordinary modelling, I would recommend just using the PC-prior; there, sigma and 1/range are given Exponential priors (for dim=2; range^{-dim/2}~Exp for other dimensions).

Finn

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/584e28b7-a021-4803-9870-a13985039970n%40googlegroups.com.

[Finn Lindgren]

Pressed send too soon...

Meant to add, that it's not clear it's worth the extra work to implement the special "independent bivariate (log)gamma" prior.

Finn

[dan]

Dear Finn

Thank you very much for you answer. I'll have to digest this and see if and how I can adapt this to my case.

Best, Daniel