Linear Regression Method Card Template (Statistics, R and Python)

Simple and multiple OLS regression: least squares, coefficient tests and CIs, R² and F test, dummies, interactions, diagnostics, HC3 robust SEs, prediction vs confidence intervals, with R and Python

AkademischVollständigLeast squares (Legendre 1805); Cook (1977); White (1980); MacKinnon and White (1985) HC3; NIST/SEMATECH e-Handbook ch. 4; APA JARS-Quant (name)stats linear regression mathAktualisiert 25.10.2026

Im Editor öffnen.md herunterladen

Anleitung: Use when a continuous outcome is modelled from one or more predictors. Define the outcome, predictors and reference levels first, check linearity, equal variance, normality, influence and VIF from the residual plots, then report estimates with 95% CIs and diagnostics. The main pitfall is interpreting a coefficient without 'holding the other predictors constant', or as a causal effect in observational data.

Vorschau

Linear Regression (OLS): [Your study question or course]

A method card for simple and multiple linear regression fitted by ordinary least squares (OLS): when to use it, the assumptions and how to check each one, the least squares formulas, coefficient tests and intervals, R² and the F test, dummy variables, interactions, diagnostics, robust standard errors and the difference between confidence and prediction intervals. It is for students preparing for a course or exam, thesis writers analysing their own data and analysts who need a reliable checklist. It speeds up learning by making you decide in a fixed order, check every assumption, solve a faded worked example and test yourself before you look at the answers. The worked examples use fictional data; replace them with your own. Based on standard regression theory (least squares, Gauss-Markov), the NIST/SEMATECH e-Handbook of Statistical Methods (public domain) and original method papers listed at the end.

Document control

Field	Value
Course or project	[STAT 210 Regression / thesis chapter 4]
Learner or analyst	[Your name]
Data set	[Example: Lakemoor University exam study, n = 20, fictional]
Software and versions	[Example: R 4.5.1 with car 3.1.5; Python 3.14 with numpy 2.5.4, pandas 3.1.0, scipy 1.18.1, statsmodels 0.15.0]
Date	YYYY-MM-DD
Version	[0.1]

When to use it

One-sentence rule: use linear regression when you want to describe, explain or predict a continuous outcome from one or more predictors, and the mean of the outcome changes roughly linearly with each predictor while the others are held constant.

Use when	Do not use when
The outcome is continuous (exam score, loaf volume, waiting time)	The outcome is binary (pass or fail): use logistic regression
You want the change in the mean outcome per unit of a predictor, adjusted for others	The outcome is a count of events (visits, defects): use Poisson or negative binomial regression
Observations are independent (one row per student, one per batch)	Rows are clustered or repeated (several exams per student, classes within schools): use a mixed model
The relationship is roughly linear, or can be made linear with a transformation or a polynomial term	The relationship is strongly curved and you have no theory for its shape: consider splines or a different model
You have enough rows: a common rule of thumb is at least 10 to 20 observations per estimated coefficient	You have more predictors than the data can support: simplify or use penalised regression

Alternatives: Pearson correlation if you only need the strength of a straight-line association between two variables; ANCOVA when the focus is a group comparison adjusted for a covariate (it is the same model with a different emphasis); robust regression or quantile regression when outliers or heavy tails dominate; generalised linear models (logistic, Poisson, negative binomial) for binary and count outcomes; mixed models for clustered data.

The question and the data

Example question: at Lakemoor University (fictional), how much does the final exam score in an introductory statistics course change with each extra weekly hour of study, after adjusting for average nightly sleep and for whether the student took the in-person or the online section?

Variable	Role	Measurement level	Unit or coding
`score`	Outcome	Continuous (interval)	Exam points, 0 to 100
`hours`	Predictor	Continuous (ratio)	Mean weekly study hours
`sleep`	Predictor (covariate)	Continuous (ratio)	Mean nightly sleep, hours
`section`	Predictor	Nominal, 2 levels	`in_person` (reference), `online`

Design: cross-sectional observational study, one row per student, sample size n = 20. Because students chose their own study hours, the coefficients describe adjusted associations, not causal effects.

