Linear Regression in a Prediction Context

STA 6543: Predictive Modeling

Overview

We begin by shifting our focus from statistical inference to predictive modeling. Our primary concern changes from “Is variable \(X\) associated with variable \(Y\)?” to “How accurately can we estimate future values of \(Y\) given an influx of new, unseen data \(X\)?”

By the end of this hour, you will be able to:

1. Split data into training and testing partitions using modern tidyverse-adjacent frameworks (rsample).

2. Construct, interpret, and evaluate a multiple linear regression model using lm().

3. Evaluate a model’s predictive performance using quantitative metrics (RMSE, MAE, \(R^2\)).

4. Critically assess core statistical assumptions to protect our model from overfitting and operational failure.

What is R?

Before we build complex predictive frameworks, let’s demystify the environment we are working in. If you strip away the advanced machine learning algorithms, R operates on a surprisingly simple three-tier hierarchy.

At its core, R is just a high-powered calculator

When you type an expression into the console, R evaluates it and immediately returns an output. It follows standard mathematical orders of operation.

# R evaluates basic arithmetic operations instantly
2 + 2

[1] 4

# More complex scaling calculation
(300 * 0.045) + (50 * 0.18) - 5

[1] 17.5

R is organized around Functions

To scale beyond basic arithmetic, we use Functions. Think of a function as a reusable machine: you give it an input (arguments), it executes an underlying mathematical formula, and it returns an output.

Instead of manually calculating a mean by adding elements and dividing, we invoke a named function:

# Defining a quick vector of ad spends
sample_spend <- c(100, 150, 200, 300)

# Passing the vector into a built-in statistical function
mean(sample_spend)

[1] 187.5

Functions are organized into Packages

R comes with a “base” suite of tools, but the true power of the language lies in its open-source ecosystem.

Developers wrap collections of specialized functions, datasets, and documentation into cohesive bundles called Packages.

Think of R as the smartphone operating system, and packages as the apps you download to give it new capabilities.

For example, standard R can fit a linear model, but we load the tidyverse package to get advanced data manipulation (mutate(), filter()) and elegant data visualization (ggplot2).

We’ll leverage these packages to build a predictive data pipeline.

The Scenario & Data Architecture

Imagine you are the Lead Data Scientist at a consumer goods firm. The CMO has provided an advertising budget dataset (marketing_data) outlining expenditures across three media channels: YouTube, Facebook, and Print Media (expressed in thousands of dollars), along with the resulting Sales (in thousands of units).

Our goal is to build a model that predicts sales based on a prospective advertising budget allocation.

Generating the Sandbox Data

To ensure total reproducibility, we will simulate a realistic corporate dataset modeling these non-linear corporate dynamics.

# Set seed for reproducibility
set.seed(42)

n <- 200
marketing_data <- tibble(
  youtube   = runif(n, 10, 300),
  facebook  = runif(n, 5, 50),
  newspaper = runif(n, 1, 40)
) %>% 
  mutate(
    # Base sales = 5,000 units. Synergistic/linear impacts added with random Gaussian noise
    sales = 5 + (0.045 * youtube) + (0.18 * facebook) + (0.005 * newspaper) + rnorm(n, 0, 2)
  )

# Preview the architecture
glimpse(marketing_data)

Rows: 200
Columns: 4
$ youtube   <dbl> 275.29375, 281.75187, 92.98047, 250.82981, 196.10620, 160.53783, 22…
$ facebook  <dbl> 44.830296, 28.269997, 43.336894, 24.925832, 12.104605, 24.904609, 4…
$ newspaper <dbl> 1.885300, 21.016342, 25.598320, 17.332093, 35.291372, 5.211496, 39.…
$ sales     <dbl> 25.457857, 24.393000, 17.190735, 22.330796, 15.887119, 16.617315, 2…

Exploratory Data Analysis (EDA)

Before fitting any predictive model, we must visualize our distributions and relationships. We want to identify raw trends and potential collinearity early.

