How to report variance estimates of random effects using gtsummary::tbl_regression()

user_11587586 0
Viewed 135

I am trying to figure out how to create nice output for a mixed effects model in R using the gtsummary function tbl_regression, but am getting error messages. I am simulating some data using the simstudy package, and then estimating a model.

library(simstudy)
library(lme4)
library(gtsummary)

### data definition

defC <- defData(varname = "c", formula = 0, variance = "1", dist="normal", id = "site")
  
defS <- defDataAdd(varname = "y", 
    formula = "1 * rx  + c", 
    variance ="10", dist = "normal")

### data generation 

set.seed(3231)
  
dc <- genData(50, defC, id = "site")
dc <- trtAssign(dc, grpName = "rx")

ds <- genCluster(dc, "site", 30, "id")
ds <- addColumns(defS, ds)

### model estimation

fit <- lmer(y ~ rx +  (1 | site), data = ds)
summary(fit)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ rx + (1 | site)
#>    Data: ds
#> 
#> REML criterion at convergence: 7767.8
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -3.3891 -0.6497  0.0002  0.6803  3.1299 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  site     (Intercept) 1.038    1.019   
#>  Residual             9.907    3.147   
#> Number of obs: 1500, groups:  site, 50
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)  0.05941    0.23393   0.254
#> rx           0.73473    0.33082   2.221
#> 
#> Correlation of Fixed Effects:
#>    (Intr)
#> rx -0.707

I can tidy the results directly using broom.mixed ...

broom.mixed::tidy(fit, scales = c("vcov", "sdcor"))
#> # A tibble: 4 × 6
#>   effect   group    term             estimate std.error statistic
#>   <chr>    <chr>    <chr>               <dbl>     <dbl>     <dbl>
#> 1 fixed    <NA>     (Intercept)        0.0594     0.234     0.254
#> 2 fixed    <NA>     rx                 0.735      0.331     2.22 
#> 3 ran_pars site     var__(Intercept)   1.04      NA        NA    
#> 4 ran_pars Residual var__Observation   9.91      NA        NA

And if I want to output the table using tbl_regression with standard deviance of the random effect, it works fine.

tbl_regression(fit, tidy_fun = broom.mixed::tidy)
Characteristic  Beta    95% CI1
rx  0.73    0.09, 1.4
site.sd__(Intercept)    1.0 
Residual.sd__Observation    3.1

But, I get an error message when I try to report the variance of the random effect (even though I know that tidy is working.

### gt output

tbl_regression(fit, 
    tidy_fun = function(x) broom.mixed::tidy(x, scales = c("vcov", "sdcor"))
)
#> 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 in tidy_fun(model, ...): unused arguments (conf.int = TRUE, conf.level = 0.95)

Any ideas what I'm doing wrong? If you need more info, let me know.

lme4 gtsummary r
1 Answers

The code you provided returns an error before we get to the tbl_regression() portion. I have therefore used the sleepstudy dataset included in the lme4 package as the example.

The issue you're running into is how you've defined the tidy function with the scales= argument. The tidy function must accept the model object, the exponentiate argument, and the conf.int argument. If you include ... as an argument, then you'll have no issues.

library(gtsummary)
#> #Uighur

fit <- lme4::lmer(Reaction ~ Days + (1 | Subject), lme4::sleepstudy)

tbl_regression(
  fit,
  tidy_fun = 
    function(x, ...) broom.mixed::tidy(x, scales = c("vcov", "sdcor"), ...)
) %>%
  as_kable() # converting to kable to display on SO
Characteristic Beta 95% CI
Days 10 8.9, 12
Subject.var__(Intercept) 1,378
Residual.var__Observation 960

Created on 2021-10-12 by the reprex package (v2.0.1)

user_87112760
Related