Student	hours	sleep	section	score
1	2	6	in_person	61
2	3	7	online	55
3	4	5	in_person	64
4	5	8	online	66
5	5	6	in_person	72
6	6	7	online	65
7	6	6	online	66
8	7	8	in_person	81
9	7	7	online	69
10	8	5	in_person	69
11	8	7	online	75
12	9	6	in_person	82
13	9	8	online	75
14	10	7	in_person	86
15	10	6	online	77
16	11	8	in_person	89
17	12	7	online	83
18	12	6	in_person	92
19	13	8	online	93
20	14	7	in_person	94

The chart shows the raw scores against study hours, coloured by section, with the simple (one-predictor) least squares line. Notice that the online students (orange) sit below the line at most values of hours: that is the section effect the multiple regression will estimate.

Assumptions and how to check each

The usual memory aid is LINE: Linearity, Independence, Normality of the errors, Equal variance. Add two practical checks: no single observation drives the result, and predictors are not so correlated that their effects cannot be separated.

Assumption	How to check	R	Python	If violated
Linearity of the mean	Residuals vs fitted plot shows no curve; component-plus-residual plots	`plot(fit, which = 1)`; `car::crPlots(fit)`	`sm.graphics.plot_ccpr_grid(fit)`	Add a polynomial or spline term, transform the predictor, or model the curve
Independence	Study design: one row per unit, no clustering or time order	Check the design; `car::durbinWatsonTest(fit)` for time order	`statsmodels.stats.stattools.durbin_watson(fit.resid)`	Mixed model, cluster-robust SEs, or time series model
Normality of errors	Normal Q-Q plot of residuals; matters mostly for small n and for prediction intervals	`plot(fit, which = 2)`; `shapiro.test(resid(fit))`	`sm.qqplot(fit.resid, line="45", fit=True)`	Transform the outcome, use bootstrap CIs, or a GLM that fits the outcome type
Equal variance (homoscedasticity)	Scale-location plot is flat; residual spread does not fan out	`plot(fit, which = 3)`; `car::ncvTest(fit)`	`het_breuschpagan(fit.resid, fit.model.exog)`	HC3 robust SEs, transform the outcome, or weighted least squares
No overly influential points	Leverage and Cook's distance; residuals vs leverage plot	`plot(fit, which = 5)`; `cooks.distance(fit)`	`fit.get_influence().cooks_distance[0]`	Check the row for errors; report results with and without it; never delete silently
No severe multicollinearity	Variance inflation factor (VIF) for each predictor	`car::vif(fit)`	`variance_inflation_factor(X, i)`	Combine or drop redundant predictors, centre before forming interactions, collect more varied data

Rules of thumb (not laws): a VIF above 5 deserves attention and above 10 is serious; a Cook's distance above 4/n (here 4/20 = 0.20) flags a point to inspect, and a value above 1 is the classic cut-off for a highly influential point; leverage above 2p/n (p = number of coefficients including the intercept) is high.

How it works

In plain language. Linear regression draws the line (or, with several predictors, the flat surface) that makes the squared vertical distances from the points to the line as small as possible. Each slope tells you how much the predicted mean outcome changes when that predictor goes up by one unit while all other predictors in the model stay the same. The residuals are what the model cannot explain. Inference (standard errors, t tests, intervals) asks how much those slopes would wobble from sample to sample if the errors behave as assumed.

Model. For observation $i = 1, \dots, n$ with $k$ predictors:

y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik} + \varepsilon_i, \qquad \varepsilon_i \sim \text{i.i.d. } N(0, \sigma^2)

In matrix form, with $\mathbf{X}$ the $n \times p$ design matrix ( $p = k + 1$ , a first column of ones for the intercept):

\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}

For simple regression (one predictor) this reduces to

\hat{\beta}_1 = \frac{S_{xy}}{S_{xx}} = \frac{\sum_i (x_i - \bar{x})(y_i - \bar{y})}{\sum_i (x_i - \bar{x})^2}, \qquad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}

Residual variance and standard errors.

\hat{\sigma}^2 = \frac{\text{SSE}}{n - p} = \frac{\sum_i (y_i - \hat{y}_i)^2}{n - p}, \qquad \widehat{\text{Var}}(\hat{\boldsymbol{\beta}}) = \hat{\sigma}^2 (\mathbf{X}^\top \mathbf{X})^{-1}

