Fit and check multiple regression and GLMs, design randomised, blocked and factorial experiments, and turn analyses into clear, reproducible reports for the researchers you work with.

Regression & Experimental Design

## Why this matters

Almost every quantitative career this course serves asks the same question in different clothes: how does an outcome change when several inputs change at once? A chemist wants to know how yield depends on temperature *and* catalyst load. A materials scientist links hardness to carbon content and cooling rate. An operations research analyst relates delivery time to distance, traffic and driver experience. **Multiple linear regression** is the workhorse for all of them, and the ideas in this module (interpretation, diagnostics, generalised models) carry straight into the design of experiments later in the course.

## The model

With one predictor you fit a line. With `p` predictors you fit a flat surface (a hyperplane):

```text
y = b0 + b1*x1 + b2*x2 + ... + bp*xp + error
```

The coefficients are chosen by **ordinary least squares (OLS)**: they minimise the sum of squared residuals, `SSE = sum (y - fitted)^2`. In matrix form, stack a column of ones (for the intercept) next to the predictors to get the model matrix `A`. The OLS solution solves the normal equations:

```text
(Aᵀ A) b = Aᵀ y
```

In practice you never invert `AᵀA` yourself; numerical libraries use a QR or SVD factorisation, which is more stable. In NumPy that is `np.linalg.lstsq`.

## Interpreting a coefficient

The key phrase is **"holding the other predictors fixed"**. In a yield model

```text
yield = 12 + 0.9*temp - 2.5*catalyst
```

`0.9` means: among batches with the *same* catalyst concentration, each extra degree is associated with 0.9 percentage points more yield. It is not the slope you would get from plotting yield against temperature alone. If temperature and catalyst were correlated in the data (say, operators added more catalyst on hot days), the simple slope would mix the two effects together. Multiple regression separates them, as far as the data allow.

Two further points:

- **The intercept** is the predicted value when every predictor is 0. That is often far outside the data (0 °C in a reactor that runs at 60–90 °C), so it is estimated imprecisely and rarely meaningful on its own. Centring predictors (subtracting their mean) makes the intercept the prediction at an average batch.
- **Units matter.** A coefficient of 0.9 per °C is 1.62 per °F. To compare predictors on different scales, look at standardised coefficients or, better, at the change over a realistic range of each input.

## Worked example in Python

You will do this by hand with NumPy in the activity. In day-to-day work, `statsmodels` gives the full table with standard errors, t-statistics, p-values and confidence intervals:

```python
import pandas as pd
import statsmodels.formula.api as smf

df = pd.read_csv("pilot_plant.csv")        # columns: yield_pct, temp, cat
model = smf.ols("yield_pct ~ temp + cat", data=df).fit()
print(model.summary())
print(model.conf_int())                     # 95% intervals for each coefficient
```

A typical (abridged) summary for 40 runs:

```text
              coef    std err      t     P>|t|    [0.025   0.975]
Intercept   15.19      2.37      6.40    0.000    10.4     20.0
temp         0.859     0.031    27.95    0.000     0.80     0.92
cat         -2.594     0.243   -10.65    0.000    -3.09    -2.10
R-squared: 0.959     Adj. R-squared: 0.957
```

How to read it:

1. **Coefficient and interval.** Each extra °C adds about 0.86 points of yield (95% CI 0.80 to 0.92), holding catalyst fixed. The interval is narrow because temperature varied a lot.
2. **t and p.** `t = coef / std err`. A tiny p-value says the coefficient is clearly not zero *given the model*. It does not say the effect is large or causal.
3. **R² = 1 − SSE/SST** is the share of the variation in `y` the model explains. Adding any predictor, even pure noise, never lowers R². **Adjusted R²** penalises extra terms and is fairer for comparing models of different size.

## Categorical predictors and interactions

Regression handles categories through **indicator (dummy) variables**. In a formula, `C(supplier)` creates one 0/1 column per supplier except a baseline, and each coefficient is the difference from that baseline. An **interaction** such as `temp:cat` (or `temp * cat`, which adds both main effects too) lets the effect of temperature depend on the catalyst level. Interactions are exactly what factorial experiments are built to estimate, which you will see in module 2.

## Common mistakes

