Linear Models (lm) — Fit and interpret linear models with lm() (intro) - R Programming

SECTION 1

Historical Context & Motivation

The idea of fitting a straight line through data points has roots stretching back over two centuries, long before anyone wrote a line of code. Linear regression is arguably the single most important statistical technique in applied science, and its implementation in R through the lm() function has become a standard tool for data analysis across disciplines. Understanding the historical trajectory helps clarify why the function's output is structured the way it is—each number in a regression summary reflects centuries of theoretical refinement.

1805 • Legendre's Least Squares

Adrien-Marie Legendre published the first formal description of the method of least squares, proposing that the best-fit line minimizes the sum of squared residuals—the mathematical criterion that lm() still uses today.

1809 • Gauss & Normal Distribution

Carl Friedrich Gauss demonstrated that least squares is the optimal estimator when errors follow a normal distribution, providing the statistical backbone for inference in linear models.

1886 • Galton Coins "Regression"

Francis Galton, studying heredity, introduced the term regression toward the mean, giving the entire class of models the name still used today.

1993 • R Language Created

Ross Ihaka and Robert Gentleman began developing R at the University of Auckland, inheriting the formula interface and lm() function from S, making linear modeling accessible to anyone with a terminal.

2000s • R Becomes the Standard

With CRAN's explosive growth and the rise of reproducible research, lm() became the de facto gateway to statistical modeling for millions of researchers, data scientists, and CS professionals.

The core question that linear regression addresses is deceptively simple: given a set of observed data points, what is the best linear relationship between one or more predictor variables and a response variable? R's lm() function wraps the entire ordinary least-squares estimation pipeline—matrix algebra, residual computation, and hypothesis testing—into a single, expressive function call, letting you focus on the modeling question rather than the numerical plumbing.

SECTION 2

Core Principles & Definitions

Before writing any R code, it is essential to internalize the fundamental concepts that underpin every call to lm(). These principles dictate what the function computes, what the output means, and when the results can be trusted. A solid grasp of these ideas will also make debugging model issues far more systematic—an important skill in any CS workflow.

Formula Interface

R uses a formula object (e.g., y ~ x1 + x2) to declaratively specify the response and predictors, separating model specification from data.

Ordinary Least Squares (OLS)

The OLS estimator finds coefficients β that minimize the sum of squared residuals Σ(yᵢ − ŷᵢ)², yielding the Best Linear Unbiased Estimator (BLUE) under Gauss-Markov assumptions.

Coefficients & Intercept

The fitted model returns an intercept (β₀) and one or more slope coefficients (β₁, β₂, …), each measuring the expected change in y per unit change in x, holding other predictors constant.

Residuals

Residuals (eᵢ = yᵢ − ŷᵢ) capture what the model fails to explain. Examining their distribution is the primary diagnostic tool for assessing model validity.

R² and Model Fit

The coefficient of determination (R²) ranges from 0 to 1 and represents the proportion of variance in the response explained by the model. It is the metric most commonly reported alongside coefficient estimates.

✦

KEY TAKEAWAY

Think of lm() like a compiler for statistical models: the formula is your source code that declaratively describes the relationship, the data frame is your input, and the returned model object is your compiled artifact containing everything—coefficients, residuals, diagnostics—ready to be inspected, queried, or piped into downstream analyses. Just as you wouldn't ship software without reading compiler warnings, you shouldn't trust model coefficients without reading the summary diagnostics.

SECTION 3

Visual Explanation — How lm() Fits a Line

The following diagram illustrates the geometric intuition behind ordinary least squares. Each data point sits at some distance from the fitted regression line, and the OLS algorithm positions the line so that the total area of the squared residual segments is minimized. The vertical red segments represent individual residuals—the quantities that lm() stores in the $residuals component of the returned object.

The cyan line represents the fitted regression equation. Violet dots are observed data points. The translucent red squares illustrate the squared residuals whose total area OLS minimizes.

Notice how some residual squares are large and others are tiny. The OLS solution is the unique line that makes the sum of all those red squares as small as possible. This is a convex optimization problem with a closed-form solution, which is why lm() returns instantly even on large datasets—there is no iterative search involved, only a single matrix computation.

SECTION 4

Mathematical Framework

Understanding the algebra behind lm() transforms the function from a black box into a transparent tool. At its core, the function solves the normal equations derived from setting the gradient of the loss function to zero. The following equations present the simple linear regression case and its multivariate generalization.

