The Model-Fitting Paradigm in R

Vincenzo Coia

library(tidyverse)
library(broom)

YouTube video for this tutorial: https://youtu.be/qivE6exNsZI

Data wrangling and plotting can get you pretty far when it comes to drawing insight from your data. But, sometimes you need to go further. For example:

  • Investigate the relationship between two or more variables
  • Predict the outcome of a variable given information about other variables
  • Get evidence towards a hypothesis through a hypothesis test.

These all involve fitting models – or at least, involve performing more advanced computations than can be done with tidyverse packages like dplyr.

In this tutorial, we’ll cover:

  1. how a model is typically fit in R;
  2. how to extract information from the model after fitting it, using the broom package; or,
  3. extract information by other means if the broom package is not set up for your model.

This tutorial is not about the specifics of fitting a model. Even though a few references to statistical concepts are made, just take these for face value.

Example: linear model

Here are a couple tasks that model-fitting would be useful to address.

  1. A car weighs 4000 lbs. Provide a numerical prediction on its mpg (miles per gallon).
  2. Does the weight of a car influence its mpg?
  3. How many more miles per gallon can we expect of a car that weighs 1000 lbs less than another car?

A simple scatterplot will give us a general idea, but can’t give us specifics. Here, we use the mtcars dataset in the datasets package. A linear model is one example of a model that can attempt to answer all three – the line corresponding to the fitted model has been added to the scatterplot.

(By the way, you can plot the model using ggplot2’s geom_smooth() layer – but you can’t interact with the model downstream.)

We can make a prediction of mpg for a car weighing 3,000 lbs by evaluating the line at X = 3:

## [1] 21.25171

We can find the p-value testing the null hypothesis that the true slope of the regression line is 0:

## [1] 1.293959e-10

And we can calculate how many more miles per gallon a car weighing 1000 less lbs would be expected to have through the slope:

## [1] 5.344472

Fitting a model in R

Fitting a model in R typically involves using a function in the following format:

method(formula, data, options)

Method:

A function such as:

  • Linear Regression: lm
  • Generalized Linear Regression: glm
  • Local regression: loess
  • Quantile regression: quantreg::rq

Formula:

In R, takes the form y ~ x1 + x2 + ... + xp (use column names in your data frame), where y here means the outcome variable you’re interested in viewing in relation to other variables x1, x2, …

Data: The data frame or tibble.

Options: Specific to the method.

Running the code returns an object that you can then work with to extract results.

For example, let’s fit a linear regression model on a car’s mileage per gallon (mpg, “Y” variable) from the car’s weight (wt, “X” variable). Notice that no special options are needed for lm().

(my_lm <- lm(mpg ~ wt, data = mtcars))
## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Coefficients:
## (Intercept)           wt  
##      37.285       -5.344

Printing the model to the screen might lead you to incorrectly conclude that the model object my_lm only contains the above text. The reality is, my_lm contains a lot more, but special instructions have been given to R to only print out a special digested version of the object. This behaviour tends to be true with model objects in general, not just for lm().

Probing the model: Broom

Now that you have the model object, there are typically three ways in which it’s useful to probe and use the model object. The broom package has three crown functions that make it easy to extract each piece of information by converting your model into a tibble.

tidy()

Use the tidy() function for a statistical summary of each component of your model, where each component gets a row in the output tibble. For lm(), tidy() gives one row per regression coefficient (slope and intercept).

tidy(my_lm)
## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    37.3      1.88      19.9  8.24e-19
## 2 wt             -5.34     0.559     -9.56 1.29e-10

tidy() only works if it makes sense to talk about model “components”.

augment()

Use the augment() function to make predictions on a dataset by augmenting predictions as a new column to your data. By default, augment() uses the dataset that was used to fit the model.

augment(my_lm) %>% 
  print(n = 5)
