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

Improve reporting of range with anisotropic models #149

Closed seananderson closed 1 year ago

seananderson commented 1 year ago

Currently naively reports the range without accounting for the rotation matrix.

Could report the maximum and minimum range in any direction.

Could not report the range and include a note to see plot_anisotropy() (the new one)

seananderson commented 1 year ago

Example of where to grab the range values from to get started:

library(sdmTMB)
mesh <- make_mesh(pcod_2011, c("X", "Y"), n_knots = 80, type = "kmeans")
fit <- sdmTMB(
  data = pcod_2011,
  formula = density ~ 1,
  mesh = mesh,
  family = tweedie(),
  share_range = FALSE,
  time = "year",
  anisotropy = TRUE
)

plot_anisotropy(fit)

d <- plot_anisotropy(fit, return_data = TRUE)
d
#>       angle        a         b      maj1       min1   model   random_field
#> 1 0.7408169 55.93917 14.934098 -41.27849  10.078887 tweedie        spatial
#> 2 0.7408169 15.82456  4.224688 -37.75284 -11.020133 tweedie spatiotemporal
#> 3 0.7408169 55.93917 14.934098 -11.67722   2.851204 tweedie        spatial
#> 4 0.7408169 15.82456  4.224688 -10.67985  -3.117472 tweedie spatiotemporal

# d$a and d$b are the ranges
# (I think always maximum and minimum?) for a given field

# example with ellipse plotting with a and b:

library(ggforce)
#> Loading required package: ggplot2

ggplot() +
  geom_ellipse(aes(x0 = 0, y0 = 0, a = 10, b = 3, angle = 0)) +
  coord_fixed(xlim = c(-10, 10), ylim = c(-10, 10))


ggplot() +
  geom_ellipse(aes(x0 = 0, y0 = 0, a = 10, b = 3, angle = 45)) +
  coord_fixed(xlim = c(-10, 10), ylim = c(-10, 10))

Created on 2023-01-27 with reprex v2.0.2

The print method needs to be updated here: https://github.com/pbs-assess/sdmTMB/blob/eff0011e8263ac7405acf2453e4f2c3efb63521e/R/print.R#L178-L207

I imagine the output would be something like 'Matern anisotropic range: 15 (maximum) and 3 (minimum) at 75°' but that's up for grabs too.

seananderson commented 1 year ago

Thanks @jdunic! It's merged into main now. This is a really nice addition.

library(sdmTMB)
fit <- sdmTMB(
  density ~ 1,
  data = pcod_2011,
  mesh = pcod_mesh_2011,
  family = tweedie(),
  share_range = FALSE,
  time = "year",
  anisotropy = TRUE
)
fit
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: density ~ 1
#> Mesh: pcod_mesh_2011 (anisotropic covariance)
#> Time column: year
#> Data: pcod_2011
#> Family: tweedie(link = 'log')
#>  
#>             coef.est coef.se
#> (Intercept)     2.79    0.46
#> 
#> Dispersion parameter: 14.76
#> Tweedie p: 1.56
#> Matérn anisotropic range (spatial): 20.11 to 65.78 at 116.98˚
#> Matérn anisotropic range (spatiotemporal): 0.01 to 0.04 at 116.98˚
#> Spatial SD: 1.99
#> Spatiotemporal SD: 1173.70
#> ML criterion at convergence: 3008.983
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
#> See ?plot_anisotropy to plot the anisotropic range.

Created on 2023-02-18 with reprex v2.0.2