Extracting Model Outputs — Extract basic model outputs (summary, coefficients, fitted values) (intro) - R Programming

SECTION 1

Historical Context & Motivation

Statistical computing has long demanded that analysts move beyond mere visual inspection of results and instead gain programmatic access to the numerical components of fitted models. In early statistical software, extracting a single regression coefficient required parsing printed output or navigating opaque data structures, a workflow that was fragile and error-prone. The development of S and its successor R introduced a fundamentally different philosophy: every model is a rich, structured object whose internals can be queried with well-defined accessor functions. This design philosophy — treating statistical results as first-class data structures — transformed the way quantitative research is conducted and automated.

1976 • S Language at Bell Labs

John Chambers and colleagues create the S language, introducing the concept of model objects as structured lists that store coefficients, residuals, and metadata in a single return value.

1993 • R's Genesis

Ross Ihaka and Robert Gentleman begin developing R at the University of Auckland, inheriting S's object-oriented model framework and expanding the generic function system for model extraction.

2000 • R 1.0.0 Released

The first stable release ships with a mature set of extractor functions — summary(), coef(), fitted(), residuals() — codifying a standard interface for model interrogation.

2014 • The broom Package

David Robinson releases broom, introducing tidy(), glance(), and augment() to convert model outputs into tidy data frames, making extraction even more systematic for modern data-science pipelines.

2020s • Tidymodels Ecosystem

The tidymodels framework standardizes model extraction across hundreds of model types, embedding extractor functions into a unified, pipeline-friendly API used in production and research.

The central question this lesson addresses is straightforward yet essential: once you have fit a statistical model in R — for instance, with lm() — how do you programmatically extract the summary table, individual coefficients, fitted values, and other diagnostics? Mastering this skill is the bridge between running a model and actually using its results in downstream code, reports, and automated decision systems.

SECTION 2

Core Principles & Definitions

Understanding how R stores and exposes model outputs requires grasping a few foundational concepts. In R, a call to a modeling function like lm() returns an S3 object — essentially a named list with a class attribute. The class attribute tells R which version of a generic function to dispatch. When you call summary(model), R dispatches to summary.lm() because the object's class is "lm". This polymorphic design means the same extraction interface works across dozens of model types.

Model Object as a List

The return value of lm() is a named list containing elements like $coefficients, $residuals, $fitted.values, and more. You can inspect all names with names(model).

Generic Extractor Functions

Functions like coef(), fitted(), and residuals() are generics that dispatch to methods based on the model's class — a form of polymorphism central to R's S3 system.

summary() Returns an Object

Calling summary(model) does not merely print text. It returns an object of class "summary.lm" with additional computed quantities like R², adjusted R², and the coefficient significance table.

The $ Operator vs. Accessor Functions

While model$coefficients accesses the raw list element, using coef(model) is preferred because the generic function may apply transformations or provide a more stable API across model types.

Fitted Values & Residuals

Fitted values (ŷ) are the model's predictions on the training data. Residuals (e = y − ŷ) measure the discrepancy. Both are stored in the model object and retrieved via fitted() and residuals().

✦

KEY TAKEAWAY

Think of a fitted model in R like a complex API response in web development. When you hit an endpoint, you get back a structured JSON object with many fields. You wouldn't manually parse the raw HTTP body — you'd use accessor methods to pull out response.status or response.data.items. Similarly, R's extractor functions like coef() and summary() are the clean, reliable accessors for the "response" that lm() returns.

SECTION 3

Visual Explanation — The Model Object Anatomy

The diagram shows the internal structure of an lm object (top) and the extractor functions (middle) that provide a stable API. Note that summary() creates a second-level object (bottom) with additional computed statistics like R² and the full coefficient significance table.

As the diagram illustrates, the model object and the summary object are two distinct layers of information. A common source of confusion for beginners is conflating coef(model) — which returns a simple named numeric vector of estimated coefficients — with summary(model)$coefficients — which returns a matrix whose columns include the estimate, standard error, t-statistic, and p-value. Recognizing this structural difference is essential when you need to extract specific quantities for downstream computation, such as pulling a p-value for an automated significance test in a pipeline.

SECTION 4

Mathematical Framework — What the Outputs Represent

To appreciate what R is storing and extracting, it is helpful to review the linear model that lm() fits and the mathematical meaning of each output component. The Ordinary Least Squares (OLS) framework underlies the model object's structure: the coefficients minimize a specific loss, and the residuals and fitted values are algebraic consequences of that minimization.

LINEAR MODEL

y = Xβ + ε

Where y is the n × 1 response vector, X is the n × p design matrix (including the intercept column), β is the p × 1 coefficient vector, and ε is the n × 1 error vector.