SIMPLE LINEAR MODEL

yᵢ = β₀ + β₁xᵢ + εᵢ , i = 1, …, n

where yᵢ is the observed response, xᵢ is the predictor, β₀ is the intercept, β₁ is the slope, and εᵢ ~ N(0, σ²) is the error term assumed to be independently and identically distributed.

OLS LOSS FUNCTION

S(β₀, β₁) = Σᵢ₌₁ⁿ (yᵢ − β₀ − β₁xᵢ)²

The OLS criterion minimizes S with respect to β₀ and β₁. Setting ∂S/∂β₀ = 0 and ∂S/∂β₁ = 0 yields the normal equations.

CLOSED-FORM SOLUTION (SIMPLE CASE)

β̂₁ = Σ(xᵢ − x̄)(yᵢ − ȳ) / Σ(xᵢ − x̄)² , β̂₀ = ȳ − β̂₁x̄

Here x̄ and ȳ are the sample means. The slope β̂₁ is the ratio of the sample covariance of x and y to the sample variance of x.

MATRIX FORM (GENERAL CASE)

β̂ = (XᵀX)⁻¹Xᵀy

In the multiple regression setting, X is the n × p design matrix (with a column of ones for the intercept), y is the n × 1 response vector, and β̂ is the p × 1 vector of estimated coefficients. This is the exact computation R performs internally using a QR decomposition for numerical stability.

ℹ

Implementation Note

R does not literally compute (XᵀX)⁻¹ because matrix inversion is numerically unstable. Instead, lm() uses QR decomposition (accessible via qr() on the model object) to solve the normal equations. This is the same numerical strategy taught in a numerical linear algebra course and is O(np²) in time complexity.

SECTION 5

Anatomy of an lm() Call and Its Output

When you call model <- lm(y ~ x, data = df), R returns an object of class "lm"—a named list containing everything needed to inspect, predict from, and diagnose the model. The diagram below maps the most important components of this object and how they relate to one another.

The lm() call returns a list object with components for coefficients, residuals, fitted values, and more. Passing this object to summary() extracts a richer set of diagnostics including standard errors, t-statistics, p-values, R², and the F-statistic.

Accessor	Returns	Use Case
`coef(model)`	Named numeric vector of β̂	Extract slope and intercept for manual prediction
`resid(model)`	Numeric vector of residuals eᵢ	Diagnostic plots, normality checks
`fitted(model)`	Numeric vector of ŷᵢ	Compare predicted vs. actual values
`summary(model)`	Summary object with R², F-stat, coefficient table	Primary tool for model interpretation
`predict(model, newdata)`	Predictions for new observations	Out-of-sample prediction, inference on new data
`confint(model)`	Matrix of 95% confidence intervals for coefficients	Uncertainty quantification for each β̂

Key accessor functions for lm() model objects

SECTION 6

Worked Example — Predicting CPU Performance

Suppose you have a data frame cpu_data with 30 observations measuring processor clock speed (GHz) and benchmark score. We want to fit a simple linear model predicting benchmark score from clock speed, interpret the output, and make a prediction.

Fitting and Interpreting lm() — CPU Benchmark

1
Step 1 — Fit the Model

We call model <- lm(score ~ clock_ghz, data = cpu_data). The formula score ~ clock_ghz tells R that score is the response and clock_ghz is the predictor. R automatically includes an intercept.

model <- lm(score ~ clock_ghz, data = cpu_data)
2
Step 2 — Read summary()

Running summary(model) produces: Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 120.5 18.3 6.58 3.2e-07 *** clock_ghz 245.8 5.7 43.12 < 2e-16 *** R² = 0.9853, Adjusted R² = 0.9848, F-statistic: 1859 on 1 and 28 df The (Intercept) estimate 120.5 is the predicted score when clock speed is 0 GHz (an extrapolation baseline). The clock_ghz coefficient 245.8 means each additional GHz is associated with a 245.8-point increase in benchmark score.

β̂₀ = 120.5, β̂₁ = 245.8, R² = 0.9853
3
Step 3 — Interpret R²

An R² of 0.9853 means that approximately 98.5% of the variance in benchmark scores is explained by clock speed alone. The adjusted R² (0.9848) penalizes for the number of predictors and is nearly identical here because we have only one predictor.

Clock speed explains ~98.5% of score variability
4
Step 4 — Assess Significance