- **Reading coefficients causally from observational data.** "Holding fixed" is a statement about the model, not about an intervention. Unmeasured variables can still drive both `x` and `y` (lesson 9 covers confounding).
- **Ignoring collinearity.** When two predictors are highly correlated, their individual coefficients become unstable and their standard errors balloon, even though predictions may be fine. Check correlations or variance inflation factors (VIFs; values above about 5–10 are a warning).
- **Extrapolating.** The fitted plane is only supported where you have data. Predicting yield at 120 °C from runs between 60 and 90 °C is a guess.
- **Chasing R².** A high R² does not mean the model is right; a low R² does not mean the effect is unimportant. Look at residuals (next lesson) and at whether the coefficients make scientific sense.
- **Throwing in every available column.** More predictors means more variance, more chances of a spurious "significant" result and harder interpretation. Choose predictors with the subject-matter expert.

## Quick recap

- Multiple regression fits `y = b0 + b1*x1 + ... + bp*xp` by least squares.
- Each slope is the expected change in `y` per unit of that predictor, holding the others fixed.
- Use confidence intervals, not just p-values; use adjusted R² to compare model sizes.
- Watch for collinearity, extrapolation and causal over-reading.

## Try it

Open the **Do** tab. You will build the model matrix, solve least squares with NumPy and compute R² for 40 simulated pilot-plant runs, then compare your estimates with the true values used to generate the data.


Fit a regression with several predictors and interpret each coefficient, its interval and R-squared.

Multiple regression

Multiple Regression, Clearly Explained!!!

Understanding Multiple Linear Regression with Python’s Statsmodels Library

## Why this matters

Software will fit a regression to anything and print confident-looking p-values. Whether those numbers mean anything depends on assumptions the software does not check for you. **Diagnostics** are how you find out whether the model fits, whether a few points are steering the result, and whether the standard errors can be trusted. In a lab or plant, a diagnostic plot is often where you discover a mislabelled sample, a drifting instrument or a second process you did not know about.

## What OLS assumes

For the coefficient estimates *and* their standard errors to be trustworthy, the classic linear model assumes:

1. **Linearity:** the mean of `y` really is a linear function of the predictors (as written in the model, including any transformed or squared terms).
2. **Independence:** errors from different observations are unrelated (no time drift, no repeated measurements of the same unit treated as separate).
3. **Constant variance (homoscedasticity):** the spread of the errors is the same everywhere.
4. **Normality of errors:** mainly matters for small samples and for prediction intervals.

None of these is about the distribution of `x` or of `y` by itself. They are about the **residuals**, `e = y - fitted`.

## The four standard plots

```python
import statsmodels.formula.api as smf
import statsmodels.api as sm
import matplotlib.pyplot as plt

model = smf.ols("strength ~ hours", data=df).fit()
infl = model.get_influence()

plt.scatter(model.fittedvalues, model.resid)          # 1. residuals vs fitted
sm.qqplot(model.resid, line="45", fit=True)           # 2. normal Q-Q
plt.scatter(model.fittedvalues, abs(infl.resid_studentized_internal) ** 0.5)  # 3. scale-location
sm.graphics.influence_plot(model)                      # 4. leverage vs residuals
```

How to read them:

- **Residuals vs fitted.** You want a shapeless band around zero. A curve means a missing nonlinear term (try adding `x**2` or a log). A funnel that widens to the right means variance grows with the mean.
- **Normal Q-Q.** Points should follow the line. S-shapes indicate heavy or light tails; a few points far off the line are candidate outliers.
- **Scale-location.** A flat trend means constant variance. An upward trend confirms the funnel.
- **Residuals vs leverage.** Points in the corners (high leverage *and* large residual) are the influential ones.

Also plot residuals **against time or run order** and against any variable you left out. A trend over run order is a classic sign of instrument drift or operator learning, and it violates independence.

## Leverage, outliers and influence

These three words mean different things:

- An **outlier** has a large residual: the model predicts it badly.
- A point with high **leverage** has unusual predictor values. Leverage is the diagonal of the hat matrix `H = A (AᵀA)⁻¹ Aᵀ`, written `h_i`. Each `h_i` lies between 0 and 1, and they sum to `p`, the number of coefficients, so the average is `p/n`. A common flag is `h_i > 2p/n`.
- An **influential** point changes the fit noticeably when removed. It usually needs both leverage and a sizeable residual.

**Cook's distance** combines the two:

```text
D_i = e_i^2 / (p * s^2) * h_i / (1 - h_i)^2        where s^2 = SSE / (n - p)
```

It equals the (scaled) total change in all fitted values when observation `i` is deleted. Values above 1, or far above the rest, deserve a look. Some analysts use `4/n` as a gentler screen.