OLS COEFFICIENT ESTIMATES

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

This is the vector returned by coef(model). It minimizes the sum of squared residuals ∑eᵢ². In simple linear regression (one predictor), β̂₁ = slope and β̂₀ = intercept.

FITTED VALUES

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

The vector ŷ is returned by fitted(model). The matrix H = X(XᵀX)⁻¹Xᵀ is the 'hat matrix' that projects y onto the column space of X.

RESIDUALS

e = y − ŷ = (I − H)y

Returned by residuals(model). The residuals measure the gap between observed and fitted values. Their sum is zero when an intercept is included, and they are orthogonal to the fitted values: ŷᵀe = 0.

ℹ

CS Perspective: Computational Cost

R does not literally compute (XᵀX)⁻¹. Instead, lm() uses a QR decomposition of X, which is numerically more stable and runs in O(np²) time. The coefficients, fitted values, and residuals are all derived from this single factorization, so extraction via coef(), fitted(), and residuals() are O(1) lookups — the results are precomputed and cached in the model object.

SECTION 5

Detailed Extraction Map — Functions and Their Returns

This section provides a comprehensive reference mapping each extractor function to the data it returns, the R type of the return value, and a typical use case. Understanding these mappings is the key to writing clean, maintainable statistical code. Rather than memorizing internal list names, you should rely on the generic accessor functions because they abstract over implementation details that may differ between model classes.

This workflow diagram traces the path from fitting a model, through extraction, to downstream use cases. Each extractor function is color-coded and connected to a practical application. The dashed boxes at the bottom show typical next steps after extraction.

Function	Return Type	What It Contains	Typical Use
`coef(model)`	Named numeric vector	Estimated β̂ values (intercept, slopes)	Build prediction equations, compare models
`fitted(model)`	Named numeric vector	ŷᵢ for each training observation	Actual-vs-predicted plots, R² computation
`residuals(model)`	Named numeric vector	eᵢ = yᵢ − ŷᵢ for each observation	Residual diagnostics, normality checks
`summary(model)`	summary.lm object (list)	R², adj R², F-stat, σ̂, coef matrix with SE/t/p	Full model assessment, significance testing
`confint(model)`	Matrix (p × 2)	95% confidence intervals for each β̂	Uncertainty quantification
`vcov(model)`	Matrix (p × p)	Variance-covariance matrix of β̂	Standard errors, hypothesis tests

Common extractor functions for lm objects

SECTION 6

Worked Example — Extracting Outputs from a Simple Linear Model

Let's walk through a complete example using R's built-in mtcars dataset. We will fit a simple linear regression predicting miles per gallon (mpg) from car weight (wt), and then extract every key output.

Extracting Outputs from lm(mpg ~ wt, data = mtcars)

1
Step 1 — Fit the Model

We begin by calling lm() and storing the result. The code is: model <- lm(mpg ~ wt, data = mtcars). This creates an S3 object of class "lm". We can verify by running class(model), which returns "lm". Calling names(model) reveals the 12 internal elements: coefficients, residuals, effects, rank, fitted.values, assign, qr, df.residual, xlevels, call, terms, model.

Model object created and stored in model
2
Step 2 — Extract Coefficients

Running coef(model) returns a named numeric vector: (Intercept) = 37.2851, wt = -5.3445. This tells us that the estimated regression equation is ŷ = 37.29 − 5.34 × wt. Each additional 1000 lbs of weight is associated with a decrease of about 5.34 mpg. You can access individual coefficients by name: coef(model)["wt"] returns -5.344472.

β̂₀ = 37.29 (intercept), β̂₁ = −5.34 (slope for wt)
3
Step 3 — Extract Fitted Values

Calling fitted(model) returns a named numeric vector of length 32 (one per car). For instance, the Mazda RX4 (wt = 2.620) has a fitted value of 37.29 − 5.34 × 2.620 ≈ 23.28 mpg. You can verify: fitted(model)["Mazda RX4"] returns approximately 23.28. These are the model's in-sample predictions and form the regression line when plotted against the predictor.

32 fitted values; Mazda RX4 ŷ ≈ 23.28 mpg
4
Step 4 — Extract Residuals

The command residuals(model) (or its alias resid(model)) returns eᵢ = yᵢ − ŷᵢ for each observation. For the Mazda RX4, the actual mpg is 21.0 and the fitted value is 23.28, so the residual is 21.0 − 23.28 = −2.28. We can confirm: sum(residuals(model)) returns a value essentially zero (within machine precision), as expected when an intercept is present.

Mazda RX4 residual ≈ −2.28; ∑eᵢ ≈ 0
5
Step 5 — Extract Summary Statistics

