Mesh building on a global scale

[mike.bl...@gmail.com]

Hi all,

 

I am new to INLA, but I have spent the last couple of weeks reading and playing with the various tutorials and the documentation. I feel that I have got to grips with the basics of the software. Thanks to you all for producing so much awesome documentation for this amazing package!

 

I would like to ask you experts and other users for some guidance and advice on my model. I apologise that this is such a long post! Please just respond to whatever you can/have time to.

 

My model is a species distribution model predicting the probability of occurrence of a species of insect as a function of various bioclimatic and demographic variables at a global scale. The data are presence/absence (~15,000 observations). The covariates are mostly raster data with a 5km*5km global resolution (WorldClim). The model has a binomial likelihood, some fixed effects covariates and a barrier.model spde spatial random effect:

 

y ~ b0 + BIOC_1 + BIOC_2 + ... + BIOC_p + f(s, model = barrier.spde) - 1

 

 

I have quite a few questions(!), but first I would like to ask for some advice on mesh building at a global scale. 

 

To start with I took a shapefile of global coastlines. I simplified the polygons (using gSimplify from rgeos) and removed all the small islands below a specific area size. The aim of this was to try to avoid too complex a shape, so I can achieve nice uniform triangulation in my sampling domain. I then used these spatial polygons to define the inner boundary of the mesh, and used the observation points as initial vertices for the triangles. (Currently I am using decimal degrees, which I appreciate is inappropriate at a global scale, but I just wanted to get my head around the mesh building before implementing the spherical mesh).

 

## Domain boundary

  bndry <- inla.sp2segment(

  shape,

  crs = crs(shape)

)

## Mesh

mesh <- inla.mesh.2d(

  loc = coords.est,  

  boundary = bndry,

  cutoff = max.edge.length/5,   

  max.edge = c(1, 5) * max.edge.length,

  min.angle = 21  

)

 

However, the observations have a pretty extreme non-random distribution. They are very clumped in certain areas. Using these locations as initial vertices seemed to be problematic, resulting in very high number of triangles and general bad mesh patterns (even when using a relatively high cutoff).

 

Next, after reading Finn's recommendation in this post, I decided to try a uniform inner triangulation (not using the observation points). This reduced the number of triangles and I can produce what I think seems like an OK mesh. After first running the model with a relatively coarse mesh to get an idea of the spatial range, the final mesh$n is ~10,000. (with an estimated spatial range of 10 degrees and max.edge = ~2 degrees).

Then I ran my model using int.strategy = "eb", strategy = "gaussian"  and diagonal = 1000 to get some computationally cheap initial values to feed the proper model (using auto and simplified lapace) and every seems good, but then when I run the proper model, it seems to struggle. I get the 'stupid search' and my DIC is inf. This made me wonder that my mesh is too poor and it is causing numerical issues? 

 

So, worried that my mesh is too complex. I read this post, and tried defining an arbitrary nonconvex hull as my boundary, rather than the shapefile itself.  

 

## Non-convex hull domain boundary

nc.hull <- inla.nonconvex.hull(

            coords.est,

            convex = 3, 

            concave = 3,

            resolution = 150)

## Mesh.nc

mesh.nc <- inla.mesh.2d(

  boundary = nc.hull,

  cutoff = max.edge.length/5,  # make max.edge/5 for full resolution 

  max.edge = c(1, 3) * max.edge.length,

  # offset = 1, #c(max.edge.length, outer.boundary),

  min.angle = 20  

)     

 

This results in a mesh$n of ~5000 and a much more uniform 'good' look. But I am worried that altering the geographic structure makes the barrier method somewhat pointless?

 

 

I would really appreciate some guidance or insight on the best way to approach this mesh building. I would like to be able to find the good trade-off between mesh complexity and computational expense! Apologies if I have made any glaring errors or oversights...

and apologies for rambling... 

 

Thanks for taking the time to read this far :)

 

Cheers,

 

Mike

 

           

 

 

[Finn Lindgren]

Hi,

the main effect of the barrier model is to eliminate the Neumann
boundary effects at the edge of a mesh by allowing the parameters to
change outside a domain of interest instead of just extending the mesh
beyond the domain of interest (which can't be done when the domain
extension would intrude on the domain of interest, such as for
islands, peninsulas, etc). The number, sizes, and shapes of triangles
instead controls the numerical approximation properties of the
resulting discretised model, in much the same way for both
alternatives for handling the boundary, so they can largely be treated
as unrelated issues.

As long as the geometry used to define the barrier in the barrier
model is close enough to your real geometry you should be fine. It is
relatively safe to move a _convex_ boundary outwards for the barrier
model, but it's
concave parts (like islands) where you would get noticeable
differences when over-simplifying the boundary/barrier.

