Programming for Applications

Chapter 20: Regression Models

Yu-You Liou (NTU)

Shih Chien University

2026-08-02

Predicting Continuous Values

What Regression Is For

A regression model relates a continuous response (dependent variable) to a set of predictors (independent/stimulus variables). Three uses:

  • prediction — a warfarin dose from a patient’s weight; a credit customer’s balance;
  • explanation — how weight relates to food intake, adjusting for age, exercise, genetics;
  • inference — visualizing trends, ANOVA, testing variable significance.

This chapter starts with the simplest model — linear regression by ordinary least squares (OLS) — then builds, evaluates, refines it, and surveys richer model types. Classification waits for Chapter 21.

A Simple Linear Model

The Data and the Question

Which batting metrics predict a team’s runs? Team batting 2000–2008 (extracted by SQL in Chapter 16, packaged as team.batting.00to08):

temp <- tempfile(fileext = ".rda")
download.file("https://raw.githubusercontent.com/cran/nutshell/master/data/team.batting.00to08.rda",
              temp, mode = "wb")
load(temp)

A lattice scatter of runs against each variable hints at importance — home runs and walks both track runs strongly:

library(lattice)
attach(team.batting.00to08)
forplot <- make.groups(
  singles=data.frame(value=singles, runs), doubles=data.frame(value=doubles, runs),
  triples=data.frame(value=triples, runs), homeruns=data.frame(value=homeruns, runs),
  walks=data.frame(value=walks, runs))
detach(team.batting.00to08)
xyplot(runs~value|which, data=forplot, scales=list(relation="free"),
       pch=19, cex=.2, layout=c(5,1))

lm: Fitting

lm fits a linear model by OLS. Note: in formulas + * - ^ are special; wrap literal arithmetic in I(), and generate polynomial terms with poly(x, degree=).

runs.mdl <- lm(
  formula=runs~singles+doubles+triples+homeruns+
            walks+hitbypitch+sacrificeflies+
            stolenbases+caughtstealing,
  data=team.batting.00to08)
runs.mdl

Call:
lm(formula = runs ~ singles + doubles + triples + homeruns + 
    walks + hitbypitch + sacrificeflies + stolenbases + caughtstealing, 
    data = team.batting.00to08)

Coefficients:
   (Intercept)         singles         doubles         triples        homeruns  
    -507.16020         0.56705         0.69110         1.15836         1.47439  
         walks      hitbypitch  sacrificeflies     stolenbases  caughtstealing  
       0.30118         0.37750         0.87218         0.04369        -0.01533

R prints almost nothing on fitting — information comes from helper functions (a deliberate contrast with SAS/SPSS).

Inspecting a Model: summary

summary(runs.mdl)

Call:
lm(formula = runs ~ singles + doubles + triples + homeruns + 
    walks + hitbypitch + sacrificeflies + stolenbases + caughtstealing, 
    data = team.batting.00to08)

Residuals:
    Min      1Q  Median      3Q     Max 
-71.902 -11.828  -0.419  14.658  61.874 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -507.16020   32.34834 -15.678  < 2e-16 ***
singles           0.56705    0.02601  21.801  < 2e-16 ***
doubles           0.69110    0.05922  11.670  < 2e-16 ***
triples           1.15836    0.17309   6.692 1.34e-10 ***
homeruns          1.47439    0.05081  29.015  < 2e-16 ***
walks             0.30118    0.02309  13.041  < 2e-16 ***
hitbypitch        0.37750    0.11006   3.430 0.000702 ***
sacrificeflies    0.87218    0.19179   4.548 8.33e-06 ***
stolenbases       0.04369    0.05951   0.734 0.463487    
caughtstealing   -0.01533    0.15550  -0.099 0.921530    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 23.21 on 260 degrees of freedom
Multiple R-squared:  0.9144,    Adjusted R-squared:  0.9114 
F-statistic: 308.6 on 9 and 260 DF,  p-value: < 2.2e-16

The report: residual distribution; a coefficient table (estimate, std. error, t-value, p-value, significance stars); and fit statistics — R² = 0.914, adjusted R² = 0.911, an overwhelmingly significant F-statistic. Singles, doubles, triples, home runs, walks, HBP, sac flies all matter; stolen bases and caught-stealing do not.

Model Accessors

Beyond summary, a family of extractors:

formula(runs.mdl)          # the formula used
runs ~ singles + doubles + triples + homeruns + walks + hitbypitch + 
    sacrificeflies + stolenbases + caughtstealing
coef(runs.mdl)[1:4]        # coefficients (alias: coefficients)
 (Intercept)      singles      doubles      triples 
-507.1601976    0.5670487    0.6911042    1.1583609 
head(residuals(runs.mdl))  # residuals; fitted() gives fitted values
         1          2          3          4          5          6 
-35.572347 -10.526473  -5.920721  57.798284  -4.279294  -4.668995 
confint(runs.mdl)[1:4, ]   # confidence intervals for coefficients
                   2.5 %       97.5 %
(Intercept) -570.8582801 -443.4621151
singles        0.5158302    0.6182671
doubles        0.5744958    0.8077126
triples        0.8175297    1.4991921

Others: effects (orthogonal effects), vcov (variance-covariance matrix), deviance, influence / influence.measures, anova (ANOVA table).

Predicting New Values