We store the summary: s <- summary(model). Now s$r.squared returns 0.7528 (about 75.3% of variance in mpg is explained by weight). The coefficient table s$coefficients is a 2×4 matrix with columns Estimate, Std. Error, t value, and Pr(>|t|). To extract the p-value for the slope: s$coefficients["wt", "Pr(>|t|)"] returns approximately 1.29 × 10⁻¹⁰, indicating the relationship is highly statistically significant.

R² = 0.7528, p-value for wt ≈ 1.29 × 10⁻¹⁰

💡

Pro Tip: Combining Extractors

You can chain extractions into a single pipeline. For example, to create a data frame with actual, fitted, and residual values: data.frame(actual = mtcars$mpg, fitted = fitted(model), residual = residuals(model)). This is exactly the kind of tidy data structure that feeds into ggplot2 or downstream analysis.

SECTION 7

Strengths & Limitations — $ vs. Extractor Functions

A natural question arises: why use coef(model) when model$coefficients seems to do the same thing? The distinction is subtle but architecturally significant, especially for CS students accustomed to thinking about interface contracts and encapsulation. The extractor functions provide an abstraction layer that shields your code from implementation changes, much like using getter methods in object-oriented programming rather than directly accessing fields.

Criterion	$ Direct Access	Generic Extractor Functions
Polymorphism	Tied to one model class's internal naming convention	Works across lm, glm, nls, lme4, and hundreds of model types
Stability	Internal names may change between R or package versions	The generic interface is part of R's public API and rarely changes
Readability	`model$fitted.values` is verbose	`fitted(model)` is concise and expressive
Potential Pitfalls	Partial matching: `model$res` may return residuals or something else	No partial matching risk; function name is explicit
Customization	Returns raw internal data only	Methods can apply transformations (e.g., deviance residuals for glm)

Comparing direct list access with generic extractor functions

🏗️

DESIGN PRINCIPLE

Using generic extractors like coef() and fitted() is analogous to programming to an interface rather than an implementation in Java or TypeScript. If you later switch your model from lm() to glm() or a random forest, your extraction code may still work without modification because the generics dispatch to the appropriate method. This is the open/closed principle in action — the system is open for extension (new model types) but closed for modification (the extraction interface stays the same).

SECTION 8

Connection to Advanced Extraction — broom and Tidymodels

The base R extraction functions we have covered are powerful, but modern data science workflows often demand model outputs in tidy data frames rather than named vectors or matrices. The broom package bridges this gap by providing three functions that convert model outputs into tibbles: tidy() for coefficient-level statistics, glance() for model-level summaries, and augment() for observation-level data (fitted values, residuals, influence measures). These tidy outputs integrate seamlessly with dplyr pipelines and ggplot2 visualization.

Aspect	Base R Extraction	broom / Tidymodels
Coefficient info	`summary(m)$coefficients` → matrix	`tidy(m)` → tibble with term, estimate, std.error, statistic, p.value columns
Model-level stats	Extract individually: `s$r.squared`, `s$sigma`	`glance(m)` → single-row tibble with R², adj R², σ, F-stat, p, df, AIC, BIC
Observation-level	Must manually combine `fitted()`, `residuals()`, etc.	`augment(m)` → tibble with .fitted, .resid, .hat, .cooksd, .std.resid
Pipeline friendly	Requires manual wrangling into data frames	Outputs are immediately pipe-ready for dplyr and ggplot2
Model agnostic	Each model type may store results differently	Consistent column names across 100+ model types

Base R vs. broom extraction paradigms

Understanding base R extraction first is essential because broom internally calls many of the same functions and because you will frequently encounter legacy code and packages that rely exclusively on the base approach. Once you are comfortable with coef(), fitted(), residuals(), and summary(), moving to the broom ecosystem will feel like a natural extension — the same information, but delivered in a format optimized for modern, reproducible analysis pipelines.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Explain why coef(model) and summary(model)$coefficients return different data structures, even though both relate to the model's coefficients. What does each one contain, and when would you prefer one over the other?

Show answer

The function coef(model) returns a named numeric vector containing only the point estimates (β̂₀, β̂₁, …). In contrast, summary(model)$coefficients returns a matrix with four columns: Estimate, Std. Error, t value, and Pr(>|t|). You would use coef() when you need only the numeric estimates — for instance, to construct a prediction equation or to pass coefficients to another function. You would use the summary coefficient matrix when you need to assess statistical significance, extract p-values, or compute confidence intervals manually. The key insight is that summary() computes additional quantities (standard errors, test statistics, p-values) that are not stored in the raw model object.