Test and interval for one coefficient (df = n - p):

t_j = \frac{\hat{\beta}_j - 0}{SE(\hat{\beta}_j)}, \qquad \hat{\beta}_j \pm t_{0.975,\, n-p} \cdot SE(\hat{\beta}_j)

Fit and the overall F test.

R^2 = 1 - \frac{\text{SSE}}{\text{SST}}, \qquad R^2_{\text{adj}} = 1 - (1 - R^2)\frac{n - 1}{n - p}, \qquad F = \frac{(\text{SST} - \text{SSE})/k}{\text{SSE}/(n - p)}

Robust (HC3) covariance, which does not assume equal variance, where $h_{ii}$ is the leverage of observation $i$ and $e_i$ its residual:

\widehat{\text{Var}}_{\text{HC3}}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \operatorname{diag}\!\left(\frac{e_i^2}{(1 - h_{ii})^2}\right) \mathbf{X} (\mathbf{X}^\top \mathbf{X})^{-1}

Confidence vs prediction interval at a new predictor row $\mathbf{x}_0$ :

\hat{y}_0 \pm t_{0.975,\,n-p}\, \hat{\sigma}\sqrt{\mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0} \quad \text{(mean)}, \qquad \hat{y}_0 \pm t_{0.975,\,n-p}\, \hat{\sigma}\sqrt{1 + \mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0} \quad \text{(new individual)}

Symbol	Meaning
$y_i$ , $\hat{y}_i$	Observed and fitted outcome for observation $i$
$x_{ij}$	Value of predictor $j$ for observation $i$
$\beta_0$	Intercept: mean outcome when every predictor is 0 (and categorical predictors are at their reference level)
$\beta_j$	Slope: change in mean outcome per one-unit increase in $x_j$ , holding the other predictors constant
$\varepsilon_i$ , $e_i$	Error (unobserved) and residual $e_i = y_i - \hat{y}_i$ (observed)
$n$ , $k$ , $p$	Sample size, number of predictors, number of coefficients ( $p = k + 1$ )
SSE, SST	Sum of squared residuals; total sum of squares $\sum (y_i - \bar{y})^2$
$\hat{\sigma}$	Residual standard error, the typical size of a residual
$h_{ii}$	Leverage: how unusual the predictor values of row $i$ are

Dummy variables and reference levels. A categorical predictor with $m$ levels enters the model as $m - 1$ indicator (0/1) columns. The level without a column is the reference: each dummy coefficient is the difference in mean outcome between that level and the reference, holding the other predictors constant. Here section has two levels, so one column sectiononline equals 1 for online students. R and statsmodels both pick the first level alphabetically unless you set it (factor(..., levels = ) in R, C(section, Treatment('in_person')) in Python).

Interactions and centring. An interaction term $x_1 \times x_2$ lets the slope of $x_1$ depend on $x_2$ . Once the interaction is in the model, the coefficient of $x_1$ alone is the slope when $x_2 = 0$ , which is often a meaningless value (zero hours of sleep, for example). Centring a continuous predictor (subtracting its mean, here 8.05 hours) makes the main effect read as the slope at the average value and reduces the correlation between the main effect and the interaction column. Centring does not change the interaction coefficient, the fitted values or R².

The flowchart is the order of work for every regression you fit.

Worked example A (full): multiple regression of exam score

Hypotheses for each slope. $H_0: \beta_j = 0$ (no linear association with the score once the other predictors are held constant) against $H_1: \beta_j \neq 0$ , two-sided, α = 0.05. The overall F test has $H_0: \beta_1 = \beta_2 = \beta_3 = 0$ .

R code.

# install.packages("car")  # run once, for vif() and hccm()
library(car)

d <- data.frame(
  hours   = c(2, 3, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 12, 12, 13, 14),
  sleep   = c(6, 7, 5, 8, 6, 7, 6, 8, 7, 5, 7, 6, 8, 7, 6, 8, 7, 6, 8, 7),
  section = factor(c("in_person", "online", "in_person", "online", "in_person",
                     "online", "online", "in_person", "online", "in_person",
                     "online", "in_person", "online", "in_person", "online",
                     "in_person", "online", "in_person", "online", "in_person"),
                   levels = c("in_person", "online")),
  score   = c(61, 55, 64, 66, 72, 65, 66, 81, 69, 69,
              75, 82, 75, 86, 77, 89, 83, 92, 93, 94))