## Worked example: a mis-recorded specimen

Twenty-five tensile specimens were cured for 2–8 hours, except one cured for 20 hours whose strength was written down as 52 MPa instead of roughly 100. With `p = 2` and `n = 25`, the leverage cut-off is `2p/n = 0.16`.

```text
specimen   hours   leverage   Cook's D
  0..23     2–8    0.04–0.08   < 0.1
   24       20.0     0.79      38.4
```

Specimen 24 sits far to the right in `x`, so it has huge leverage; because its `y` is wrong, it drags the line down and gets a Cook's distance of 38. The fitted slope with it is about 0.54 MPa per hour; without it, about 2.96 (close to the true 3). One bad record changed the conclusion.

What should you do?

1. **Check the record**, not the statistics. Is there a lab notebook entry, a calibration log, a photo of the broken specimen?
2. If it is a **recording error**, correct it, document the correction and refit.
3. If the value is **real**, it is telling you something: maybe strength levels off after long cures, so a straight line is the wrong model. Do not delete it just because it is inconvenient.
4. Report results **with and without** the point when you cannot resolve it.

## Fixing common problems

| Symptom | Likely remedy |
|---|---|
| Curved residuals | Add a squared term, transform `x`, or use a spline |
| Funnel shape | Log-transform `y`, use weighted least squares, or a GLM (next lesson) |
| Unequal variance only | Heteroscedasticity-robust standard errors (`fit(cov_type="HC3")`) |
| Trend over run order | Add time as a predictor; randomise run order in future experiments |
| Repeated measures of the same unit | Mixed model (module 3) |

## Common mistakes

- **Testing normality of `y` instead of the residuals.** A bimodal `y` can come from a binary predictor and be perfectly fine.
- **Deleting points with large residuals until the fit looks good.** That is data dredging; it makes p-values meaningless.
- **Relying on formal tests alone.** With thousands of rows, a normality test rejects trivial departures; with ten rows, it misses serious ones. Look at the plots.
- **Forgetting that leverage depends only on `x`.** A high-leverage point that lies on the line is harmless, and it actually improves precision.

## Quick recap

- Diagnose with residual plots: vs fitted, Q-Q, scale-location, leverage, and run order.
- Leverage `h_i` measures unusual `x`; Cook's distance measures how much one point moves the fit.
- Investigate influential points; correct errors, keep real data and say what you did.

## Try it

In the **Do** tab you will compute leverages from the hat matrix and Cook's distance for the tensile data, and confirm that one specimen dominates the fit.


Check a regression's assumptions with residual plots and find high-leverage and influential points.

Diagnostics

Simple Linear Regression:  Checking Assumptions with Residual Plots

Leverage and Influential Points in Simple Linear Regression

## Why this matters

Ordinary regression assumes a continuous response with constant, roughly normal noise. Many responses are nothing like that. Did the assay pass or fail? Did the customer churn? How many defects are on this board, how many colonies on this plate, how many calls arrived in this hour? Forcing these into OLS gives predictions outside the possible range (probabilities above 1, negative counts) and wrong standard errors. **Generalised linear models (GLMs)** keep the familiar linear predictor but connect it to the response in a way that respects what the response is.

## Three ingredients

Every GLM has:

1. **A distribution** for `y` from the exponential family: normal, binomial (yes/no or successes out of `n`), Poisson (counts), gamma (positive, skewed amounts), and others.
2. **A linear predictor**, `η = b0 + b1*x1 + ... + bp*xp`, exactly as in multiple regression.
3. **A link function** `g` that connects the mean to the linear predictor: `g(mean of y) = η`.

| Response | Distribution | Usual link | Mean as a function of η |
|---|---|---|---|
| Continuous | Normal | identity | `μ = η` (this is OLS) |
| Yes/no | Binomial | logit | `p = 1 / (1 + exp(-η))` |
| Count | Poisson | log | `μ = exp(η)` |
| Positive, skewed | Gamma | log | `μ = exp(η)` |

The link guarantees the mean stays in range: a logistic model can never predict a probability of 1.3, and a Poisson model can never predict −2 defects.

GLMs are fitted by **maximum likelihood**, usually with an algorithm called **iteratively reweighted least squares (IRLS)**, which is Newton's method in disguise: repeatedly fit a weighted regression, update the weights, and stop when the coefficients settle. You will write it in the activity.