predict(object, newdata, interval=c("none","confidence","prediction"), level=0.95, ...) applies the model to fresh data — optionally with intervals:

newdata <- team.batting.00to08[1, ]
predict(runs.mdl, newdata, interval="confidence")
       fit      lwr      upr
1 899.5723 890.1239 909.0208

anova(model1, model2, ...) compares nested models; update(model, formula) refits after a formula change — far cheaper than building from scratch on large data.

Diagnostic Plots

plot(lm.object) draws the standard diagnostics — residuals vs. fitted, a normal Q-Q of residuals, scale-location, and residuals vs. leverage:

par(mfrow=c(2,2))
plot(runs.mdl)
par(mfrow=c(1,1))

Look for: structure in residuals (nonlinearity), departures from the Q-Q line (non-normality), a fanning scale-location (heteroscedasticity), and high-leverage outliers.

lm in Detail

lm(formula, data, subset, weights, na.action,
   method = "qr", model = TRUE, x = FALSE, y = FALSE,
   qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, ...)

The model is \(y = c_0 + c_1 x_1 + \cdots + c_n x_n + \varepsilon\); OLS chooses the \(c_i\) minimizing the sum of squared residuals. Key arguments: formula (required); data, subset, weights, na.action; contrasts controls factor coding (treatment, Helmert, …); weights gives weighted least squares. Companions glm, nls, loess, rlm share many of these argument names.

When OLS Isn’t Enough

Assumptions and Their Tests

OLS is technically valid only under assumptions — yet often predicts well even when they’re bent:

  • linearity of the response in the predictors;
  • constant error variance (homoscedasticity) — test with car::ncvTest (non-constant variance score test);
  • predictors not collinear (singular fits);
  • normally distributed, independent errors.

When they fail badly, reach for the alternatives below.

Robust and Resistant Regression

Outliers wreck OLS. Two defenses (package MASS):

  • Resistant regression — lqs (least trimmed squares, least median of squares): fits using a subset of points, ignoring extreme ones;
  • Robust regression — rlm: M-estimation, down-weighting (rather than discarding) large residuals iteratively.
library(MASS)
lqs(runs~homeruns+walks, data=team.batting.00to08)   # resistant
rlm(runs~homeruns+walks, data=team.batting.00to08)   # robust

Subset Selection and Shrinkage

Too many predictors? Three remedies:

  • Stepwise selection — step (and MASS::stepAIC, wider model range) add/drop terms by AIC.
  • Ridge regression — MASS::lm.ridge(formula, data, lambda=) shrinks coefficients toward zero (L2 penalty), taming collinearity.
  • Lasso / least angle regression — lars::lars(x, y, type=c("lasso","lar","forward.stagewise","stepwise")) (L1 penalty) shrinks and selects, zeroing some coefficients; it computes the entire path at once.

Ridge and lasso are both penalized regression; the modern unified tool is glmnet.

step(runs.mdl)                          # stepwise by AIC
MASS::lm.ridge(runs~., team.batting.00to08, lambda=seq(0,1,0.1))

Principal Components and PLS Regression

When predictors are many and correlated, regress on components instead — package pls, both via mvr (aliases pcr, plsr):

  • Principal components regression (pcr) — regress on the leading principal components of the predictors;
  • Partial least squares (plsr) — components chosen to also correlate with the response.

Generalized and Nonlinear Models

glm: Generalized Linear Models

glm extends linear models to non-normal responses via a family and a link function — predicting the link of the expected response:

glm(formula, family = gaussian, data, weights, subset,
    na.action, start = NULL, ...)
Family Typical link Use
gaussian identity ordinary linear regression
binomial logit binary / proportion outcomes (Chapter 21)
poisson log counts
Gamma inverse positive continuous, skewed
inverse.gaussian 1/μ²

(glmnet also fits penalized GLMs — family, lambda, alpha choosing ridge↔︎lasso.)

Nonlinear Models

When the relationship isn’t linear in the parameters:

  • nls(formula, data, start, ...) — nonlinear least squares; you supply the functional form and starting values.
  • Survival models — survival::survreg (parametric) and coxph (Cox proportional hazards), reporting Wald tests on coefficients.
nls(y ~ a * exp(b * x), data=d, start=list(a=1, b=0.1))
library(survival)
coxph(Surv(time, status) ~ age + sex, data=lung)

Smoothing: Splines, Surfaces, Kernels

For flexible, nonparametric trends:

  • Splines — smooth.spline, and splines::ns/bs for use inside formulas;
  • Polynomial surfaces — loess (local regression) fits a smooth surface from local polynomial fits;
  • Kernel smoothing — ksmooth, density.
plot(cars$speed, cars$dist)
lines(smooth.spline(cars$speed, cars$dist), col="blue")
lines(loess.smooth(cars$speed, cars$dist), col="red")

Machine-Learning Regression

The chapter closes by previewing tree-based regression (full treatment: Chapters 21–22):

  • rpart — recursive partitioning (regression trees);
  • bagging, boosting (gbm), random forests (randomForest) — ensembles of trees, each highly parallelizable (Chapter 26);
  • PRIM (patient rule induction method) — bump hunting.

Tip

A workflow, not a function. Fit with lm, read summary, look at plot, then iterate with update — adding, dropping, or transforming terms. Reach for robust/penalized/nonlinear tools only when the diagnostics demand it.