Question 2: BASIC CALCULATION

Given the output coef(model) returns (Intercept) = 12.5, x1 = 3.2, x2 = -1.8, write the R code to extract only the coefficient for x1, and manually compute the fitted value for an observation where x1 = 4 and x2 = 2.

Show answer

To extract the coefficient for x1: coef(model)["x1"] returns 3.2. To compute the fitted value: ŷ = 12.5 + 3.2 × 4 + (−1.8) × 2 = 12.5 + 12.8 − 3.6 = 21.7. In R code, you could compute this as sum(coef(model) * c(1, 4, 2)) which multiplies the intercept by 1, x1's coefficient by 4, and x2's coefficient by 2, then sums the products. Alternatively, predict(model, newdata = data.frame(x1 = 4, x2 = 2)) achieves the same result more robustly.

Question 3: INTERMEDIATE

You fit model <- lm(y ~ x, data = df) and store s <- summary(model). Write R code that: (a) extracts the R² value, (b) extracts the p-value for the slope coefficient, and (c) creates a logical variable indicating whether the slope is statistically significant at α = 0.01.

Show answer

(a) r_sq <- s$r.squared extracts the coefficient of determination. (b) p_val <- s$coefficients["x", "Pr(>|t|)"] extracts the p-value for the slope by indexing the coefficient matrix by row name "x" and column name "Pr(>|t|)". Alternatively, s$coefficients[2, 4] uses numeric indices (row 2 = slope, column 4 = p-value), though named indexing is more robust. (c) is_significant <- p_val < 0.01 creates a logical. This entire sequence demonstrates extracting nested elements from the summary object and using them programmatically for automated decision-making in a pipeline.

Question 4: APPLIED

You are building an automated model-validation report. Write an R function called model_report(model) that takes an lm object and returns a named list with elements: r_squared, rmse (root mean squared error of residuals), max_abs_residual, and significant_predictors (a character vector of predictor names with p < 0.05).

Show answer

Here is the function:

model_report <- function(model) { s <- summary(model); resids <- residuals(model); coef_mat <- s$coefficients; sig_rows <- coef_mat[-1, "Pr(>|t|)"] < 0.05; list(r_squared = s$r.squared, rmse = sqrt(mean(resids^2)), max_abs_residual = max(abs(resids)), significant_predictors = names(which(sig_rows))) }

. Key details: we use coef_mat[-1, ] to exclude the intercept row (since the intercept is typically not a 'predictor' in the practical sense). The RMSE is computed from the raw residuals as √(mean(e²)). The function names(which(sig_rows)) converts the logical mask into a character vector of significant predictor names. This demonstrates combining multiple extractor calls into a reusable, production-ready function.

Question 5: CRITICAL THINKING

Consider the following scenario: you switch from lm() to glm(family = poisson). Discuss which of the following extraction calls will still work unchanged, which will return different quantities, and which may break: (a) coef(model), (b) fitted(model), (c) model$fitted.values, (d) summary(model)$r.squared. What does this reveal about the value of generic extractor functions versus direct list access?

Show answer

(a) coef(model) works for both lm and glm objects — the generic dispatches to coef.glm() — but for Poisson GLMs the coefficients are on the log scale, so their interpretation changes. (b) fitted(model) also works: for Poisson GLMs it returns predicted counts (on the response scale), not log-counts, because the method applies the inverse link function automatically. (c) model$fitted.values still works here because glm objects happen to store this element with the same name, but this is an implementation detail, not a guarantee. (d) summary(model)$r.squared will return NULL for a GLM because R² is not defined the same way for generalized linear models; the summary.glm object does not contain this element. This demonstrates the fragility of direct list access and the subtlety of model-type-specific outputs. Generic extractors provide a more portable and semantically correct interface, though the analyst must still understand what the extracted quantities mean in the new model context.

SUMMARY

Summary — Extracting Model Outputs in R

Summary

R's modeling functions return S3 objects — structured named lists with a class attribute that enables polymorphic dispatch. The four fundamental extraction tools are coef() for coefficient estimates (β̂), fitted() for in-sample predictions (ŷ), residuals() for the differences between observed and fitted values (e = y − ŷ), and summary() which returns a richer object containing R², the full coefficient significance matrix, and overall model diagnostics. Always prefer these generic extractor functions over direct $ access to ensure portability across model types.

A critical distinction is that coef(model) returns a simple vector while summary(model)$coefficients returns a matrix with standard errors, t-statistics, and p-values. For modern pipelines, the broom package extends this paradigm by converting all model outputs into tidy data frames via tidy(), glance(), and augment(). Mastering base extraction first provides the foundation for understanding what these higher-level tools compute and return.