R/R Markdown Table Options for Multilevel Models

The Models

For this example, I am going to use Allison Horst’s palmerpenguins package for the data and produce two simple and non-substantive multi-level models. I will use lme4::lmer and lmerTest::lme4, which produces p-values and a few different estimates. Note: I will not be interpreting these models.

library(palmerpenguins)

# unconditional model
model_0 <- lme4::lmer(body_mass_g ~ 1 + (1 | island), REML = F, data=penguins)

model_0_lmertest <- lmerTest::lmer(body_mass_g ~ 1 + (1 | island), REML = F, data=penguins)

model_1 <- lme4::lmer(body_mass_g ~ sex + (1 | island), REML = F, data=penguins)

model_1_lmertest <- lmerTest::lmer(body_mass_g ~ sex + (1 | island), REML = F, data=penguins)

The Packages

The most popular packages that can take a model object and produce a neat table summary include modelsummary, gtsummary, huxtable, and sjPlot. I will generate default tables with each of these packages with the goal of being able to insert them into a Word document via knitr. I will try to keep customization to a minimum.

modelsummary

modelsummary offers the msummary function and allows me to create side-by-side tables. This will work for lme4 models, but will not work for lmerTest models unless you load the package broom.mixed first.

modelsummary::msummary(list(
  "Null Model" = model_0, 
  "Model 1" = model_1))

	Null Model	Model 1
(Intercept)	4048.143	3713.505
	(275.116)	(273.909)
sd_(Intercept).island	471.900	468.149
sd_Observation.Residual	626.335	532.804
sexmale		674.365
		(58.402)
AIC	5393.7	5147.3
BIC	5405.2	5162.5
Log.Lik.	-2693.833	-2569.644

Here, I have loaded broom.mixed and added a call for significance stars with stars=T.

library(broom.mixed)

## Registered S3 methods overwritten by 'broom.mixed':
##   method         from 
##   augment.lme    broom
##   augment.merMod broom
##   glance.lme     broom
##   glance.merMod  broom
##   glance.stanreg broom
##   tidy.brmsfit   broom
##   tidy.gamlss    broom
##   tidy.lme       broom
##   tidy.merMod    broom
##   tidy.rjags     broom
##   tidy.stanfit   broom
##   tidy.stanreg   broom

modelsummary::msummary(list(
  "Null Model" = model_0_lmertest, 
  "Model 1" = model_1_lmertest), stars = T)

	Null Model	Model 1
(Intercept)	4048.143***	3713.505***
	(275.116)	(273.909)
sd__(Intercept)	471.900NA	468.149NA
sd__Observation	626.335NA	532.804NA
sexmale		674.365***
		(58.402)
AIC	5393.7	5147.3
BIC	5405.2	5162.5
Log.Lik.	-2693.833	-2569.644
* p < 0.1, ** p < 0.05, *** p < 0.01

Screenshot from Word

Screenshot from Word for msummary

Notes

• msummary produces a decent side-by-side table in Word.
• The default style looks like APA, which is something I would want.
• However, it splits the fixed effects of intercept and sex up so that intercept is above the random effects and sex is below.
• Finally, the default width is a bit narrow, meaning additional customization is needed. It can easily be customized using flextable, which works well in Word. The following code coverts it to flextable:

modelsummary::msummary(list("Null Model" = model_0, "Model 1" = model_1), output = 'flextable')

gtsummary

gtsummary is based on the gt package, which allows easy and advanced customization of tables in R.

I could not get gtsummary to make a table for my null model, model_0. It produced the error: Error: must be a character vector, not a logical vector.. However, it worked for model_1 and model_1_lmertest. Here is the no frills default:

gtsummary::tbl_regression(model_1_lmertest)

That table is clearly not a good table. I did some googling, and found I had to add the following tidy function to produce a table with all estimates:

gtsummary::tbl_regression(model_1_lmertest, tidy_fun = broom.mixed::tidy)

Screenshot from Word

Screenshot from Word of gtsummary

huxtable

huxtable works with both lme4 and lmerTest. It uses the huxreg function and can easily make side-by-side tables.

huxtable::huxreg(list(
                 "Null Model" = model_0_lmertest, 
                 "Model 1" = model_1_lmertest))

## Warning in huxtable::huxreg(list(`Null Model` = model_0_lmertest, `Model 1` = model_1_lmertest)): Unrecognized statistics: r.squared
## Try setting `statistics` explicitly in the call to `huxreg()`

Screenshot from Word

Screenshot from Word of huxtable

sjPlot

sjPlot is a table and visualization package meant for a number of social science methods (linear models, generalized linear models, PCA, robust standard errors, etc.).

Below is the HTML table. This table can be styled with CSS.

sjPlot::tab_model(model_0, model_1, 
                  dv.labels = c("Null Model", "Model 1"))

Final Thoughts

In a perfect world, we would be able to knit the sjPlot table to Word and be done with it. Unfortunately, this does not work for now. The best and easiest option would be using modelsummary::msummary. If you need p-values, use lmertest and broom.mixed, and be sure to call stars = T inside msummary. Then, make sure to also include output = 'flextable' to add more customization and allow your table to knit to Word.

The following example creates a simple function to automatically customize the table (useful if you will have multiple tables), and adds ICC and R2. I did not have time to organize the fixed and random effects, nor add the Greek symbols (ideally, in parentheses), but I think it’s easy to see how well modelsummary and flextable work together for your multilevel models!

library(tidyverse)
library(flextable)

#table formatting function
flex <- function(data, title=NULL) {
  set_table_properties(data, layout = "autofit", width = 1) %>%
  fontsize(size=12, part="all") %>%
      set_caption(title, 
                  autonum = officer::run_autonum(seq_id = "tab", 
                                                 pre_label = "Table ", 
                                                 post_label = "\n", 
                                                 bkm = "anytable"))
}

#adding ICC and R2 - note you can also write a function with glance_custom
# https://vincentarelbundock.github.io/modelsummary/articles/newmodels.html
mlm_added <- tribble(~term, ~"Null Model", ~"Model 1",
                     "ICC", 
                        performance::icc(model_0)$ICC_conditional, 
                        performance::icc(model_1)$ICC_conditional,
                     "pseudo R2",
          performance::r2(model_0)[["R2_conditional"]][["Conditional R2"]],
          performance::r2(model_1)[["R2_conditional"]][["Conditional R2"]])

#changing names of coefficients
cm <- c("sexmale"="Male",
        "sd__(Intercept)"="SD Intercept (τ00)",
        "sd__Observation" = "SD Observation (σ2)",
        "(Intercept)" = "Intercept")

#the final table
modelsummary::msummary(list("Null Model" = model_0, 
                            "Model 1" = model_1), 
                       output = 'flextable', 
                       add_rows = mlm_added,
                       coef_map = cm) %>%
  flex("Multilevel Model of Palmer Penguins")

Screenshot from Word

Screenshot from Word of Customized modelsummary

---

Do you know a better way to make great MLM tables in R Markdown, or have you figured out a sjPlot hack!? Please let me know.

---

Thanks to Vincet Arel-Bundock (https://twitter.com/VincentAB; author of modelsummary) for the tip about loading broom.mixed

Thanks to Daniel (https://twitter.com/strengejacke; author of sjPlot) for pointing out to use dv.labels