Supported models

modelsummary automatically supports all the models supported by the broom package. Specifically, modelsummary uses the glance and tidy methods from the broom package to extract information from model objects. To see if a package is supported out-of-the-box, load the broom package, and see if the tidy and glance methods produce valid output:

library(broom)

tidy(mod)
glance(mod)

If the tidy and glance functions do not work on your model type, you may consider opening an issue on the Github website of the broom package: https://github.com/tidymodels/broom/issues

Note that several packages supply additional tidy and glance functions. Loading these packages can increase the range of models supported by modelsummary. For example, loading the broom.mixed package allows modelsummary to support a wide variety of mixed (hierarchical, multi-level) models (e.g., lme4, brms, rstanarm).

library(broom.mixed)
library(modelsummary)

modelsummary(mod)

Customizing existing models

The information that modelsummary can display about a model is determined by the output of the broom::glance and broom::tidy functions. glance is used the extract the information at the the bottom of the table, where goodness-of-fit measures are displayed. tidy is used to extract the information in the body of the table: parameter and uncertainty estimates.

Sometimes users will want to include information that is not supplied by those functions. A pretty easy way to include extra information is to define new glance_custom and tidy_custom methods. To illustrate, we estimate two linear regression models using the lm function:

library(modelsummary)

mod <- list()
mod[[1]] <- lm(hp ~ mpg + drat, mtcars)
mod[[2]] <- lm(wt ~ mpg + drat + am, mtcars)

In R, the lm function produces models of class “lm”:

class(mod[[1]])
#> [1] "lm"

Let’s say you would like to print the dependent variable for each model of this particular class. All you need to do is define a new method called glance_custom.lm. This method should return a data.frame (or tibble) with 1 row, and 1 column per piece of information you want to display. For example:

glance_custom.lm <- function(x) {
    dv <- as.character(formula(x)[2])
    out <- tibble::tibble(`Dependent variable:` = dv)
    return(out)
}

; ow, let’s customize the body of the table. The statistic_override argument already allows users to customize uncertainty estimates. But imagine you want to override the coefficient estimates of your “lm” models. Easy! All you need to do is define a tidy_custom.lm method which returns a data.frame (or tibble) with one column called “term” and one column called “estimate”.

Here, we’ll substitute estimates by an up/down arrow which represents their signs:

tidy_custom.lm <- function(x) {
    s <- summary(x)$coefficients
    out <- tibble::tibble(term = row.names(s),
                          estimate = ifelse(s[,1] > 0, '⭡', '⭣'))
    return(out)
}

After you define the glance_custom and tidy_custom methods, modelsummary will automatically display your customized model information:

modelsummary(mod)

	Model 1	Model 2
(Intercept)	⭡	⭡
	(55.415)	(0.728)
mpg	⭣	⭣
	(1.792)	(0.019)
drat	⭡	⭣
	(20.198)	(0.245)
am		⭣
		(0.240)
Num.Obs.	32	32
R2	0.614	0.803
R2 Adj.	0.588	0.782
AIC	337.9	46.4
BIC	343.7	53.7
Log.Lik.	-164.940	-18.201
Dependent variable:	hp	wt
F	23.100	38.066

Note that you can define a std.error column in tidy_custom.lm to replace the uncertainty estimates instead of the coefficients.

Fixed effects regression

One common use-case for glance_custom is to display the list of “fixed effects” included in a linear regression model.

The feols function from the fixest package offers an extremely fast way to estimate fixed effects regression models. The feols function produces an object of class “fixest”, which stores the fixed effects variables in a character vector called “fixef_vars”. To display those fixed effects, we create a glance_custom.fixest and call modelsummary:

glance_custom.fixest <- function(x) {
    out <- tibble::tibble(.rows = 1)
    for (n in x$fixef_vars) {
        out[[paste('FE: ', n)]] <- 'X'
    }
    return(out)
}

library(fixest)
library(modelsummary)

url <- 'https://vincentarelbundock.github.io/Rdatasets/csv/plm/EmplUK.csv'
dat <- read.csv(url)
mod <- list()
mod[[1]] <- feols(emp ~ wage, dat)
mod[[2]] <- feols(emp ~ wage | firm, dat)
mod[[3]] <- feols(emp ~ wage | firm + year, dat)

modelsummary(mod)

	Model 1	Model 2	Model 3
(Intercept)	15.228
	(2.149)
wage	-0.307	-0.137	-0.062
	(0.087)	(0.061)	(0.050)