Both p-values are far below 0.05, so both the intercept and slope are statistically significant. The F-statistic of 1859 with a p-value < 2 × 10⁻¹⁶ confirms that the model as a whole is highly significant—the predictor contributes meaningful explanatory power.

p < 2 × 10⁻¹⁶ — model is highly significant
5
Step 5 — Predict a New Value

To predict the score for a 3.5 GHz processor: predict(model, newdata = data.frame(clock_ghz = 3.5)). Manually: ŷ = 120.5 + 245.8 × 3.5 = 120.5 + 860.3 = 980.8. The predicted benchmark score is approximately 981 points.

ŷ = 980.8

SECTION 7

Strengths, Limitations & Common Pitfalls

No model is universally appropriate. Understanding when lm() is the right tool—and when it will mislead you—is as important as knowing how to call it. The table below contrasts its strengths and limitations to help you make informed modeling decisions.

Strengths	Limitations
Closed-form solution — no tuning, no convergence issues, O(np²) complexity	Assumes a linear relationship; cannot capture non-linear patterns without feature engineering
Highly interpretable — each coefficient has a direct unit-change meaning	Sensitive to outliers — a single extreme point can heavily influence β̂
Well-understood theory — inference (p-values, CIs) is exact under assumptions	Requires assumptions (linearity, independence, homoscedasticity, normality of errors)
Rich ecosystem — works seamlessly with predict(), anova(), broom::tidy(), ggplot2	Multicollinearity inflates standard errors and makes individual coefficient estimates unstable
Serves as baseline — essential benchmark before trying complex models	No built-in regularization — prone to overfitting with many predictors relative to n

Strengths vs. limitations of lm() for linear regression

✦

KEY TAKEAWAY

In software engineering terms, lm() is like a well-documented, stable API with clearly stated preconditions (the Gauss-Markov assumptions). If you call the API with inputs that violate those preconditions, it will still return an answer—but the guarantees on correctness are void. Always run diagnostics (plot(model)) just as you would run unit tests before deploying code.

SECTION 8

Connection to Advanced Modeling in R

The lm() function is the entry point in a large family of modeling tools in R. Once you are comfortable with its interface and interpretation, extending to more powerful frameworks requires surprisingly little additional syntax—the formula interface is shared across nearly all modeling functions. The table below contrasts lm() with its natural successors.

Feature	lm()	glm()	glmnet / lasso
Response type	Continuous (numeric)	Continuous, binary, count, etc.	Same as glm (with regularization)
Link function	Identity (y = Xβ)	Logit, log, inverse, etc.	Same as glm
Regularization	None	None	L1 (Lasso), L2 (Ridge), Elastic Net
Feature selection	Manual (stepAIC, p-values)	Manual	Automatic (Lasso shrinks coefficients to zero)
Estimation	Closed-form (QR decomposition)	Iteratively reweighted least squares	Coordinate descent
When to use	Baseline, continuous outcome, interpretable model	Non-normal response (logistic, Poisson)	High-dimensional data, p ≈ n or p > n

lm() in context: progression toward more flexible models

As you advance, you will also encounter nlme::lme() and lme4::lmer() for mixed-effects models, mgcv::gam() for generalized additive models, and machine learning frameworks like tidymodels that wrap lm() in a consistent pipeline. Every one of these tools shares the same formula syntax you learned here, which means mastering lm() gives you transferable fluency across the entire R modeling ecosystem.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Explain what the formula y ~ x1 + x2 passed to lm() specifies about the statistical model, and describe how R interprets the tilde (~) and the plus sign (+) in this context. Why does R automatically include an intercept term even though it is not written explicitly?

Show answer

The tilde (~) separates the response variable (left side) from the predictors (right side), expressing the relationship yᵢ = β₀ + β₁x₁ᵢ + β₂x₂ᵢ + εᵢ. The plus sign (+) indicates additive (not multiplicative) inclusion of predictors in the linear predictor. R includes an intercept by default because omitting it would force the regression surface through the origin, which is rarely justified on substantive grounds. If you want to suppress it, you must write y ~ x1 + x2 - 1 or equivalently y ~ 0 + x1 + x2. The formula interface is a domain-specific language embedded in R that maps directly to the design matrix X—each term on the right side generates one or more columns in X.

Question 2: BASIC CALCULATION

You fit model <- lm(time ~ size, data = jobs) and coef(model) returns (Intercept) = 2.3, size = 0.045. Predict the execution time for a job of size 500. What does the slope of 0.045 mean in plain language?

