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

Real data is clustered. Several measurements come from the same patient, plate, batch, field plot, school or machine. Observations within a cluster are more alike than observations from different clusters, which breaks the independence assumption behind OLS and ANOVA. Ignore it and your standard errors are too small, your p-values too optimistic, and your "significant" result may not replicate. **Mixed-effects models** (also called multilevel or hierarchical models) handle clustering directly, and they are now the default analysis for repeated measures, split-plot experiments and multi-site studies.

## Fixed effects and random effects

A mixed model contains both:

- **Fixed effects** are the effects you want to estimate and report: the treatment, the dose, the temperature. Their levels are chosen deliberately, and you would use the same levels again if you repeated the study.
- **Random effects** describe variation between clusters that are a sample from a larger population: these 12 batches out of all batches the plant will ever make, these 30 patients out of all patients. You are not interested in batch 7 for its own sake; you want to know **how much batches vary** and to account for it.

A useful test: if you repeated the study, would you use the **same** levels (fixed) or **new** ones drawn from the same population (random)? Treatments are fixed; batches, subjects and days are usually random.

## The random-intercept model

Suppose 20 cell-culture flasks are split between two media, and each flask is measured on 5 days. A random-intercept model is:

```text
y_ij = b0 + b1*medium_i + b2*day_ij + u_i + e_ij
u_i  ~ Normal(0, σ_u²)     flask-to-flask variation
e_ij ~ Normal(0, σ_e²)     day-to-day noise within a flask
```

Each flask gets its own baseline `b0 + u_i`, drawn from a normal distribution. Two variance components are estimated: `σ_u²` (between flasks) and `σ_e²` (within flasks). Their ratio gives the **intraclass correlation**:

```text
ICC = σ_u² / (σ_u² + σ_e²)
```

An ICC of 0.6 means 60% of the variance is between flasks: measurements from the same flask are strongly correlated, and the effective sample size for the medium comparison is much closer to 20 flasks than to 100 measurements.

## Random slopes

If the growth rate itself varies by flask, add a **random slope**: each flask gets its own deviation from the average slope for day, with its own variance (and usually a correlation with the random intercept). Random slopes are important when the fixed effect of interest varies across clusters; leaving them out can make tests of that effect anti-conservative.

## Fitting one in Python

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

# Random intercept per flask
m1 = smf.mixedlm("growth ~ C(medium) + day", data=df, groups=df["flask"]).fit()
print(m1.summary())       # fixed effects plus "Group Var" (σ_u²) and residual scale (σ_e²)

# Random intercept and random slope for day
m2 = smf.mixedlm("growth ~ C(medium) + day", data=df,
                 groups=df["flask"], re_formula="~day").fit()
