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:
- how a model is typically fit in R;
- how to extract information from the model after fitting it, using the
- extract information by other means if the
broompackage 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.
- A car weighs 4000 lbs. Provide a numerical prediction on its mpg (miles per gallon).
- Does the weight of a car influence its mpg?
- 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
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:
##  21.25171
We can find the p-value testing the null hypothesis that the true slope of the regression line is 0:
##  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:
##  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)
A function such as:
- Linear Regression:
- Generalized Linear Regression:
- Local regression:
- Quantile regression:
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
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
(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
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() function for a statistical summary of each component of your model, where each component gets a row in the output tibble. For
tidy() gives one row per regression coefficient (slope and intercept).
## # 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() 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() 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.
## # 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
Here’s how to work with these model objects in that case.
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)
##  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()
##  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
- For a model fit with
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
## (Intercept) wt ## 37.285126 -5.344472
Standard error of the residuals with
##  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
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:
##  "coefficients" "residuals" "effects" "rank" ##  "fitted.values" "assign" "qr" "df.residual" ##  "xlevels" "call" "terms" "model"
Want the degrees of freedom of the residuals? Just extract that entry:
##  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
## ## 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:
##  "call" "terms" "residuals" "coefficients" ##  "aliased" "sigma" "df" "r.squared" ##  "adj.r.squared" "fstatistic" "cov.unscaled"
Now I can get more pieces of information, like the adjusted R-squared value:
##  0.7445939