Boundary constraints in spatial models

[michael]

Hi!

I was experimenting with employing a barrier model for spatial interpolation in a rather complex coastal domain.

However, I want to implement a model with a boundary cosntraint (Dirichlet boundary condition) on land boundaries to set the value there to 0.

From what I gather, the Dirichlet boundary condition is supported for 1D meshes in the INLA library.

There is an old question from 2015 (https://groups.google.com/g/r-inla-discussion-group/c/WKl2PTp_prA/m/42vFiqbm5SMJ) which proposes a few ways of doing this.

Is the solution of removing the rows and columns in the precission matrix and A matrix corresponding to the boundary nodes proposed there still the easiest way to achieve this? Or has further work been done on this?

Best regards,
Michael

[INLA help]

you can also bypass this issue by introducing fake observations, that observe the spatial to be zero on land, with a fixed high (enough) precision, and then everything else will be automatic. 

H

--
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/8082b87e-f1e5-4071-9a82-c6154068e36bn%40googlegroups.com.

[Finn Lindgren]

Hi,

yes, removing the appropriate rows&columns is the most efficient approach, and we've unfortunately not automated that yet.

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/8082b87e-f1e5-4071-9a82-c6154068e36bn%40googlegroups.com.

--

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

[michael]

Thank you both for the fast replies. 
As far as I understood this is what was done in the example code uploaded in the 2015 question by you Finn so I will study that closer. 

I don't have a clear intuition for what removing the rows and columns actually mean in the model other than there will not be any predictions made for these nodes.

With adding observations as proposed in the first response the effects of these values would propagate through the field by the spatial autocorrelation. 
Would I still get this effect by removing the rows and columns?

Kind regards, 
Michael

[Finn Lindgren]

Using the "zero observations" approach will in theory give the same result as Dirichlet conditions; A dirichlet condition it precisely that; an imposed zero value at/outside the boundary.

The "remove rows/columns" approach achieves the same thing by noting that in the basis function representation

u(s) =  \sum_{j=1}^N \psi_j(s) w_j,

setting u(s)=0 on the boundary is the same thing as setting all w_j=0 for nodes on the boundary, which in turn is the same thing as not including those basis functions in the sum.

For the precision, the argument is similar; the precision for the interior node weights is the submatrix of the precision belonging to the interior nodes.

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/0db2c6bd-ba35-4657-828e-132713d3f2d1n%40googlegroups.com.

[michael]

That makes sense.

Thank you both!

Michael

[michael]

Hi!

To do this, I need to specify different precissions for different observations as I do not want the same precision for my normal observations. It is not clear to me how I can do this from the documentation. 
I would have assumed it is done in the control.family argument to INLA, but can't find any mention on how to differentiate between observations there. 

Is this something that can be done easily? 

Best regards,
Michael 

[Finn Lindgren]

With control.inla you would need a separate likelihood for the special
"observations".
But there's also a "scale" parameter to inla() itself, that _might_ be
useful, in that you can set a large scaling factor for the precision
for those particular "observations".
However, since they will still be involved in the parameter
estimation, I would worry that they would influence the precision
estimates too much, so the multi-likelihood options seems much safer.

You just need to make the observation variable a matrix.
Instead of y = c(ordinary_obs, special_obs), use
y=rbind(cbind(ordinary_obs, NA), cbind(NA, special_obs)) .
(that's for plain inla; in inlabru, that matrix is constructed
automatically when building multi-likelihood models, so you don't need
the NA stuff there).

Then
family = c("gaussian", "gaussian")
and
control.family = list(list(...), list(...))
containing the two separate hyper-definitions.

Finn

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/e3a7793f-579d-46ab-8ff2-f2c776223dadn%40googlegroups.com.

[michael sederlin]

Hi Finn, 

Thank you, that helps. 

Would this also work for a survival model where we have both a dependant variable "y" and a censoring variable "cens"? I would imagine that i follow the procedure and create the observation matrices for both the "y" and "cens"? Or is it sufficient to only do this for y and let "cens" be a 1 x (n+m) matrix where n is the normal data and m is the boundary observations? 

Best regards, 
Michael

You received this message because you are subscribed to a topic in the Google Groups "R-inla discussion group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/r-inla-discussion-group/GHJInPRd3T0/unsubscribe.
To unsubscribe from this group and all its topics, 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/CAHC4gzZiOpeZi%3DqHMF1i5igF3cQDHxGDhyCbEm79ty4BudvJsw%40mail.gmail.com.

--

//Michael Sederlin

[michael]

Hi again, 
Just to add to the following message. 

Naturally, the boundary observations are not censored. So a regular Gaussian likelihood could be used for these as well, if this is possible. 

Best regards, 
Michael

[Finn Lindgren]

I'm not sure about the details when involving survival models, but
something along those line should work, yes.

I hope to add more automatic handling of special likelihoods to
inla.stack and/or inlabru soon but haven't had time to look into the
details yet.

Finn

[michael]

Hi!

That would be great!

I have been testing this but been unable to get the observation approach to work in a survival model. It seems to interpret the censoring indicators as well as the observation matrix as observations and complain that the "the number of families 2 does not match number of response variables 6" where my y-variable has dim 2x(N+M) and cens  has size 1 x (N + M). 

But. 
If I may ask a question about the masking implementation of the Dirichlet condition. 
In the example code provided earlier (appended below)  we make a call to inla.spde2.generic directly. It would be very useful to be able to fix the range and standard deviation to better understand the model behaviour. 
In the newer implementation inla.spde.pcmatern this would be done by providing prior.range = c(prior.val, NA) and prior.sigma = c(prior.val, NA). 

I have not managed to figure out how to achieve this in the implementation below, but would it be straight-forward? 

Best regards, 
Michael 


```
library(INLA)

dirichlet.mask <- function(mesh)
{
  mask <- rep(TRUE, mesh$n)
  mask[unique(mesh$segm$bnd$idx)] <- FALSE
  mask
}

dirichlet.fem <- function(mesh, order=2)
{
  output <- c("c0", "c1",
              paste("k", seq_len(order), sep = ""))
  fem <- inla.fmesher.smorg(mesh$loc, mesh$graph$tv, fem = order,
                            output = output)
  mask <- dirichlet.mask(mesh)
  fem$c0 <- fem$c0[mask,mask,drop=FALSE]
  fem$c1 <- fem$c1[mask,mask,drop=FALSE]
  for (idx in seq_len(order)) {
    fem[[paste("g", idx, sep = "")]] <-
      fem[[paste("k", idx, sep = "")]][mask,mask,drop=FALSE]
  }
  fem
}

matern.dirichlet <- function(mesh, B.tau, B.kappa,
                             theta.mu=rep(0, ncol(B.tau)-1),
                             theta.Q=diag(ncol(B.tau)-1))
{
  fem <- dirichlet.fem(mesh)
  inla.spde2.generic(M0=fem$c0, M1=fem$g1, M2=fem$g2,
                     B0=B.tau, B1=2*B.kappa,
                     B2=cbind(1, B.tau[,-1,drop=FALSE]*0.0),
                     theta.mu=theta.mu, theta.Q=theta.Q,
                     transform="identity")
}

dirichlet.A <- function(mesh, ...)
{
  A <- inla.spde.make.A(mesh, ...)
  A[, dirichlet.mask(mesh), drop=FALSE]
}

########################################################################

prior.sd <- 1 ## Change this to something reasonable for your problem
prior.range <- 5 ## Change this to something reasonable for your problem

## For parameterisation details, see Eq (4) and onwards in
## Lindgren & Rue, "Bayesian Spatial Modelling with R-INLA,
## Journal of Statistical Software, 2015
## (also as the "read this first" spde tutorial on r-inla.org)
log.kappa0 <- log(8)/2 - log(prior.range)
log.tau0 <- -log(4*pi)/2 - log(prior.sd) - log.kappa0
B.tau <- cbind(log.tau0, -1, 1)
B.kappa <- cbind(log.kappa0, 0, -1)
theta.mu <- c(0, 0)
theta.prec <- c(0.1, 1)
theta.Q <- diag(theta.prec, length(theta.mu))

lattice <- inla.mesh.lattice(1:30, 1:20)
mesh <- inla.mesh.create(lattice=lattice, extend=FALSE)

## Create a Neumann model
spdeN <- inla.spde2.matern(mesh, B.tau=B.tau, B.kappa=B.kappa,
                           theta.prior.mean=theta.mu,
                           theta.prior.prec=theta.prec)
## Create a Dirichlet model
spdeD <- matern.dirichlet(mesh, B.tau=B.tau, B.kappa=B.kappa,
                          theta.mu=theta.mu,
                          theta.Q=theta.Q)
```

[Finn Lindgren]

Yes, to fix tau and kappa you just need to let B.tau and B.kappa only
have a single column, which gives the log(tau) and log(kappa) values.
For other parameterisations (sigma and range), you just need to work
out the transformation.
When you already have B.tau and B.kappa for a given parameterisation,
and you want to fix the theta vector to some theta0, then you need
B.tau.fixed <- B.tau[, 1, drop = FALSE] + B.tau[, -1, drop = FALSE] %*% theta0
B.kappa.fixed <- B.kappa[, 1, drop = FALSE] + B.kappa[, -1, drop =
FALSE] %*% theta0


Finn

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/a8744170-f4e9-48e4-929a-4d8cceb369f2n%40googlegroups.com.

[michael sederlin]

Hi Finn, 

I see! It is not trivial (to me at least: ) ) to get a thorough understanding of the inner workings on R-INLA.

So my understanding so far is that we implement the boundary condition through masking the mesh nodes in the fem - objects created by fmesher.smorg. Normally this is called in inla.spde2.matern which then calls inla.spde2.generic. 
I would like to use the pc-matern prior with the parametrisation prior.range = c(rho0, prho) and prior.sigma = c(sigma0, psigma)  and have the ability to fix these parameters for debugging and understanding the model (as I mentioned before). 

So, as far as I now understand it would be possible to create modified versions of inla.spde2.pcmatern and inla.spde2.matern which perform the masking and output the model object to achieve this? 

I have appended code of a small example with the modified functions inla.dirichlet.spde2.pcmatern and inla.dirichlet.spde2.matern in dirichlet.R and a simple test case in dirichlet_model_test.R 
The only change is that the masking vector is provided as an input to pcmatern which is further passed to matern() which performs the masking before creating the spde model object with the default inla.spde2.generic() 

From the testing I have done, this seems to work ok. But it would also be very comforting with a confirmation that I have not misunderstood. 

It uses INLA 22.05.03

"This is INLA_22.05.03 built 2022-05-03 07:58:22 UTC. 

 - See www.r-inla.org/contact-us for how to get help."

All the best, 
Michael 

--

//Michael Sederlin