Show answer

Using the fitted equation ŷ = 2.3 + 0.045 × size, for size = 500: ŷ = 2.3 + 0.045 × 500 = 2.3 + 22.5 = 24.8 time units. The slope of 0.045 means that for each one-unit increase in job size, the model predicts an additional 0.045 time units of execution time, holding everything else constant. In R, you would compute this with predict(model, newdata = data.frame(size = 500)). Note that this prediction is only reliable within the range of observed size values (interpolation), not for extrapolation far beyond the training data.

Question 3: INTERMEDIATE

After fitting model <- lm(latency ~ threads + memory_gb, data = perf), the summary shows R² = 0.72, Adjusted R² = 0.68, and the p-value for memory_gb is 0.34. What do you conclude about the memory_gb predictor? What single command would you use to compare the full model against a reduced model without memory_gb?

Show answer

A p-value of 0.34 for memory_gb means we fail to reject the null hypothesis H₀: β_memory = 0 at any conventional significance level (0.05 or 0.10). This suggests that, conditional on threads being in the model, memory_gb does not contribute statistically significant explanatory power to latency. The gap between R² (0.72) and Adjusted R² (0.68) further hints at a predictor that adds complexity without proportional improvement. To formally compare models, fit the reduced model reduced <- lm(latency ~ threads, data = perf) and then run anova(reduced, model) which performs a partial F-test for the nested models.

Question 4: APPLIED

You are building a model to predict HTTP response time from request payload size and number of concurrent connections. Write the complete R code to: (a) fit the model, (b) print the coefficient table, (c) generate diagnostic plots, and (d) predict response time for a 2048-byte payload with 100 concurrent connections. Assume the data frame is called web_logs with columns resp_ms, payload_bytes, and connections.

Show answer

(a) model <- lm(resp_ms ~ payload_bytes + connections, data = web_logs). (b) summary(model) displays the coefficient table with estimates, standard errors, t-values, and p-values. (c) par(mfrow = c(2, 2)); plot(model) produces four diagnostic plots: Residuals vs. Fitted, Q-Q, Scale-Location, and Residuals vs. Leverage. (d) predict(model, newdata = data.frame(payload_bytes = 2048, connections = 100)). For prediction intervals, add interval = "prediction" to the predict call. Always check the diagnostic plots before trusting any prediction—look for non-random residual patterns, heavy tails on the Q-Q plot, or high-leverage points.

Question 5: CRITICAL THINKING

A colleague fits lm(bugs ~ lines_of_code + num_devs + commits, data = projects) and obtains R² = 0.95 on the training data. They claim this proves that more lines of code cause more bugs. Identify at least three flaws in this reasoning, referring to specific statistical concepts. What additional analyses would you recommend?

Show answer

Flaw 1: Correlation does not imply causation. A high R² indicates association, not a causal mechanism. Lines of code and bugs may both be driven by project complexity (a confounding variable). Flaw 2: A training R² of 0.95 with three predictors may reflect overfitting, especially if n is small. The colleague should report adjusted R², perform cross-validation, or use a holdout test set to assess generalization. Flaw 3: Multicollinearity is likely—lines_of_code, num_devs, and commits are plausible correlated. This inflates standard errors and makes individual coefficient estimates unreliable. Check VIF values using car::vif(model). Recommended additional analyses: run diagnostic plots to check residual assumptions, compute VIF for collinearity, use cross-validation for honest predictive assessment, and consider a causal inference framework (e.g., DAGs or instrumental variables) if the goal is genuinely to establish causality rather than prediction.

SUMMARY

Summary

R's lm() function fits ordinary least-squares linear models using R's declarative formula interface (e.g., y ~ x1 + x2). Internally, it solves β̂ = (XᵀX)⁻¹Xᵀy via QR decomposition, returning a rich model object containing coefficients, residuals, fitted values, and metadata. The summary() function extracts the coefficient table with standard errors, t-statistics, and p-values, along with R² and Adjusted R² for goodness-of-fit assessment and the F-statistic for overall model significance.

Always validate models with diagnostic plots (plot(model)) to check the Gauss-Markov assumptions: linearity, independence, homoscedasticity, and normality of residuals. Use predict() for out-of-sample predictions and confint() for confidence intervals. Mastering lm() provides the foundation for R's entire modeling ecosystem, from glm() for generalized linear models to regularized regression and beyond.