# Pivot longer to generate clean, faceted ggplot visualizations
marketing_data %>%
  pivot_longer(cols = c(youtube, facebook, newspaper), 
               names_to = "channel", 
               values_to = "budget") %>%
  ggplot(aes(x = budget, y = sales)) +
  geom_point(alpha = 0.6, color = "#636363") +
  geom_smooth(method = "lm", se = TRUE, color = "#de2d26", linetype = "dashed") +
  facet_wrap(~ channel, scales = "free_x") +
  theme_minimal(base_size = 14) +
  labs(
    title = "Sales Volume vs. Advertising Spend by Channel",
    subtitle = "Initial linear trend extraction for predictive viability",
    x = "Ad Budget (Thousands of $)",
    y = "Sales (Thousands of Units)"
  )

Say Something Smart: I see strong, clean linear trends for youtube and facebook. newspaper appears highly scattered around the trendline, signaling that it might act as statistical noise in our predictive model.

Fitting the Model

Let’s start with what you know. We will use 100% of our data to build the model, and then use that exact same data to assess how “accurate” we are.

\[ \text{Sales} = \beta_0 + \beta_1(\text{YouTube}) + \beta_2(\text{Facebook}) + \beta_3(\text{Newspaper}) + \epsilon \]

# Model fitting using the entire dataset
full_model <- lm(sales ~ youtube + facebook + newspaper, data = marketing_data)

# Tidy the model for later
model_coefficients <- broom::tidy(full_model)

tidy(full_model)

Whoa that’s ugly, we can clean it up

# Tidy, format, and add custom significance stars
full_model %>%
  tidy() %>%
  mutate(
    # Round statistical estimates for readability
    estimate   = round(estimate, 3),
    std.error  = round(std.error, 3),
    statistic  = round(statistic, 2),
    
    # Create custom significance stars based on p-value thresholds
    stars = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.1   ~ ".",
      TRUE            ~ ""
    ),
    
    # Format p-value column cleanly (handling incredibly small values)
    p.value = if_else(p.value < 0.001, "< 0.001", as.character(round(p.value, 4))),
    
    # Combine p-value and stars into a single clean column
    `p-value` = paste(p.value, stars)
  ) %>%
  select(Term = term, Estimate = estimate, `Std. Error` = std.error, `t-statistic` = statistic, `p-value`)

gt_summary has some nice functions for displaying regression results

tbl_regression(full_model)

Characteristic	Beta	95% CI	p-value
youtube	0.04	0.04, 0.05	<0.001
facebook	0.20	0.18, 0.22	<0.001
newspaper	0.01	-0.01, 0.04	0.2
Abbreviation: CI = Confidence Interval

Deconstructing the Output

• YouTube (\(\beta_1\) = 0.042): Holding all other ad spend constant, every $1,000 increase in the YouTube budget yields an estimated increase of 42 units sold.

• Facebook (\(\beta_2\) = 0.202: Every $1,000 invested drives an estimated 202 units sold, assuming other variables remain static.

Notice that Facebook exhibits a massive marginal return compared to YouTube. In a business context, if the CMO has an unallocated budget, capital should pivot aggressively toward Facebook to maximize rapid acquisition scaling.

Evaluating Training Performance Metrics

Now, let’s calculate our Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) using this full dataset.

# 1. Generate predictions using modelr
full_predictions <- marketing_data %>%
  add_predictions(full_model, var = "pred_sales")

# 2. Create a unified yardstick metric set
marketing_metrics <- yardstick::metric_set(yardstick::rmse, yardstick::mae, yardstick::rsq)

# 3. Calculate all metrics simultaneously into a clean data frame
training_metrics_table <- full_predictions %>% 
  marketing_metrics(truth = sales, estimate = pred_sales)

# 4. Extract the raw numeric values for our callout text below
# (Using pull() safely extracts the number from the .estimate column)
train_rmse_num <- training_metrics_table %>% filter(.metric == "rmse") %>% pull(.estimate)
train_mae_num  <- training_metrics_table %>% filter(.metric == "mae")  %>% pull(.estimate)
train_r2_num   <- training_metrics_table %>% filter(.metric == "rsq")  %>% pull(.estimate)

