\(~\)

Fitting, Interpreting, and Evaluating Linear Regression Models

As seen in our lecture, the lm() function is the first step in a linear regression analysis, and we can then use the summary() function on the fitted model to obtain additional information, such as model’s slope and intercept, hypothesis testing results on the slope and intercept, and \(R^2\). The code below briefly reviews this workflow:

## Load some data
grinnell_peers = read.csv("https://remiller1450.github.io/data/grinnell_peers.csv")

## Fit and store the lm model
my_model = lm(Med_Grad_Earnings_4Y ~ Ret_Rate, data = grinnell_peers)

## Summarize the fitted model
summary(my_model)
## 
## Call:
## lm(formula = Med_Grad_Earnings_4Y ~ Ret_Rate, data = grinnell_peers)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -7042  -3586  -1369   3596  11328 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -121877      37607  -3.241 0.005119 ** 
## Ret_Rate      189541      40408   4.691 0.000246 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5306 on 16 degrees of freedom
## Multiple R-squared:  0.579,  Adjusted R-squared:  0.5527 
## F-statistic:    22 on 1 and 16 DF,  p-value: 0.0002455

A few details to clarify here:

  1. The model is specified using a formula, which you can read as “outcome variable ~ explanatory variable”
  2. We will focus on “Multiple R-squared” in this lab. Adjusted R-squared is primarily useful when comparing models that use different numbers of explanatory variables. Since simple linear regression only uses one explanatory variable, there’s no reason to worry about Adjusted R-squared.
  3. The slope of this model is 189541, which is difficult to interpret because it reflects a 1-unit change in Ret_Rate, ranges from 0.8656 to 0.9632. We can find the expected effect of a 0.01-unit change in Ret_Rate, which is more realistic in this application, by dividing the slope by 100 (so 189541/100 = 1895.41)
  4. The intercept of this model is -121877, which is meaningless in this application because no colleges will have a retention rate of exactly 0.

Question #1: For this question you should use the Happy Planet dataset introduced in our previous lab.

