Calculating index values for subsets of the domain

[Cole-Monnahan-NOAA]

Say I've got three areas, A, B, and C which are essentially non-overlapping subsets of the large spatial domain. Each observation can be categorized as in an area.

Is it possible to get index values and SEs for each area, apart from the total?

I'm thinking of a couple possible ways to do this. First is to generating a polygon in R and generating a ton of points and then integrate it by pushing it through as newdata. Or to create separate meshes for each and somehow get those to work internally to sdmTMB.

Any quick thoughts on whether this is easy to do or not? I could also make a reprex is that's helpful.

[seananderson]

If you want an index for one or more subsets of the full survey domain (which you've fit your model to), you can just subset your prediction grid (or triangles or whatever assuming you have their areas) to match your domain subsets. I would overlay the polygon on the prediction grid to find those within the grid (and possibly bother to deal with partial areas in the outer grid cells) and pass that subsetted grid to the newdata argument in predict. The same fitted model can be used for all subsets. Is this what you're looking to do or am I misunderstanding?

Also, nice to see you here!

[Cole-Monnahan-NOAA]

Hi Sean,

Yeah it's nice to be here. In full disclosure I steered a student here (away from VAST) and haven't dug properly yet. But he's got it working and making progress which is great.

Yes I understand the response and I think that'll work. The only part I don't know how to do is pull out the prediction grid filtering. Are you saying put the polygon over the prediction mesh? Or do you mean the points that are generated and used for prediction (if this is even what you do?).

[Lewis-Barnett-NOAA]

The idea is to simply remove the locations from the prediction grid/extrapolation grid that you want to exclude. You could do this using spatial objects, e.g., using a polygon defining the boundary of the subset of the domain to clip/crop a grid or raster you may have made for the entire domain. Or, if the boundaries are simple, like all areas west of 140 degrees longitude, you could just do a simple filter operation on a data frame where each row is a location representing the centroid of the grid. Either way, you just specify the resulting object to be the newdata.

[Cole-Monnahan-NOAA]

Right, but how do I extract the extrapolation grid? I couldn't find that in the documentation.

[Lewis-Barnett-NOAA]

There is no prediction/extrapolation grid created internally by the package at this time. You have to do some external steps to create this, which then defines both the resolution and extent of predictions.

If the fitted model includes covariates, values of those will also need to be in this prediction grid.

[seananderson]

See this example I just added to the bottom of the index vignette: https://github.com/pbs-assess/sdmTMB/blob/2edaf9891987fcae9435064e0b46e0d923c4fb24/vignettes/index-standardization.Rmd#L138-L156

The extrapolation grid is whatever you supply to newdata in predict.sdmTMB(). In this case the extrapolation grid is based on internal package data qcs_grid, but any data frame will do.

[seananderson]

There's an example of using sf::st_intersects() internally in the barrier mesh function: https://github.com/pbs-assess/sdmTMB/blob/2edaf9891987fcae9435064e0b46e0d923c4fb24/R/mesh.R#L419-L423

That's for intersecting a barrier on the mesh. Here you'd be intersecting the survey domain on the prediction grid. The prediction grid could be a triangularized mesh if you want, as long as you get the centroids and their area, but from some experimentation I did a while ago, it didn't seem to make much of a difference.

[Cole-Monnahan-NOAA]

Thanks Sean, I understand now. I came at it from a VAST perspective but the predict method really simplifies and makes it flexible. I like that. We'll try out this method on our analysis and let you know if we have any further questions or issues. Thanks for the quick help!