# 5. Format the table columns nicely for the lecture layout
training_metrics_formatted <- training_metrics_table %>% 
  mutate(
    Metric = c("RMSE (Standard Deviation of Error)", 
               "MAE (Average Absolute Error)", 
               "R² (Coefficient of Determination)"),
    `Training Performance Value` = round(.estimate, 3)
  ) %>% 
  select(Metric, `Training Performance Value`)

# 6. Render the table beautifully using kable
datatable(training_metrics_formatted)

The Illusion of Success

Our model boasts an \(R^2\) of 0.855, meaning it explains 85.5% of the variance in our data. On average, our predictions only miss by 1.45 thousand units.

If you show this to the CMO, they may be thrilled. But mathematically, we have cheated. We gave the model the answers to the test before it took it. We have no objective baseline to know if this model generalizes to next quarter’s unseen marketing campaigns.

In a traditional econometrics class, you use 100% of the data to estimate our model parameters. In predictive modeling, this is a cardinal sin. It leads to overfitting—capturing the random noise of the historical dataset instead of the true underlying signal.

We must evaluate our model’s performance on unseen data. We will partition our data into an 80% Training Set and a 20% Testing Set.

The Remedy: Train vs. Test Partitioning

To protect our firm from operational failure, we must evaluate our model on unseen data. We use the rsample package to partition our original dataset into an 80% Training Set (to build the model) and a 20% Testing Set (held back completely as our proxy for the future).

# Using rsample for rigorous partitioning
set.seed(123) 
data_split <- initial_split(marketing_data, prop = 0.80)

train_data <- training(data_split)
test_data  <- testing(data_split)

# Re-fit the model ONLY using the training slice
split_model <- lm(sales ~ youtube + facebook + newspaper, data = train_data)

The Ultimate Comparison: Train vs. Test Error

Let’s look at how the model performs on the data it was allowed to see (Train) versus how it performs on the data we hid from it (Test).

# 1. Predictions and Metrics on the Training Data (Pulling just the raw number)
train_pred <- train_data %>% add_predictions(split_model, var = "pred_sales")

metrics_train_rmse <- train_pred %>% 
  yardstick::rmse(truth = sales, estimate = pred_sales) %>% 
  pull(.estimate)

metrics_train_mae  <- train_pred %>% 
  yardstick::mae(truth = sales, estimate = pred_sales) %>% 
  pull(.estimate)

metrics_train_rsq  <- train_pred %>% 
  yardstick::rsq(truth = sales, estimate = pred_sales) %>% 
  pull(.estimate)

# 2. Predictions and Metrics on the Test Data (Pulling just the raw number)
test_pred <- test_data %>% add_predictions(split_model, var = "pred_sales")

metrics_test_rmse  <- test_pred %>% 
  yardstick::rmse(truth = sales, estimate = pred_sales) %>% 
  pull(.estimate)

metrics_test_mae   <- test_pred %>% 
  yardstick::mae(truth = sales, estimate = pred_sales) %>% 
  pull(.estimate)

metrics_test_rsq   <- test_pred %>% 
  yardstick::rsq(truth = sales, estimate = pred_sales) %>% 
  pull(.estimate)

# 3. Build and format the comparison table
comparison_table <- tibble(
  Metric = c("RMSE (Standard Deviation of Error)", "MAE (Average Absolute Error)", "Rsq (Coefficient of Determination)"),
  `Training Data (Seen)` = c(metrics_train_rmse, metrics_train_mae, metrics_train_rsq),
  `Testing Data (Unseen)` = c(metrics_test_rmse, metrics_test_mae, metrics_test_rsq)
) %>% 
  mutate(across(where(is.numeric), ~ round(.x, 3)))

# 4. Render the table beautifully for the lecture slides/HTML document
datatable(comparison_table, caption = "Out-of-Sample Performance Comparison")

While this result runs counter to the standard textbook rule that “training error is always lower than test error,” seeing a slightly lower test error is a very common phenomenon in practice. It does not mean our model is broken or that our code is wrong.