## Logistic regression: reading the coefficients

For a yes/no outcome, the model is linear in the **log-odds**:

```text
log( p / (1 - p) ) = b0 + b1*x1 + ...
```

So `exp(b1)` is an **odds ratio**: the multiplicative change in the odds of "yes" for a one-unit increase in `x1`, holding the rest fixed. If a logistic model for assay failure has `b_humidity = 0.12` per percentage point of humidity, then `exp(0.12) ≈ 1.13`: each extra point multiplies the odds of failure by about 1.13 (13% higher odds). Odds ratios are not risk ratios; when the outcome is rare they are close, but when it is common they are not.

## Poisson regression: reading the coefficients

For counts, the model is linear in the log of the mean:

```text
log(μ) = b0 + b1*x1 + ...      so      μ = exp(b0) * exp(b1)^x1 * ...
```

`exp(b1)` is a **rate ratio**. In the activity's solder example the slope is about 0.8 per unit of temperature offset, and `exp(0.8) ≈ 2.2`: moving one unit hotter roughly doubles the expected defects per board.

When counts come from different amounts of exposure (boards of different area, observation windows of different length), add `log(exposure)` as an **offset**, a term with its coefficient fixed at 1, so the model describes a rate per unit of exposure.

## Worked example with statsmodels

```python
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Logistic: pass/fail by humidity and operator
logit = smf.glm("failed ~ humidity + C(operator)", data=assays,
                family=sm.families.Binomial()).fit()
print(np.exp(logit.params))          # odds ratios
print(np.exp(logit.conf_int()))      # their 95% intervals

# Poisson with an exposure offset: defects per square cm of board
pois = smf.glm("defects ~ offset_temp + C(line)", data=boards,
               family=sm.families.Poisson(),
               offset=np.log(boards["area_cm2"])).fit()
print(np.exp(pois.params))           # rate ratios
print(pois.pearson_chi2 / pois.df_resid)   # dispersion check
```

Exponentiate both the coefficients and the interval ends; never exponentiate a standard error.

## Overdispersion

The Poisson distribution has variance equal to its mean. Real counts are often more variable (some boards come from a bad reel, some plates were contaminated). The quick check is the **dispersion ratio**: Pearson chi-square divided by residual degrees of freedom. Near 1 is fine. Well above 1 (say 1.5 or more) means the Poisson standard errors are too small and p-values too optimistic. Remedies:

- **Quasi-Poisson** (scale the standard errors by the dispersion), or
- **Negative binomial** regression (`sm.families.NegativeBinomial()`), which adds a parameter for extra variance, or
- a **mixed model** with a random effect for the clustering source (module 3).

Logistic models on grouped data (successes out of `n`) can be overdispersed in the same way.

## Model checking for GLMs

Ordinary residuals are not very useful because the variance changes with the mean. Use **deviance** or **Pearson residuals** instead (`model.resid_deviance`, `model.resid_pearson`) and plot them against fitted values and each predictor, as in the previous lesson. Compare nested models with a **likelihood-ratio test** (difference in deviance, chi-square with the difference in parameters) and non-nested ones with **AIC**.

## Common mistakes

- **Reading logistic coefficients as changes in probability.** They are changes in log-odds. To talk about probabilities, compute predicted probabilities at meaningful predictor values.
- **Ignoring exposure** in count models, so a line that simply ran longer looks worse.
- **Skipping the dispersion check** and reporting tiny p-values from an overdispersed Poisson fit.
- **Perfect separation** in logistic regression: if a predictor splits yes and no perfectly, the estimates run off to infinity. Software warns about it; consider penalised (Firth) regression or collect more varied data.

## Quick recap

- A GLM = a distribution + a linear predictor + a link function.
- Logistic (binomial, logit): `exp(b)` is an odds ratio. Poisson (log): `exp(b)` is a rate ratio.
- Use offsets for exposure; check overdispersion; use deviance residuals.

## Try it

In the **Do** tab you will turn a worked logistic IRLS loop into a Poisson regression, fit it to 500 simulated circuit boards and compute the rate ratio.


Model yes/no outcomes and counts with logistic and Poisson GLMs and interpret odds and rate ratios.

Generalised linear models

Understanding Generalized Linear Models (Logistic, Poisson, etc.)

Poisson regression - clearly explained

Fit and interpret multiple regression, check residuals, leverage and influence, and model yes/no outcomes and counts with logistic and Poisson GLMs.