Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add names of random effects in tidy() #337

Open
seananderson opened this issue Apr 22, 2024 · 1 comment
Open

Add names of random effects in tidy() #337

seananderson opened this issue Apr 22, 2024 · 1 comment

Comments

@seananderson
Copy link
Member

         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 #41 (comment)

@seananderson
Copy link
Member Author

@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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant