Add regions / sub-areas to get_index?

[ericward-noaa]

Question from a user: what's a good way to generate regional indices in a survey domain? E.g. in addition to generating a total estimate of density with CIs, could we do the same if we were able to map each grid cell to a survey area ("North", "Central", "South") or similar. Also referenced over here:

#166

[seananderson]

We could possibly add something in get_index(), but you'd also have to iterate over the predictions first, and by then it's fairly simple to iterate over the index calculations. Unless we want to make something that wraps both the prediction and index calculations.

Below is an example. You could use lapply() and do.call("rbind", ...) if you wanted to skip purrr.

library(sdmTMB)
library(ggplot2)

fit <- sdmTMB(
  density ~ as.factor(year),
  data = pcod_2011,
  mesh = pcod_mesh_2011, 
  time = "year", 
  spatiotemporal = "off", 
  family = tweedie()
)
fit
#> Spatial model fit by ML ['sdmTMB']
#> Formula: density ~ as.factor(year)
#> Mesh: pcod_mesh_2011 (isotropic covariance)
#> Time column: year
#> Data: pcod_2011
#> Family: tweedie(link = 'log')
#>  
#>                     coef.est coef.se
#> (Intercept)             2.96    0.44
#> as.factor(year)2013     0.15    0.18
#> as.factor(year)2015     0.32    0.19
#> as.factor(year)2017    -0.36    0.20
#> 
#> Dispersion parameter: 15.44
#> Tweedie p: 1.59
#> Matern range: 32.64
#> Spatial SD: 2.02
#> ML criterion at convergence: 3017.303
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

# assuming a very recent version of sdmTMB:
grids <- sdmTMB::replicate_df(qcs_grid, "year", unique(pcod_2011$year))
grids$region <- ifelse(grids$Y > 5725, "north", "south")

preds <- grids |>
  split(grids$region) |> 
  lapply(function(x) predict(fit, newdata = x, return_tmb_object = TRUE))

inds <- purrr::map_dfr(preds, get_index, area = 4, .id = "region")

head(inds)
#>   region year      est      lwr       upr  log_est        se
#> 1  north 2011 654038.2 489383.8  874091.0 13.39092 0.1479715
#> 2  north 2013 762516.6 560039.1 1038198.1 13.54438 0.1574608
#> 3  north 2015 905095.2 675079.8 1213482.4 13.71580 0.1495993
#> 4  north 2017 456086.0 329941.0  630459.5 13.03044 0.1651905
#> 5  south 2011 349478.2 257196.8  474869.7 12.76420 0.1564312
#> 6  south 2013 407442.4 301364.4  550859.0 12.91765 0.1538698

ggplot(inds, aes(year, est, ymin = lwr, ymax = upr, fill = region)) +
  geom_ribbon(alpha = 0.3) +
  geom_line(aes(colour = region))

^{Created on 2023-01-23 with reprex v2.0.2}

[ericward-noaa]

Thanks, this seems like a fine solution. I think we can close this -- and we can come back to it if we wanted to wrap predictions / indices together