R - Multiple Regression

R - Multiple Linear Regression

Multiple Linear Regression is an extension of simple linear regression used to model the relationship between a single continuous response variable and two or more predictor variables. It allows us to assess the joint effect of several predictors on an outcome.

The general mathematical equation for multiple regression is:

y = β₀ + β₁x₁ + β₂x₂ + ... + βₚxₚ + ε

y is the response variable.
β₀ is the intercept.
β₁, β₂, … βₚ are the coefficients for the predictor variables.
x₁, x₂, … xₚ are the predictor variables.
ε is the random error term.

In R, we create the regression model using the lm() function. The model estimates the values of the coefficients (the β values) from the input data. Once the model is built, we can use it to predict the response for new data and evaluate the significance of the predictors.

The `lm()` Function

This function creates a linear model based on the relationship between predictor(s) and a response variable.

Syntax

lm(formula, data)

formula: An object of class "formula" specifying the model. For multiple regression, it looks like response ~ predictor1 + predictor2 + ....
data: A data frame containing the variables in the model.

Comprehensive Example: Predicting Car Mileage

We’ll use the built-in mtcars dataset to model a car’s mileage (mpg) based on its displacement (disp), horsepower (hp), and weight (wt).

Step 1: Load Packages and Prepare Data

Using the dplyr package from the tidyverse makes data selection clean and readable.

# Install tidyverse if you haven't already: install.packages("tidyverse")
library(dplyr)

# Select the relevant columns from the mtcars dataset
input_data <- mtcars %>%
  select(mpg, disp, hp, wt)

# Display the first few rows of our dataset
print(head(input_data))

The prepared data looks like this:

                   mpg disp  hp    wt
Mazda RX4         21.0  160 110 2.620
Mazda RX4 Wag     21.0  160 110 2.875
Datsun 710        22.8  108  93 2.320
Hornet 4 Drive    21.4  258 110 3.215
Hornet Sportabout 18.7  360 175 3.440
Valiant           18.1  225 105 3.460

Step 2: Create and Evaluate the Model

We build the model with lm() and then use summary() to get a detailed analysis, which is the standard practice.

# Create the relationship model
model <- lm(mpg ~ disp + hp + wt, data = input_data)

# Get a comprehensive summary of the model
summary(model)

The summary() output is rich with information:

Call:
lm(formula = mpg ~ disp + hp + wt, data = input_data)

Residuals:
    Min      1Q  Median      3Q     Max
-3.8924 -1.6421 -0.1785  1.0954  5.5494

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
disp        -0.000937   0.010350  -0.091  0.92851
hp          -0.031157   0.011355  -2.744  0.01049 *
wt          -3.800891   1.066178  -3.565  0.00132 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.61 on 28 degrees of freedom
Multiple R-squared:  0.8268,  Adjusted R-squared:  0.8083
F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

Interpreting the Summary:

Coefficients (Estimate): These are the β values. For example, for every 1 unit increase in weight (wt), mpg is predicted to decrease by 3.8, holding other variables constant.
Pr(>|t|) (p-value): This tells us if a predictor is statistically significant. A small p-value (e.g., < 0.05) suggests the predictor has a meaningful relationship with the response. Here, hp and wt are significant, while disp is not in the presence of the other two.
Adjusted R-squared: This indicates that our model explains approximately 80.8% of the variance in mpg, which is a strong fit.

Step 3: Tidying Model Output with `broom`

The broom package is a modern tool for converting model objects into tidy data frames, making them easier to plot and report.

# Install broom if you haven't already: install.packages("broom")
library(broom)

# Get a tidy data frame of model coefficients
tidy_model <- tidy(model)
print(tidy_model)

# A tibble: 4 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept) 37.1      2.11        17.6    1.78e-16
2 disp        -0.000937 0.0104      -0.0905 9.29e- 1
3 hp          -0.0312   0.0114      -2.74   1.05e- 2
4 wt          -3.80     1.07        -3.57   1.32e- 3

Step 4: Predicting New Values

The standard and safest way to predict is with the predict() function, which applies the model equation to new data.

# Create a data frame with new car specifications
new_cars <- data.frame(
  disp = c(221, 250),
  hp = c(102, 110),
  wt = c(2.91, 3.20)
)

# Predict the MPG for the new cars
predicted_mpg <- predict(model, newdata = new_cars)

# Show the predictions
print(predicted_mpg)

The model predicts the mileage for the new car data:

       1        2
22.71021 21.43889