Thx @bgoodri, as.matrix(), as.array(), or as.data.frame() will only be one of many intermediate steps to achieve an output consistent with lme4 or brms, though.
Here is a workaround. Workaround because it depends on several packages. Ideally, there will be an rstanarm::ranef(..., summary = FALSE) implementation for consistency with lme4 and brms. For package developers this consistency is important to make a package available for many interfaces. I tested this with cross-classified and nested models with several varying intercepts, including varying slopes and interactions, so it should be fairly general. May still require some adjustment for specific cases, though.
library(magrittr)
library(dplyr)
library(tidyr)
library(stringr)
library(purrr)
library(tidybayes)
get_rstanarm_sims <- function(model) {
constcomp_sims <- model %>%
magrittr::use_series(coefficients) %>%
names %>%
magrittr::extract(!str_detect(., "^b\\[")) %>%
as.matrix(model, pars = .)
varcomp_terms <- model %>%
rstanarm::VarCorr() %>%
purrr::map(dimnames) %>%
purrr::map(1)
varcomp_sims <- model %>%
tidybayes::spread_draws(b[term, groups], sep = " ") %>%
dplyr:: mutate(group = stringr::str_extract(groups, paste0(str_c("^", names(varcomp_terms)),
collapse = "|")),
level = stringr::str_remove(groups, paste0(str_c("^", names(varcomp_terms), ":"),
collapse = "|"))) %>%
dplyr::group_by(group) %>%
dplyr::group_split() %>%
purrr::imap(~ arrange(., match(.x$term, varcomp_terms[[.y]]))) %>%
purrr::map(~ array(data = .x$b,
dim = c(max(.x$.draw),
length(unique(.x$level)),
length(unique(.x$term))),
dimnames = list(NULL,
unique(.x$level),
unique(.x$term)))) %>%
magrittr::set_names(names(rstanarm::ngrps(model)))
list(constcomp = constcomp_sims,
varcomp = varcomp_sims)
}
sims <- get_rstanarm_sims(model = rstanarm_model)
str(sims)
This returns output consistent with arm:sim():
List of 2
$ constcomp: num [1:50, 1:6] -1.78 -1.75 -1.64 -1.49 -1.56 ...
..- attr(*, "dimnames")=List of 2
.. ..$ iterations: NULL
.. ..$ parameters: chr [1:6] "(Intercept)" "party2" "party3" "sex2" ...
$ varcomp :List of 2
..$ age : num [1:50, 1:4, 1:6] 0.442 0.372 0.438 0.799 0.191 ...
.. ..- attr(*, "dimnames")=List of 3
.. .. ..$ : NULL
.. .. ..$ : chr [1:4] "age1" "age2" "age3" "age4"
.. .. ..$ : chr [1:6] "(Intercept)" "party2" "party3" "sex2" ...
..$ race: num [1:50, 1:4, 1] -0.035 0.1601 0.2817 0.0299 0.2457 ...
.. ..- attr(*, "dimnames")=List of 3
.. .. ..$ : NULL
.. .. ..$ : chr [1:4] "race1" "race2" "race3" "race4"
.. .. ..$ : chr "(Intercept)"
whereby constcomp is fixef and varcomp is ranef in lme4.
I noticed that brms::ranef(..., summary = FALSE) is also only partially consistent with arm::sim(). For nested models, e.g., something like state/county, brms returns the parameters for these groups always at the end of the list and names them differently. The ordering is no problem if we index by name, however. For instance, lme4 would expand ``state/county` to two varying components, “state” and “county:state”, but brms returns “state” and “state:county”.
Using this code to get the draws in an almost arm::sim() consistent output from a brms model:
simulations <- list(constcomp = brms::fixef(model, summary = FALSE),
varcomp = brms::ranef(model, summary = FALSE))
str(simulations)
List of 2
$ constcomp: num [1:50, 1] 0.0171 -0.0114 -0.0664 -0.0912 -0.1037 ...
..- attr(*, "dimnames")=List of 2
.. ..$ iterations: NULL
.. ..$ parameters: chr "Intercept"
$ varcomp :List of 2
..$ state : num [1:50, 1:15, 1] -0.0881 -0.0372 -0.05 -0.0797 -0.0706 ...
.. ..- attr(*, "dimnames")=List of 3
.. .. ..$ : NULL
.. .. ..$ : chr [1:15] " DE – Baden-Württemberg" " DE – Bayern" " DE – Berlin" " DE – Brandenburg" ...
.. .. ..$ : chr "Intercept"
..$ state:county: num [1:50, 1:152, 1] -0.63 -0.72 -0.982 -0.673 -0.61 ...
.. ..- attr(*, "dimnames")=List of 3
.. .. ..$ : NULL
.. .. ..$ : chr [1:152] " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Breisgau-Hochschwarzwald" " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Konstanz" " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Rottweil" " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Schwarzwald-Baar-Kreis" ...
.. .. ..$ : chr "Intercept"
This code can be used to get the names right:
names_new <- simulations %>%
magrittr::use_series(varcomp) %>%
names %>%
stringr::str_split(":") %>%
purrr::map(rev) %>%
purrr::map(paste0, collapse = ":")
names(simulations$varcomp) <- names_new
str(simulations)
List of 2
$ constcomp: num [1:50, 1] 0.0171 -0.0114 -0.0664 -0.0912 -0.1037 ...
..- attr(*, "dimnames")=List of 2
.. ..$ iterations: NULL
.. ..$ parameters: chr "Intercept"
$ varcomp :List of 2
..$ state : num [1:50, 1:15, 1] -0.0881 -0.0372 -0.05 -0.0797 -0.0706 ...
.. ..- attr(*, "dimnames")=List of 3
.. .. ..$ : NULL
.. .. ..$ : chr [1:15] " DE – Baden-Württemberg" " DE – Bayern" " DE – Berlin" " DE – Brandenburg" ...
.. .. ..$ : chr "Intercept"
..$ county:state: num [1:50, 1:152, 1] -0.63 -0.72 -0.982 -0.673 -0.61 ...
.. ..- attr(*, "dimnames")=List of 3
.. .. ..$ : NULL
.. .. ..$ : chr [1:152] " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Breisgau-Hochschwarzwald" " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Konstanz" " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Rottweil" " DE – Baden-Württemberg_DE – Baden-Württemberg – Freiburg – Schwarzwald-Baar-Kreis" ...
.. .. ..$ : chr "Intercept"