Is it necessary to use the exact same data to make the mesh as to run the model?

[dmkaplan2000]

Hi,

I am a bit confused about the extent to which the mesh used to run a model should have been made with the exact same data used in the model. I thought that one could make the mesh with a given set of data and then run sdmTMB with that mesh, but using another set of data, but at some point I got an error message saying that the data does not seem to fit the mesh, making me wonder to what extent the data is linked with the mesh (the mesh object does save the coordinates of the data used to create it, so I am wondering if this information is used when running the model). For example, can you rearrange the rows in the data after having created the mesh? Can you add columns to the data frame? Can you filter out rows after creating the mesh?

The reason I am asking is that I first create basically the dataset I want to use for models, but then as a last step, depending on what exact model I plan on running, I may filter out some rows and potentially rearrange the order of the rows. I could make the mesh with this final final data just before running the models, but it would be kind of nice for all models to use the same mesh. More importantly, with the final filtered data, make_mesh sometimes produces a mesh that does not cover the original data, leading models to fail. For example, see the mesh below that for some reason I can't quite figure out does not cover the data:

In this case, if I make the mesh before filtering, the mesh covers the original data better (though it would probably be best for me to figure out why make_mesh isn't producing a mesh that covers all the data).

Thanks,
David

[dmkaplan2000]

FYI, I found a solution to my mesh not covering the data problem. I unioned the data and then buffered it by a reasonable amount and then supplied that to make_mesh as the boundary argument. This custom boundary forced the mesh to cover the data.

But I would still like to know if the data and the mesh need to be formally linked together or not.

[ericward-noaa]

Hi David -- yes, the data is directly linked to the mesh, so dropping / adding a row will mess things up without a new call to make_mesh(). One of the elements of the mesh is the projection matrix A, which maps observations on to mesh vertices (letting you predict from the smaller number of vertices to larger number of locations of observations). More on the projector matrix here: https://becarioprecario.bitbucket.io/spde-gitbook/ch-intro.html

So if the mesh is constructed on the data after your filtering, I think it should cover it just fine.

[dmkaplan2000]

Thanks, this explains things to me.

[seananderson]

You can recalculate the bilinear projection matrix by supplying the knots/vertices from a previously generated mesh. E.g.:

library(sdmTMB)

mesh1 <- make_mesh(pcod, c("X", "Y"), cutoff = 10)
df2 <- pcod[1:100,]
mesh2 <- make_mesh(df2, c("X", "Y"), mesh = mesh1$mesh)

plot(mesh1)

plot(mesh2)

^{Created on 2024-10-10 with reprex v2.1.1}

Regarding points being outside the mesh. I too have seen this sometimes. Unfortunately, this is caused by the use of the fairly simple fmesher::fm_rcdt_2d_inla() internally and is out of my control (it's under the control of fmesher). In retrospect, I wish I had used what is now fmesher::fm_mesh_2d_inla(), although you can easily swap that out in the fmesher_func argument.

There's an example of that in the latest preprint version.

You can also generate your own mesh using the fmesher (or INLA) functions following any of the tutorials you'll find online and supply that to the mesh argument. In retrospect, it might have been a better idea to have the user due that, supply that to sdmTMB() and skip the make_mesh() function all together. I was trying to shield the user from the complexities of mesh creation for simple use cases.

Example:

library(sdmTMB)
inla_mesh <- fmesher::fm_mesh_2d_inla(
  loc = cbind(pcod$X, pcod$Y), # coordinates
  max.edge = c(25, 50), # max triangle edge length; inner and outer meshes
  offset = c(5, 25),  # inner and outer border widths
  cutoff = 5 # minimum triangle edge length
)
mesh <- make_mesh(pcod, c("X", "Y"), mesh = inla_mesh)
plot(mesh)

^{Created on 2024-10-10 with reprex v2.1.1}

[dmkaplan2000]

Thanks, these are all good suggestions. As I said above, the boundary argument seems to have fixed things a bit, but recalculating the projection matrix before fitting was the step I was missing. I will do this.

Overall, I think you did a good thing in trying to shield the user from the complexities of mesh generation. It is just that for me the links between the mesh and the data weren't 100% clear (perhaps my fault for not reading enough of the documentation). One thing is that my models using a mesh disconnected from the data generally ran without returning an error that clearly indicated that something was off with the mesh (they generally failed due to singularities, which may or may not be related to the mesh). One time I did get a "mesh is bad" error, but I think this is because I randomly used a mesh made with data that had fewer lines than the data supplied to the model, whereas in all other cases I was using data in the model that had fewer lines than the data used to create the mesh, which did not seem to generate a bad mesh error.

Perhaps to assure that this doesn't happen, additional checks on the match between the model and the data should be performed before running the model? Perhaps make_mesh could include a hash on something like paste0(x,y,sep=",",collapse=";") that would be check on input data in sdmTMB?

[dmkaplan2000]

A related question - is it bad if the mesh is much bigger than the data area? To assure that the mesh covers the data, I am putting a good sized cushion around the edge of the data. Is this a problem?

[seananderson]

is it bad if the mesh is much bigger than the data area?

You want the mesh to be a bit bigger than your data to avoid edge effects. At least the size of the correlation range is often suggested and maybe more than that. The only problem with going bigger is the additional random effects (every vertex is a random effect) meaning slower model fitting.

[dmkaplan2000]

Thanks for the help!