happy_planet <- read.csv("https://remiller1450.github.io/data/HappyPlanet.csv")
  • Part A: Fit a linear regression model using a country’s life expectancy to predict its happiness score, storing your fitted model in an object named life_model. Use the summary function to print a summary of this fitted model.
  • Part B: In statistical modeling we work with a population-level model and a fitted model. Consider the two model equations, A and B, found below. Which of these is the population-level model equation, and which is the fitted model equation? In which equation do the slope and intercept from the summary you found in Part A belong?
    • A: \(\text{Happiness} = \beta_0 + \beta_1 \cdot \text{Life Expectancy} + \epsilon\)
    • B: \(\widehat{\text{Happiness}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \text{Life Expectancy}\)
  • Part C: Use the hypothesis test of \(H_0: \beta_1 = 0\) and the estimated slope to write a 1-sentence conclusion about how life expectancy relates to happiness scores across the countries found in this data set. Your conclusion should address the strength of evidence for these variables being related, and if there is evidence of relationship it should describe the relationship that was found in the sample data.
  • Part D: Use the fitted model to predict the happiness score of a country with a life expectancy of 75 years. Does every country with a life expectancy of 75 years have this happiness score? Briefly explain.
  • Part E: In many behavioral/social sciences, \(R^2\) values tend to be low, with values exceeding 0.09 being described as “medium effects”, and values exceeding 0.25 being described as “large effects” (reference). Based upon these guidelines, does life expectancy appear to explain a large amount of variability in country happiness scores?

\(~\)

Residuals and Model Diagnostics

The \(T\)-test results you see in the output produced by the summary() function are only valid if the test statistic actually follows a \(t\)-distribution under the null hypothesis, which will only occur when a certain conditions are met. These conditions can be summarized as:

  • Model errors are independent, and follow identical Normal distributions (same standard deviation/spread, regardless of the \(X\) value)

In practice, this will occur when there are no problematic patterns in the residuals (the model errors we observe in the sample data). There are two main types of residual plots we should check:

## Fitted model
my_model = lm(Med_Grad_Earnings_4Y ~ Ret_Rate, data = grinnell_peers)

## Residual vs. fitted and QQ-plot
plot(my_model, which = c(1,2))

  • The Residuals vs. Fitted plot helps us assess the independence and constant spread of the observed errors. We are looking for no patterns, and similar amounts of variability everywhere along the X-axis.
  • The QQ-plot (shorthand for quantile-quantile plot) shows how each percentile of the observed residuals compares to where it theoretically should be if the residuals followed a Normal distribution.
    • We would like this to be a perfect 45-degree line, suggesting the quantiles are perfectly aligned
    • If the residuals are right-skewed, the points will bend upwards towards the right side of the plot (suggesting a longer right tail than a Normal distribution should have)
    • If the residuals are left-skewed, the points will bend downward towards the left side of the plot (suggesting a longer left tail than a Normal distribution should have)

In the example shown above, there are no problematic findings in either of these diagnostic plots, as the Residuals vs. Fitted plot shows seemingly random scatter and the QQ-plot stays reasonably close to the 45-degree line.

Question #2:

  • Part A: Using the fitted model from Question #1. Display the two diagnostic plots shown in this section’s example. Note that these diagnostic plots suggest the conditions for statistical testing using the model are reasonable. With real world data, we do not ever expect to see perfection in these plots, and the deviations in the upper right of the QQ-plot are relatively small for real data.
  • Part B: Fit a model using GDPperCapita to predict Happiness and display the two diagnostic plots covered in this section. What is problem(s) do you see in the Residuals vs. Fitted plot? Briefly explain.

\(~\)

Accommodating Non-linear Relationships

One potential cause of a pattern in the Residuals vs. Fitted plot is a non-linear relationship between an explanatory variable and the outcome. For the model using GDPperCapita to predict Happiness we can see signs of non-linearity in a scatterplot of these variables:

library(ggplot2)
ggplot(happy_planet, aes(x = GDPperCapita, y = Happiness)) + geom_point()

There are many ways to introduce non-linear effects into a regression model, but in this class we will only consider quadratic effects:

Feature Simple Linear Regression Quadratic Regression
Population Model \(Y = \beta_0 + \beta_1X + \epsilon\) \(Y = \beta_0 + \beta_1X + \beta_2X^2 + \epsilon\)
Fitted Model \(\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1X\) \(\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1X + \hat{\beta}_2X^2\)
R Syntax lm(Y ~ X, data = my_data) lm(Y ~ X + I(X^2), data = my_data)

The quadratic model will find the values of \(\hat{\beta}_0\) (intercept), \(\hat{\beta}_1\) (linear coefficient), and \(\hat{\beta}_2\) (quadratic coefficient) that minimize the squared residual error for the sample data. Or, in other words, the best fitting quadratic curve for the data.

## Example quadratic model
quad_model = lm(Med_Grad_Earnings_4Y ~ Ret_Rate + I(Ret_Rate^2), data = grinnell_peers)
summary(quad_model)
## 
## Call:
## lm(formula = Med_Grad_Earnings_4Y ~ Ret_Rate + I(Ret_Rate^2), 
##     data = grinnell_peers)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -6482  -3530  -1626   4169  11264 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)     396616    1353641   0.293    0.774
## Ret_Rate       -940530    2949393  -0.319    0.754
## I(Ret_Rate^2)   614964    1604846   0.383    0.707
## 
## Residual standard error: 5453 on 15 degrees of freedom
## Multiple R-squared:  0.5831, Adjusted R-squared:  0.5275 
## F-statistic: 10.49 on 2 and 15 DF,  p-value: 0.001415

One challenge posed by this model is that we no longer have a simple interpretation of the effect of Ret_Rate on Med_Grad_Earnings_4Y. We can use the library visreg to visualize the fitted model:

# install.packages('visreg')
library(visreg)

## Visualize the model predictions using Ret_Rate on the x-axis 
visreg(quad_model, xvar = "Ret_Rate", gg = TRUE)

Here we can see that the prediction line gradually gets steeper as Ret_Rate increases. If we wanted to quantitatively understand the relationship between Ret_Rate and Med_Grad_Earnings_4Y in the fitted model, we have two options:

  1. Find the instantaneous slope at certain values of Ret_Rate that are interesting
  2. Find the minimum (or maximum) of the quadratic curve as a “turning point” where the effect switches from positive to negative

For the first approach (instantaneous slope), we can rely on the fact that the derivative of \(\hat{\beta}_0 + \hat{\beta}_1X + \hat{\beta}_2X^2\) with respect to \(X\) is linear:

\[\text{Slope at } X = a \text{ is given by }\hat{\beta}_1 + 2\hat{\beta}_2a\]

So in our example model the instantaneous slope at retention rate of 0.96 is \(-940530 + 2\cdot614964\cdot0.96 = 240200.9\), which is slightly steeper than the slope of \(189541\) in the simple linear regression model (Med_Grad_Earnings_4Y ~ Ret_Rate).

For the second approach, we can also rely on calculus to give us the minimum/maximum of a quadratic function, which will occur at: \[x_{min/max} = \frac{-\hat{\beta}_1}{2\hat{\beta}_2}\]

For our example quadratic model this is \(\frac{940530}{2 \cdot 614964} = 0.765\). So, while such a retention rate is possible, it is also outside the range of our observed data, so we wouldn’t want to report it as a point of interest.

Question #3:

  • Part A: Fit a quadratic model using GDPperCapita (and its squared term) to predict Happiness. Use the summary() function to compare the \(R^2\) of this model to simple linear regression model you fit using these variables in Part B of Question 2. Why model explains more variability in the outcome?
  • Part B: Use visreg() to visualize the quadratic model. Does the quadratic curve show an increasing or decreasing slope for higher values of GDPperCapita?
  • Part C: Find the instantaneous slope at GDPperCapita values of 5000 (representing lower income countries) and 40000 (representing higher income countries). How do these values compare? What do they tell you about the relationship between GDPperCapita and Happiness in the fitted model?
  • Part D: Find the value of GDPperCapita where the estimated effect is expected to turn negative. Is this value meaningful in this application? Or is it a problematic extrapolation?
  • Part E: Using the fitted model from Part A. Display the two diagnostic plots discussed in the previous section of the lab. Do the modeling assumptions seem better satisfied for the quadratic model used in this analysis than for the simple linear regression you looked at in Part B of Question 2? Briefly explain.

\(~\)

Application

For this application you’ll use the Iowa City home sales dataset provided below. These data were assembled by a University of Iowa professor looking to buy a home in 2009, and they include all homes sold between January 2005 and December 2008.

ic_homes = read.csv("https://remiller1450.github.io/data/IowaCityHomeSales.csv")

Question #4: Our analysis of these data will use sale.amount as the outcome, and consider bedrooms and area.living as two possible predictors of sale price.

  • Part A: Use the select() function from the dplyr library to create a version of these data including only the three variables mentioned above. Then, use the pairs.panels() function from the psych library to construct a scatterplot/correlation matrix (see our previous lab for an example). Which predictor, bedrooms or area.living, seems more strongly associated with sale.amount? Does the stronger predictor appear to have an approximately linear relationship with sale.amount?
  • Part B: Fit a simple linear regression model using the variable you identified in Part A to predict the outcome, sale.amount. Then use the summary() function to print a summary of the fitted model.
  • Part C: How much variation in sale.amount is explained by the fitted model from Part B? Does this suggest the variable you used is a good predictor of `sale.amount?
  • Part D: Provide a 1-sentence conclusion for the hypothesis test for the slope coefficient provided by the summary() function. You should include an interpretation of the slope coefficient itself.
  • Part E: Using your fitted model, display the two diagnostic plots discussed in the previous section of the lab. Do the assumptions of the model seem appropriate? Briefly explain what you see.
  • Part F: Now consider a quadratic model involving the predictor you identified in Part A. Fit this model and report the \(R^2\) value. Does it seem to explain more variation in sale.amount than the simple linear regression model you fit earlier? Is the difference large enough for you to favor this model? Briefly explain.