Load Data

In this dataset, each row is a diamond!

library(tidyverse) # for data wrangling
library(viridis) # for colors
library(broom) # for regression
library(moderndive) # for regression

# Save the diamonds dataset as an object in my environment
mydiamonds <- read_csv("mydiamonds.csv")

View Data

# View first 3 rows of dataset
mydiamonds %>% head(3)

Codebook

In this dataset, our variables mean:

• price: price of diamond in US dollars (from $326 to $18,823!)

• carat: weight of the diamond (0.2 to 5.01 carats)

• cut: quality of the cut (Fair, Good, Very Good, Premium, Ideal)

Predicted Values

All models have a lot of information stored inside them, which you can access using the $ sign. (Note: select() can’t extract them, because model objects are not data.frames)

# This code gives you the original values put into your model (y ~ x...)
# We're using head() here to show just the first few rows.
m$model %>% head(3)

This code gives you the predicted values for your outcome, given each row of data (\(Y_{Predicted}\)).

m$fitted.values %>% head(3)

##        1        2        3 
## 6908.982 2524.950 6455.461

This code gives you the residuals for your outcome, given each row of data (\(Y_{Observed} − Y_{Predicted}\)).

m$residuals %>% head(3)

##          1          2          3 
## -2312.9817  -590.9498 -1615.4611

Last, we can use the moderndive package to combine all these into a nice data.frame, where price_hat means \(Y_{predicted}\) / fitted.values from above.

library(moderndive)
mdat <- m %>%
  get_regression_points()
# Let's view just the first few rows
mdat %>% head(3)

mdat <- data.frame(
  # Put the model variables in columns with their own variable names
  m$model,
  # Put the fitted values in a column called predicted
  price_hat = m$fitted.values,
  # Put the residuals in a column called residuals
  residuals = m$residuals)
# Look at first few rwos
mdat %>% head(3)

Residuals

Residuals represent the difference between the observed outcome and the outcome predicted by the model. A regression model finds the line which minimizes these residuals, thus getting the best “model” for the overall trend in the data.

The animation below below visualizes residuals as lines stemming from the best fit line. The alpha and beta coefficients are varied to show changing residuals for each of these different lines, compared to the blue line, which is the actual line of best fit.

• Cases well predicted by the model are tiny, and close to the best fit line.

• Cases poorly predicted by the model are BIG, and far from the best fit line.

# Visualize residuals, pulling from our new data.frame mdat
mdat %>%
  ggplot(mapping = aes(
    # Make x axis diamond carats,
    x = carat,
    # Make y axis price of diamonds (outcome)
    y = price,
    # Make the size of points how far away they are from predicted value
    size = abs(residuals),
    # Draw a line whose starting point is the predicted value (line)
    ymin = price_hat,
    # And whose end point is the observed value (dot)
    ymax = price)) +
  # Now draw the points
  geom_point(color = "darkgrey", alpha = 0.5) +
  # Draw the line of best fit
  geom_smooth(method = "lm", se = FALSE) +
  # Draw the 'residuals', the distance between predicted and observed
  geom_linerange(color = "steelblue", size = 0.5, alpha = 0.5) +
  # Add a classic theme
  theme_classic() +
  # Remove the size legend
  guides(size = "none")

R-squared

R-squared (\(R^{2}\)) is a simple statistic that ranges from 0 - 1. R2 measures the percentage of variation in the outcome that is explained by the model.

# Using our data.frame of model inputs and outputs,
# let's create two summary statistics
mdat %>%
  summarize(
    # get sum of each squared difference between observed and mean price
    tss = sum((price - mean(price))^2),
    # get sum of each squared difference between observed and predicted price
    rss = sum((price - price_hat)^2)) %>%
  # calculate R2 based on these stats
  mutate(R2 = 1 - rss / tss)

For example, the following animation shows how each of the possible lines plotted above produces a different residual sum of squares, leading to a different \(R^{2}\). The sweet spot, where \(R^{2}\) is highest (and therefore the residuals are minimized) is when the slope is closest to our actual observed beta value, $7559 per carat. Otherwise, both higher and lower slopes lead to a lower \(R^{2}\).

F-statistic

Finally, we also want to measure how useful this model is, compared to the intercept. Models are imperfect approximations of trends in data. The simplest possible model of data is the intercept line, the amount of the outcome you would have if X were zero. (Imagine a flat line through your data at the intercept.) If we have a good model, it had better explain more variation than the amount explained by the intercept line, right? To do this, we can calculate an F-statistic for our model.

• F statistics (just like Chi-squared) range from 0 to infinity.

• If your F statistic is small, the model ain’t much better than the intercept.

• If your F statistic is large, the model explains much, much more variation than the intercept.

• We can compare your F statistic to a null distribution of scores we’d get due to chance, and find the probability we got this statistic due to chance, our p-value. The broom package’s glance function lets us do this below, giving us a statistic and p.value.

m %>% 
  glance()

ingredients <- m %>%
  get_regression_points() %>%
  summarize(
    # Calculate residual sum of squares
    residual_sum_of_squares = sum( (price - price_hat) ^2),
    # Calculate total sum of squares
    total_sum_of_squares = sum( (price - mean(price))^2),
    # Get sample size
    n = n(),
    # Calculate number of variables in your model
    # Here, I grabbed the model, asked it to give us the column names,
    # and calculated the length of that vector (2 names, so p = 2)
    p = length(c("price", "carat") )) %>%
  # Calculate explained sum of squares
  mutate(explained_sum_of_squares = total_sum_of_squares - residual_sum_of_squares)

# View our result!
ingredients %>% glimpse()

## Rows: 1
## Columns: 5
## $ residual_sum_of_squares  <dbl> 2321500044
## $ total_sum_of_squares     <dbl> 14222520095
## $ n                        <int> 1000
## $ p                        <int> 2
## $ explained_sum_of_squares <dbl> 11901020051

ingredients %>%
  mutate(
  # Mean Squares due to Regression, given the no. of predictors
  mean_squares_due_to_regression = explained_sum_of_squares / (p - 1),
  # Mean Squared Error
  # How much variation was NOT explained, given the sample size and no. of variables
  mean_squared_error = residual_sum_of_squares / (n - p)) %>%
  # Compute the F-statistic, which is a ratio of explained vs. unexplained variation
  mutate(f_statistic = mean_squares_due_to_regression / mean_squared_error) %>%
  # Finally, throw it into this pf function,
  # which plots a theoretical null distribution 
  # based on the number of variables and sample size,
  # and identifies how extreme our F-statistic is compared to one we'd get by chance
  # Computes p-value
  mutate(p.value = pf(f_statistic, df1 = p - 1, df2 = n - p, lower.tail = FALSE)) %>%
  
  # View output
  glimpse()

## Rows: 1
## Columns: 9
## $ residual_sum_of_squares        <dbl> 2321500044
## $ total_sum_of_squares           <dbl> 14222520095
## $ n                              <int> 1000
## $ p                              <int> 2
## $ explained_sum_of_squares       <dbl> 11901020051
## $ mean_squares_due_to_regression <dbl> 11901020051
## $ mean_squared_error             <dbl> 2326152
## $ f_statistic                    <dbl> 5116.183
## $ p.value                        <dbl> 0