fit <- lm(score ~ hours + sleep + section, data = d)
summary(fit)
confint(fit)

# Diagnostics
par(mfrow = c(2, 2)); plot(fit)          # residual vs fitted, Q-Q, scale-location, leverage
vif(fit)
ncvTest(fit)                              # Breusch-Pagan style test
shapiro.test(resid(fit))                  # weak check of residual normality
round(cooks.distance(fit), 3)

# Robust (HC3) standard errors with t tests on the residual df
se_hc3 <- sqrt(diag(hccm(fit, type = "hc3")))
t_hc3  <- coef(fit) / se_hc3
round(cbind(se_hc3, t_hc3, p = 2 * pt(-abs(t_hc3), df.residual(fit))), 4)

# Confidence vs prediction interval for an online student, 10 hours, 7 hours sleep
new <- data.frame(hours = 10, sleep = 7,
                  section = factor("online", levels = levels(d$section)))
predict(fit, new, interval = "confidence")
predict(fit, new, interval = "prediction")

R output (R 4.5.1, car 3.1.5; output abridged: R prints all 20 Cook's distances in row order, shown here as the four largest, and the bracketed interval labels are added for reading).

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    39.9427     4.5958   8.691 1.86e-07 ***
hours           2.8941     0.1976  14.644 1.09e-10 ***
sleep           2.3946     0.7434   3.221  0.00534 ** 
sectiononline  -7.4080     1.3469  -5.500 4.85e-05 ***

Residual standard error: 2.757 on 16 degrees of freedom
Multiple R-squared:  0.9512,	Adjusted R-squared:  0.9421 
F-statistic:   104 on 3 and 16 DF,  p-value: 1.044e-10

                    2.5 %    97.5 %
(Intercept)    30.2000561 49.685300
hours           2.4751664  3.313090
sleep           0.8186002  3.970610
sectiononline -10.2633362 -4.552634

   hours    sleep  section 
1.114784 1.290612 1.193360 

Non-constant Variance Score Test 
Chisquare = 0.4455617, Df = 1, p = 0.50445

	Shapiro-Wilk normality test
W = 0.96077, p-value = 0.5593

# four largest of round(cooks.distance(fit), 3)
   10    19    20     2 
0.501 0.181 0.158 0.104 

              se_hc3   t_hc3      p
(Intercept)   6.5764  6.0736 0.0000
hours         0.2229 12.9813 0.0000
sleep         0.9813  2.4402 0.0267
sectiononline 1.2072 -6.1367 0.0000

       fit      lwr     upr          (confidence)
1 78.23821 76.16382 80.3126
       fit      lwr      upr         (prediction)
1 78.23821 72.03635 84.44007

Python code.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import variance_inflation_factor

d = pd.DataFrame({
    "hours": [2, 3, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 12, 12, 13, 14],
    "sleep": [6, 7, 5, 8, 6, 7, 6, 8, 7, 5, 7, 6, 8, 7, 6, 8, 7, 6, 8, 7],
    "section": ["in_person", "online", "in_person", "online", "in_person",
                "online", "online", "in_person", "online", "in_person",
                "online", "in_person", "online", "in_person", "online",
                "in_person", "online", "in_person", "online", "in_person"],
    "score": [61, 55, 64, 66, 72, 65, 66, 81, 69, 69,
              75, 82, 75, 86, 77, 89, 83, 92, 93, 94]})

fit = smf.ols("score ~ hours + sleep + C(section, Treatment('in_person'))", data=d).fit()
print(fit.summary().tables[1])
print("R2 =", round(fit.rsquared, 4), " adj R2 =", round(fit.rsquared_adj, 4),
      " F =", round(fit.fvalue, 2), " p =", f"{fit.f_pvalue:.3g}")

X = fit.model.exog
print("VIF:", [round(float(variance_inflation_factor(X, i)), 3) for i in range(1, X.shape[1])])
print("Max Cook's D:", round(fit.get_influence().cooks_distance[0].max(), 3))

# HC3: fit(cov_type="HC3") uses z (normal) tests by default; use_t=True gives t tests like R
hc3 = fit.get_robustcov_results(cov_type="HC3", use_t=True)
print("HC3 SE:", np.round(hc3.bse, 4), " p:", np.round(hc3.pvalues, 4))

new = pd.DataFrame({"hours": [10], "sleep": [7], "section": ["online"]})
print(fit.get_prediction(new).summary_frame(alpha=0.05).round(3))

Python output (Python 3.14, statsmodels 0.15.0; output abridged: the long label C(section, Treatment('in_person'))[T.online] is shortened to section[T.online]).

                       coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------
Intercept           39.9427      4.596      8.691      0.000      30.200      49.685
section[T.online]   -7.4080      1.347     -5.500      0.000     -10.263      -4.553
hours                2.8941      0.198     14.644      0.000       2.475       3.313
sleep                2.3946      0.743      3.221      0.005       0.819       3.971
R2 = 0.9512  adj R2 = 0.9421  F = 104.04  p = 1.04e-10
VIF: [1.193, 1.115, 1.291]
Max Cook's D: 0.501
HC3 SE: [6.5764 1.2072 0.2229 0.9813]  p: [0.     0.     0.     0.0267]
     mean  mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  obs_ci_upper
0  78.238    0.979         76.164         80.313        72.036         84.44

R and Python agree on every estimate. Two differences are by design and worth remembering. First, statsmodels puts the categorical term before the numeric ones, so its VIF list is in the order section, hours, sleep (1.193, 1.115, 1.291), while car::vif prints hours, sleep, section. Second, robust standard errors requested with fit(cov_type="HC3") use normal (z) tests unless you add use_t=True (get_robustcov_results inherits t tests from the plain OLS fit); with z tests the HC3 p-value for sleep is 0.015 instead of the t-based 0.027 that R gives. The t version is the safer choice in small samples.

Checks. The residual vs fitted plot shows no curve; the score test for non-constant variance gives p = .504, and the scale-location plot is flat; the Q-Q plot is close to the line (Shapiro-Wilk W = 0.961, p = .559, a weak check at n = 20). All VIFs are below 1.3, so collinearity is not a problem. Student 10 (8 hours, only 5 hours of sleep, score 69) has Cook's distance 0.501, above the 4/n = 0.20 flag but below 1. Refitting without that student moves the sleep slope from 2.39 to 1.59 points (95% CI [0.21, 2.97]), so the sleep estimate is sensitive to one row, while the hours and section estimates barely change. That is reported, not hidden.

Calculation by hand for the key statistics. Using the output values:

$t$ for hours $= 2.8941 / 0.1976 = 14.65$ (software 14.644 from unrounded numbers), df = 20 - 4 = 16.
95% CI for hours: $t_{0.975,16} = 2.120$ , so $2.8941 \pm 2.120 \times 0.1976 = 2.8941 \pm 0.419 = [2.475, 3.313]$ .
SST = 2494.20 and SSE = 121.62 (from anova() against the null model), so $R^2 = 1 - 121.62 / 2494.20 = 0.9512$ and $R^2_{\text{adj}} = 1 - 0.0488 \times 19/16 = 0.9421$ .
$F = \dfrac{(2494.20 - 121.62)/3}{121.62/16} = \dfrac{790.86}{7.601} = 104.0$ .

Effect size with CI. The unstandardised slope is the primary effect size: each extra weekly study hour goes with 2.89 more exam points, 95% CI [2.48, 3.31], holding sleep and section constant. For the whole model, $R^2 = .95$ (adjusted .94). The online section scored 7.41 points lower, 95% CI [4.55, 10.26] lower, which is about 2.7 residual standard deviations (7.41 / 2.757), a large difference for this course.

Interpretation. With sleep and section held constant, every additional hour of weekly study is associated with about 2.9 more exam points; an extra hour of nightly sleep with about 2.4 more points (an estimate that depends noticeably on one student); and taking the online section with about 7.4 fewer points. These are associations in observational data: students who study more may differ in other ways. For an online student who studies 10 hours a week and sleeps 7 hours, the model predicts 78.2 points. We are 95% confident the mean score of all such students lies in [76.2, 80.3], but a single such student's score will fall in the much wider prediction interval [72.0, 84.4], because it includes individual variation around the mean.

APA-style report. A multiple linear regression predicted exam score from weekly study hours, nightly sleep and course section (reference: in-person), F(3, 16) = 104.04, p < .001, R² = .95, adjusted R² = .94. Study hours (b = 2.89, 95% CI [2.48, 3.31], t(16) = 14.64, p < .001) and sleep (b = 2.39, 95% CI [0.82, 3.97], t(16) = 3.22, p = .005) were positively associated with score, and online students scored lower than in-person students (b = -7.41, 95% CI [-10.26, -4.55], t(16) = -5.50, p < .001). Residual diagnostics showed no clear departures from linearity, normality or equal variance (all VIF < 1.3); HC3 robust standard errors gave the same conclusions (sleep: p = .027), and the sleep coefficient fell to 1.59 when one influential student (Cook's D = 0.50) was removed.

Worked example B (half solved): simple regression

Same students, one predictor only: score ~ hours. You are given these summary values: $\bar{x} = 8.05$ , $\bar{y} = 75.70$ , $S_{xx} = 216.95$ , $S_{xy} = 680.30$ , and this R output.

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   50.457      2.644   19.08 2.17e-13 ***
hours          3.136      0.304   10.31 5.54e-09 ***

Residual standard error: 4.478 on 18 degrees of freedom
Multiple R-squared:  0.8553,	Adjusted R-squared:  0.8472 
F-statistic: 106.4 on 1 and 18 DF,  p-value: 5.537e-09

Slope by hand: $\hat{\beta}_1 = S_{xy}/S_{xx} = 680.30 / 216.95 =$ [ ].
Intercept by hand: $\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} = 75.70 - [\ ] \times 8.05 =$ [ ].
95% CI for the slope with $t_{0.975,18} = 2.101$ : $3.136 \pm 2.101 \times 0.304 =$ [ , ].
In simple regression, $F = t^2$ : check $10.31^2 =$ [ ].
The simple slope (3.14) is larger than the adjusted slope in Example A (2.89). One reason: [ ].
Report sentence: Study hours predicted exam score, b = [ ], 95% CI [ ], t([ ]) = [ ], p [ ], R² = [ ].

Exercise C (unsolved)

Brightwell Bakery (fictional) tests 12 batches of sourdough. For each batch it records the proofing time in minutes, the dough hydration in percent and the loaf volume in millilitres.

Batch	proof (min)	hydration (%)	volume (mL)
1	30	62	1510
2	35	65	1555
3	40	68	1610
4	45	64	1590
5	50	70	1680
6	55	66	1650
7	60	72	1745
8	65	68	1710
9	70	63	1660
10	75	70	1760
11	80	66	1720
12	85	72	1800

Tasks: (a) fit volume ~ proof + hydration in R and Python; (b) interpret both slopes with "holding the other constant"; (c) give the 95% CI for the hydration slope; (d) give R² and adjusted R²; (e) give the 95% prediction interval for a batch proofed 60 minutes at 68% hydration and explain why it is wider than the confidence interval for the mean.

Reading software output

Output line (R `summary(lm)`)	statsmodels equivalent	Meaning	Report it?
`Estimate`	`coef`	$\hat{\beta}_j$ : change in mean outcome per unit of the predictor, others held constant	Yes, with unit
`Std. Error`	`std err`	Sampling uncertainty of the estimate	Usually in a table
`t value`, `Pr(>\|t\|)`	`t`, `P>\|t\|`	Test of $\beta_j = 0$ given the other predictors	Yes, t(df) and p
`confint(fit)`	`[0.025 0.975]`	95% CI for $\beta_j$	Yes, always
`Residual standard error`	`np.sqrt(fit.scale)`	$\hat{\sigma}$ , typical prediction error in outcome units	Often useful
`Multiple R-squared` / `Adjusted R-squared`	`R-squared` / `Adj. R-squared`	Share of variance explained; adjusted penalises extra predictors	Yes
`F-statistic ... p-value`	`F-statistic`, `Prob (F-statistic)`	Test that all slopes are zero together	Yes, F(k, n - p)
`(Intercept)`	`Intercept`	Mean outcome at all predictors 0 and reference levels	Only if 0 is meaningful
AIC from `AIC(fit)`	`fit.aic`	Model comparison on the same data	When comparing models

AIC trap. For the same model, R reports AIC = 102.86 and statsmodels reports 100.86. R counts the residual variance $\sigma^2$ as an estimated parameter and statsmodels does not, so the two differ by exactly 2. Compare AIC values only within one software package.

Reporting template

Fill-in sentence. A multiple linear regression predicted [outcome] from [predictors] (reference level: [level]), F([k], [n - p]) = [F], p [=/<] [p], R² = [R²], adjusted R² = [adj R²]. Holding the other predictors constant, [predictor 1] was associated with a [b]-[unit] [increase/decrease] in [outcome] per [unit] (95% CI [lower, upper], t([df]) = [t], p [=/<] [p]). Diagnostics: [linearity, variance, normality, influence, VIF result].

Predictor	b	SE	95% CI	t(df)	p
Example: Intercept	39.94	4.60	[30.20, 49.69]	8.69	< .001
Example: Study hours (per hour)	2.89	0.20	[2.48, 3.31]	14.64	< .001
Example: Sleep (per hour)	2.39	0.74	[0.82, 3.97]	3.22	.005
Example: Online (vs in-person)	-7.41	1.35	[-10.26, -4.55]	-5.50	< .001
[Predictor]	[ ]	[ ]	[ ]	[ ]	[ ]

Note under the table: n = [20]; R² = [.95]; residual SE = [2.76]; [robust SE type if used].

Results checklist.

Common mistakes and misconceptions

"The slope is the effect of hours." In observational data it is an adjusted association. Causal language needs a design (randomisation) or a defended causal model.
Forgetting "holding the others constant". The hours slope changed from 3.14 (alone) to 2.89 (with sleep and section). A coefficient only has meaning together with the list of other predictors in the model.
Checking normality of the raw outcome or predictors. The assumption is about the errors. A skewed predictor is fine; skewed residuals are the issue.
Reading a high R² as proof the model is right. R² = .95 here, yet one student still moves the sleep slope by a third. A curved relationship can also have a high R² with a straight line.
Comparing raw coefficients to say which predictor "matters more". Hours and sleep are both in hours but have different spreads. Use standardised coefficients or, better, a substantive comparison of realistic changes.
Interpreting a main effect alone when an interaction is in the model. With hours * section, the hours coefficient is the slope for the reference section only. Centre continuous predictors and report simple slopes.
Treating a prediction interval and a confidence interval as the same. The CI is for the mean of many similar students, the PI is for one student, and the PI is always wider.
Deleting influential points to "clean" the data. Check them for data entry errors; if they are real, keep them and report a sensitivity analysis.
Coding a categorical variable with more than two levels as one number (1, 2, 3). That forces equal steps between categories. Use dummies.
p = 0.005 for sleep does not mean a 0.5% chance that sleep has no association. It is the probability of data at least this extreme if the true slope were zero and the model were right.

Self-test questions

In Example A, interpret the coefficient -7.41 for sectiononline in one sentence, including the reference level and the phrase about the other predictors.
Why is the 95% prediction interval for one student [72.0, 84.4] so much wider than the confidence interval for the mean [76.2, 80.3]?
A colleague adds 8 random noise predictors to a model and R² goes up. What happens to adjusted R², and why?
The scale-location plot fans out to the right. Name two remedies and say which one leaves the coefficients unchanged.
In score ~ hours_c * section + sleep (with hours_c = hours - 8.05), the output gives hours_c = 2.64 and hours_c:sectiononline = 0.58. What is the estimated hours slope for online students?
Your software reports VIF = 12 for one predictor. What does that mean, and what is one thing you could do?

My notes and next steps

Item	Entry
Confidence with this method (1-5)	[3]
What I can now do	[Fit, check and report a multiple regression]
Still unclear	[Interactions between two continuous predictors]
Mistakes to log in academic/stats-error-log	[Forgot to set the reference level]
Next review date (spacing)	YYYY-MM-DD
Next template	[academic/stats-logistic-regression]

Sources and further reading

Structure aligned with the method card used across this statistics library; formulas are standard regression results written in our own notation.

Legendre, A.-M. (1805). Nouvelles méthodes pour la détermination des orbites des comètes (first published account of least squares).
Cook, R. D. (1977). Detection of influential observation in linear regression. Technometrics, 19(1), 15-18.
Breusch, T. S., and Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287-1294.
White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817-838.
MacKinnon, J. G., and White, H. (1985). Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties. Journal of Econometrics, 29(3), 305-325.
Long, J. S., and Ervin, L. H. (2000). Using heteroscedasticity consistent standard errors in the linear regression model. The American Statistician, 54(3), 217-224.
Aiken, L. S., and West, S. G. (1991). Multiple Regression: Testing and Interpreting Interactions. Sage.
Fox, J., and Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Sage (the car package).
NIST/SEMATECH e-Handbook of Statistical Methods, chapter 4 Process Modeling, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
Further reading by name: de Smith, M. J., Statistical Analysis Handbook (StatsRef, 2024 edition); OpenIntro Statistics (Diez, Çetinkaya-Rundel and Barr); APA Publication Manual (7th ed.) and APA JARS-Quant; Wasserstein, R. L., and Lazar, N. A. (2016), The ASA statement on p-values.

Answer key

Example B.

$680.30 / 216.95 = 3.136$ .
$75.70 - 3.136 \times 8.05 = 75.70 - 25.24 = 50.46$ .
$3.136 \pm 0.639 = [2.497, 3.775]$ (R confint: [2.497, 3.774]).
$10.31^2 = 106.3$ , matching F = 106.4 (rounding).
Omitted-variable bias: students who study more also sleep a little more in this sample (correlation r = 0.28 between hours and sleep), and sleep is positively related to score. Without sleep in the model, part of the sleep association is credited to hours. Section is almost uncorrelated with hours here (r = -0.05), so it mainly reduces the residual error.
b = 3.14, 95% CI [2.50, 3.77], t(18) = 10.31, p < .001, R² = .86.

Exercise C (R 4.5.1 output):

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 651.8469    59.1728   11.02 1.59e-06 ***
proof         3.2775     0.1766   18.55 1.76e-08 ***
hydration    12.2908     0.9422   13.04 3.77e-07 ***
Residual standard error: 9.281 on 9 degrees of freedom
Multiple R-squared:  0.9908,	Adjusted R-squared:  0.9887 
                 2.5 %     97.5 %
hydration    10.159418  14.422183
       fit      lwr      upr       (prediction, proof 60, hydration 68)
1 1684.269 1662.362 1706.177

(b) Each extra minute of proofing goes with 3.28 mL more volume at a fixed hydration; each extra percentage point of hydration goes with 12.29 mL more volume at a fixed proofing time. (c) 95% CI [10.16, 14.42] mL. (d) R² = .991, adjusted .989. (e) Predicted 1684 mL, 95% PI [1662, 1706]. The PI adds the batch-to-batch scatter ( $\hat{\sigma}$ = 9.3 mL) to the uncertainty in the fitted mean, so it is always wider than the CI for the mean.

Self-test.

Holding study hours and sleep constant, online students scored on average 7.41 points lower than in-person students (the reference level), 95% CI [4.55, 10.26] points lower.
The CI only reflects uncertainty in the estimated mean; the PI adds the variation of one student around that mean ( $\hat{\sigma}$ = 2.76 points), the "1 +" term under the square root.
Adjusted R² will usually go down (or rise only slightly by chance), because it charges for each extra coefficient through the factor $(n - 1)/(n - p)$ . R² can never decrease when predictors are added.
HC3 robust standard errors (coefficients unchanged, only SEs, tests and CIs change), or a transformation of the outcome such as log (changes the coefficients and their meaning); weighted least squares is a third option.
$2.64 + 0.58 = 3.22$ points per hour for online students (the interaction was not significant, p = .136, so do not over-interpret the difference).
The variance of that coefficient is 12 times larger than it would be if the predictor were uncorrelated with the others. Drop or combine a redundant predictor, centre before forming interactions or polynomials, or accept wide CIs if the predictor is needed for adjustment only.

Revision history

Version	Date	Author	Change
0.1	YYYY-MM-DD	[Name]	First draft from the template
[0.2]	YYYY-MM-DD	[Name]	[Added HC3 sensitivity analysis]