---
bibliography: regression_output.bib
biblio-style: apalike
link-citations: yes
output:  
  bookdown::html_document2:
    theme: default
    highlight: textmate
    code_folding: show
    mathjax: local
    self_contained: FALSE
    number_sections: FALSE
    css: regression_output.css
---

# Introduction
  
The general linear model (GLM), more commonly called linear regression, is the most common statistical modeling tool a statistician or data scientist will use.  As such, it is crucial to know how to present results from a GLM in a way that is understandable to your audience.

The GLM assumes that the conditional distribution of the $i^{th}$ response is normal with expectation that is a linear function of a set of $p$ predictors with uncorrelated errors that represent the aggregate omitted causes of Y (equation \@ref(eq:ols)) [@fox2008applied].

\begin{equation} 
  Y_{i} \sim N(\beta_0 + X_{i1}\beta_{1} + + X_{i2}\beta_{2} + ... + X_{ip}\beta_{p}, \sigma_{\epsilon}^2)
  (\#eq:ols)
\end{equation}

This presentation uses the R programming language and assumes the end user is taking advantage of RStudio IDE to compile their R markdown files into HTML [@Rlang2019; @Rstudio].  All of the files needed to reproduce these results can be downloaded from the Git repository <code>git clone <span>https</span>://git.waderstats.com/regression_output/</code>.

# Required Libraries

The libraries <i>knitr</i>, <i>bookdown</i>, and <i>kableExtra</i> are used to generate the HTML output [@knitr; @bookdown; @kableExtra].  The <i>ggplot2</i> library is loaded for the example data set that is used in this presentation [@ggplot2].  The <i>Hmisc</i> library provides functionality needed to create variable labels [@Hmisc].  Regression diagnostics are plotted using the <i>ggfortify</i> package [@ggfortify].  The <i>stargazer</i> library will be used create nice looking regression output [@stargazer].

```{r libraries, eval = TRUE, echo = TRUE, results = FALSE, warning = FALSE, message = FALSE}
package_loader <- function(x, ...) {
  if (x %in% rownames(installed.packages()) == FALSE) install.packages(x)
  library(x, ...)
}

packages <- c("knitr", "bookdown", "kableExtra", "ggplot2", "Hmisc", "ggfortify", "stargazer")

invisible(sapply(X = packages, FUN = package_loader, character.only = TRUE))
```

# Example Data Setup

The data set used in this presentation is <b>mpg</b> from the <i>ggplot2</i> package.  From the description in the manual:

> This dataset contains a subset of the fuel economy data that the EPA makes available <a href="http://fueleconomy.gov" target="_blank">here</a>. It contains only models which had a new release every year between 1999 and 2008 - this was used as a proxy for the popularity of the car.

```{r ExampleDataSetup, eval = TRUE, echo = TRUE, results = TRUE, warning = FALSE, message = FALSE}
set.seed(123)
data(mpg)
mpg <- data.frame(mpg)

colnames(mpg)[which(colnames(mpg) == "manufacturer")] <- "manu"

mpg$manu <- factor(mpg$manu)
mpg$model <- factor(mpg$model)
mpg$displ <- as.numeric(mpg$displ)
mpg$year <- factor(mpg$year, levels = c("1999", "2008"), ordered = TRUE)

mpg$dp <- as.Date(NA, origin = "1970-01-01")
mpg$dp[which(mpg$year == "1999")] <- sample(seq(as.Date('1999-01-01', format = "%Y-%m-%d", origin = "1970-01-01"), as.Date('1999-12-25', format = "%Y-%m-%d", origin = "1970-01-01"), by="day"), dim(mpg)[1]/2)
mpg$dp[which(mpg$year == "2008")] <- sample(seq(as.Date('2008-01-01', format = "%Y-%m-%d", origin = "1970-01-01"), as.Date('2008-12-25', format = "%Y-%m-%d", origin = "1970-01-01"), by="day"), dim(mpg)[1]/2)

mpg$cyl <- factor(mpg$cyl, levels = c(4, 5, 6, 8), ordered = TRUE)
mpg$trans <- factor(mpg$trans)
mpg$drv <- factor(mpg$drv, levels = c("f", "r", "4"), labels = c("front-wheel drive", "rear wheel drive", "4wd"))
mpg$fl <- factor(mpg$fl)
mpg$class <- factor(mpg$class)

mpg$rn <- rnorm(dim(mpg)[1], mean = 10, sd = 5)
mpg$rn[sample(1:length(mpg$rn), size = 50)] <- NA

mpg$party <- factor(sample(c("republican", "democrat", "independent", NA), dim(mpg)[1], replace = TRUE), levels = c("republican", "democrat", "independent"))
mpg$comments <- sample(c("I like this car!", "Meh.", "This is the worst car ever!", "Does it come in green?", "want cheese flavoured cars.", "Does it also fly?", "Blah, Blah, Blah, Blah, Blah, Blah, Blah, Blah", NA), dim(mpg)[1], replace = TRUE)

label(mpg$manu) <- "manufacturer"
label(mpg$model) <- "model name"
label(mpg$displ) <- "engine displacement, in litres"
label(mpg$year) <- "year of manufacture"
label(mpg$dp) <- "date of purchase"
label(mpg$cyl) <- "number of cylinders"
label(mpg$trans) <- "type of transmission"
label(mpg$drv) <- "drive type"
label(mpg$cty) <- "city miles per gallon"
label(mpg$hwy) <- "highway miles per gallon"
label(mpg$fl) <- "fuel type"
label(mpg$class) <- "type of car"
label(mpg$rn) <- "some random numbers that are generated from a normal distrubtion with mean = 10 and sd = 5"
label(mpg$party) <- "some random political parties"
label(mpg$comments) <- "some random comments"

kable(head(mpg), caption = "Header of <b>mpg</b>.", booktabs = TRUE, escape = FALSE) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
```

# City MPG Regression

Model building starts with what predictors we hypothesize that the outcome depends on.  In this case, I hypothesize that miles per gallon in the city and highway depends on a constant, engine displacement, year of manufacture, the number of cylinders, drive type, fuel type, and type of car, some random numbers, and political party affiliation (equations \@ref(eq:cty) and \@ref(eq:hwy)).

\begin{equation} 
  cty_{i} = \beta_0 + dipl_{i}\beta_{1} + year_{i}\beta_{2} + cyl_{i}\beta_{3} + drv_{i}\beta_{4} + \\class_{i}\beta_{5} + rn_{i}\beta_{6} + party_{i}\beta_{7} + \epsilon_{i}
  (\#eq:cty)
\end{equation}

\begin{equation} 
  hwy_{i} = \beta_0 + dipl_{i}\beta_{1} + year_{i}\beta_{2} + cyl_{i}\beta_{3} + drv_{i}\beta_{4} + \\class_{i}\beta_{5} + rn_{i}\beta_{6} + party_{i}\beta_{7} + \epsilon_{i}
  (\#eq:hwy)
\end{equation}

## Regression Diagnostics

The regression diagnostics show a few problems.  The residuals vs. fitted plot shows a few points that have large residuals.  These points translate to a normal quantile plot with a large discrepancy between the standardized residuals and the theoretical quantiles on the tail.  From the scale location plots, we see relative homoscedasticity except for the outliers.  The plot of Cook's distance, a measure of influence on the fit, shows that the outliers have a disproportionate influence on the estimates in the final fit.

```{r ctyRegression, eval = TRUE, echo = TRUE, results = TRUE, warning = FALSE, message = FALSE, fig.width = 8,  fig.asp = 0.7, fig.align = 'center', fig.cap = "City MPG Regression Diagnostics"}
ctyFit <- lm(cty ~ displ + year + cyl + drv + class + rn + party, data = mpg)

autoplot(ctyFit, which = c(1, 2, 3, 4), ncol = 2, label.size = 3)
```

```{r hwyRegression, eval = TRUE, echo = TRUE, results = TRUE, warning = FALSE, message = FALSE, fig.width = 8,  fig.asp = 0.7, fig.align = 'center', fig.cap = "Highway MPG Regression Diagnostics"}
hwyFit <- lm(hwy ~ displ + year + cyl + drv + class + rn + party, data = mpg)

autoplot(hwyFit, which = c(1, 2, 3, 4), ncol = 2, label.size = 3)
```

## Regression Output {#regression}

The regression starts with re-fitting the models with the problem points removed.  Finally, we use the <i>stargazer</i> function to print the results. To customzie the disply the results, we also need a bit of css which can be downloaded <a href="https://waderstats.com/rstudio-files/regression_output/regression_output.css" target="_blank">here</a>.  The <i>id</i> for the CSS is generated by adding <i>{#regression}</i> to the header of this section.  Check the <a href="https://waderstats.com/rstudio-files/regression_output/regression_output.Rmd" target="_blank">markdown</a> to see it in context.

```{r regressionOutput, eval = TRUE, echo = TRUE, results = 'asis', warning = FALSE, message = FALSE}
ctyFit <- lm(cty ~ displ + year + cyl + drv + class + rn + party, data = mpg[c(-213, -223), ])

hwyFit <- lm(hwy ~ displ + year + cyl + drv + class + rn + party, data = mpg[c(-213, -223), ])

depLabels <- c(label(mpg$cty), label(mpg$hwy))

indepLabels <- c(
  "intercept",
  label(mpg$displ), 
  paste(label(mpg$year), ": ", levels(mpg$year)[2], " (ref = ", levels(mpg$year)[1], ")", sep = ""),
  paste(label(mpg$cyl), ": ", levels(mpg$cyl)[2], " (ref = ", levels(mpg$cyl)[1], ")", sep = ""),
  paste(label(mpg$cyl), ": ", levels(mpg$cyl)[3], " (ref = ", levels(mpg$cyl)[1], ")", sep = ""),
  paste(label(mpg$cyl), ": ", levels(mpg$cyl)[4], " (ref = ", levels(mpg$cyl)[1], ")", sep = ""),
  paste(label(mpg$drv), ": ", levels(mpg$drv)[2], " (ref = ", levels(mpg$drv)[1], ")", sep = ""), 
  paste(label(mpg$drv), ": ", levels(mpg$drv)[3], " (ref = ", levels(mpg$drv)[1], ")", sep = ""), 
  paste(label(mpg$class), ": ", levels(mpg$class)[2], " (ref = ", levels(mpg$class)[1], ")", sep = ""), 
  paste(label(mpg$class), ": ", levels(mpg$class)[3], " (ref = ", levels(mpg$class)[1], ")", sep = ""), 
  paste(label(mpg$class), ": ", levels(mpg$class)[4], " (ref = ", levels(mpg$class)[1], ")", sep = ""), 
  paste(label(mpg$class), ": ", levels(mpg$class)[5], " (ref = ", levels(mpg$class)[1], ")", sep = ""), 
  paste(label(mpg$class), ": ", levels(mpg$class)[6], " (ref = ", levels(mpg$class)[1], ")", sep = ""), 
  paste(label(mpg$class), ": ", levels(mpg$class)[7], " (ref = ", levels(mpg$class)[1], ")", sep = ""), 
  label(mpg$rn),
  paste(label(mpg$party), ": ", levels(mpg$party)[2], " (ref = ", levels(mpg$party)[1], ")", sep = ""),
  paste(label(mpg$party), ": ", levels(mpg$party)[3], " (ref = ", levels(mpg$party)[1], ")", sep = "")
)

stargazer(ctyFit, hwyFit, type = "html",
          title="Regression Results", 
          dep.var.labels = depLabels,
          covariate.labels = indepLabels,
          ci = TRUE,
          ci.level = 0.95,
          single.row = TRUE,
          align = TRUE,
          intercept.bottom = FALSE,
          intercept.top = TRUE
)
```

## Regression Summary

Using copy and paste to summarize results leaves a ton of room to make errors.  To prevent this, we can summarize the results using inline R code.  Check out the <a href="https://waderstats.com/rstudio-files/regression_output/regression_output.Rmd" target="_blank">markdown</a> file for how I did it in this case.

```{r regressionSummary, eval = TRUE, echo = TRUE, results = TRUE, warning = FALSE, message = FALSE}
ctyCoef <- coef(ctyFit)
ctyConf <- confint(ctyFit)
ctyRes <- cbind(ctyCoef, ctyConf)

hwyCoef <- coef(hwyFit)
hwyConf <- confint(hwyFit)
hwyRes <- cbind(hwyCoef, hwyConf)

inlineReg <- function(x, lab) {
  r1 <- round(x[lab, 1], 1)
  r2 <- round(x[lab, 2], 1)
  r3 <- round(x[lab, 3], 1)

  res <- paste(r1, " (95% CI = [", r2, ",", r3, "])", sep = "")
  
  return(res)
}
```

The results show that highway MPG is much better than city MPG, with `r inlineReg(x = ctyRes, lab = "(Intercept)")` and `r inlineReg(x = hwyRes, lab = "(Intercept)")` respectively.  We found that as compared to 1999, cars built in 2008 had an increase of  `r inlineReg(x = ctyRes, lab = "year.L")` in city MPG and `r inlineReg(x = hwyRes, lab = "year.L")` in highway MPG.  We found that minivans have the worst MPG as compared to the referent category, `r inlineReg(x = ctyRes, lab = "classminivan")` MPG lower in the city and `r inlineReg(x = hwyRes, lab = "classminivan")` MPG lower on the highway.

# Session Info

```{r sessionInfo, eval = TRUE, echo = TRUE, results = TRUE, warning = FALSE, message = FALSE}
sessionInfo()
```

# References

<!-- This is needed to resize iframe content. Please ignore if duplicationg this R markdown file --> 
<script type="text/javascript" src="https://www.waderstats.com/wp-includes/js/iframe-resizer/js/iframeResizer.contentWindow.min.js" defer></script>