inclusion of brms and rstanarm as vetted models in `tbl_regression()`

[storopoli]

I've been using a lot of gtsummary for my frequentist work, but I miss a lot gtsummary's functionality for bayesian models.

brms(brmsfit object) and rstanarm (stanreg objects) models are already implemented in broom.mixed for tidy, glance, augment and several other methods. See broom.mixed vignette (section "Capabilities").

[ddsjoberg]

Please add code examples, and I'll see what I can do. Please test on the most recent version of gtsummary as well. I think I added support for some of them?

[larmarange]

Dear @storopoli

could you please also test such models with broom.helpers https://larmarange.github.io/broom.helpers/ and in particular to check that tidy_plus_plus() returns the expected results.

As tbl_regression() is now based on broom.helpers, the first thing to check is that broom.helpers works as expected with such models.

Regards

[storopoli]

Dear @ddsjoberg, yes you are corrected you have already implemented them from rstanarm. But they are missing from brms models.

library(brms)
library(rstanarm)

rstanarm

model_rstanarm <- stan_glm(mpg ~ hp + vs, data = mtcars)
gtsummary::tbl_regression(model_rstanarm)

It plots a gtsummary table object for the model

broom.helpers::tidy_plus_plus(model_rstanarm)
 # A tibble: 2 x 14
  term  variable var_label var_class var_type var_nlevels contrasts contrasts_type reference_row label estimate std.error
  <chr> <chr>    <chr>     <chr>     <chr>          <int> <chr>     <chr>          <lgl>         <chr>    <dbl>     <dbl>
1 hp    hp       hp        numeric   continu…          NA NA        NA             NA            hp     -0.0542    0.0145
2 vs    vs       vs        numeric   continu…          NA NA        NA             NA            vs      2.57      1.96  
# … with 2 more variables: conf.low <dbl>, conf.high <dbl>

brms

model_brms <- brm(mpg ~ hp + vs, data = mtcars)
gtsummary::tbl_regression(model_brms)
x There was an error calling `tidy_fun()`. Most likely, this is because
  the function supplied in `tidy_fun=` was misspelled, does not exist,
  is not compatible with your object, or was missing necessary arguments (e.g.
  `conf.level=` or `conf.int=`). See error message below.
Error: Error: No tidy method for objects of class brmsfit

broom.helpers::tidy_plus_plus(model_brms)
x There was an error calling `tidy_fun()`. Most likely, this is because
  the function supplied in `tidy_fun=` was misspelled, does not exist,
  is not compatible with your object, or was missing necessary arguments (e.g.
  `conf.level=` or `conf.int=`). See error message below.
Error: Error: The model does not contain group-level effects.

[ddsjoberg]

@storopoli can you try adding tidy_fun = broom.mixed::tidy to both tbl_regression and broom.helpers? (for some reason the brms models are not compiling on my machine)

[storopoli]

Here it is

> broom.helpers::tidy_plus_plus(model_brms, tidy_fun = broom.mixed::tidy)
Registered S3 method overwritten by 'broom.mixed':
  method      from 
  tidy.gamlss broom
x There was an error calling `tidy_fun()`. Most likely, this is because
  the function supplied in `tidy_fun=` was misspelled, does not exist,
  is not compatible with your object, or was missing necessary arguments (e.g.
  `conf.level=` or `conf.int=`). See error message below.
Error: Error: The model does not contain group-level effects.
> 

> gtsummary::tbl_regression(model_brms, tidy_fun = broom.mixed::tidy)
x There was an error calling `tidy_fun()`. Most likely, this is because
  the function supplied in `tidy_fun=` was misspelled, does not exist,
  is not compatible with your object, or was missing necessary arguments (e.g.
  `conf.level=` or `conf.int=`). See error message below.
Error: Error: The model does not contain group-level effects.
> 

[storopoli]

@ddsjoberg the models also do not compile on my machine (I think it's a rstan bug). You can try with cmdstanr.

remotes::install_github("stan-dev/cmdstanr")
library(cmdstanr)
install_cmdstan()
model_brms <- brm(mpg ~ hp + vs, data = mtcars, backend = "cmdstanr")

[ddsjoberg]

@storopoli is there is tidier you know of that works on these models already?

[storopoli]

Yes, tidybayes and effectsize::standardize_parameters().

[ddsjoberg]

@storopoli if you write a tidier using the effectsize package, I can add it to gtsummary.

I think effectsize has a function to convert their model summary table to the broom tidy format

[storopoli]

@ddsjoberg I can sure try to help. Could you point me to an example or template to follow?

[ddsjoberg]

@storopoli I just looked into the details, and I remember why I didn't add brmsfit models: the model errors when we pass it to broom.helpers::tidy_plus_plus()

@larmarange I'll try to transfer this issue to broom.helpers to look into adding at some point if it's possible.

library(broom.mixed)
#> Registered S3 method overwritten by 'broom.mixed':
#>   method      from 
#>   tidy.gamlss broom

# load previously fit model
load(system.file("extdata", "brms_example.rda", package="broom.mixed"))
broom.mixed::tidy(brms_crossedRE)
#> # A tibble: 7 x 8
#>   effect   component group   term          estimate std.error conf.low conf.high
#>   <chr>    <chr>     <chr>   <chr>            <dbl>     <dbl>    <dbl>     <dbl>
#> 1 fixed    cond      <NA>    (Intercept)     33.9       3.81   26.3       40.5  
#> 2 fixed    cond      <NA>    wt              -3.91      1.39   -6.27      -1.17 
#> 3 ran_pars cond      cyl     sd__(Interce~    4.57      3.64    0.350     14.6  
#> 4 ran_pars cond      gear    sd__(Interce~    4.57      3.21    0.0637    12.8  
#> 5 ran_pars cond      gear    sd__wt           1.71      1.32    0.121      4.92 
#> 6 ran_pars cond      gear    cor__(Interc~   -0.287     0.583  -0.963      0.886
#> 7 ran_pars cond      Residu~ sd__Observat~    2.59      0.414   2.03       3.63

broom.helpers::tidy_plus_plus(brms_crossedRE)
#> x Unable to identify the list of variables.
#>   
#>   This is usually due to an error calling `stats::model.frame(x)`or `stats::model.matrix(x)`.
#>   It could be the case if that type of model does not implement these methods.
#>   Rarely, this error may occur if the model object was created within
#>   a functional programming framework (e.g. using `lappy()`, `purrr::map()`, etc.).
#> Error in terms.default(model): no terms component nor attribute

^{Created on 2021-01-26 by the reprex package (v0.3.0)}

[ddsjoberg]

@storopoli once the issue is solved on broom.helpers, please re-open this issue and we'll finalize the gtsummary portion (very easy to do)

[larmarange]

@ddsjoberg brms models are now supported in broom.helpers. Please note that duplicated terms are not disambiguated for mixed models. I let you see how it is managed by tbl_regression().

Please let me know if at some point you need a new release of broom.helpers on CRAN.

[storopoli]

@ddsjoberg broom.helpers 1.2.0 was released 3 days ago. How do I reopen the issue to finalize the gtsummary portion?

[larmarange]

To be noted that now in broom.helpers 1.2.0 are, by default, disambiguated.