The support jtools
provides for helping to understand and report the results of regression
models falls into a few broad categories:
summ)effect_plot; see other vignette)plot_coefs, plot_summs)export_summs)summWhen 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:
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:
lm, glm, svyglm
(survey), merMod (lme4), and
rq (quantreg) models.lm and glm
plus quantreg’s built-in options for rq
models)lm only) can optionally be included in the outputlm, 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.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:
You can force knitr to give the console style of output
by setting the chunk option render = 'normal_print'.
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.
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.
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.
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 has been a proponent of dividing by 2 standard deviations; if
you want to do things that way, give the argument
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.
In many cases, you’ll learn more by looking at confidence intervals
than p-values. You can request them from summ.
You can adjust the width of the confidence intervals, which are by default 95% CIs.
You might also want to drop the p-values altogether.
Remember that you can omit p-values regardless of whether you have requested confidence intervals.
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.
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!
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.
Note that the summary of linear mixed models will omit p-values by
default unless the 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 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.
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.
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:
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.
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.
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:
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.
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.
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:
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:
summ() methodplot_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().
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:
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:
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):
You can likewise export to PDF ("PDF"), HTML
("HTML"), or Excel format ("xlsx").
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.
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.
You can pre-set the number of digits you want printed for all
jtools functions with the jtools-digits
option.
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.
summYou 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:
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:
digitsmodel.infomodel.fitpvalsrobustconfintci.widthvifsconf.method (merMod models only)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.
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).