```

In R, the usual tool is `lme4`: `lmer(growth ~ medium + day + (1 | flask))` for a random intercept and `(1 + day | flask)` for a random slope. For non-normal responses there are **generalised linear mixed models** (`glmer` in R, `BinomialBayesMixedGLM` or other packages in Python).

## Shrinkage: why random effects are "partially pooled"

There are three ways to handle the flasks:

1. **Complete pooling:** ignore flasks (plain OLS). Wrong standard errors.
2. **No pooling:** a separate fixed intercept for every flask. Uses many parameters, and flasks with few measurements get noisy estimates.
3. **Partial pooling (random effects):** each flask's estimated intercept is pulled toward the overall mean, more strongly when the flask has few or noisy measurements. This **shrinkage** is a feature: it borrows strength across clusters and usually predicts new data better.

## Where mixed models appear in experiments

- **Blocked designs with many blocks** (days, lots, sites): blocks as random effects.
- **Split-plot designs:** whole plots (an oven run, a field strip) are a random effect, so the hard-to-change factor is tested against whole-plot variation, not the smaller subplot error.
- **Repeated measures** over time on the same subjects or units.
- **Multi-site trials and inter-laboratory studies** (sites or labs as random effects; the variance components give repeatability and reproducibility).
- **Nested sampling** such as several aliquots per sample, several samples per batch.

## Practical cautions

- **You need enough clusters.** Estimating a variance from 3 batches is very imprecise; 5–6 is a bare minimum and more is much better. With very few clusters, treat them as fixed.
- **REML vs ML.** Variance components are usually estimated by restricted maximum likelihood (REML, the default). To compare models that differ in their **fixed** effects with a likelihood-ratio test, refit with ML.
- **Degrees of freedom** for fixed-effect tests are not obvious in mixed models. Methods such as Satterthwaite or Kenward–Roger (in R's `lmerTest` / `pbkrtest`) give better small-sample p-values than a naive z-test.
- **Convergence warnings** often mean the random-effects structure is too complex for the data. Simplify (drop a correlation or a random slope) and say so.

## Common mistakes

- Treating repeated measurements as independent replicates (pseudo-replication again).
- Making the treatment a random effect because it has "levels".
- Reporting only fixed effects and hiding the variance components, which are often the most useful result for process improvement.
- Using a random effect for a factor with two or three levels.

## Quick recap

- Fixed effects: levels of interest. Random effects: clusters sampled from a population.
- A random intercept models cluster-to-cluster baseline variation; the ICC says how strong it is.
- Mixed models give honest standard errors for clustered, repeated and split-plot data and shrink cluster estimates toward the mean.

## Try it

This lesson's **Do** tab is a short quiz on choosing fixed versus random effects and reading variance components.


Explain fixed and random effects and when clustered or repeated data need a mixed model.

Mixed models

Simple Explanation of Mixed Models (Hierarchical Linear Models, Multilevel Models)

Linear mixed effects models - the basics

## Why this matters

A two-level factorial tells you which factors matter and in which direction. It cannot tell you where the best setting is if the optimum lies **inside** the region, because two levels only see straight lines. **Response surface methodology (RSM)**, developed by Box and Wilson in 1951 for chemical processes, is the standard toolkit for the next step: move toward better conditions, then fit a curved model and find the optimum. It is used to tune reaction yields, formulation recipes, machining parameters, fermentation conditions and many other processes.

## The sequential strategy

RSM is a campaign, not a single experiment:

1. **Screen** many factors with a fractional factorial (module 2) and keep the few that matter.
2. **Climb.** Far from the optimum, the surface is roughly a tilted plane. Fit a first-order model and move along the **path of steepest ascent** (step sizes proportional to the coefficients) until the response stops improving.
3. **Check for curvature** with centre points. If the mean at the centre differs clearly from the mean of the factorial corners, you are near a peak (or valley).
4. **Fit a second-order model** with a design that supports it, and locate the optimum.
5. **Confirm** with new runs at the predicted best settings.

## The second-order model

For two coded factors `x1` and `x2`:

```text
y = b0 + b1*x1 + b2*x2 + b11*x1^2 + b22*x2^2 + b12*x1*x2 + error
```

Six coefficients. The squared terms give curvature, and the cross term tilts and rotates the surface. You fit it by ordinary least squares exactly as in module 1, with the extra columns added to the model matrix.

## Designs for a second-order model

A 2^k factorial cannot estimate the squared terms (every `x²` column is all ones). Two standard designs can:

- **Central composite design (CCD):** a 2^k factorial (the corners) plus `2k` **axial** (star) points at distance `α` along each axis, plus several **centre** runs. For two factors with `α = √2 ≈ 1.414`, that is 4 + 4 + (3 to 5) runs. The axial points add the levels needed for curvature; the centre runs give pure error and stabilise the prediction variance.
- **Box–Behnken design:** three levels per factor, no corner runs. Useful when the extreme corners are unsafe or impossible (for example, maximum temperature with maximum pressure).

```text
CCD, two factors, coded units
 x1      x2
-1      -1        factorial
+1      -1
-1      +1
+1      +1
-1.414   0        axial
+1.414   0
 0     -1.414
 0     +1.414
 0       0        centre (x3)
