mesh size for continental-scale data?

[jasonwi...@u.boisestate.edu]

 I'm looking at environmental predictors of nest initiation date in a raptor species that breeds across North America. I have a dataset of >1500 nests in the US and Canada, but the observations are pretty clumpy (much of the data comes from local nest box studies where boxes are placed ~0.5 km from each other).

I created a leaflet map of the nest sites here so you have a better idea of the distribution of nests: http://rpubs.com/jmwiniarski/nestmap2

 Is this distribution (over a huge area, but clumped together tightly in certain places) ok when it comes to creating a mesh? Does anyone have recommendations for building a mesh in this situation? I've attached a copy of my nest locations if that is helpful.

Thanks,

Jay

[Finn Lindgren]

Hi,

If I understand correctly you are modelling a quantity that is measured at the nest locations (the initiation time), and not modelling the nest locations themselves.

In this situation, I would recommend one of two approaches:

1) don’t use the nest locations for building the mesh, but just a boundary covering the whole domain, with a uniform triangle size

2) using the nest locations to get adaptive triangle sizes, but with a “cutoff” parameter large enough avoid small clusters.

(In both cases, use a second outer boundary layer with larger triangles to move spde edge effects)

The new INLA::meshbuilder() tool can be used to assess the approximation properties of meshes in a “shiny” way. A separate inla.mesh.assessment() function is also available for command line use. (In the development INLA at least, not sure if it was added before or after the latest testing version)

See the 03_meshes and 03_mesh_assessment practicals linked here for more info:

https://sites.google.com/inlabru.org/inlabru/workshops/26-30-march-2018-st-andrews

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.

<nests.csv>

[Finn Lindgren]

PS. The meshassessment.R file in the same folder as the practical actually contains a version of the assessment code that is now in the INLA package.

See the 03_meshes slides for some more info on mesh assessment. The key is to make the “sd deviation” flat; if triangles are too big compared to the spatial correlation range, the bad linear approximation can be seen by the std.dev. inside the triangle being clearly lower than at the vertices. In a “good” mesh, the std.dev. is (almost) constant.

Finn

On 13 Apr 2018, at 23:16, jasonwi...@u.boisestate.edu wrote:

--

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.

<nests.csv>

[jasonwi...@u.boisestate.edu]

Hi Finn,

Thanks for those resources, I will take a look at the mesh assessment function.

The response variable is the day of the year that a nest is initiated. The data comes from 2012-2017, so there are sometimes multiple nest records for the same nest box location. I would like to look at the environmental predictors of laydate (e.g., temperature, precipitation, vegetation green-up, aritificial-light-at-night) that I extracted for each nest record from its nest box location. 

Would you still recommend the approach you describe above? If so, could I just use the border of the continental US or an arbitrary polygon that encapsulates all of the nest box locations? Or use inla.nonconvex.hull to create a buffer around all of the points?

Thanks,

Jay 

On Friday, April 13, 2018 at 4:35:15 PM UTC-6, Finn Lindgren wrote:

PS. The meshassessment.R file in the same folder as the practical actually contains a version of the assessment code that is now in the INLA package.

See the 03_meshes slides for some more info on mesh assessment. The key is to make the “sd deviation” flat; if triangles are too big compared to the spatial correlation range, the bad linear approximation can be seen by the std.dev. inside the triangle being clearly lower than at the vertices. In a “good” mesh, the std.dev. is (almost) constant.

Finn

On 13 Apr 2018, at 23:16, jasonwi...@u.boisestate.edu wrote:

 I'm looking at environmental predictors of nest initiation date in a raptor species that breeds across North America. I have a dataset of >1500 nests in the US and Canada, but the observations are pretty clumpy (much of the data comes from local nest box studies where boxes are placed ~0.5 km from each other).

I created a leaflet map of the nest sites here so you have a better idea of the distribution of nests: http://rpubs.com/jmwiniarski/nestmap2

 Is this distribution (over a huge area, but clumped together tightly in certain places) ok when it comes to creating a mesh? Does anyone have recommendations for building a mesh in this situation? I've attached a copy of my nest locations if that is helpful.

Thanks,

Jay

--
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-group+unsub...@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.

<nests.csv>

[Finn Lindgren]

On 13 Apr 2018, at 23:50, jasonwi...@u.boisestate.edu wrote:

