pbs-assess / sdmTMB

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

Add names of random effects in tidy() #337

Open seananderson opened 2 months ago

seananderson commented 2 months ago
         Is there a way to extract the names of the random effects (e.g. right now, they just return Sigma_G, and if you have multiple random effects, you don't know which one is which)?

Originally posted by @sschooler in https://github.com/pbs-assess/sdmTMB/issues/41#issuecomment-2071050294

seananderson commented 2 months ago

@sschooler This is now a separate issue and we'll get to this eventually. At the same time, we should add better names for the ranges (@ecophilina). In the meantime, here's an example:

library(sdmTMB)
x <- runif(500, -1, 1)
y <- runif(500, -1, 1)
loc <- data.frame(x = x, y = y)
mesh <- make_mesh(loc, c("x", "y"), n_knots = 50, type = "kmeans")

s <- sdmTMB_simulate(
  ~1,
  data = loc,
  mesh = mesh,
  range = 1.4,
  phi = 0.1,
  sigma_O = 0.2,
  seed = 1,
  B = 0
)

g <- rep(gl(30, 10), 999)
set.seed(134)
RE_vals <- rnorm(30, 0, 0.4)
h <- rep(gl(40, 10), 999)
set.seed(1283)
RE_vals2 <- rnorm(40, 0, 0.2)
s$g <- g[seq_len(nrow(s))]
s$h <- h[seq_len(nrow(s))]
s$observed <- s$observed + RE_vals[s$g] + RE_vals2[s$h]

fit <- sdmTMB(
  data = s, time = NULL,
  formula = observed ~ 1 + (1 | g) + (1 | h), mesh = mesh
)

# you can see them here with names to match up:
fit
#> Spatial model fit by ML ['sdmTMB']
#> Formula: observed ~ 1 + (1 | g) + (1 | h)
#> Mesh: mesh (isotropic covariance)
#> Data: s
#> Family: gaussian(link = 'identity')
#>  
#>             coef.est coef.se
#> (Intercept)     0.22    0.15
#> 
#> Random intercepts:
#>   Std. Dev.
#> g      0.44
#> h      0.17
#> 
#> Dispersion parameter: 0.10
#> Matérn range: 1.25
#> Spatial SD: 0.18
#> ML criterion at convergence: -253.583
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

# yes, it would be helpful if these were named better:
tidy(fit, 'ran_pars')
#> # A tibble: 5 × 3
#>   term    estimate std.error
#>   <chr>      <dbl>     <dbl>
#> 1 range      1.25    0.364  
#> 2 phi        0.103   0.00369
#> 3 sigma_O    0.180   0.0349 
#> 4 sigma_G    0.436   0.0596 
#> 5 sigma_G    0.175   0.0275

# you can currently grab them from:
fit$split_formula[[1]]$barnames  # this internal structure may someday change
#> [1] "g" "h"

# or access just the relevant printing function (which uses the above):
sdmTMB:::print_iid_re(fit)
#>   Std. Dev.
#> g      0.44
#> h      0.17

# here (the values themselves) they are named:
tidy(fit, 'ran_vals')
#> # A tibble: 70 × 3
#>    term  estimate std.error
#>    <chr>    <dbl>     <dbl>
#>  1 g_1    -0.295      0.126
#>  2 g_2     0.691      0.126
#>  3 g_3     0.0805     0.126
#>  4 g_4    -0.307      0.126
#>  5 g_5    -0.604      0.126
#>  6 g_6     0.0498     0.126
#>  7 g_7     0.556      0.126
#>  8 g_8    -0.0244     0.126
#>  9 g_9     0.0985     0.126
#> 10 g_10    0.660      0.126
#> # ℹ 60 more rows

Created on 2024-04-22 with reprex v2.1.0