```

Coded units make the coefficients comparable and the design symmetric. Convert back to natural units at the end: `temp = 75 + 10 * x1` if 65–85 °C is coded −1 to +1.

## Finding the optimum

At the top of the surface both partial derivatives are zero:

```text
∂y/∂x1 = b1 + 2*b11*x1 + b12*x2 = 0
∂y/∂x2 = b2 + 2*b22*x2 + b12*x1 = 0
```

That is a 2 × 2 linear system:

```text
[ 2*b11   b12  ] [x1]     [ b1 ]
[  b12   2*b22 ] [x2] = − [ b2 ]
```

Its solution is the **stationary point**. Whether it is a maximum, minimum or saddle depends on the eigenvalues of the matrix `B = [[b11, b12/2], [b12/2, b22]]`:

- all negative: a **maximum**;
- all positive: a **minimum**;
- mixed signs: a **saddle** (the best point is on the edge of the region; use ridge analysis or constrained optimisation).

## Worked example

Suppose the fitted model in coded units is

```text
y = 80 + 4*x1 + 2*x2 - 3*x1^2 - 2*x2^2 + 1*x1*x2
```

The system is `[[-6, 1], [1, -4]] x = -[4, 2]`. Solving gives `x1 = 18/23 ≈ 0.78`, `x2 = 16/23 ≈ 0.70`. The eigenvalues of `B = [[-3, 0.5], [0.5, -2]]` are about −3.21 and −1.79, both negative, so it is a maximum. The predicted response there is `80 + 0.5 × (4 × 0.78 + 2 × 0.70) ≈ 82.3` (at a stationary point, `ŷ = b0 + ½ bᵀx`). The point lies inside the design region (|x| < 1.414), so the prediction is an interpolation and worth confirming with a few runs.

```python
import numpy as np
b1, b2, b11, b22, b12 = 4, 2, -3, -2, 1
x_s = np.linalg.solve([[2 * b11, b12], [b12, 2 * b22]], [-b1, -b2])
print(x_s, np.linalg.eigvalsh([[b11, b12 / 2], [b12 / 2, b22]]))
```

For contour plots, evaluate the fitted model on a grid (`np.meshgrid`) and draw it with `matplotlib.pyplot.contourf`; for several responses at once (yield high, impurity low), overlay the contours or use a **desirability function**.

## Common mistakes

- **Fitting a quadratic to a 2^k design without axial points**: the squared terms are not estimable.
- **Trusting an optimum outside the design region.** Extrapolated quadratics are unreliable; run another experiment centred nearer the new point.
- **Skipping confirmation runs.** The stationary point is a prediction with uncertainty.
- **Optimising one response and ignoring others** (cost, safety, impurity). Set constraints with the process owner first.
- **Forgetting lack-of-fit.** With replicated centre points, test whether the quadratic is adequate; if not, transform the response or narrow the region.

## Quick recap

- RSM is sequential: screen, climb by steepest ascent, check curvature, fit a quadratic, confirm.
- CCD and Box–Behnken designs support second-order models.
- The stationary point solves a linear system; eigenvalues tell maximum, minimum or saddle.

## Try it

In the **Do** tab you will fit a second-order model to an 11-run central composite design with NumPy and solve for its stationary point.


Fit a second-order response surface from a central composite design and locate its optimum.

Response surfaces

RSM: Introduction to Response Surface Methodology

Factorial vs fractional vs response surface designs | when to use what?

## Why this matters

The whole point of randomised experiments is that they support **causal** conclusions: the treatment caused the difference. But many questions cannot be randomised. You cannot assign smoking, a customer's region, or which supplier a plant chose ten years ago. Analysts still have to say something about cause and effect from observational data, and most of the famous statistical blunders come from doing it carelessly. **Causal inference** gives you a disciplined way to state your assumptions and decide which variables to adjust for. The main tool is a simple picture: the **directed acyclic graph (DAG)**.

## Potential outcomes in one paragraph

For each unit, imagine two outcomes: `Y(1)` if treated and `Y(0)` if not. The **causal effect** is the difference, but you only ever see one of them; the other is a counterfactual. Randomisation solves this on average: treated and untreated groups are comparable, so the difference in group means estimates the **average treatment effect**. Without randomisation, the groups may differ in other ways, and those differences, not the treatment, may explain the outcome.

## DAGs

A DAG is a diagram with an arrow from each cause to its effects, and no loops. You draw it from subject-matter knowledge **before** looking at the results.

```text
          Age
         /    \
        v      v
   Exercise --> Blood pressure