Would you still recommend the approach you describe above?

Yes! You’re describing a common type of model.

If so, could I just use the border of the continental US or an arbitrary polygon that encapsulates all of the nest box locations? Or use inla.nonconvex.hull to create a buffer around all of the points?

I’d use inla.nonconvex.hull, and ignore the ocean when looking the the results.

_If_ you want the geography to be part of the mesh, please don’t use the US border; I’m pretty sure the nests don’t really care about politics enough to have boundary effects there. ;-)

Finn

Thanks,

Jay 

On Friday, April 13, 2018 at 4:35:15 PM UTC-6, Finn Lindgren wrote:

PS. The meshassessment.R file in the same folder as the practical actually contains a version of the assessment code that is now in the INLA package.

See the 03_meshes slides for some more info on mesh assessment. The key is to make the “sd deviation” flat; if triangles are too big compared to the spatial correlation range, the bad linear approximation can be seen by the std.dev. inside the triangle being clearly lower than at the vertices. In a “good” mesh, the std.dev. is (almost) constant.

Finn

On 13 Apr 2018, at 23:16, jasonwi...@u.boisestate.edu wrote:

 I'm looking at environmental predictors of nest initiation date in a raptor species that breeds across North America. I have a dataset of >1500 nests in the US and Canada, but the observations are pretty clumpy (much of the data comes from local nest box studies where boxes are placed ~0.5 km from each other).

I created a leaflet map of the nest sites here so you have a better idea of the distribution of nests: http://rpubs.com/jmwiniarski/nestmap2

 Is this distribution (over a huge area, but clumped together tightly in certain places) ok when it comes to creating a mesh? Does anyone have recommendations for building a mesh in this situation? I've attached a copy of my nest locations if that is helpful.

Thanks,

Jay

--
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-group+unsub...@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.

<nests.csv>

--
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.

[Vaughn Bodden]

Should the outer mesh also be relatively constant? 

On Friday, April 13, 2018 at 11:35:15 PM UTC+1, Finn Lindgren wrote:

PS. The meshassessment.R file in the same folder as the practical actually contains a version of the assessment code that is now in the INLA package.

See the 03_meshes slides for some more info on mesh assessment. The key is to make the “sd deviation” flat; if triangles are too big compared to the spatial correlation range, the bad linear approximation can be seen by the std.dev. inside the triangle being clearly lower than at the vertices. In a “good” mesh, the std.dev. is (almost) constant.

Finn

On 13 Apr 2018, at 23:16, jasonwi...@u.boisestate.edu wrote:

 I'm looking at environmental predictors of nest initiation date in a raptor species that breeds across North America. I have a dataset of >1500 nests in the US and Canada, but the observations are pretty clumpy (much of the data comes from local nest box studies where boxes are placed ~0.5 km from each other).

I created a leaflet map of the nest sites here so you have a better idea of the distribution of nests: http://rpubs.com/jmwiniarski/nestmap2

 Is this distribution (over a huge area, but clumped together tightly in certain places) ok when it comes to creating a mesh? Does anyone have recommendations for building a mesh in this situation? I've attached a copy of my nest locations if that is helpful.

Thanks,

Jay

--
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-group+unsub...@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.

<nests.csv>

[Finn Lindgren]

No, the standard deviation in the extension doesn’t need to be constant; mesh extensions are used to move the non-stationary boundary effect away from the domain of interest.

Finn

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/bdb1a9b4-4109-4a2f-a20d-d594735eb8c5%40googlegroups.com.

[Jason Winiarski]

Hi Finn,

I am finally returning to this analysis again, and in the meantime I have received more data points across a broader geographic area. Nest box locations are scatted across the continental US, as well as some locations in southern Canada and Alaska. Nest boxes in this data set are as close together as several meters, and the farthest distance between boxes is ~6000 km:

I've started the process of creating a non-convex hull to make the mesh:

library(tidyverse)

library(sf)

library(INLA)

library(inlabru)

library(rnaturalearth)

theme_set(theme_minimal())

# for plotting and converted to a meter projection

n_america <- ne_countries(country = c('united states of america', 'canada'), returnclass = 'sf')

north_american_albers_equal_area_conic_projection <- '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs' 

# read in data

