Tutorial · Intermediate
Regression Models in R: Tests, Diagnostics, and Interpretation
A complete applied linear-regression workflow, from inspecting the data to reporting statistically defensible results.
Fuel economy in the built-in mtcars dataset
The working model is mpg ~ wt + hp. Every section explains what the procedure tests, how to run it, and what its result does—and does not—justify.
- Observations
- 32
- Outcome
- MPG
- Predictors
- Weight + HP
Visual guide
See the key ideas before the tests
These schematic figures show what to look for in the plots produced by R.Foundation
Prepare R and inspect the data
The tutorial models fuel economy (mpg) using vehicle weight (wt, in 1,000 lb) and horsepower (hp) from R’s built-in mtcars dataset. Begin by checking variable types, missing values, summaries, and impossible values.
data(mtcars)
head(mtcars[c("mpg", "wt", "hp")])
str(mtcars[c("mpg", "wt", "hp")])
summary(mtcars[c("mpg", "wt", "hp")])
colSums(is.na(mtcars[c("mpg", "wt", "hp")]))> str(mtcars[c("mpg", "wt", "hp")])
'data.frame': 32 obs. of 3 variables:
$ mpg: num 21 21 22.8 21.4 18.7 ...
$ wt : num 2.62 2.88 2.32 3.21 3.44 ...
$ hp : num 110 110 93 110 175 ...
> colSums(is.na(mtcars[c("mpg", "wt", "hp")]))
mpg wt hp
0 0 0The three model variables are numeric and contain no missing observations. A clean summary does not prove that the data are error-free; always compare ranges and units with the study context.
Foundation
Explore relationships before estimating
Scatterplots and correlations reveal direction, approximate form, unusual observations, and strong relationships among predictors before a model imposes a linear structure.
pairs(mtcars[c("mpg", "wt", "hp")],
pch = 19, col = adjustcolor("#2563eb", 0.55))
cor(mtcars[c("mpg", "wt", "hp")])
plot(mpg ~ wt, data = mtcars, pch = 19, col = "#2563eb")
abline(lm(mpg ~ wt, data = mtcars), col = "#dc2626", lwd = 2)> round(cor(mtcars[c("mpg", "wt", "hp")]), 3)
mpg wt hp
mpg 1.000 -0.868 -0.776
wt -0.868 1.000 0.659
hp -0.776 0.659 1.000
Graphical output: a scatterplot matrix and an mpg-versus-weight
scatterplot with a downward-sloping fitted line.In these data, mpg is strongly negatively correlated with weight (−0.868) and horsepower (−0.776). Weight and horsepower are positively correlated (0.659), so multicollinearity should later be assessed.
Foundation
Estimate the multiple linear regression
The model estimates the conditional mean of mpg as a linear function of weight and horsepower: mpgᵢ = β₀ + β₁wtᵢ + β₂hpᵢ + εᵢ.
model <- lm(mpg ~ wt + hp, data = mtcars)
summary(model)
coef(model)
fitted(model)
residuals(model)Call:
lm(formula = mpg ~ wt + hp, data = mtcars)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.22727 1.59879 23.285 < 2e-16
wt -3.87783 0.63273 -6.129 1.12e-06
hp -0.03177 0.00903 -3.519 0.00145
Residual standard error: 2.593 on 29 degrees of freedom
Multiple R-squared: 0.8268; Adjusted R-squared: 0.8148
F-statistic: 69.21 on 2 and 29 DF, p-value: 9.109e-12The fitted equation is mpg = 37.227 − 3.878(wt) − 0.0318(hp). Holding horsepower constant, an additional 1,000 lb is associated with about 3.88 fewer miles per gallon. Holding weight constant, one additional horsepower is associated with about 0.032 fewer mpg.
Statistical inference
Test individual coefficients with t-tests
Each coefficient t-test evaluates whether one partial regression coefficient differs from a hypothesized value while the other included predictors remain in the model.
For predictor j: H₀: βⱼ = 0 versus H₁: βⱼ ≠ 0. At α = 0.05, reject H₀ when p < 0.05, while also examining effect size and confidence intervals.
summary(model)$coefficients
# Test a non-zero value, for example H0: beta_wt = -3
b <- coef(model)["wt"]
se <- summary(model)$coefficients["wt", "Std. Error"]
t_value <- (b - (-3)) / se
p_value <- 2 * pt(abs(t_value), df.residual(model), lower.tail = FALSE)
c(t = t_value, p = p_value)> summary(model)$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.22727012 1.59878754 23.28469 2.565e-20
wt -3.87783074 0.63273349 -6.12870 1.120e-06
hp -0.03177295 0.00902971 -3.51871 1.451e-03
> c(t = t_value, p = p_value)
t p
-1.38736 0.17587For weight, t = −6.129 and p = 1.12×10⁻⁶; for horsepower, t = −3.519 and p = 0.00145. Both p-values are below 0.05, providing evidence that each partial association differs from zero in this model.
Statistical inference
Calculate and interpret confidence intervals
A confidence interval communicates a range of coefficient values compatible with the model and sampling uncertainty.
confint(model, level = 0.95)
# Manual 95% interval for the weight coefficient
b <- coef(model)["wt"]
se <- summary(model)$coefficients["wt", "Std. Error"]
b + c(-1, 1) * qt(0.975, df.residual(model)) * se> confint(model, level = 0.95)
2.5 % 97.5 %
(Intercept) 33.95738245 40.49715778
wt -5.17191604 -2.58374544
hp -0.05024078 -0.01330512The 95% interval is [−5.172, −2.584] for weight and [−0.0502, −0.0133] for horsepower. Neither includes zero, agreeing with the two-sided tests at the 5% level. For weight, plausible conditional associations range from roughly 2.58 to 5.17 fewer mpg per additional 1,000 lb.
Statistical inference
Evaluate overall model significance with the F-test
The overall F-test compares the fitted model with an intercept-only model and asks whether the predictors jointly improve explained variation.
H₀: βwt = βhp = 0 versus H₁: at least one slope is not zero.
summary(model)$fstatistic
null_model <- lm(mpg ~ 1, data = mtcars)
anova(null_model, model)Analysis of Variance Table
Model 1: mpg ~ 1
Model 2: mpg ~ wt + hp
Res.Df RSS Df Sum of Sq F Pr(>F)
1 31 1126.05
2 29 195.05 2 931 69.21 9.109e-12 ***F(2, 29) = 69.21 with p = 9.11×10⁻¹². Reject the joint null: weight and horsepower collectively improve the explanation of mpg relative to an intercept-only model.
Statistical inference
Conduct partial and joint F-tests
Nested-model tests assess whether a block of variables adds explanatory value after variables in a smaller model are already included.
For the example below, H₀: the coefficients on hp and am are jointly zero after controlling for wt.
restricted <- lm(mpg ~ wt, data = mtcars)
unrestricted <- lm(mpg ~ wt + hp + am, data = mtcars)
anova(restricted, unrestricted)
# Equivalent general linear hypotheses
install.packages("car") # Run once
library(car)
linearHypothesis(unrestricted, c("hp = 0", "am = 0"))Analysis of Variance Table
Model 1: mpg ~ wt
Model 2: mpg ~ wt + hp + am
Res.Df RSS Df Sum of Sq F Pr(>F)
1 30 278.32
2 28 180.29 2 98.034 7.613 0.00225 **
Linear hypothesis test:
hp = 0
am = 0
F = 7.613, Df = (2, 28), p-value = 0.00225If the nested-model p-value is below the chosen significance level, reject the restrictions and retain evidence that the tested block contributes jointly. If it is large, the sample does not provide enough evidence that the larger model improves fit.
Statistical inference
Assess model fit without overinterpreting R²
R² measures the sample proportion of outcome variation explained by the fitted values. Adjusted R² penalizes adding predictors that contribute little.
summary(model)$r.squared
summary(model)$adj.r.squared
sigma(model)
actual_vs_fitted <- data.frame(
actual = mtcars$mpg,
fitted = fitted(model),
residual = residuals(model)
)> summary(model)$r.squared
[1] 0.8267855
> summary(model)$adj.r.squared
[1] 0.8148396
> sigma(model)
[1] 2.593412R² = 0.8268 and adjusted R² = 0.8148, so the model explains about 82.7% of the sample variation in mpg. The residual standard error is 2.593 mpg, an estimate of the typical conditional error scale.
Model diagnostics
Check linearity and functional form
Residual-versus-fitted plots should show an approximately patternless cloud around zero. Curvature suggests omitted nonlinear terms or an unsuitable functional form.
Ramsey RESET: H₀: the linear specification is adequate against an alternative involving powers of the fitted values.
plot(model, which = 1)
abline(h = 0, lty = 2, col = "#dc2626")
install.packages("lmtest") # Run once
library(lmtest)
resettest(model, power = 2:3, type = "fitted")
# Example nonlinear extension
quadratic_model <- lm(mpg ~ wt + I(wt^2) + hp, data = mtcars)Graphical output: the Residuals vs Fitted plot displays residuals
around the horizontal zero line. A visible curved smooth would suggest
that the linear form is incomplete.
> resettest(model, power = 2:3, type = "fitted")
RESET test
RESET = 0.834, df1 = 2, df2 = 27, p-value = 0.445For the plot, look for systematic curves rather than isolated points. For RESET, a small p-value is evidence against the current functional form; a large p-value means the test did not detect that particular misspecification.
Model diagnostics
Examine residual normality
Normality is mainly relevant for exact small-sample t and F inference. The conditional mean can still be linear when residuals are non-normal.
Shapiro–Wilk: H₀: residuals follow a normal distribution versus H₁: they do not.
plot(model, which = 2) # Normal QQ plot
hist(residuals(model), breaks = 8, col = "#93c5fd",
main = "Model residuals", xlab = "Residual")
shapiro.test(residuals(model))> shapiro.test(residuals(model))
Shapiro-Wilk normality test
W = 0.9279, p-value = 0.0337
Graphical output: a residual histogram and a normal QQ plot.
The small p-value signals a detectable departure from normality.In the QQ plot, broad alignment with the line supports approximate normality; tail departures reveal skewness or heavy tails. For Shapiro–Wilk, p < 0.05 is evidence against normality, while p ≥ 0.05 is not proof of normality.
Model diagnostics
Test constant error variance
Heteroskedasticity means the conditional residual variance changes with predictors or fitted values. OLS coefficients can remain unbiased under exogeneity, but conventional standard errors may be wrong.
Breusch–Pagan: H₀: constant error variance versus H₁: variance depends on the regressors.
plot(model, which = 1)
plot(model, which = 3) # Scale-location plot
library(lmtest)
bptest(model)
# A more flexible White-style auxiliary test
bptest(model, ~ fitted(model) + I(fitted(model)^2))> bptest(model)
studentized Breusch-Pagan test
BP = 0.8807, df = 2, p-value = 0.6438
Graphical output: Residuals vs Fitted and Scale-Location plots.
Here, the test does not reject constant variance at the 5% level.A fan or funnel shape is visual evidence of changing variance. For Breusch–Pagan, p < 0.05 supports heteroskedasticity; p ≥ 0.05 means the test did not find sufficient evidence against constant variance.
Model diagnostics
Test residual autocorrelation when observations are ordered
Serial correlation is primarily a concern for time-ordered, spatially ordered, or clustered observations. It violates independence and can distort conventional standard errors.
Durbin–Watson: H₀: no first-order residual autocorrelation. The common two-sided alternative is non-zero first-order autocorrelation.
library(lmtest)
dwtest(model)
# Visual check when an ordering is substantively meaningful
plot(residuals(model), type = "o", pch = 19,
ylab = "Residual", xlab = "Observation order")
abline(h = 0, lty = 2)> dwtest(model)
Durbin-Watson test
DW = 1.3624, p-value = 0.0714
alternative hypothesis: true autocorrelation is greater than 0
Graphical output: residuals connected in the current row order.
That ordering is illustrative rather than temporal for mtcars.A Durbin–Watson statistic near 2 is consistent with little first-order autocorrelation. Values substantially below 2 suggest positive correlation; values above 2 suggest negative correlation. Use the p-value for the stated alternative.
Model diagnostics
Diagnose multicollinearity with VIF
Multicollinearity occurs when predictors contain overlapping information. It inflates coefficient uncertainty and can make estimates unstable without necessarily harming overall prediction.
install.packages("car") # Run once
library(car)
vif(model)
# Manual VIF for wt
auxiliary <- lm(wt ~ hp, data = mtcars)
1 / (1 - summary(auxiliary)$r.squared)> vif(model)
wt hp
1.766625 1.766625
> 1 / (1 - summary(auxiliary)$r.squared)
[1] 1.766625
Both VIF values are well below common screening thresholds of 5 or 10.VIF equals 1 when a predictor is unrelated to the others. Values above 5 deserve attention and values above 10 are often treated as serious, but these are diagnostic conventions rather than universal hypothesis-test cutoffs.
Model diagnostics
Identify leverage, outliers, and influential observations
Leverage measures unusual predictor combinations, studentized residuals flag unusual outcomes conditional on predictors, and Cook’s distance summarizes influence on the fitted model.
influence_table <- data.frame(
car = rownames(mtcars),
leverage = hatvalues(model),
studentized = rstudent(model),
cooks_d = cooks.distance(model)
)
influence_table[order(-influence_table$cooks_d), ][1:5, ]
plot(model, which = 4) # Cook's distance
plot(model, which = 5) # Residuals versus leverage> influence_table[order(-influence_table$cooks_d), ][1:5, ]
car leverage studentized cooks_d
17 Chrysler Imperial 0.0764 2.572 0.4236
20 Toyota Corolla 0.1067 2.129 0.2087
31 Maserati Bora 0.3942 0.929 0.1574
28 Lotus Europa 0.1942 1.413 0.1295
16 Lincoln Continental 0.1097 -1.514 0.0830
Graphical output: Cook's-distance and residuals-versus-leverage plots.In this model, Chrysler Imperial has the largest Cook’s distance (about 0.424). Rules such as Cook’s D > 4/n or leverage > 2p/n are screening aids; they identify observations worth investigating, not automatic deletion targets.
Application
Use heteroskedasticity-robust standard errors
Robust covariance estimators keep the OLS coefficients unchanged but adjust estimated uncertainty when constant variance is doubtful.
install.packages(c("sandwich", "lmtest")) # Run once
library(sandwich)
library(lmtest)
coeftest(model, vcov. = vcovHC(model, type = "HC3"))
# Robust confidence intervals
robust_vcov <- vcovHC(model, type = "HC3")
robust_se <- sqrt(diag(robust_vcov))
critical <- qt(0.975, df.residual(model))
cbind(estimate = coef(model),
lower = coef(model) - critical * robust_se,
upper = coef(model) + critical * robust_se)> coeftest(model, vcov. = vcovHC(model, type = "HC3"))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.227270 2.229805 16.6953 < 2.2e-16
wt -3.877831 0.768519 -5.0458 2.23e-05
hp -0.031773 0.009385 -3.3855 0.00206
The coefficients are unchanged; only their estimated uncertainty changes.Compare robust and conventional standard errors. Material differences indicate that the homoskedastic formula was influential. Interpret estimates in the same units, but base tests and intervals on the selected robust covariance matrix.
Application
Compare candidate models responsibly
Nested F-tests evaluate restrictions, while adjusted R², AIC, BIC, and validation address different aspects of model adequacy and parsimony.
model_small <- lm(mpg ~ wt, data = mtcars)
model_full <- lm(mpg ~ wt + hp, data = mtcars)
anova(model_small, model_full) # Nested F-test
AIC(model_small, model_full)
BIC(model_small, model_full)
cbind(
adjusted_R2 = c(summary(model_small)$adj.r.squared,
summary(model_full)$adj.r.squared),
RMSE = c(sigma(model_small), sigma(model_full))
)> anova(model_small, model_full)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 30 278.32
2 29 195.05 1 83.27 12.381 0.001451 **
> AIC(model_small, model_full)
df AIC
model_small 3 166.029
model_full 4 156.652
> BIC(model_small, model_full)
df BIC
model_small 3 170.426
model_full 4 162.515For nested models, a small ANOVA p-value supports the added terms. Smaller AIC or BIC indicates a preferred trade-off between fit and complexity within the compared set. These criteria do not prove causal validity or guarantee external performance.
Application
Generate predictions and distinguish interval types
A confidence interval estimates the mean response for cars with specified predictors. A prediction interval targets one new car and is wider because it includes individual residual variation.
new_car <- data.frame(wt = 3, hp = 120)
predict(model, newdata = new_car,
interval = "confidence", level = 0.95)
predict(model, newdata = new_car,
interval = "prediction", level = 0.95)> predict(model, newdata = new_car, interval = "confidence")
fit lwr upr
1 21.78102 20.77178 22.79027
> predict(model, newdata = new_car, interval = "prediction")
fit lwr upr
1 21.78102 16.38174 27.18031For wt = 3 and hp = 120, predicted fuel economy is 21.78 mpg. The 95% confidence interval for the conditional mean is [20.77, 22.79], while the 95% prediction interval for one new car is much wider: [16.38, 27.18].
Application
Report results with statistical and substantive meaning
A complete report states the model, sample, units, estimates, uncertainty, diagnostic evidence, limitations, and the exact interpretation supported by the research design.
# A compact publication-style table
install.packages("modelsummary") # Run once
library(modelsummary)
modelsummary(model,
statistic = "conf.int",
stars = TRUE,
gof_map = c("nobs", "r.squared", "adj.r.squared", "rmse"))Model summary
────────────────────────────────────────
Estimate 95% CI
(Intercept) 37.227 [33.957, 40.497]
wt -3.878 [ -5.172, -2.584]
hp -0.032 [ -0.050, -0.013]
────────────────────────────────────────
Num.Obs. 32
R2 0.827
R2 Adj. 0.815
RMSE 2.593
────────────────────────────────────────Example: ‘Controlling for horsepower, a 1,000-lb increase in weight was associated with 3.88 fewer mpg (95% CI [−5.17, −2.58], p < 0.001). The model explained 82.7% of sample variation (adjusted R² = 0.815).’ Add diagnostic and design limitations before making broader claims.