pbs-assess / sdmTMB

:earth_americas: An R package for spatial and spatiotemporal GLMMs with TMB
https://pbs-assess.github.io/sdmTMB/
184 stars 26 forks source link

Add ability to predict using stars data objects #103

Closed plthompson closed 1 year ago

plthompson commented 2 years ago

@seananderson, adding this issue that we discussed. Looking to be able to make predictions using sdmTMB models based on spatiotemporal environmental data formatted using the stars package. This works for non-spatial models (e.g. lm, glm - https://r-spatial.github.io/stars/articles/stars7.html) but it doesn't know how to deal with the spatial coordinates specified in the mesh within sdmTMB. When I try predicting with a sdmTMB model it gives the following error:

Error in `predict()`:
! `xy_cols` (the column names for the x and y coordinates) are not in `newdata`.
Did you miss specifying the argument `xy_cols` to match your data?
The newer `make_mesh()` (vs. `make_spde()`) takes care of this for you.

The following code should allow you to reproduce this (I'll send you the data separately):

mesh <- make_mesh(data = Pcod, xy_cols = c("X", "Y"), cutoff = 15)

glm_m <- glm(data = Pcod, formula = catch_kgpkm2 ~ poly(temperature, 2, raw = TRUE))
temperature_stars$cpue <- predict(glm_m, temperature_stars)
ggplot()+
  geom_stars(data = temperature_stars[2])+
  facet_wrap(~year(time))

sdm_m <- sdmTMB(data = Pcod, formula = catch_kgpkm2 ~ poly(temperature, 2, raw = TRUE), mesh = mesh, spatiotemporal = "IID", time = "time")
smd_prediction <- predict(sdm_m, temperature_stars)

If you do decide to make this possible, one issue will be how to save the output. You could put all of the prediction vectors into a stars object but this seems a bit clunky. It probably makes more sense to convert the stars data to a dataframe and output that with new columns for the model predictions.

plthompson commented 2 years ago

After working with the stars package a bit more, I don't think this is a priority. You can convert a stars object to an sf dataframe using st_as_sf(x, points = TRUE, long = TRUE). I guess if you wanted to, you could code this into the predict function for sdmTMB objects but it probably isn't necessary.