Num.Obs.	1031	1031	1031
R2	0.012	0.981	0.983
R2 Adj.	0.011	0.978	0.980
R2 Pseudo
R2 Within		0.016	0.003
AIC	8625.3	4825.0	4751.9
BIC	8635.1	5521.2	5487.7
Log.Lik.	-4310.633	-2271.475	-2226.949
FE: firm		X	X
FE: year			X
Std. errors	Standard	Clustered (firm)	Clustered (firm)

Adding new models

modelsummary relies on two functions from the broom package to extract model information: tidy and glance. If broom doesn’t support the type of model you are trying to summarize, modelsummary won’t support it out of the box. Thankfully, it is extremely easy to add support for most models using custom methods.

For example, models produced by the MCMCglmm package are not currently supported by broom. To add support, you simply need to create a tidy and a glance method:

# load packages and data
library(modelsummary)
library(MCMCglmm)
data(PlodiaPO)

# add custom functions to extract estimates (tidy) and goodness-of-fit (glance) information
tidy.MCMCglmm <- function(object, ...) {
    s <- summary(object, ...)
    ret <- tibble::tibble(term = row.names(s$solutions),
                          estimate = s$solutions[, 1],
                          conf.low = s$solutions[, 2],
                          conf.high = s$solutions[, 3])
    ret
}
glance.MCMCglmm <- function(object, ...) {
    ret <- tibble::tibble(dic = object$DIC,
                          n = nrow(object$X))
    ret
}

# estimate a simple model
model <- MCMCglmm(PO ~ 1 + plate, random = ~ FSfamily, data = PlodiaPO, verbose=FALSE, pr=TRUE)

# summarize the model
modelsummary(model, statistic = 'conf.int')

Two important things to note. First, the methods are named tidy.MCMCglmm and glance.MCMCglmm because the model object I am trying to summarize is of class MCMCglmm. You can find the class of a model by running: class(model).

Second, in the example above, we used the statistic = 'conf.int' argument. This is because the tidy method produces conf.low and conf.high columns. In most cases, users will define std.error column in their custom tidy methods, so the statistic argument will need to be adjusted.

If you create new tidy and glance methods, please consider contributing them to broom so that the rest of the community can benefit from your work: https://github.com/tidymodels/broom

2-stage instrumental variable estimation

2-stage instrumental variable (IV) estimators are popular in many disciplines. Analysts often want to draw tables that include:

The 2nd stage estimates
The 1st stage estimates, but only for variables not included in the 2nd stage (i.e., the “instruments”).

The problem is that many IV regression packages do not supply results for the 1st stage. Even when these results are accessible, the tidy methods supplied by broom only return the 2nd stage results.

To achieve our goal, we will:

Define a new tidy.iv method which extracts information about both stages of the IV estimation.
Change the class of our model objects to iv. This ensure that the new tidy method is triggered by modelsummary.
Sort the coefficients using the coef_map argument, to ensure that 1st stage results are displayed below 2nd stage results.

For convenience, modelsummary includes two unexported functions which can do most this for you. Reading the annotated code for these two functions could be an instructive exercise for those who wish to exploit everything that modelsummary has to offer.

To print IV regression results from lfe::felm models, we execute this code:

library(lfe)

models <- list(
  felm(wt ~ drat | 0 | (disp ~ hp), data = mtcars),
  felm(wt ~ drat + vs | 0 | (disp ~ hp), data = mtcars)
  )

# define a new method to combine tidy information about both stages
tidy.iv <- modelsummary:::tidy_felm_iv

# change the model classes to trigger the new method
class(models[[1]]) <- class(models[[2]]) <- c('iv', class(models[[1]]))

# use helper function to sort coefficients with 1st stage last
cm <- modelsummary:::coef_map_felm_iv(models)

modelsummary(models, coef_map=cm)

	Model 1	Model 2
(Intercept)	3.243	2.611
	(1.159)	(1.543)
drat	-0.371	-0.293
	(0.254)	(0.289)
`disp(fit)`	0.006	0.007
	(0.001)	(0.002)
vs		0.237
		(0.307)
Stage 1 hp	1.069	0.844
	(0.175)	(0.226)
Num.Obs.	32	32
R2	0.801	0.811
R2 Adj.	0.787	0.790

Adding and customizing statistical models

Supported models

Customizing existing models

Fixed effects regression

Adding new models

2-stage instrumental variable estimation