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

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

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.

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?

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.

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.

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:

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

For simple regression (one predictor) this reduces to

Residual variance and standard errors.

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

Fit and the overall F test.

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

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

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

1. Slope by hand: $\hat{\beta}_1 = S_{xy}/S_{xx} = 680.30 / 216.95 =$ [ ].
2. Intercept by hand: $\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} = 75.70 - [\ ] \times 8.05 =$ [ ].
3. 95% CI for the slope with $t_{0.975,18} = 2.101$: $3.136 \pm 2.101 \times 0.304 =$ [ , ].
4. In simple regression, $F = t^2$: check $10.31^2 =$ [ ].
5. The simple slope (3.14) is larger than the adjusted slope in Example A (2.89). One reason: [ ].
6. 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.

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

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].

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

Results checklist.

• Outcome, predictors, units and reference levels stated
• Sample size n and any excluded rows (with reasons) reported
• Estimates with 95% CIs, not only p-values
• Overall F test and R² (and adjusted R²) reported
• Each assumption checked, with the plot or test named
• Influential points examined; results with and without reported if they matter
• Robust SEs used or justified as unnecessary
• Interactions interpreted with simple slopes, not main effects alone
• Causal language avoided for observational data
• Software and versions named

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

1. In Example A, interpret the coefficient -7.41 for sectiononline in one sentence, including the reference level and the phrase about the other predictors.
2. 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]?
3. A colleague adds 8 random noise predictors to a model and R² goes up. What happens to adjusted R², and why?
4. The scale-location plot fans out to the right. Name two remedies and say which one leaves the coefficients unchanged.
5. 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?
6. 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

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.

1. $680.30 / 216.95 = 3.136$.
2. $75.70 - 3.136 \times 8.05 = 75.70 - 25.24 = 50.46$.
3. $3.136 \pm 0.639 = [2.497, 3.775]$ (R confint: [2.497, 3.774]).
4. $10.31^2 = 106.3$, matching F = 106.4 (rounding).
5. 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.
6. 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.

1. 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.
2. 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.
3. 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.
4. 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.
5. $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).
6. 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