laydates_sf <- read_rds('/Users/Jay/Desktop/laydates_sf_example.rds')

# plot points over north america

laydates_sf %>%

  ggplot() +

  geom_sf(data = n_america) +

  geom_sf(pch = 1)

# convert to meter projection

laydates_sf <- laydates_sf %>%

  st_transform(., crs = north_american_albers_equal_area_conic_projection) %>%

  mutate(

    x = map_dbl(geometry, 1),

    y = map_dbl(geometry, 2)

  ) %>%

  as_tibble() %>%

  select(-geometry)

laydates_sf

# put coordinates on km scale

laydates_sf$x <- laydates_sf$x/1000

laydates_sf$y <- laydates_sf$y/1000

nest_locations <- cbind(laydates_sf$x, laydates_sf$y)

bound1 <- inla.nonconvex.hull(

  nest_locations,

  convex = 250

  )

mesh1 <- inla.mesh.2d(

  boundary = bound1,

  max.edge = c(250, 600),

  cutoff = 100

)

ggplot() + 

  gg(mesh1) + 

  geom_point(data = laydates_sf, aes(x, y)) +

  coord_equal() +

  labs(

    x = 'x (km)',

    y = 'y (km)'

  )

Which looks like this:

But I'm still a bit unclear on a few things for the mesh construction. First, in this situation where I have some disjunct data points (in Alaska) is it best to keep all of the data points within the same boundary, or is it ok for those points to have their own boundary (like the plot above)? Second, is there a good rule of thumb for choosing starting values for inner and outer max.edge and cut-off?

Thanks,

Jay

[Finn Lindgren]

Hi Jason,

disjoint regions of finer resolution, connected by regions of coarse resolution, as in your example, is perfectly ok. In your case, I would construct not only the inner boundary as a non-convex hull, but also the outer boundary; that will allow the connection between Alaska and the rest to use fewer triangles but still remain connected.

Outer max.edge can usually be set quite large, since the min.angle criterion will typically force it to add triangles anyway; Increasing min.angle from the default 21 is usually a good idea (26-27 is usually ok). Just make sure to set the max.n.strict and max.n parameters (see ?inla.mesh.2d) to prevent an attempt at making infinitely many triangles. (You can set it to some values much larger than you want your mesh to use, just to prevent an infinite loop, and adjust min.angle until the mesh is reasonably uniform)

Since you're not imposing any physical boundaries on the domain, you can set cutoff to a value on the order of the inner max.edge, e.g. max.edge/2.

For good approximation accuracy, max.edge should be a fracetion of the correlation range. range/10 is a safe value, but often too strict.

Use the INLA::meshbuilder() shiny app to experiment with the approximation accuracy visualisation tool to get a feeling for how it works.

The extension should be large enough to move the boundary effects away from the data, approx the correlation range, or 2*range to be on the safe side.

Finn

Thanks,

Jay 

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.

<nests.csv>

--
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.

--

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/7b9b0a23-7a35-43d8-8750-1d9a278cfe6d%40googlegroups.com.

--

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

[Jason Winiarski]

Hi Finn, these are really helpful starting points for me. For also constructing the outer boundary as a non-convex hull, would that be manipulated using the inla.mesh.2d() function? I guess I'm not sure which argument would be used for treating the outer boundary as a non-convex hull.

Thanks,

Jay

Thanks,

Jay 

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsub...@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.

<nests.csv>

--
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-group+unsub...@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.

--
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-group+unsub...@googlegroups.com.

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/7b9b0a23-7a35-43d8-8750-1d9a278cfe6d%40googlegroups.com.

[Finn Lindgren]

Hi Jason,

I was thinking something like this:

bound1 <- inla.nonconvex.hull(

  nest_locations,

  convex = 250

  )

bound2 <- inla.nonconvex.hull(

  nest_locations,

  convex = 1000

  )

mesh1 <- inla.mesh.2d(

  boundary = list(bound1, bound2),

  max.edge = c(250, 600),

  cutoff = 100

)

Thanks,

Jay 

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.

<nests.csv>

--
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.

--
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/7b9b0a23-7a35-43d8-8750-1d9a278cfe6d%40googlegroups.com.

--

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

--

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/0a92a9f8-f40d-4b6c-9888-cb4ef97f6231%40googlegroups.com.