Quick summary
Summarize this blog with AI
Introduction
R’s lm() fits a model linear in its coefficients. The workflow below uses supplied fictional marketing data to fit a model, interpret coefficients, inspect residuals and compare held-out error with a baseline. A high training R-squared alone does not establish predictive usefulness or a causal advertising effect.
Tested with R 4.5.1. Run the code blocks in order. Install the named packages before starting; inputs are supplied below.
Create the fixture and reserve evaluation rows
Each row is a fictional independent market observation. sales and ad_spend use consistent arbitrary units; region is categorical. The fixed split makes the example reproducible. A real split should respect time, customers or repeated units, and all preprocessing must be learned within the training partition.
dat <- data.frame(id = 1:12, ad_spend = 1:12,
region = factor(rep(c("A", "B"), 6)),
sales = c(11, 19, 17, 23, 19, 29, 23, 33, 27, 37, 31, 41))
test_ids <- c(3L, 6L, 9L, 12L)
training <- dat[!dat$id %in% test_ids, ]
testing <- dat[dat$id %in% test_ids, ]
stopifnot(!anyNA(dat), all(is.finite(dat$ad_spend)), all(is.finite(dat$sales)),
length(intersect(training$id, testing$id)) == 0,
nrow(training) == 8, nrow(testing) == 4)
print(training)
Fit numeric and categorical predictors
In lm(sales ~ ad_spend + region), + adds main effects and the intercept is included by default. With A as the reference, regionB is the fitted B-versus-A difference holding ad_spend fixed. ad_spend * region would also add an interaction; I(ad_spend^2) adds a raw squared term while remaining linear in coefficients.
fit <- lm(sales ~ ad_spend + region, data = training, na.action = na.fail)
design <- model.matrix(fit)
stopifnot(fit$rank == ncol(design), !anyNA(coef(fit)))
print(round(coef(fit), 4))
print(summary(fit)$coefficients)
stopifnot(max(abs(as.numeric(design %*% coef(fit)) - fitted(fit))) < 1e-10)
The fitted intercept is 8.2174, ad_spend slope is 2.1304 and regionB coefficient is 7.0000. The intercept is the fitted regionA mean at zero ad spend, which lies outside this training range; it need not have a useful standalone interpretation. The slope is an association conditional on region. Causal interpretation requires additional design and identification assumptions.
Separate training fit from statistical inference
Training R-squared compares residual squared error with variation around the training mean for this intercept model. Adding predictors cannot reduce ordinary training R-squared on the same rows, so it is not a model-selection guarantee. Standard coefficient tests assume an appropriate conditional-mean specification and independent, constant-variance errors; normal errors provide exact small-sample t inference.
r_squared <- summary(fit)$r.squared
manual_r2 <- 1 - sum(residuals(fit)^2) / sum((training$sales - mean(training$sales))^2)
print(round(c(training_R2 = r_squared), 4))
stopifnot(abs(r_squared - manual_r2) < 1e-12)
A small coefficient p-value does not prove reliability, causality or out-of-sample accuracy. Strongly correlated predictors can make individual coefficient estimates unstable even when fitted values are stable; exact collinearity makes some coefficients unidentifiable.
Inspect residual pattern, spread and tails
Residuals versus fitted values can reveal curvature or changing spread. A Q-Q plot checks error-shape plausibility; it does not assess independence. The sample is intentionally tiny, so these are teaching diagnostics, not proof that assumptions hold. Investigate influential observations and collection errors without automatically deleting them.
draw_lm_checks <- function() {
old <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
on.exit(par(old))
plot(fitted(fit), residuals(fit), pch = 19, xlab = "Fitted sales", ylab = "Residual",
main = "Residuals versus fitted")
abline(h = 0, col = "#a63d28")
qqnorm(residuals(fit), pch = 19, main = "Residual Q-Q plot")
qqline(residuals(fit), col = "#a63d28")
}
draw_lm_checks()
Compare held-out error with a training-mean baseline
predict.lm() uses the fitted coefficients on new rows. Calculate errors on the same untouched test observations for both model and baseline. The baseline must use the training mean, not the test mean.
predicted <- predict(fit, newdata = testing)
baseline <- rep(mean(training$sales), nrow(testing))
rmse <- function(actual, predicted) sqrt(mean((actual - predicted)^2))
scores <- c(model_RMSE = rmse(testing$sales, predicted),
baseline_RMSE = rmse(testing$sales, baseline),
model_MAE = mean(abs(testing$sales - predicted)))
print(round(data.frame(id = testing$id, actual = testing$sales, predicted = predicted), 3))
print(round(scores, 3))
stopifnot(all(is.finite(predicted)), scores["model_RMSE"] < scores["baseline_RMSE"])
ci <- predict(fit, newdata = testing, interval = "confidence")
pi <- predict(fit, newdata = testing, interval = "prediction")
stopifnot(all(pi[, "upr"] - pi[, "lwr"] > ci[, "upr"] - ci[, "lwr"]))
Test RMSE is 1.315 for the model versus 9.421 for the baseline, and model MAE is 1.000. The model beats the baseline on this four-row fixture, which is insufficient evidence of dependable real-world performance. A confidence interval describes the conditional mean; a prediction interval includes individual-observation variation and is wider under the model. Neither guarantees coverage if the assumptions or deployment population change.
Use training-only folds for development
The following four-fold example leaves the final test rows untouched. It refits the entire lm model in each training fold and records one out-of-fold prediction per training row. If selecting features, tuning transformations or imputing missing values, perform those steps inside the folds as well. Reusing the same folds for many selection decisions can still overfit the development process.
fold <- c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L)
out_of_fold <- rep(NA_real_, nrow(training))
for (k in sort(unique(fold))) {
inner <- lm(sales ~ ad_spend + region, data = training[fold != k, ], na.action = na.fail)
out_of_fold[fold == k] <- predict(inner, training[fold == k, ])
}
cat("Training-only cross-validation RMSE:", round(rmse(training$sales, out_of_fold), 3), "
")
stopifnot(!anyNA(out_of_fold))
fails <- function(expr) inherits(tryCatch(expr, error = identity), "error")
missing_training <- training; missing_training$sales[1] <- NA
stopifnot(fails(lm(sales ~ ad_spend + region, missing_training, na.action = na.fail)))
unknown <- data.frame(ad_spend = 5, region = "C")
stopifnot(fails(predict(fit, unknown)),
length(predict(fit, testing[FALSE, ])) == 0)
Training-only cross-validation RMSE is 0.817 for this fixed fold assignment. New categorical levels require an explicit modeling policy. Predictions far beyond observed ad-spend values are extrapolations. For time-series data use time-respecting validation and model temporal dependence; regressing stock price on date alone does not establish a useful forecasting strategy.