I'd recommend taking a look at the meshbuilder() tool for
interactively exploring mesh quality in relation to computational
approximation properties.
The important parts can also be accessed directly for any mesh,
through the inla.mesh.assessment() function, that evaluates the
pointwise standard deviations
for a specified correlation range in a basic SPDE/GMRF/Matern model.
A numerically accurate model would have a constant variance throughout
the domain of interest.
The barrier models are a bit more forgiving at the barriers, but the
tool should still give a good idea of the required triangle sizes.
The usually appropriate quantity to plot is "sd.dev", which is an
estimate of the ratio between the pointwise variances of the discrete
model and the desired continuous domain model.
Currently, the function returns "sd", "sd.dev", and "edge.len", which
is the average triangle edge length near each location (mostly just
informational).

The example below uses the ggplot2 features from inlabru (see
http://github.com/fbachl/inlabru for installation instructions; either
CRAN or the more bug-fixed development version) to illustrate this.
(The interactive meshbuilder tool has more options, but doesn't handle
spherical meshes; inla.mesh.assessment also doesn't fully handle that
case, but will currently just show you one side of the sphere...)

bnd <- inla.mesh.segment(cbind(c(0, 10, 10, 0, 0),
c(0, 0, 10, 10, 0)), bnd = TRUE)
mesh <- inla.mesh.2d(boundary = bnd, max.edge = 1)
out <- inla.mesh.assessment(mesh, spatial.range = 3, alpha = 2)

ggplot2::ggplot() + inlabru::gg(out, mapping=aes_(col=~sd.dev))

To build a spherical mesh, you just need to convert the
longitude-latitude cutoff&max.edge parameters to radians (multiply by
pi/180), and supply
crs = inla.CRS("sphere")
to inla.mesh.2d(), which will automatically transform inputs with
proper crs info (so you'll need to specify a crs=inla.CRS("longlat")
when you construct the nonconvex hull).

One can in principle compute the same assessment directly for the
chosen barrier model, but the barrier models weren't around when the
meshbuilder assessment code was written, so it needs to go on the TODO
list.

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 post to this group, send email to
> r-inla-disc...@googlegroups.com.
> Visit this group at https://groups.google.com/group/r-inla-discussion-group.
> For more options, visit https://groups.google.com/d/optout.



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

[mike.bl...@gmail.com]

Hi Finn,

Thanks so much for your detailed response. That's really helpful.

I had been playing with the `inla.mesh.assessment` tool, but I wasn't exactly sure whether I should  be monitoring at `sd` or `sd.dev`. It's much clearer now, especially with the `inlabru` ggplot tools. Really nice, thanks!

My reasoning was that I was worried that my boundary was too wiggly and this may have been causing numerical issues related to the triangles shape and size in difficult areas. I thought a more simplified non-convex hull could fix this, but I was concerned that over simplification would modify the geographic boundaries and negate the benefits of the barrier model by allowing smoothing and borrowing of strength from areas that are actually ocean, or discontiguous land areas. But like you say, as long as the geometry of the hull is close enough to the actual geometry it should hopefully be ok. 

Using the `inla.mesh.assessment` tool both meshes (using either the shapefile as a boundary or the hull) look to be OK, they both exhibit constant variance throughout 

the domain of interest. I think I will go ahead with the hull mesh, because it is a bit less complex and has more uniform triangles, so will hopefully be more numerically stable/accurate.

Thanks again for your help. I'll be back with more questions!

Cheers,

Mike

> email to r-inla-discussion-group+unsub...@googlegroups.com.

[Nicola Criscuolo]

Hello Mike and Finn, and thank you both for your detailed posts.

I am also facing the challenge of building a mesh, and fitting an LGCP model with inlabru(), on a point process at the global scale. In my case, my points are on the land. More specifically, such points should be were there are human settlements, since I am trying to model the distribution of persons belonging to a specific working category.

So far, I've collected more points in Europe than in the other parts of the world. Below there is the code I used to build the mesh, where the study_area_shp, and all my raster covariates, have the following CRS: +proj=longlat +datum=WGS84 +no_defs:

mesh <- inla.mesh.2d(boundary = study_area_shp,
                                        crs = proj4string(obj = study_area_shp),
                                        offset = c(0.25,
                                                          1.5),
                                        cutoff = 0.001, 
                                        max.edge = c(0.6,
                                                                 2))

You can see the mesh built for this continent in the image attached. Before building the mesh, I chose a simplified version of the European shapefile in order to avoid having to many edges, especially along the coastlines, and I identified and got rid of the points falling outside of the "boundary" area, following Bakka's tutorial. In addition, for efficiency, I created bigger triangles in the ocean and smaller on the land. Now, my doubts:

- Number of triangles. For an LGCP, what should be the correct size of triangles? With a very fine mesh, the computational time increases a lot.

- Measure units. I have created this mesh with numbers that "work" with my projection system, but I don't get the measure units. For example, the first value of the "max.edge" argument is 0.6, but this 0.6 refers to 0.6° since I have a "longlat" projection system?

- Preliminary evaluation. What is, in your opinion, the best way to evaluate that a mesh at the global level is "OK" before fitting a model on it? Should I first run a separate model with different parameters or using directly the inla.mesh.assessment() function (and, if yes, with which range?).

Every advice is very appreciated, I can't find many posts of people working at the global scale and I feel this really increases the complexity of several aspects of the analysis. Thank you really a lot in advance for your help!

Nico