## # A tibble: 32 x 9
##   .rownames           mpg    wt .fitted .resid .std.resid   .hat .sigma  .cooksd
##   <chr>             <dbl> <dbl>   <dbl>  <dbl>      <dbl>  <dbl>  <dbl>    <dbl>
## 1 Mazda RX4          21    2.62    23.3 -2.28     -0.766  0.0433   3.07  1.33e-2
## 2 Mazda RX4 Wag      21    2.88    21.9 -0.920    -0.307  0.0352   3.09  1.72e-3
## 3 Datsun 710         22.8  2.32    24.9 -2.09     -0.706  0.0584   3.07  1.54e-2
## 4 Hornet 4 Drive     21.4  3.22    20.1  1.30      0.433  0.0313   3.09  3.02e-3
## 5 Hornet Sportabout  18.7  3.44    18.9 -0.200    -0.0668 0.0329   3.10  7.60e-5
## # … with 27 more rows

We can also predict on new datasets. Here are predictions of mpg for cars weighing 3, 4, and 5 thousand lbs.

augment(my_lm, newdata = tibble(wt = 3:5))
## # A tibble: 3 x 2
##      wt .fitted
##   <int>   <dbl>
## 1     3    21.3
## 2     4    15.9
## 3     5    10.6

Notice that only the “X” variables are needed in the input tibble (wt), and that since the “Y” variable (mpg) wasn’t provided, augment() couldn’t calculate anything besides a prediction.

glance()

Use the glance() function to extract a summary of the model as a whole, in the form of a one-row tibble. This will give you information related to the model fit.

glance(my_lm)
## # A tibble: 1 x 12
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.753         0.745  3.05      91.4 1.29e-10     1  -80.0  166.  170.
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Probing the Model: Traditional Method

In order for a model to work with the broom package, someone has to go out of their way to contribute to the broom package for that model. While this has happened for many common models, many others remain without broom compatibility.

Here’s how to work with these model objects in that case.

Prediction

If broom::augment() doesn’t work, then the developer of the model almost surely made it so that the predict() function works (not part of the broom package). The predict() function typically takes the same format of the augment() function, but usually doesn’t return a tibble.

Here are the first 5 predictions of mpg on the my_lm object, defaulting to predictions made on the original data:

predict(my_lm) %>% 
  unname() %>% 
  head(5)
## [1] 23.28261 21.91977 24.88595 20.10265 18.90014

Here are the predictions of mpg made for cars with a weight of 3, 4, and 5 thousand lbs:

predict(my_lm, newdata = tibble(wt = 3:5)) %>% 
  unname()
## [1] 21.25171 15.90724 10.56277

Checking the documentation of the predict() function for your model isn’t obvious, because the predict() function is a “generic” function. Your best bet is to append the model name after a dot. For example:

  • For a model fit with lm(), try ?predict.lm
  • For a model fit with rq(), try ?predict.rq (from the quantreg package)

If that doesn’t work, just google it: "Predict function for rq"

Functions designed for the model

The developer of the model probably made a suite of other functions to help you probe the model object. Check the documentation to find these.

For example, I can find the regression coefficients of a lm() object with coef():

coef(my_lm)
## (Intercept)          wt 
##   37.285126   -5.344472

Standard error of the residuals with sigma():

sigma(my_lm)
## [1] 3.045882

Don’t be surprised if there is no function to extract what you want. If that’s the case, read on…

Manually Extracting Information

If you can’t extract model information from built-in functions like coef() or sigma(), you can manually dig in to the model object. In most cases, a model object is just a list in disguise! (Lists are discussed in STAT 545B). You can therefore extract information like you would with any other list.

To help figure out what’s in the list, use the names() function. Here are all the entries in the lm() object we fit earlier:

names(my_lm)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

Want the degrees of freedom of the residuals? Just extract that entry:

my_lm$df.residual
## [1] 30

To complicate things further this might only be part of the information made available in the model object! More info might be stored in yet another list after applying the summary() function:

summary(my_lm)
## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Like the original model object, summary() looks like it returns a bunch of text – but it’s actually another list! Let’s again use names() to get the components of this list:

names(summary(my_lm))
##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

Now I can get more pieces of information, like the adjusted R-squared value:

summary(my_lm)$adj.r.squared
## [1] 0.7445939