```

Three basic structures cover almost everything:

1. **Fork (confounder):** `X <- Z -> Y`. Z causes both. Example: age affects both exercise and blood pressure. The association between X and Y is partly spurious. **Adjust for Z** (include it in the model, stratify, or match on it).
2. **Chain (mediator):** `X -> M -> Y`. The effect of X flows through M. Example: a training programme raises skill, which raises output. If you want the total effect of X, **do not adjust for M**; adjusting removes the part of the effect you care about.
3. **Collider:** `X -> C <- Y`. C is caused by both. Example: a journal publishes papers that are either very novel or very rigorous. Among published papers (conditioning on C), novelty and rigour look negatively related even if they are unrelated in general. **Do not adjust for C**; doing so creates a spurious association (collider bias, also called selection bias).

## The backdoor criterion

A **backdoor path** is any path from X to Y that starts with an arrow **into** X. These paths carry non-causal association. To estimate the effect of X on Y, adjust for a set of variables that:

- blocks every backdoor path, and
- contains no descendant of X (no mediators, no consequences of treatment).

A path is blocked by adjusting for a fork or chain variable on it, and it is blocked *without* adjustment if it passes through a collider that you leave alone.

## Worked example

A materials lab finds that specimens annealed longer have **lower** hardness, and wonders whether annealing time should be cut. The DAG drawn with the process engineer:

```text
Alloy grade --> Anneal time
Alloy grade --> Hardness
Anneal time --> Grain size --> Hardness
Anneal time --> Hardness
Hardness --> Passed QC
Anneal time --> Passed QC
```

- **Alloy grade** is a confounder: softer grades are annealed longer by procedure. Adjust for it.
- **Grain size** is a mediator. For the total effect of annealing time, leave it out. (Include it only if you want the direct effect not through grain size, which needs extra assumptions.)
- **Passed QC** is a collider (caused by both hardness and time). Analysing only specimens that passed QC would bias the result. Use all specimens.

The adjustment model is `hardness ~ anneal_time + C(alloy_grade)`, not `hardness ~ anneal_time + grain_size + ...` on the QC-passed subset.

## Simpson's paradox

A trend can reverse when you split by a confounder. Classic case: a treatment looks worse overall but better within every severity group, because the sickest patients were more likely to receive it. The DAG tells you which answer is right: if severity causes both treatment choice and outcome, the stratified (adjusted) comparison answers the causal question. If the splitting variable were a collider or mediator, the overall comparison would be the right one. Data alone cannot decide; the causal structure does.

## Tools beyond adjustment

- **Randomised experiments** remain the gold standard when possible, including A/B tests.
- **Matching and propensity scores** balance measured confounders between groups.
- **Difference-in-differences**, **regression discontinuity** and **instrumental variables** use features of how treatment was assigned to handle some unmeasured confounding, under their own assumptions.
- Software: **DAGitty** (dagitty.net, free, in the browser) finds adjustment sets from a drawn DAG; in Python, `DoWhy` implements the same workflow.

All of these rely on assumptions you cannot fully test. State them, and do a **sensitivity analysis**: how strong would an unmeasured confounder need to be to explain away the result?

## Common mistakes

- **"Adjust for everything."** Adding mediators or colliders introduces bias rather than removing it.
- **Choosing covariates by p-value** or by which ones change the estimate the most.
- **Drawing the DAG after seeing the results** to justify a preferred model.
- **Using causal language** ("increases", "leads to") for an analysis that only supports association, without stating the assumptions.

## Quick recap

- Causal effects compare potential outcomes; randomisation makes groups comparable.
- DAGs encode assumptions: forks (adjust), chains (do not adjust for total effects), colliders (never adjust).
- Block all backdoor paths, avoid descendants of the treatment, and test sensitivity to unmeasured confounding.

## Try it

This lesson's **Do** tab is a quiz: you will read small DAGs and pick the right adjustment set.


Draw a causal DAG, spot confounders, mediators and colliders, and choose what to adjust for.

Causal inference intro

Introduction to Causal Graphs

Confounding Graphically - Causal Inference

Model clustered and repeated data with random effects, find optimal settings with response surface designs, and reason about cause and effect with DAGs.