R - Linear Regression

Modern Linear Regression in R

Linear regression is a foundational statistical method used to model the relationship between a dependent variable (response) and one or more independent variables (predictors). It aims to find the best-fitting straight line through the data points.

The general equation for a simple linear regression is: y = β₀ + β₁x + ε, where y is the response, x is the predictor, β₀ is the intercept, β₁ is the slope, and ε is the error term. R’s core function lm() is the standard for this, but we’ll enhance our workflow with the tidyverse and broom packages for a modern approach.

The Workflow for Regression Analysis

A typical regression analysis project follows these steps:

• 1. Data Preparation: Collect and structure your data into a data frame.
• 2. Exploration: Visualize the relationship between variables (e.g., with a scatter plot).
• 3. Model Fitting: Use the lm() function to create the regression model.
• 4. Model Interpretation: Analyze the model’s summary to understand coefficients and model fit.
• 5. Model Diagnostics: Check the model’s assumptions to ensure its validity.
• 6. Prediction: Use the trained model to predict outcomes for new data.

Example: Predicting Weight from Height

Let’s model the relationship between a person’s height (predictor) and weight (response).

# Load the necessary modern packages
# install.packages(c("tidyverse", "broom"))
library(tidyverse)
library(broom)

# Step 1: Prepare the data in a tibble (a modern data frame)
height_weight_data <- tibble(
  height_cm = c(151, 174, 138, 186, 128, 136, 179, 163, 152, 131),
  weight_kg = c(63, 81, 56, 91, 47, 57, 76, 72, 62, 48)
)

# Step 2: Explore the data with a scatter plot using ggplot2
ggplot(height_weight_data, aes(x = height_cm, y = weight_kg)) +
  geom_point() +
  labs(
    title = "Relationship between Height and Weight",
    x = "Height (cm)",
    y = "Weight (kg)"
  ) +
  theme_minimal()

The scatter plot should show a positive linear trend, suggesting that as height increases, weight tends to increase as well.

Step 3 & 4: Fit and Interpret the Model

We use lm() to fit the model. The formula weight_kg ~ height_cm reads as “model weight as a function of height.”

# Step 3: Fit the linear model
model <- lm(weight_kg ~ height_cm, data = height_weight_data)

# Step 4: Interpret the model with summary()
summary(model)

The summary() output is rich with information:

Call:
lm(formula = weight_kg ~ height_cm, data = height_weight_data)

Residuals:
    Min      1Q  Median      3Q     Max
-3.9775 -1.8944 -0.0412  1.6629  6.3002

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -38.45509    8.04901  -4.778  0.00139 **
height_cm     0.67461    0.05191  12.997 1.16e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.253 on 8 degrees of freedom
Multiple R-squared:  0.9548,    Adjusted R-squared:  0.9491
F-statistic: 168.9 on 1 and 8 DF,  p-value: 1.164e-06

Key Takeaways from the Summary:

• Coefficients: The intercept (-38.45) and height_cm (0.67) give us the regression equation: weight_kg = -38.45 + 0.67 * height_cm.
• P-value (Pr(>|t|)): The very small p-value for height_cm (1.16e-06) indicates it is a statistically significant predictor of weight.
• Multiple R-squared: 0.9548 means that approximately 95.5% of the variability in weight can be explained by height in our model, which is a very strong fit.

Modern Interpretation with broom

The broom package can tidy the model output into a data frame, which is excellent for reporting and further analysis.

# Get a clean tibble of model coefficients
tidy(model)

# Get a tibble of model-level statistics
glance(model)

Step 5: Visualize the Model Fit

ggplot2 makes it easy to overlay the regression line on our scatter plot.

ggplot(height_weight_data, aes(x = height_cm, y = weight_kg)) +
  geom_point(color = "blue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "red") + # Add the regression line
  labs(
    title = "Height & Weight Regression",
    x = "Height in cm",
    y = "Weight in Kg"
  ) +
  theme_bw()

Step 6: Predict New Values

Use the predict() function to estimate the weight for a person with a new height. It’s best practice to provide the new data in a data frame with the same column name as the predictor.

# Create a data frame for the new observation
new_height <- tibble(height_cm = 170)

# Predict the weight
predicted_weight <- predict(model, newdata = new_height)
print(predicted_weight)

This will output the predicted weight for a person of 170cm height.

       1
76.22869

Best Practices and Next Steps

• Model Diagnostics: A critical step is to check the assumptions of linear regression (linearity, independence, normality of residuals, and equal variance). You can do this by plotting the model object: par(mfrow = c(2, 2)); plot(model). The augment() function from broom is also excellent for creating data frames for diagnostic plots.
• Multiple Regression: You can add more predictors to the model using the + sign in the formula, e.g., lm(y ~ x1 + x2, data = df).
• Team Collaboration: Use version control like Git to track changes to your R scripts. Use a tool like renv to manage project-specific package dependencies, ensuring your analysis is reproducible by others.