# R/R Markdown Table Options for Multilevel Models

I am currently taking a multilevel modeling class in my PhD program. It is week 3-ish and I am learning a lot. The course is software-agnostic, meaning we can use any software we want (SPSS, SAS, Stata, R, HLM). The first assignment comes as a data set and Word document. I am working completely in R / R Markdown to generate both the data and model analyses and to complete the Word-based homework. Because I am working in an R -> Word workflow, I need to know what options I have for creating nice tables of my models.

So, I am taking a break from my first assignment to compare the default output of the most popular R model summary packages, look at some extended features, and make a cross-package comparison for my own future reference.

Let’s begin!

## 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

#### 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)
Characteristic Beta 95% CI1 p-value
sex
female
male 674 560, 789 <0.001

1 CI = Confidence Interval

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)
Characteristic Beta 95% CI1 p-value
sex
female
male 674 560, 789 <0.001
sd__(Intercept) 468
sd__Observation 533

1 CI = Confidence Interval

#### Notes

• This table is created in gt, but as of right now, gt doesn’t work with Word so it was automatically converted to a kable table.
• You can do a lot of customization using gt and then convert it to kable or flextable (using as_flex_table()).
• While adding confidence intervals are nice, the lack of fit statistics makes the default output inadequate. Looking through the gtsummary reference, I do not see a way to add these.
• Side-by-side tables requires saving each table as an object and then using tbl_merge.

## 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()
 Null Model Model 1 (Intercept) 4048.143 *** 3713.505 *** (275.116) (273.909) sd__(Intercept) 471.900 468.149 (NA) (NA) sd__Observation 626.335 532.804 (NA) (NA) sexmale 674.365 *** (58.402) N 342 333 logLik -2693.833 -2569.644 AIC 5393.665 5147.289 *** p < 0.001; ** p < 0.01; * p < 0.05.

#### Notes

• By default huxreg produces a nice side-by-side table.
• The default looks like simple ASCII/Markdown table, but it knits to Word nicely without any need to make changes.
• Output looks like modelsummary::msummary but lacks ICC values.
• huxtable also separates fixed effects like modelsummary
• Can be customized via huxtable or flextable via as_flextable - though I could not get this to work.

## 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"))
Null Model Model 1
Predictors Estimates CI p Estimates CI p
(Intercept) 4048.14 3508.92 – 4587.36 <0.001 3713.50 3176.65 – 4250.36 <0.001
sex [male] 674.37 559.90 – 788.83 <0.001
Random Effects
σ2 392294.95 283879.82
τ00 222689.64 island 219163.06 island
ICC 0.36 0.44
N 3 island 3 island
Observations 342 333
Marginal R2 / Conditional R2 0.000 / 0.362 0.185 / 0.540

NONE

#### Notes

• It works with lmerTest. However, if you only use lme4, sjPlot automatically adds p-values.
• sjPlot creates a very nice table with random effects listed with their Greek values, $$\sigma^2$$ and $$\tau$$, plus ICC and $$R^2$$!
• Unfortunately, it produces only an HTML table and this table cannot be automatically added to Word.
• Styling the table for HTML requires CSS.
• It cannot be converted to flextable or kable

## 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',
coef_map = cm) %>%
flex("Multilevel Model of Palmer Penguins")
 Null Model Model 1 Male 674.365 (58.402) SD Intercept (t00) 471.900 468.149 SD Observation (s2) 626.335 532.804 Intercept 4048.143 3713.505 (275.116) (273.909) AIC 5393.7 5147.3 BIC 5405.2 5162.5 Log.Lik. -2693.833 -2569.644 ICC 0.362 0.355 pseudo R2 0.362 0.540

#### Screenshot from Word

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.