Load the Cars Dataset

First, let’s load the built-in cars dataset to a native R dataframe

df <- as.data.frame(cars)

head(df)
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10

Model Creation

We can use R’s built-in linear model (lm) to create

model <- lm(cars$dist ~ cars$speed)

summary(model)
## 
## Call:
## lm(formula = cars$dist ~ cars$speed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## cars$speed    3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Model Analysis

There’s 3 main criteria we need to consider when creating (and evaluating) a linear model:

  1. Linear relationship - the dependent variable should have a roughly linear relationship to an independent variable
  2. Homoscedasticity - the variance of our residuals as a function of our residuals as a function of the independent variable is constant
  3. Independence of Errors - there’s no pattern/correlation between residuals

Let’s plot our linear model using the geom_abline functino from ggplot2

# Plotting our datapoints overlaid with our linear model
c <- coef(model)
ggplot(df, aes(x=speed, y=dist)) + geom_point() + 
  geom_abline(slope = c[["cars$speed"]], 
              intercept =c[["(Intercept)"]])

This plot looks roughly linear, so our condition of a linear relationship is met.

Next, let’s plot the distribution of residuals of our model. It’s important that our residuals are normally distributed, as that’s a main criterion for simple linear regression, as it helps us check the condition of homoscedasticity

res <- resid(model) #alternatively, model$residuals
hist(res)

The above histogram looks roughly normal by the eye test. We can also plot our residuals as a function of our model’s predicted values (\(\hat{y}\))

plot(fitted(model), res)

There’s no clear pattern in this plot, though we should verify further that our residuals are normally distributed.

We should also generate a Q-Q plot to help in determining if our residuals follow a normal distribution. In this framework, the plotted residuals should be linear, which appears to be the case with the exception of the right tail of the distribution.

# Create Q-Q plot of our residuals
qqnorm(res)
qqline(res)

We can also produce a panel plot of diagnostic charts as in our textbook example

par(mfrow=c(2,2))
plot(model)