Lab 8: Introduction to linear regression

Getting Started

library(tidyverse)
library(openintro)

The data

Exercise 1

What are the dimensions of the dataset?

dim(hfi)

## [1] 1458  123

The dataset has 1458 rows and 123 columns.

Exercise 2

What type of plot would you use to display the relationship between the personal freedom score, pf_score, and one of the other numerical variables? Plot this relationship using the variable pf_expression_control as the predictor. Does the relationship look linear? If you knew a country’s pf_expression_control, or its score out of 10, with 0 being the most, of political pressures and controls on media content, would you be comfortable using a linear model to predict the personal freedom score?

ggplot(hfi, aes(pf_expression_control, pf_score)) +
  geom_point() +
  labs(title = "PF Score Predicted by PF Expression Control")

## Warning: Removed 80 rows containing missing values or values outside the scale range
## (`geom_point()`).

A scatter plot effectively shows the relationship between two numerical variables, helping us assess whether the relationship is linear or nonlinear.

If the relationship looks linear, we can quantify the strength of the relationship with the correlation coefficient.

hfi %>%
  summarise(cor(pf_expression_control, pf_score, use = "complete.obs"))

## # A tibble: 1 × 1
##   `cor(pf_expression_control, pf_score, use = "complete.obs")`
##                                                          <dbl>
## 1                                                        0.796

Exercise 3

Looking at your plot from the previous exercise, describe the relationship between these two variables. Make sure to discuss the form, direction, and strength of the relationship as well as any unusual observations.

This is a positive linear relationship: as the expression control value increases, the score also tends to rise. The relationship is linear in form, with a positive direction and a strong correlation of 0.796.

# This will only work interactively (i.e. will not show in the knitted document)
hfi <- hfi %>% filter(complete.cases(pf_expression_control, pf_score))

Exercise 4

Using plot_ss, choose a line that does a good job of minimizing the sum of squares. Run the function several times. What was the smallest sum of squares that you got? How does it compare to your neighbors?

df <- hfi %>% select(pf_expression_control, pf_score) %>% drop_na()

m1 <- lm(pf_score ~ pf_expression_control, data = hfi)

The first argument in the function lm is a formula that takes the form y ~ x. Here it can be read that we want to make a linear model of pf_score as a function of pf_expression_control. The second argument specifies that R should look in the hfi data frame to find the two variables.

The output of lm is an object that contains all of the information we need about the linear model that was just fit. We can access this information using the summary function.

summary(m1)

## 
## Call:
## lm(formula = pf_score ~ pf_expression_control, data = hfi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8467 -0.5704  0.1452  0.6066  3.2060 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.61707    0.05745   80.36   <2e-16 ***
## pf_expression_control  0.49143    0.01006   48.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8318 on 1376 degrees of freedom
## Multiple R-squared:  0.6342, Adjusted R-squared:  0.634 
## F-statistic:  2386 on 1 and 1376 DF,  p-value: < 2.2e-16

Exercise 5

Fit a new model that uses pf_expression_control to predict hf_score, or the total human freedom score. Using the estimates from the R output, write the equation of the regression line. What does the slope tell us in the context of the relationship between human freedom and the amount of political pressure on media content?

m1 <- lm(pf_score ~ pf_expression_control, data = hfi)
summary(m1)

## 
## Call:
## lm(formula = pf_score ~ pf_expression_control, data = hfi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8467 -0.5704  0.1452  0.6066  3.2060 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.61707    0.05745   80.36   <2e-16 ***
## pf_expression_control  0.49143    0.01006   48.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8318 on 1376 degrees of freedom
## Multiple R-squared:  0.6342, Adjusted R-squared:  0.634 
## F-statistic:  2386 on 1 and 1376 DF,  p-value: < 2.2e-16

m2 <- lm(hf_score ~ pf_expression_control, data = hfi)
summary(m2)

## 
## Call:
## lm(formula = hf_score ~ pf_expression_control, data = hfi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6198 -0.4908  0.1031  0.4703  2.2933 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           5.153687   0.046070  111.87   <2e-16 ***
## pf_expression_control 0.349862   0.008067   43.37   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.667 on 1376 degrees of freedom
## Multiple R-squared:  0.5775, Adjusted R-squared:  0.5772 
## F-statistic:  1881 on 1 and 1376 DF,  p-value: < 2.2e-16

Prediction and prediction errors

ggplot(data = hfi, aes(x = pf_expression_control, y = pf_score)) +
  geom_point() +
  stat_smooth(method = "lm", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

Exercise 6

If someone saw the least squares regression line and not the actual data, how would they predict a country’s personal freedom school for one with a 6.7 rating for pf_expression_control? Is this an overestimate or an underestimate, and by how much? In other words, what is the residual for this prediction?

The personal freedom score is 7.9

Model diagnostics

ggplot(data = m1, aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Fitted values") +
  ylab("Residuals")

Exercise 7

Is there any apparent pattern in the residuals plot? What does this indicate about the linearity of the relationship between the two variables?

ggplot(data = m1, aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Fitted values") +
  ylab("Residuals")

The residual plot shows a random distribution of points, suggesting that a linear regression model would be appropriate for this dataset.

Nearly normal residuals

ggplot(data = m1, aes(x = .resid)) +
  geom_histogram(binwidth = 25) +
  xlab("Residuals")

normal probability plot of the residuals.

ggplot(data = m1, aes(sample = .resid)) +
  stat_qq()

Exercise 8

Based on the histogram and the normal probability plot, does the nearly normal residuals condition appear to be met?

The histogram and probability plot of the residuals display a normal distribution, indicating that the normality condition is satisfied.

Exercise 9

Based on the residuals vs. fitted plot, does the constant variability condition appear to be met?

Yes, the constant variance condition is satisfied. The variability remains consistent in both the positive and negative directions along the x-axis.

More practice

Choose another freedom variable and a variable you think would strongly correlate with it. Produce a scatterplot of the two variables and fit a linear model. At a glance, does there seem to be a linear relationship?

hfi %>% ggplot(aes(x = pf_religion_restrictions , y = pf_religion_harassment)) +
  geom_point()

## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

Based on a preliminary look, there doesn’t appear to be a linear relationship, as the points are dispersed around the 9 and 10 range.

lm_model <- lm(pf_religion_harassment ~ pf_religion_restrictions, data = hfi)
summary(lm_model)

## 
## Call:
## lm(formula = pf_religion_harassment ~ pf_religion_restrictions, 
##     data = hfi)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.76723 -0.41159  0.06432  0.47713  2.03376 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               6.28916    0.08628   72.89   <2e-16 ***
## pf_religion_restrictions  0.34151    0.01151   29.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7272 on 1362 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.3926, Adjusted R-squared:  0.3921 
## F-statistic: 880.2 on 1 and 1362 DF,  p-value: < 2.2e-16

From the above, it appears that the first model has a significantly stronger correlation.