---
title: "Tools for summarizing and visualizing regression models"
author: "Jacob Long"
date: "`r Sys.Date()`"
output:
html_document:
toc: true
toc_float: true
theme: "spacelab"
vignette: >
%\VignetteIndexEntry{Tools for summarizing and visualizing regression models}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r echo=FALSE}
required <- c("survey", "huxtable", "broom", "lme4", "quantreg")
if (!all(sapply(required, requireNamespace, quietly = TRUE)))
knitr::opts_chunk$set(eval = FALSE)
knitr::opts_chunk$set(message = FALSE, warning = FALSE, fig.width = 6,
fig.height = 4, dpi = 125, render = knitr::normal_print)
library(jtools)
```
The support `jtools` provides for helping to understand and report the results
of regression models falls into a few broad categories:
* Generating flexible table output in the console that includes multiple
standard error specifications, toggles for confidence intervals, VIFs, p values,
and so on (`summ`)
* Plotting predicted data from models to aid in substantive interpretation
and understanding model fit, including models with interactions (`effect_plot`;
see other vignette)
* Plotting regression coefficients and their uncertainty in a visually appealing
way (`plot_coefs`, `plot_summs`)
* Exporting regression summaries as tables in PDF/LaTeX and Word formats for
publication (`export_summs`)
# `summ`
When sharing analyses with colleagues unfamiliar with R, I found that the output
generally was not clear to them. Things were even worse if I wanted to give
them information that is not included in the `summary()` like robust
standard errors, scaled coefficients, and VIFs since the functions for
estimating these don't append them to a typical regression table. After creating
output tables "by hand" on multiple occasions, I thought it would be best to
pack things into a reusable function: It became `summ()`.
For example purposes, we'll create a model using the `movies` data from this
package. These data comprise information about over 800 movies across several
decades. We will be predicting the Metacritic `metascore`, which ranges from
0 to 100 (where higher numbers reflect more positive reviews) using
the gross revenue in the United States (`us_gross`), the fan rating at IMDB
(`imdb_rating`), and a categorical variable reflecting the genre (`genre5`) with
"Action" as the reference level.
With no user-specified arguments except a fitted model, the output of `summ()`
looks like this:
```{r}
library(jtools) # Load jtools
data(movies) # Telling R we want to use this data
fit <- lm(metascore ~ imdb_rating + log(us_gross) + genre5, data = movies)
summ(fit)
```
Like any output, this one is somewhat opinionated — some information is shown
that perhaps not everyone would be interested in, some may be missing. That,
of course, was the motivation behind the creation of the function; I didn't
like the choices made by R's core team with `summary()`!
Here's a quick (not comprehensive) list of functionality supported by `summ`:
* Summaries for `lm`, `glm`, `svyglm` (`survey`), `merMod` (`lme4`),
and `rq` (`quantreg`) models.
* Variable scaling and centering (i.e., calculating standardized coefficients)
* Robust standard errors (for `lm` and `glm` plus `quantreg`'s built-in
options for `rq` models)
* Confidence intervals, VIFs, and partial correlations (`lm` only) can
optionally be included in the output
* p-values can be dropped from the output
* R^2 (`lm`, linear `svyglm`), pseudo-R^2 (`glm`, `svyglm`, `merMod`),
R^1 (`rq`), and other model fit statistics are calculated and reported.
These can also be suppressed if you don't want them.
* Ability to choose defaults for many options using `set_summ_defaults()` to
reduce the need to do redundant typing in interactive use.
Model types supported are `lm`, `glm`, `svyglm`, `merMod`, and `rq`, though
not all will be reviewed in detail here.
**Note:** The output in this vignette will mimic how it looks in the R console,
but if you are generating your own RMarkdown documents and have `kableExtra`
installed, you'll instead get some prettier looking tables like this:
```{r render = 'knit_print'}
summ(fit)
```
You can force `knitr` to give the console style of output by setting the
chunk option `render = 'normal_print'`.
## Report robust standard errors
One of the problems that originally motivated the creation of `summ()` was
the desire to efficiently report robust standard errors — while it is easy
enough for an experienced R user to calculate robust standard errors, there are
not many simple ways to include the results in a regression table as is common
with the likes of Stata, SPSS, etc.
Robust standard errors require the user to have the `sandwich`
package installed. It does not need to be loaded.
There are multiple types of robust standard errors that you may use, ranging
from "HC0" to "HC5". Per the recommendation of the authors of the `sandwich`
package, the default is "HC3" so this is what you get if you set
`robust = TRUE`. Stata's default is "HC1", so you may want to use
that if your goal is to replicate Stata analyses. To toggle the type of
robust errors, provide the desired type as the argument to `robust`.
```{r}
summ(fit, robust = "HC1")
```
Robust standard errors can also be calculated for generalized linear models
(i.e., `glm` objects) though there is some debate whether they should be used
for models fit iteratively with non-normal errors. In the case of `svyglm`, the
standard errors that package calculates are already robust to
heteroskedasticity, so any argument to `robust` will be ignored with a
warning.
You may also specify with `cluster` argument the name of a variable in the input
data or a vector of clusters to get cluster-robust standard errors.
## Standardized/scaled coefficients
Some prefer to use scaled coefficients in order to avoid dismissing an
effect as "small" when it is just the units of measure that are small.
scaled betas are used instead when `scale = TRUE`. To be clear,
since the meaning of "standardized beta" can vary depending on who you talk to,
this option mean-centers the predictors as well but does not alter the dependent
variable whatsoever. If you want to scale the dependent variable too,
just add the `transform.response = TRUE` argument. This argument does not
do anything to factor variables and doesn't scale binary variables by default.
```{r}
summ(fit, scale = TRUE)
```
If you have transformed
variables (e.g., `log(us_gross)`), the function will scale the
already-transformed variable. In other words, it is similar to the result
you would get if you did `scale(log(us_gross))` rather than
`log(scale(us_gross))` which would cause an error since you cannot apply
`log()` to numbers <== 0.
You can also choose a different number of standard deviations to divide by for
standardization.
[Andrew Gelman](http://www.stat.columbia.edu/~gelman/research/published/standardizing7.pdf)
has been a proponent of dividing by 2 standard
deviations; if you want to do things that way, give the argument `n.sd = 2`.
```{r}
summ(fit, scale = TRUE, n.sd = 2)
```
Note that this is achieved by refitting the model. If the model took a long time
to fit initially, expect a similarly long time to refit it.
### Mean-centered variables
In the same vein as the standardization feature, you can keep the original scale
while still mean-centering the predictors with the `center = TRUE` argument.
As with `scale`, this is not applied to the response variable unless
`transform.response = TRUE`.
```{r}
summ(fit, center = TRUE)
```
## Confidence intervals
In many cases, you'll learn more by looking at confidence intervals than
p-values. You can request them from `summ`.
```{r}
summ(fit, confint = TRUE, digits = 3)
```
You can adjust the width of the confidence intervals, which are by default
95% CIs.
```{r}
summ(fit, confint = TRUE, ci.width = .5)
```
## Removing p values
You might also want to drop the p-values altogether.
```{r}
summ(fit, confint = TRUE, pvals = FALSE)
```
Remember that you can omit p-values regardless of whether you have requested
confidence intervals.
## Generalized and Mixed models
`summ` has been expanding its range of supported model types. `glm` models are
a straightforward choice. Here we can take our previous model, but make it
a probit model to reflect the fact that `metascore` is bound at 0 and 100
(analogous to 0 and 1). We'll use the `quasibinomial` family since this is
a percentage rather a binary outcome.
```{r}
fitg <- glm(metascore/100 ~ imdb_rating + log(us_gross) + genre5, data = movies,
family = quasibinomial())
summ(fitg)
```
For exponential family models, especially logit and Poisson, you may be
interested in getting the exponentiated coefficients rather than the linear
estimates. `summ` can handle that!
```{r}
summ(fitg, exp = TRUE)
```
Standard errors are omitted for odds ratio estimates since the confidence
intervals are not symmetrical.
You can also get summaries of `merMod` objects, the mixed models from the
`lme4` package.
```{r message = FALSE, warning = FALSE}
library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
summ(fm1)
```
Note that the summary of linear mixed models will omit p-values by default
unless the \code{pbkrtest} package is installed for linear models.
There's no clear-cut way to derive p-values with linear mixed models and
treating the t-values as you would for OLS models can lead to
inflated Type 1 error rates. Confidence intervals are
better, but not perfect. Kenward-Roger calculated degrees of freedom are
fairly good under many circumstances and those are used by default when
\code{pbkrtest} package is installed. Be aware that for larger datasets, this
procedure can take a long time.
See the documentation (`?summ.merMod`) for more info.
You also get an estimated model R-squared for mixed models using the
Nakagawa & Schielzeth (2013) procedure with code adapted from the `piecewiseSEM`
package.
### svyglm
I won't run through any examples here, but `svyglm` models are supported and
provide near-equivalent output to what you see here depending on whether they
are linear models or generalized linear models. See `?summ.svyglm` for details.
### Quantile regression
`summ()` also supports quantile regression models, as estimated by the `rq`
package. See `?summ.rq` for details.
# `effect_plot()`
Sometimes to really understand what your model is telling you, you need to
see the kind of predictions it will give you. For that, you can use
`effect_plot()`, which does what it sounds like. There is a separate vignette
available to explore all it can offer, but here's a basic example with our
`glm` model:
```{r}
effect_plot(fitg, pred = imdb_rating, interval = TRUE, plot.points = TRUE,
jitter = 0.05)
```
Now we're really learning something about our model---and you can see the close
but imperfect agreement between fans and critics.
# `plot_summs()` and `plot_coefs()`
When it comes time to share your findings, especially in talks, tables are
often not the best way to capture people's attention and quickly convey the
results. Variants on what are known by some as "forest plots" have been
gaining popularity for presenting regression results.
For that, `jtools` provides `plot_summs()` and `plot_coefs()`. `plot_summs()`
gives you a plotting interface to `summ()` and allows you to do so with
multiple models simultaneously (assuming you want to apply the same
arguments to each model).
Here's a basic, single-model use case.
```{r}
plot_summs(fit)
```
Note that the intercept is omitted by default because it often distorts the
scale and generally isn't of theoretical interest. You can change this behavior
or omit other coefficients with the `omit.coefs` argument.
We may still want to use other features of `summ()`, like having robust standard
errors. No problem.
```{r}
plot_summs(fit, robust = TRUE)
```
Note that by default the width of the confidence
interval is .95, but this can be changed with the `ci_level` argument. You
can also add a thicker band to convey a narrow interval using the
`inner_ci_level` argument:
```{r}
plot_summs(fit, inner_ci_level = .9)
```
### Plot coefficient uncertainty as normal distributions
Most of our commonly used regression models make an assumption that the
coefficient estimates are asymptotically normally distributed, which is how we
derive our confidence intervals, p values, and so on. Using the
`plot.distributions = TRUE` argument, you can plot a normal distribution
along the width of your specified interval to convey the uncertainty. This
is also great for didactic purposes.
While the common OLS model assumes a *t* distribution, I decided that they are
visually sufficiently close that I have opted not to try to plot the points
along a *t* distribution.
```{r}
plot_summs(fit, plot.distributions = TRUE, inner_ci_level = .9)
```
### Comparing model coefficients visually
Comparison of multiple models simultaneously is another benefit of plotting.
This is especially true when the models are nested. Let's fit a second model
and compare.
```{r}
fit2 <- lm(metascore ~ imdb_rating + log(us_gross) + log(budget) + genre5,
data = movies)
plot_summs(fit, fit2)
```
Doing this with `plot.distributions = TRUE` creates a nice effect:
```{r}
plot_summs(fit, fit2, plot.distributions = TRUE)
```
By providing a list of `summ()` arguments to `plot_summs()`, you can compare
results with different `summ()` arguments (each item in the list corresponds to
one model; the second list item to the second model, etc.). For instance,
we can look at how the standard errors differ with different `robust` arguments:
```{r}
plot_summs(fit, fit, fit, robust = list(FALSE, "HC0", "HC5"),
model.names = c("OLS", "HC0", "HC5"))
```
### Support for models with no `summ()` method
`plot_coefs()` is very similar to `plot_summs()`, but does not offer the
features that `summ()` does. The tradeoff, though, is that it allows for model
types that `summ()` does not — any model supported by `tidy()` from the
`broom` or `broom.mixed` packages should work.
Note: If you provide unsupported model
types to `plot_summs()`, it just passes them to `plot_coefs()`.
# Table output for Word and RMarkdown documents
Sometimes you really do want a table, but it can't be standard R output.
For that, you can use `export_summs()`. It is a wrapper around `huxtable`'s
`huxreg()` function that will give you nice looking output if used in
RMarkdown documents or, if requested, printed to a Word file. In the latter
case, complicated models often need more fine-tuning in Word, but it gets
you started.
Like `plot_summs()`, `export_summs()` is designed to give you the features
available in `summ()`, so you can request things like robust standard errors
and variable scaling.
Here's an example of what to expect in a document like this one:
```{r eval = FALSE}
export_summs(fit, fit2, scale = TRUE)
```
```{r echo = FALSE, results = 'asis'}
huxtable::print_html(export_summs(fit, fit2, scale = TRUE))
```
When using RMarkdown, set `results = 'asis'` for the chunk with `export_summs()`
to get the right formatting for whatever type of output document (HTML, PDF,
etc.)
To format the error statistics, simply put the statistics desired in curly
braces wherever you want them in a character string. For example, if you want
the standard error in parentheses, the argument would be `"({std.error})"`,
which is the default. Some other ideas:
* `"({statistic})"`, which gives you the test statistic in parentheses.
* `"({statistic}, p = {p.value})"`, which gives the test statistic followed by
a "p =" p value all in parentheses. Note that you'll have to pay special
attention to rounding if you do this to keep cells sufficiently narrow.
* `"[{conf.low}, {conf.high}]"`, which gives the confidence interval in the
standard bracket notation. You could also explicitly write the confidence
level, e.g., `"95% CI [{conf.low}, {conf.high}]"`.
Here's an example with confidence intervals instead of standard errors:
```{r eval = FALSE}
export_summs(fit, fit2, scale = TRUE,
error_format = "[{conf.low}, {conf.high}]")
```
```{r echo = FALSE, results = 'asis'}
huxtable::print_html(export_summs(fit, fit2, scale = TRUE,
error_format = "[{conf.low}, {conf.high}]"))
```
There's a lot more customization that I'm not covering here: Renaming the
columns, renaming/excluding coefficients, realigning the errors, and so on.
If you want to save to a Word doc, use the `to.file` argument (requires
the `officer` and `flextable` packages):
```{r eval = FALSE}
export_summs(fit, fit2, scale = TRUE, to.file = "docx", file.name = "test.docx")
```
You can likewise export to PDF (`"PDF"`), HTML (`"HTML"`), or Excel format
(`"xlsx"`).
# Other options
## Adding and removing written output
Much of the output with `summ` can be removed while there are several other
pieces of information under the hood that users can ask for.
To remove the written output at the beginning, set `model.info = FALSE` and/or
`model.fit = FALSE`.
```{r}
summ(fit, model.info = FALSE, model.fit = FALSE)
```
## Choose how many digits past the decimal to round to
With the `digits =` argument, you can decide how precise you want the outputted
numbers to be. It is often inappropriate or distracting to report quantities with
many digits past the decimal due to the inability to measure them so precisely or
interpret them in applied settings. In other cases, it may be necessary to use
more digits due to the way measures are calculated.
The default argument is `digits = 2`.
```{r}
summ(fit, model.info = FALSE, digits = 5)
```
```{r}
summ(fit, model.info = FALSE, digits = 1)
```
You can pre-set the number of digits you want printed for all `jtools` functions
with the `jtools-digits` option.
```{r}
options("jtools-digits" = 2)
summ(fit, model.info = FALSE)
```
```{r echo = F}
options("jtools-digits" = NULL)
```
Note that the `summ` object contains the non-rounded values if you want to use
them later. The digits option just affects the printed output.
```{r}
j <- summ(fit, digits = 3)
j$coeftable
```
## Set default arguments to `summ`
You may like some of the options afforded to you by `summ()` but may not
like the inconvenience of typing them over and over. To streamline your
sessions, you can use the `set_summ_defaults()` function to avoid redundant
typing.
It works like this:
```{r eval = F}
set_summ_defaults(digits = 2, pvals = FALSE, robust = "HC3")
```
If you do that, you will have 2 digits in your output, no p values displayed,
and "HC3" sandwich robust standard errors in your `summ` output for the rest
of the R session. You can also use this in a .RProfile, but remember that it
should be included in scripts so that your code runs the same on every computer
and every session.
Here are all the options that can be toggled via `set_summ_defaults`:
* `digits`
* `model.info`
* `model.fit`
* `pvals`
* `robust`
* `confint`
* `ci.width`
* `vifs`
* `conf.method` (merMod models only)
## Calculate and report variance inflation factors (VIF)
When multicollinearity is a concern, it can be useful to have VIFs reported
alongside each variable. This can be particularly helpful for model comparison
and checking for the impact of newly-added variables. To get VIFs reported in
the output table, just set `vifs = TRUE`.
```{r}
summ(fit, vifs = TRUE)
```
There are many standards researchers apply for deciding whether a VIF is too
large. In some domains, a VIF over 2 is worthy of suspicion. Others set the bar
higher, at 5 or 10. Others still will say you shouldn't pay attention to these
at all. Ultimately, the main thing to consider is that small effects are more
likely to be "drowned out" by higher VIFs, but this may just be a natural,
unavoidable fact with your model (e.g., there is no problem with high VIFs when
you have an interaction effect).