---
title: "Tools for summarizing and visualizing regression models"
author: "Jacob Long"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tools for summarizing and visualizing regression models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r echo=FALSE}
required <- c("survey", "huxtable", "broom")
if (!all(sapply(required, requireNamespace, quietly = TRUE)))
  knitr::opts_chunk$set(eval = FALSE)
knitr::opts_chunk$set(message = F, warning = F, fig.width = 6, fig.height = 5)
library(jtools)
```

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 best to pack things
into a reusable function: It became `summ`.

With no user-specified arguments except a fitted model, the output of `summ`
looks like this:

```{r}
# Fit model
fit <- lm(Income ~ Frost + Illiteracy + Murder, data = as.data.frame(state.x77))
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`!

## Report robust standard errors

One of the problems that originally motivated the creation of this function 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". Stata's default is "HC1", so you may want to use
that if your goal is to replicate Stata analyses.

```{r}
summ(fit, robust = TRUE, robust.type = "HC3")
```

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 a `robust = TRUE` parameter 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 beta 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 `scale.response = TRUE` argument.

```{r}
summ(fit, scale = TRUE)
```

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`. 

```{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.

```{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 = 2)
```

You can adjust the width of the confidence intervals, which are by default
95% CIs.

```{r}
summ(fit, confint = TRUE, ci.width = .5, digits = 2)
```

You might also want to drop the p-values altogether.

```{r}
summ(fit, confint = TRUE, pvals = FALSE, digits = 2)
```

Note 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` was a
natural extension and will cover most use cases.

```{r}
fitg <- glm(vs ~ drat + mpg, data = mtcars, family = binomial)

summ(fitg)
```

For exponential family models, especially logit and Poisson, you may be
interested in getting odds ratios rather than the linear beta estimates. 
`summ` can handle that!

```{r}
summ(fitg, odds.ratio = 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 omits 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 like you would for OLS models will 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 can also get an estimated model R-squared for mixed models using the 
Nakagawa & Schielzeth (2013) procedure with code adapted from the `MuMIn`
package. To get those, simply set `r.squared = TRUE` but note that it sometimes
takes a long time to run and is a bit error-prone. 

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

# Plotting coefficients and predicted lines

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

In the above example, the differing scales of the variables makes it kind of
difficult to make a quick assessment. No problem, we can just use the tools
built into `summ` for that.

```{r}
plot_summs(fit, scale = TRUE)
```

See? Now we have a better idea of how the uncertainty and magnitude of effect
differs for these variables. Note that by default the width of the confidence
interval is .95, but this can be changed with the `ci_level` argument.

Another compelling use case for `plot_summs` is robust standard errors.

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(Income ~ Frost + Illiteracy + Murder + `HS Grad`,
           data = as.data.frame(state.x77))
plot_summs(fit, fit2, scale = TRUE)
```

This is a classic case in which adding a new predictor causes another one's
estimate to get much closer to zero. 

`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` package 
should work. Additionally, if you want to compare summ outputs with 
different arguments, you can use `plot_coefs` for that because `plot_coefs`
treats `summ` objects as a model type.

Let's look at how we can use `plot_coefs` to look at what happens when we
use non-robust standard errors and two different robust standard error 
specifications.

First, we save the `summ` objects:

```{r}
s1 <- summ(fit, scale = TRUE)
s2 <- summ(fit, scale = TRUE, robust = TRUE, robust.type = "HC0")
s3 <- summ(fit, scale = TRUE, robust = TRUE, robust.type = "HC3")
```

Now we give the `summ` objects to `plot_coefs`. Note the use of the 
`model.names` argument to make which model is which clearer.

```{r}
plot_coefs(s1, s2, s3, model.names = c("OLS","HC0","HC3"))
```

So in this case we can see that for `Frost` and `Murder`, each subsequent 
standard error specification increases the standard errors. The pattern doesn't
quite hold for `Illiteracy`.

## `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 is similar to `interact_plot` (see other vignette) but
for main effects in particular.

Here's the most basic use with our linear model:

```{r}
effect_plot(fit, pred = Illiteracy)
```

Okay, not so illuminating. Let's see the uncertainty around this line.

```{r}
effect_plot(fit, pred = Illiteracy, interval = TRUE)
```

Now we're getting somewhere.

How does this compare to the observed data? Let's see!

```{r}
effect_plot(fit, pred = Illiteracy, interval = TRUE, plot.points = TRUE)
```

Now we're really learning something about our model---and it looks like the
linear fit is basically correct. 

Another time that `effect_plot` is useful is for non-linear (generalized)
models. Let's plot our logit model from earlier.

```{r}
effect_plot(fitg, pred = mpg)
```

Aha! A lovely logistic curve. Now you know that when `mpg` is from about 30
and higher, it no longer is giving any appreciable increase in probability
because it has already reached 1. 

Next, we'll plot the observed values. Note the `jitter = 0`; by default, a 
small "jitter" is applied to the plotted
points in case they overlap, but with a logit model it would make it appear 
that some of the points are something other than 0 or 1.

```{r}
effect_plot(fitg, pred = mpg, plot.points = TRUE, jitter = 0)
```

Now you can see how at the low and high ends of `mpg`, there are clusters of
0 and 1 values while near the middle of the range, there are values of both 0 
and 1. 


# 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'}
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'}
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.word` argument (requires
the `officer` and `flextable` packages):

```{r eval = FALSE}
export_summs(fit, fit2, scale = TRUE, to.word = TRUE, word.file = "test.docx")
```

# 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 return object has non-rounded values if you wish to use them later.

```{r}
j <- summ(fit, digits = 3)

j$coeftable
```

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

Note that the `car` package is needed to calculate VIFs.

```{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. Ultimately, the main thing to consider is that small
effects are more likely to be "drowned out" by higher VIFs.