Random Sampling Luck (The “Easy Test” Effect)

When we split a dataset into an 80/20 partition, rsample::initial_split() assigns observations completely at random.

Because our testing set only contains 20% of the data (40 observations out of 200), it is highly susceptible to sampling variance. By pure random chance, the 40 observations that landed in our test set happen to cluster slightly closer to the true regression line, containing fewer extreme outliers or anomalous “noise” points than the training set.

Think of it like a professor writing an exam: the training set was a comprehensive, rugged study guide with a few incredibly hard trick questions (outliers). The test set happened to randomly select 40 of the more straightforward questions. The model didn’t suddenly get smarter; it just took a slightly “easier” test.

What This Tells You About our Model

Seeing our test metrics match or beat our training metrics tells you three highly positive things about our predictive architecture:

• Zero Overfitting: If our model had overfit, it would have memorized the training data’s noise, causing the training error to be tiny (e.g., 0.5) and the test error to explode (e.g., 4.5). The fact that they are tightly bounded (1.907 vs 1.661) proves our model captured the true underlying signal without over-indexing on sample quirks.

• High Structural Stability: our model is incredibly stable. The parameters (\(\beta\) weights) estimated on our training data generalize beautifully to independent data.

• The Variance is Small: The difference between an RMSE of 1.91 and 1.66 on a sample size of this scale is statistically negligible. It is well within the expected bounds of standard sampling error.

---

How to Turn This Into a Brilliant Lecture Takeaway

This is the ultimate setup for our take-home exercise or live demonstration. You can show our students how to prove that this is just a quirk of random sampling variance by introducing Cross-Validation or Re-sampling.

Tell our students: “Look at how the metrics shift if we change the random seed that slices the data.”

If you rerun the split with a few different seeds, watch how the relationship bounces around naturally:

# Seed 123 (current split)  : Test Error < Training Error (Lucky Test Split)
# Seed 456                  : Test Error > Training Error (Standard Textbook Split)
# Seed 789                  : Test Error ≈ Training Error (Perfect Equilibrium Split)

The test error isn’t a fixed mathematical law but a random variable with its own distribution, you elevate their understanding from basic formula-following to true predictive data scientists. our model is robust, our code is working perfectly, and here’s a real-world data talking point for today’s class

Analyzing the Results in a Business Context

Usually when lm() optimizes parameters, it minimizes errors specifically for the training sample’s unique quirks and random noise. When confronted with the Testing Data, that specific noise is gone. Only the true underlying signal remains. If our testing error is close to our training error (as it is here), you have built a robust model that can be safely used to project future sales for the corporate marketing strategy. Generally, Testing Error is higher than the Training Error. This is the fundamental reality of predictive modeling.

Diagnostic Audits: Visualizing the Errors (10 Minutes)

Let’s visualize our testing errors using ggplot2 to ensure our linear assumptions hold up out-of-sample.

ggplot(test_pred, aes(x = pred_sales, y = sales)) +
  geom_point(color = "#636363", size = 2, alpha = 0.8) +
  geom_abline(intercept = 0, slope = 1, color = "#de2d26", linetype = "dashed", size = 1) +
  theme_minimal(base_size = 14) +
  labs(
    title = "Actual Sales vs. Out-of-Sample Predictions",
    subtitle = "The dashed red line represents perfect prediction",
    x = "Predicted Sales (Unseen Data)",
    y = "Actual Sales (Unseen Data)"
  )

Summary for the C-Suite

By establishing foundational programming elements, making a structural evaluation mistake, and correcting it through out-of-sample validation, we proved our model’s practical utility. We can confidently tell executive leadership that when deploying this model out into the wild, our average forecasting error will be approximately 1.31 thousand units.

Mini-Exercise

Add a highly complex, unnecessary polynomial term to the model (e.g., I(youtube^3)).

Re-calculate the Full Training Error and the Test Error. Watch how the Training Error drops even lower, while the Test Error spikes wildly. This will showcase extreme overfitting in action!