# 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. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "description": "Exam score against weekly study hours with the simple least squares line", "width": 380, "height": 260, "data": {"values": [ {"hours": 2, "score": 61, "section": "in_person"}, {"hours": 3, "score": 55, "section": "online"}, {"hours": 4, "score": 64, "section": "in_person"}, {"hours": 5, "score": 66, "section": "online"}, {"hours": 5, "score": 72, "section": "in_person"}, {"hours": 6, "score": 65, "section": "online"}, {"hours": 6, "score": 66, "section": "online"}, {"hours": 7, "score": 81, "section": "in_person"}, {"hours": 7, "score": 69, "section": "online"}, {"hours": 8, "score": 69, "section": "in_person"}, {"hours": 8, "score": 75, "section": "online"}, {"hours": 9, "score": 82, "section": "in_person"}, {"hours": 9, "score": 75, "section": "online"}, {"hours": 10, "score": 86, "section": "in_person"}, {"hours": 10, "score": 77, "section": "online"}, {"hours": 11, "score": 89, "section": "in_person"}, {"hours": 12, "score": 83, "section": "online"}, {"hours": 12, "score": 92, "section": "in_person"}, {"hours": 13, "score": 93, "section": "online"}, {"hours": 14, "score": 94, "section": "in_person"} ]}, "layer": [ { "mark": {"type": "point", "filled": true, "size": 70}, "encoding": { "x": {"field": "hours", "type": "quantitative", "title": "Weekly study hours"}, "y": {"field": "score", "type": "quantitative", "title": "Exam score", "scale": {"domain": [50, 100]}}, "color": {"field": "section", "type": "nominal", "title": "Section"} } }, { "mark": {"type": "line", "color": "#555555"}, "transform": [{"regression": "score", "on": "hours"}], "encoding": { "x": {"field": "hours", "type": "quantitative"}, "y": {"field": "score", "type": "quantitative"} } } ] } ``` ## 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. ```mermaid flowchart TD A["Define outcome, predictors and reference levels"] --> B["Plot the data"] B --> C["Fit the OLS model"] C --> D["Residual vs fitted: curve?"] D -->|Yes| E["Add a term or transform, then refit"] E --> C D -->|No| F["Scale-location: fanning?"] F -->|Yes| G["Use HC3 robust SEs or transform y"] F -->|No| H["Q-Q plot and leverage and Cook's distance"] G --> H H --> I["VIF for collinearity"] I --> J["Interpret coefficients holding others constant"] J --> K["Report estimates, CIs, R2 and diagnostics"] ``` ## 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.** ```r # 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). ```text 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.** ```python 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]`). ```text 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. ```text 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. | 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.** - [ ] 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 | 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.** 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): ```text 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 | 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] |