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.
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.
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
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
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
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
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() and residuals().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.Visual Explanation — The Model Object Anatomy
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.
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.
coef(model). It minimizes the sum of squared residuals ∑eᵢ². In simple linear regression (one predictor), β̂₁ = slope and β̂₀ = intercept.fitted(model). The matrix H = X(XᵀX)⁻¹Xᵀ is the 'hat matrix' that projects y onto the column space of X.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.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.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.
| 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 |
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.
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.modelcoef(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.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.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.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.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.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) |
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).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 |
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.
Practice Problems
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?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.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.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).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?Summary — Extracting Model Outputs in R
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.