Refresher

STAT 6160 • Experimental Design • Fall 2026

Siddharth Vishwanath

2026-08-27

Agenda

\[ \newcommand{\pa}[1]{\left(#1\right)} \newcommand{\pb}[1]{\left\{#1\right\}} \newcommand{\pc}[1]{\left[#1\right]} \newcommand{\pA}[1]{\Big(#1\Big)} \newcommand{\pB}[1]{\Big\{#1\Big\}} \newcommand{\pC}[1]{\Big[#1\Big]} \newcommand{\eq}[1]{\begin{align}{#1}\end{align}} \newcommand{\eref}[1]{Eq.~\eqref{#1}} \newcommand{\pref}[1]{part~(\ref{#1})} \newcommand{\nn}{\nonumber} \newcommand{\n}{\nonumber\\} \newcommand{\f}{\frac} \newcommand{\defeq}{\stackrel{\small\bullet}{=}} \newcommand{\ra}{\rangle} \newcommand{\la}{\langle} \newcommand{\norm}[1]{\left\|#1\right\|} \newcommand{\abs}[1]{\left\lvert#1\right\rvert} \newcommand{\Abs}[1]{\Bigl\lvert#1\Bigr\rvert} \newcommand{\pr}{{\mathbb P}} \newcommand{\qr}{{\mathbb Q}} \newcommand{\xv}{{\boldsymbol{x}}} \newcommand{\av}{{\boldsymbol{a}}} \newcommand{\bv}{{\boldsymbol{b}}} \newcommand{\cv}{{\boldsymbol{c}}} \newcommand{\dv}{{\boldsymbol{d}}} \newcommand{\ev}{{\boldsymbol{e}}} \newcommand{\fv}{{\boldsymbol{f}}} \newcommand{\gv}{{\boldsymbol{g}}} \newcommand{\hv}{{\boldsymbol{h}}} \newcommand{\nv}{{\boldsymbol{n}}} \newcommand{\sv}{{\boldsymbol{s}}} \newcommand{\tv}{{\boldsymbol{t}}} \newcommand{\uv}{{\boldsymbol{u}}} \newcommand{\vv}{{\boldsymbol{v}}} \newcommand{\wv}{{\boldsymbol{w}}} \newcommand{\zerov}{{\mathbf{0}}} \newcommand{\onev}{{\mathbf{1}}} \newcommand{\phiv}{{\boldsymbol{\phi}}} \newcommand{\cc}{{\check{C}}} \newcommand{\xv}{{\boldsymbol{x}}} \newcommand{\Xv}{{\boldsymbol{X}\!}} \newcommand{\yv}{{\boldsymbol{y}}} \newcommand{\Yv}{{\boldsymbol{Y}}} \newcommand{\zv}{{\boldsymbol{z}}} \newcommand{\Zv}{{\boldsymbol{Z}}} \newcommand{\Iv}{{\boldsymbol{I}}} \newcommand{\Jv}{{\boldsymbol{J}}} \newcommand{\Cv}{{\boldsymbol{C}}} \newcommand{\Ev}{{\boldsymbol{E}}} \newcommand{\Fv}{{\boldsymbol{F}}} \newcommand{\Gv}{{\boldsymbol{G}}} \newcommand{\Hv}{{\boldsymbol{H}}} \newcommand{\alphav}{{\boldsymbol{\alpha}}} \newcommand{\epsilonv}{{\boldsymbol{\epsilon}}} \newcommand{\betav}{{\boldsymbol{\beta}}} \newcommand{\deltav}{{\boldsymbol{\delta}}} \newcommand{\gammav}{{\boldsymbol{\gamma}}} \newcommand{\etav}{{\boldsymbol{\eta}}} \newcommand{\piv}{{\boldsymbol{\pi}}} \newcommand{\thetav}{{\boldsymbol{\theta}}} \newcommand{\tauv}{{\boldsymbol{\tau}}} \newcommand{\muv}{{\boldsymbol{\mu}}} \newcommand{\argmin}{\mathop{\text{arg}\ \!\text{min}}} \newcommand{\arginf}{\mathop{\text{arg}\ \!\text{inf}}} \newcommand{\argmax}{\mathop{\text{arg}\ \!\text{max}}} \newcommand{\argsup}{\mathop{\text{arg}\ \!\text{sup}}} \newcommand{\f}{\frac} \newcommand{\e}{\epsilon} \newcommand{\b}{\bullet} \newcommand{\inv}{^{-1}} \newcommand{\A}{{\mathcal A}} \newcommand{\C}{{\mathbb C}} \newcommand{\E}{{\mathbb E}} \newcommand{\F}{{\mathcal{F}}} \newcommand{\N}{{\mathbb N}} \newcommand{\R}{{\mathbb R}} \newcommand{\Q}{{\mathbb Q}} \newcommand{\Z}{{\mathbb Z}} \newcommand{\X}{{\mathbb{X}}} \newcommand{\Y}{{\mathbb{Y}}} \newcommand{\G}{{\mathcal{G}}} \newcommand{\M}{{\mathcal{M}}} \newcommand{\Xn}{X_1, \dots, X_n} \newcommand{\Xm}{X_1, \dots, X_m} \newcommand{\Yn}{Y_1, \dots, Y_n} \newcommand{\Ym}{Y_1, \dots, Y_m} \newcommand{\Xb}{{\mathbb{X}}} \newcommand{\Yb}{{\mathbb{Y}}} \newcommand{\s}{{{\sigma}}} %%%%%%%%%%%%%%%%%%%%%%%%%%% \newcommand{\c}{^c} \newcommand{\bar}{\overline} \newcommand{\indep}{\perp\kern-6pt \perp} \newcommand{\pr}{{\mathbb P}} \newcommand{\R}{{\mathbb R}} \newcommand{\E}{{\mathbb E}} \newcommand{\var}{\text{Var}} \newcommand{\sd}{\text{SD}} \newcommand{\se}{\text{SE}} \newcommand{\ssres}{{\mathrm{SS}_{\mathrm{Res}}}} \newcommand{\sstrt}{{\mathrm{SS}_{\mathrm{Trt}}}} \newcommand{\sstot}{{\mathrm{SS}_{\mathrm{Tot}}}} \newcommand{\rank}{\mathop{\mathrm{rank}}} \newcommand{\supp}{\text{supp}} \newcommand{\Unif}{{\text{Unif}}} \newcommand{\Ber}{{\text{Ber}}} \newcommand{\Bin}{{\text{Bin}}} \newcommand{\Geo}{{\text{Geo}}} \newcommand{\Poi}{{\text{Poi}}} \newcommand{\Exp}{{\text{Exp}}} \newcommand{\Chisq}{{\chi^2}} \newcommand{\iid}{{\stackrel{iid}{\sim}}} \newcommand{\px}{p_{X}} \newcommand{\py}{p_{Y}} \newcommand{\pxy}{p_{X,Y}} \newcommand{\fx}{f_{X}} \newcommand{\fz}{f_{Z}} \newcommand{\fy}{f_{Y}} \newcommand{\fxy}{f_{X,Y}} \newcommand{\Fx}{F_{X}} \newcommand{\Fz}{F_{Z}} \newcommand{\Fy}{F_{Y}} \newcommand{\Fxy}{F_{X,Y}} \newcommand{\Xbar}{\overline{X}} \newcommand{\Ybar}{\overline{Y}} \newcommand{\Zbar}{\overline{Z}} \newcommand{\nXbar}{n \cdot \overline{X}} \newcommand{\nYbar}{n \cdot \overline{Y}} \newcommand{\nZbar}{n \cdot \overline{Z}} \newcommand{\Xn}{X_1, X_2, \dots, X_n} \newcommand{\Xm}{X_1, X_2, \dots, X_m} \newcommand{\Ym}{Y_1, Y_2, \dots, Y_m} \newcommand{\sumXn}{X_1 + X_2 + \dots + X_n} \newcommand{\la}{\ell_\alpha} \newcommand{\ua}{u_\alpha} \newcommand{\at}{{\alpha/2}} %%\renewcommand{\hat}{\widehat} \newcommand{\mux}{\mu_{X}} \newcommand{\muy}{\mu_{Y}} \newcommand{\sx}{\sigma_{X}} \newcommand{\sy}{\sigma_{Y}} \newcommand{\pvalue}{$p$-value} \newcommand{\Ho}{H_{0}} \newcommand{\Ha}{H_{a}} \]

§1. Setting up an inference problem

§2. Sampling distributions

§3. Confidence intervals

§4. Hypothesis testing

§5. Simple linear regression

§6. Multiple regression in matrix form

§7. Estimating the parameters

§8. ANOVA and general linear hypotheses

[1:00, 1.5 min] §1 to §4 run close to an hour once we stop for the two worked examples in §4; §5 to §8 fill the rest. That’s snug for 90 minutes with this much discussion built in. A few slides are flagged in their notes as compressible. Use those if we’re running long rather than rushing §7 and §8, which is where the term actually starts.

Say what is not covered: no measure theory, no asymptotic theory beyond the CLT, no Cramér-Rao. Say that §2 is the section to pay attention to if they only pay attention to one.

§1. Setting up an inference problem

The data collection process

Every inference problem starts with a description of how the data came to exist.

Assumption

\(X_1, \ldots, X_n \sim D(\theta)\), drawn independently from a distribution \(D\) with unknown parameter \(\theta\)

\(D\)	Assumption	\(\supp(X)\)	\(\theta\)
Bernoulli	\(\Xn \sim \Ber(p)\)	{0, 1}	\(p \in (0, 1)\)
Poisson	\(\Xn \sim \Poi(\lambda)\)	{0, 1, 2, …}	\(\lambda \in (0, \infty)\)
Normal	\(\Xn \sim \mathcal{N}(\mu, \sigma^2)\)	\(\R\)	\(\mu \in \R\), \(\sigma^2 \in (0, \infty)\)
Exponential	\(\Xn \sim \mathrm{Exp}(\lambda)\)	\((0, \infty)\)	\(\lambda \in (0, \infty)\)

Two separate claims are being made here:

1. a claim about the form of \(D\), and
2. a claim about the dependence structure, here independence.

In a designed experiment, the second claim is something you arrange through randomization rather than something you assume.

[2:30, 2 min] Point out that the assumption is a model of the data collection process, not of the data.

The last line is the connection to lecture 1: randomization is a physical procedure that produces the independence, rather than an assumption about it. Blocking and split plots later in the term will break the iid part, and we will have to write a different assumption.

The support determines the assumption

The support of \(X_i\), \(\supp(X)\), is the set of values it can take. It constrains which \(D\) is available.

Type	Example support	Typical \(D\)
Binary	\(\{0, 1\}\)	\(\mathrm{Ber}(\theta)\)
Count	\(\{0, 1, 2, \ldots\}\)	\(\mathrm{Poi}(\lambda)\), \(\mathrm{Bin}(n,\theta)\)
Categorical	\(\{\text{A}, \text{B}, \text{C}\}\)	\(\mathrm{Mult}(1, \boldsymbol{\pi})\)
Ordinal	\(\{\text{low}, \text{med}, \text{high}\}\)	ordered categorical
Positive continuous	\((0, \infty)\)	\(\mathrm{Exp}(\lambda)\), \(\mathrm{Gamma}\), lognormal
Continuous	\(\R\)	\(\mathcal{N}(\mu, \sigma^2)\)

The support narrows the choice of \(D\), but it does not determine it. Two distributions with the same support can disagree about everything else.

Categorical and ordinal responses have to be encoded before they enter a linear model. We will use one-hot and dummy encodings for this in §6.

Typos on a page is the stock Poisson example. Typos cluster: a bad paragraph has several. Which Poisson assumption does that break, and which way does it push \(\var(X)\) relative to \(\E[X]\)?

[4:30, 2.5 min] Work down the table quickly.

The row worth a comment is the ordinal one: coding low/med/high as 1/2/3 imposes equal spacing, which is an assumption most people make without noticing.

The fragment is the point of the slide. Support gets you as far as “count, so Poisson”, and no further. Poisson also asserts a constant rate and independent occurrences, and clustered typos break both, which gives \(\var(X) > \E[X]\) instead of the Poisson’s \(\var(X) = \E[X]\). Overdispersion is the name for it. Take answers before saying any of that; someone usually produces the clustering observation without the vocabulary.

Flag the encoding sentence, since factor levels in a designed experiment are exactly this problem and it comes back in §6 and §8.

Parameter and statistic

Population parameter

\(\theta\) is a fixed, unknown number describing the population. It does not depend on the sample.

Statistic

\(\hat\theta = \hat\theta(X_1, \ldots, X_n)\) is a function of the data. It is our best guess for \(\theta\).

The observed data \(x_1, \ldots, x_n\) are realizations of the random variables \(X_1, \ldots, X_n\). Because \(\hat\theta\) is a function of random variables, \(\hat\theta\) is itself a random variable.

\[ \theta \ \text{is fixed and unknown}, \qquad \hat\theta \ \text{is random and observed} \]

A clinic has exactly 300 patients and you measure the outcome on all 300. Is the mean of those 300 numbers a parameter or a statistic?

[7:00, 2 min] This distinction is the one that has to be solid before §2 and §3 make sense.

Write the last line on the board. Everything in §3 follows from which of the two symbols is the random one.

The fragment has no single right answer, which is why it is worth asking. If the population is those 300 patients, the number is \(\theta\): you have run a census, there is nothing to infer, and the standard error is zero. If you regard the 300 as a sample from some wider population of patients the clinic might treat, the same number is \(\hat\theta\) and it has a sampling distribution. The arithmetic is identical either way. What changes is the population you defined, and that is a modelling decision you make, not something the data announces.

Uncertainty is variance of the statistic

For a sample of size \(n\) from \(\mathcal{N}(\mu, \sigma^2)\) with \(\hat\mu = \bar X\):

\[ \E[\hat\mu] = \mu, \qquad \var(\hat\mu) = \frac{\sigma^2}{n}, \qquad \se(\hat\mu) = \frac{\sigma}{\sqrt{n}} \]

Standard error

The standard error of a statistic is the standard deviation of its distribution.

Reporting \(\hat\theta\) alone is incomplete. What makes it an inference is the attached statement about how much \(\hat\theta\) varies, and that statement comes from the distribution of \(\hat\theta\).

That \(\sigma^2/n\) used independence. If instead every pair of observations has correlation \(\rho > 0\), what is \(\var(\hat\mu)\), and what does it tend to as \(n \to \infty\)?

[9:00, 2 min] Note that \(\mathrm{Var}(\hat\mu)\) involves \(\sigma^2\), the population variance, and \(n\), the design. Both come back in §7 when we look at \(\mathrm{Var}(\hat\betav)\).

The fragment is worth the board. With equicorrelated errors \(\var(\hat\mu) = \frac{\sigma^2}{n}\left[1 + (n-1)\rho\right]\), which tends to \(\rho\sigma^2\), not to zero. So under clustering, more observations do not drive the uncertainty down. This is why measuring the same plot twenty times is not the same as twenty plots, and it is the reason blocking and split-plot designs need their own error terms later in the term.

Then move to §2, which is where the distribution of \(\hat\theta\) gets defined.

§2. Sampling distributions

Definition

Sampling distribution

The sampling distribution of \(\hat\theta\) is the distribution of \(\hat\theta\) when \(X_1, \ldots, X_n\) are drawn from \(D(\theta)\).

It depends on three things:

• the distribution \(D\),
• the parameter \(\theta\), and
• the sample size \(n\).

It does not depend on the data you happened to observe. It is a property of the data collection process.

[11:00, 1.5 min] Read the definition. Emphasise the three dependencies, since the third one is where design enters.

Note that changing \(n\) changes the sampling distribution without changing anything about the population. That is the whole basis of sample size calculations.

How to think about it

Call the process of collecting \(X_1, \ldots, X_n \sim D(\theta)\) an experiment. Repeat the experiment \(N\) times:

\[ \begin{aligned} \text{1st experiment:} \quad & X^{(1)}_1, \ldots, X^{(1)}_n \sim D(\theta) \ \longrightarrow \ \hat\theta^{(1)} \\ \text{2nd experiment:} \quad & X^{(2)}_1, \ldots, X^{(2)}_n \sim D(\theta) \ \longrightarrow \ \hat\theta^{(2)} \\ & \ \ \vdots \\ N\text{th experiment:} \quad & X^{(N)}_1, \ldots, X^{(N)}_n \sim D(\theta) \ \longrightarrow \ \hat\theta^{(N)} \end{aligned} \]

The histogram of \(\hat\theta^{(1)}, \ldots, \hat\theta^{(N)}\) approximates the sampling distribution of \(\hat\theta\).

In practice you run the experiment once. The other \(N-1\) experiments are hypothetical, and the sampling distribution is a statement about them.

On the next slide you can change both \(n\) and \(N\). Commit now: which of the two changes the width of the histogram of \(\hat\theta\), and which does not?

[12:30, 2.5 min] This is the operational definition. Spend time here.

Two indices to keep straight: \(n\) is the size of one sample, \(N\) is the number of repetitions. \(N\) is a device for describing the distribution; \(n\) is a property of the experiment.

Make them commit out loud before you open the widget, and hold them to it. Only \(n\) changes the width; \(N\) only makes the histogram a smoother picture of a distribution whose width was already fixed. Rooms split on this, and the split is exactly the confusion between the experiment and the description of the experiment. Having predicted, they watch the next slide for a reason.

Then go to the widget and let them watch it happen.

Simulating a sampling distribution

[15:00, 4.5 min] Left panel is one experiment, right panel is the sampling distribution built from \(N\) of them.

Hand the controls to the room: literally ask someone to call out a population and someone else a statistic, rather than driving it yourself.

Things to do:

1. Exponential population, sample mean, \(n=2\). The right panel is skewed. Raise \(n\) and it becomes symmetric. That is the CLT.
2. Change \(N\) with \(n\) fixed. The right panel gets smoother but does not change shape or width. \(N\) is resolution, \(n\) is the experiment.
3. Switch the statistic to the sample maximum. The normal reference is a bad fit at every \(n\). Not all statistics are asymptotically normal.
4. Read the reported se and check it against the board work: doubling \(n\) should move it by \(1/\sqrt2\), not \(1/2\).

Some sampling distributions

For \(X_1, \ldots, X_n\) iid:

\[ \begin{aligned} X_i \sim \mathrm{Ber}(\theta), \quad \hat\theta = \bar X \quad &\Longrightarrow \quad \hat\theta \sim \tfrac{1}{n}\mathrm{Bin}(n, \theta) \\[4pt] X_i \sim \mathcal{N}(\mu, \sigma^2), \quad \hat\mu = \bar X \quad &\Longrightarrow \quad \hat\mu \sim \mathcal{N}\!\pa{\mu, \tfrac{\sigma^2}{n}} \\[4pt] X_i \sim \mathcal{N}(\mu, \sigma^2), \quad \hat\sigma^2 = \tfrac{1}{n-1}\textstyle\sum_i (X_i - \bar X)^2 \quad &\Longrightarrow \quad \hat\sigma^2 \sim \tfrac{\sigma^2}{n-1}\chi^2_{n-1} \\[4pt] X_i \sim \mathrm{Exp}(\lambda), \quad \hat\lambda = \bar X \quad &\Longrightarrow \quad \hat\lambda \sim \tfrac{1}{n}\mathrm{Gamma}(n, n\lambda) \end{aligned} \]

In the normal case \(\bar X\) and \(\hat\sigma^2\) are independent, and

\[ \frac{\hat\mu - \mu}{\sqrt{\hat\sigma^2 / n}} \ \sim \ t_{n-1} \]

[19:30, 2 min] These are results they have seen. Read them off.

The last one is the one that matters for §3 and §4: replacing \(\sigma\) by \(\hat\sigma\) changes the sampling distribution from normal to \(t\). This is the whole difference between the \(Z\)-test and \(t\)-test we work through in §4.

The independence of \(\bar X\) and \(\hat\sigma^2\) is special to the normal distribution. It comes back in §7 as the independence of \(\hat\betav\) and the residual sum of squares.

The normal family

Let \(Z, Z_1, Z_2, \ldots\) be iid \(\mathcal{N}(0,1)\).

\[ \begin{aligned} \chi^2_\nu \ &= \ \sum_{i=1}^{\nu} Z_i^2 & &\E = \nu, \quad \mathrm{Var} = 2\nu \\[6pt] t_\nu \ &= \ \frac{Z}{\sqrt{\chi^2_\nu / \nu}}, \quad Z \indep \chi^2_\nu & &\text{symmetric, heavier tails than } \mathcal{N}(0,1) \\[6pt] F_{\nu_1, \nu_2} \ &= \ \frac{\chi^2_{\nu_1} / \nu_1}{\chi^2_{\nu_2} / \nu_2}, \quad \chi^2_{\nu_1} \indep \chi^2_{\nu_2} & &\text{ratio of two variance estimates} \end{aligned} \]

Useful identities:

\[ t_\nu^2 \ \stackrel{d}{=} \ F_{1, \nu}, \qquad \nu_1 F_{\nu_1, \infty} \ \stackrel{d}{=} \ \chi^2_{\nu_1}, \qquad t_\nu \ \xrightarrow{d} \ \mathcal{N}(0,1) \ \text{ as } \nu \to \infty \]

\(\chi^2\) appears for sums of squares, \(t\) for a single coefficient or contrast, \(F\) for comparing two sums of squares.

[21:30, 3 min] Define each as a function of normals rather than by its density.

The independence conditions in the \(t\) and \(F\) definitions are doing work. Most of the distribution theory in §7 goes into establishing them.

\(\chi^2\) shows up for sums of squares, \(t\) for a single coefficient or contrast, \(F\) for comparing two sums of squares. That is the whole of the ANOVA table in §8, and the \(t_\nu \to \mathcal{N}(0,1)\) limit is why the \(Z\)-test and the \(t\)-test in §4 agree once \(n\) is large.

Degrees of freedom: in each case \(\nu\) is the dimension of a subspace. We will see this concretely in §7.

Why the normal appears so often

Law of large numbers

\(\displaystyle \frac{1}{n}\sum_{i=1}^n X_i \longrightarrow \E[X]\) as \(n \to \infty\).

Central limit theorem

\(\displaystyle \frac{\bar X - \E[X]}{\sqrt{\mathrm{Var}(\bar X)}} \Longrightarrow \mathcal{N}(0,1)\) as \(n \to \infty\).

The LLN says the statistic converges to the parameter. The CLT says the sampling distribution of the statistic converges to a normal, whatever \(D\) is.

This is why normal-based intervals and tests are used for statistics that are averages, and why they are not automatic for statistics that are not.

[24:30, 2 min] State both, note the difference: LLN is about a number, CLT is about a distribution.

Refer back to the widget: the exponential population with the sample mean was the CLT, and the sample maximum was the case where it does not apply.

§3. Confidence intervals

Definition

Confidence interval

A \(100 \times (1-\alpha)\%\) confidence interval for \(\theta\) is a pair of statistics \([L, U]\) with \[\pr\pa{L \leq \theta \leq U} = 1 - \alpha\]

\(L\) and \(U\) are computed from the data, so they are random. \(\theta\) is fixed. The probability statement is about the interval, not about \(\theta\).

For a sample mean with \(\sigma^2\) unknown,

\[ \mathrm{CI}(\mu, \alpha) \ = \ \pc{\hat\mu - t_{n-1;\,\alpha/2}\cdot \frac{\hat\sigma}{\sqrt{n}}, \ \ \hat\mu + t_{n-1;\,\alpha/2}\cdot \frac{\hat\sigma}{\sqrt{n}}} \]

[26:30, 1.5 min] The interval is built from a quantile of the sampling distribution of \(\hat\mu\), which is why §2 came first.

If \(\sigma^2\) were known, the same construction with \(z_{\alpha/2}\) and \(\sigma/\sqrt{n}\), the interval version of the \(Z\)-test coming up in §4.

Coverage

The coverage probability is the proportion of repeated experiments in which the interval contains \(\theta\).

\[ \text{coverage} \ = \ \pr\pa{\theta \in [L, U]} \]

Interpretation:

If we repeat the experiment many times and compute \([L, U]\) each time, then \(100 \times (1-\alpha)\%\) of those intervals will contain \(\theta\).

This is the same repeated-experiments picture as §2. A confidence interval is a statement about the procedure, evaluated over the sampling distribution.

Two independent experiments each produce a 95% interval for the same \(\mu\). What is the probability that both cover \(\mu\)? That at least one does?

[28:00, 2 min] State the interpretation exactly. The common error is to read \(1-\alpha\) as a probability about \(\theta\) for the one interval in hand.

The fragment is arithmetic they can do in their heads, and the arithmetic is the argument: \(0.95^2 = 0.9025\) and \(1 - 0.05^2 = 0.9975\). Multiplying probabilities across experiments only makes sense because coverage is a property of the procedure, not a statement about \(\mu\). If anyone objects that \(\mu\) either is or is not in a given interval, agree: that is the point, and it is why the 95% has to attach to the repetitions rather than to this one interval.

Then simulate it.

Coverage, simulated

[30:00, 3.5 min] Left panel: one vertical bar per experiment, red when the interval misses. Right panel: the sampling distribution the intervals come from.

Ask a student to name \(n\) and \(\alpha\), another to pick the seed.

Things to do:

1. At \(n = 10\), \(\alpha = 0.05\), the reported coverage is close to \(0.95\).
2. Change the seed. Individual intervals move, coverage does not. This is the fragment from the last slide, made visible.
3. Switch to the \(z\) quantile at \(n = 5\). Coverage drops to roughly \(0.88\). This is the practical difference between \(z\) and \(t\).
4. Increase \(n\) with \(z\) selected. The gap closes, since \(t_\nu \to \mathcal{N}(0,1)\).
5. Change \(\alpha\) to \(0.10\). Intervals narrow, more red bars.

§4. Hypothesis testing

Anatomy of a test

Assumption

\(X_1, \ldots, X_n \sim D(\theta)\).

Population parameter

\(\theta\), the quantity the hypotheses are about.

Null hypothesis \(H_0\)

A statement about \(\theta\) under the most conservative assumption.

Alternative hypothesis \(H_a\)

The statement about \(\theta\) we want to show evidence in favour of.

[33:30, 1.5 min] Same anatomy they have seen. The purpose of listing it is that we fill it in twice in a few minutes (, ), and again in §8 for the general linear hypothesis.

Anatomy of a test, continued

Sample statistic

\(\hat\theta\), our best guess for \(\theta\).

Test statistic

\(\hat T = \hat T(X_1, \ldots, X_n)\), the function of the data used to decide between \(H_0\) and \(H_a\).

Rejection region

\(R(\hat T, \alpha)\), the set of values of \(\hat T\) for which we reject \(H_0\).

	\(H_0\) true	\(H_a\) true
Fail to reject \(H_0\)	Correct	Type II error
Reject \(H_0\)	Type I error	Correct

[35:00, 1.5 min] The rejection region is chosen so that the Type I error rate equals \(\alpha\). That choice requires the distribution of \(\hat T\) under \(H_0\), which is the next slide.

Designing a test statistic

A test statistic is built so that its value separates \(H_0\) from \(H_a\). For \(H_0 : \mu = \mu_0\) against \(H_a : \mu \neq \mu_0\):

The distribution of \(\hat T\) under \(H_0\) involves neither \(\mu_0\) nor \(\sigma\). Which of the two operations removed which?

\[ \hat T \ = \ \frac{\hat\mu - \mu_0}{\hat\sigma / \sqrt{n}} \]

Two things are happening:

1. Centering at \(\mu_0\), so that \(\hat T\) is near zero under \(H_0\) and away from zero under \(H_a\).
2. Scaling by the standard error, so that the distribution of \(\hat T\) under \(H_0\) does not depend on \(\sigma\).

The scaling is what makes the distribution of \(\hat T\) under \(H_0\) tabulatable: here \(\hat T \sim t_{n-1}\) regardless of \(\mu_0\) and \(\sigma\).

[36:30, 2.5 min] Reveal the formula first, then ask the fragment. The attribution is the content: centering at \(\mu_0\) is what removes \(\mu_0\), and dividing by the standard error is what removes \(\sigma\). Students who have only memorised the formula can rarely say which piece does which.

Centering and scaling are the two design choices in every test statistic in the course, including both worked examples coming up.

Without the scaling, the null distribution would depend on \(\sigma\) and there would be no single critical value. This is why test statistics are ratios.

In §8 the same two operations appear: the \(F\) statistic centres by fitting under \(H_0\) and scales by the residual mean square.

Calibration by the null sampling distribution

The rejection region is chosen from the sampling distribution of \(\hat T\) under \(H_0\):

\[ R(\hat T, \alpha) \ = \ \pb{\, t : \abs{t} > t_{n-1;\,\alpha/2} \,}, \qquad \pr_{H_0}\pa{\hat T \in R} = \alpha \]

\(p\)-value

\(p \ = \ \pr_{H_0}\pa{\abs{\hat T} \geq \abs{\hat t_{\text{obs}}}}\), the probability under \(H_0\) of a test statistic at least as extreme as the one observed.

Reject \(H_0\) when \(p < \alpha\). This is the same rule as \(\hat T \in R\).

[39:00, 2 min] Everything is computed under \(H_0\). The alternative enters only through the choice of which tail counts as extreme, one tail for \(H_a: \mu \neq \mu_0\)’s two-sided cousin is written here, but does a one-sided version, so flag that the region shape depends on \(H_a\).

Note the equivalence of the two rules, since students often hold them as separate procedures.

Confidence intervals and tests are also connected: the interval is the set of \(\mu_0\) that would not be rejected.

Power

Power

\(\text{power}(\theta) = \pr_\theta\pa{\text{reject } H_0}\), computed under a specific alternative \(\theta\).

Power depends on:

• the size of the departure from \(H_0\),
• the standard error of \(\hat\theta\), hence on \(\sigma^2\) and \(n\), and
• the level \(\alpha\).

The first is a property of the population. The second is a property of the design. Sample size and variance reduction are the two things an experimenter controls, and both act through the standard error.

Double \(n\), halve \(\sigma\), or double the effect size. Two of these move the power by the same amount. Which two, and why?

[41:00, 1.5 min] Keep this short. Power calculations for designed experiments get a full treatment later in the term.

For the fragment: power depends on the data only through \(\abs{\mu - \mu_0} \big/ (\sigma/\sqrt{n})\). Halving \(\sigma\) and doubling the effect size both double that quantity; doubling \(n\) multiplies it by \(\sqrt 2\) only. So the two that match are halving \(\sigma\) and doubling the effect, and replication is the weakest of the three per unit of effort. Since the effect size is usually not yours to set, variance reduction is where the leverage is, and that is most of the rest of this course.

Then move to the two worked examples.

Worked example: a one-sample \(Z\)-test

A filling line is supposed to dispense \(\mu_0 = 355\) mL per can. Years of process data put the fill-volume standard deviation at \(\sigma = 4\) mL, stable enough to treat as known. After a nozzle change, a sample of \(n = 25\) cans has mean fill \(\hat\mu = 353.2\) mL. Has the line drifted low?

Assumption	\(X_1, \ldots, X_{25} \sim \mathcal{N}(\mu, 4^2)\), \(\sigma\) known
Parameter / statistic	\(\mu\) / \(\hat\mu = \bar X = 353.2\)
\(H_0\) vs \(H_a\)	\(\mu = 355\) vs \(\mu < 355\)
Test statistic	\(\hat Z = \dfrac{\hat\mu - \mu_0}{\sigma/\sqrt n} \sim \mathcal{N}(0,1)\) under \(H_0\)

\[ \hat z \ = \ \frac{353.2 - 355}{4/\sqrt{25}} \ = \ \frac{-1.8}{0.8} \ = \ -2.25 \]

Suppose \(\sigma\) were not known, you substituted \(\hat\sigma\), but kept the \(\mathcal{N}(0,1)\) critical value. Would the true Type I error rate be above or below \(\alpha\)?

[42:30, 2.5 min] The anatomy from the four slides above, filled in. The only new move is \(\sigma\) known replacing \(\hat\sigma\), which turns the \(t_{n-1}\) into \(\mathcal{N}(0,1)\). That is the whole difference between a \(Z\)-test and a \(t\)-test.

Have someone compute \(\hat z\) out loud before revealing it; the arithmetic is clean on purpose.

The fragment: above \(\alpha\). The \(t_{n-1}\) distribution has heavier tails than \(\mathcal{N}(0,1)\), so the normal critical value cuts off more than \(\alpha\) of the true null distribution and the test over-rejects. They have already seen the size of the effect: the coverage widget in §3 reported about \(0.88\) instead of \(0.95\) at \(n = 5\) when the \(z\) quantile was selected. Same mistake, viewed from the interval side.

Worth saying plainly that a genuine \(Z\)-test is rare in practice. It is here because it isolates the \(\sigma\)-known case before the \(t\)-test complicates it.

A one-sample \(Z\)-test, continued

[45:00, 4 min] The black curve is the null sampling distribution from , the orange dashed line is the observed \(\hat z = -2.25\), red is the rejection region for whichever \(\alpha\) is selected.

Things to do:

1. Start at \(\alpha = 0.05\), one-sided (less). Reject.
2. Ask the room: for what \(\alpha\) does the decision flip? Walk the dropdown down: it flips between \(0.025\) and \(0.01\), since \(z_{0.01} = -2.326\) is just past \(-2.25\). This is the \(p\)-value from made visible: \(p \approx 0.0122\) is exactly the boundary.
3. Switch to two-sided. The critical value moves out to \(\pm1.96\) at \(\alpha=0.05\); ask why the one-sided test rejects more easily for the same \(\alpha\), the same tail area split differently.
4. This is the exercise from the old course notes this deck grew out of: “if \(H_0\) were true, what would the test statistic have looked like?” The black curve is the answer, and the red region is what “surprising” means once you fix \(\alpha\).

Worked example: a two-sample \(t\)-test

Two catalysts, independent trials, reaction yield in %.

	Catalyst A	Catalyst B
\(n\)	6	6
\(\bar y\)	71.5	68.3
\(s\)	2.1	2.4

Assumption

\(Y_{A1},\ldots,Y_{A6} \sim \mathcal{N}(\mu_A,\sigma^2)\) and \(Y_{B1},\ldots,Y_{B6} \sim \mathcal{N}(\mu_B,\sigma^2)\), independent, common \(\sigma^2\).

Parameter \(\theta = \mu_A - \mu_B\). Same two moves as §4: center the difference at \(0\), scale by its standard error.

\[ \hat T \ = \ \frac{\pa{\bar y_A - \bar y_B} - 0}{S_p\sqrt{\tfrac1{n_A}+\tfrac1{n_B}}} \ \sim \ t_{n_A+n_B-2}, \qquad S_p^2 \ = \ \frac{(n_A-1)S_A^2 + (n_B-1)S_B^2}{n_A+n_B-2} \]

Compute \(\hat T\) and decide, at \(\alpha = 0.05\). \(H_0: \mu_A = \mu_B\) against \(H_a: \mu_A \neq \mu_B\).

[49:00, 4 min] Send someone to the board for this one. It is the longest hand computation of the lecture, and worth doing in full rather than summarising. In order: \(S_p^2 = \tfrac{5(2.1^2)+5(2.4^2)}{10} = 5.085\), \(S_p = 2.255\), \(\mathrm{se} = 2.255\sqrt{1/6+1/6} = 1.302\), \(\hat t = 3.2/1.302 = 2.46\), \(df=10\).

Compare to \(t_{10;\,0.025} = 2.228\): \(2.46 > 2.228\), reject. \(p \approx 0.034\).

Two things to say once the arithmetic is done:

Pooling. \(S_p^2\) is a weighted average of the two sample variances, and it only makes sense under the common-\(\sigma^2\) assumption above. If that assumption looks shaky, software defaults to Welch’s version instead, which doesn’t pool and adjusts the degrees of freedom down (R and Python’s scipy/statsmodels both default to it). The formula isn’t worth board time; knowing it exists and when it kicks in is enough.

Where this is going. This is the \(g=2\) case of the one-way ANOVA in §8, and \(\hat t^2 = \hat F\) for \(q=1\) is the identity from §2. Everything in §8 is this same test statistic generalized to more than two groups.

[53:00, natural point for a one- or two-minute breather before switching gears to linear models.]

§5. Simple linear regression

The model

For pairs \((x_1, y_1), \ldots, (x_n, y_n)\),

\[ y_i \ = \ \beta_0 \ + \ \beta_1 x_i \ + \ \e_i, \qquad i = 1, \ldots, n \]

The assumptions, in the order they will matter:

\[ \begin{aligned} \text{A1.} \quad & \E[\e_i] = 0 && \text{the mean of } y_i \text{ is } \beta_0 + \beta_1 x_i \\[3pt] \text{A2.} \quad & \mathrm{Var}(\e_i) = \sigma^2 && \text{constant variance} \\[3pt] \text{A3.} \quad & \mathrm{Cov}(\e_i, \e_j) = 0, \ i \neq j && \text{uncorrelated errors} \\[3pt] \text{A4.} \quad & \e_i \sim \mathcal{N}(0, \sigma^2) && \text{normality} \end{aligned} \]

\(x_i\) is treated as fixed. The randomness is in \(\e_i\), so \(y_i\) is the random variable and \(\beta_0, \beta_1, \sigma^2\) are the parameters.

[53:00, 2 min] Connect back to: \(\beta_0, \beta_1\) are population parameters, \(\hat\beta_0, \hat\beta_1\) will be statistics with sampling distributions.

A1 to A3 are second-moment assumptions and give unbiasedness; §7 adds a brief optimality result on top of that (least squares turns out to be the best linear unbiased choice). We’ll flag it there, not prove it. A4 gives exact \(t\) and \(F\) distributions. We will need A4 only in §7.

In a designed experiment \(x_i\) being fixed is literally true, since we set it.

Fitted values and residuals

Given estimates \(\hat\beta_0, \hat\beta_1\):

\[ \begin{aligned} \text{fitted value} \quad & \hat y_i \ = \ \hat\beta_0 + \hat\beta_1 x_i \\[4pt] \text{residual} \quad & e_i \ = \ y_i - \hat y_i \\[4pt] \text{residual sum of squares} \quad & \ssres \ = \ \sum_{i=1}^n e_i^2 \end{aligned} \]

The residual \(e_i\) is an observable quantity computed from the fit. The error \(\e_i\) is unobservable. They are not the same thing.

Least squares chooses \(\hat\beta_0, \hat\beta_1\) to minimise \(\ssres\):

\[ (\hat\beta_0, \hat\beta_1) \ = \ \argmin_{b_0, b_1} \ \sum_{i=1}^n \pa{y_i - b_0 - b_1 x_i}^2 \]

[55:00, 1.5 min] The residual/error distinction is worth stating explicitly; it comes back when we look at residual diagnostics.

Note that the criterion is a choice. Least absolute deviations would give a different estimator. The reasons we use squares are that it has a closed form and, as we’ll flag in §7, it’s provably the best linear unbiased choice under A1 to A3 alone.

Fitting a line

[56:30, 3.5 min] Green segments are the residuals, orange is the line they control, grey dashed is the least squares fit.

Things to do:

1. Move \(b_0\) and \(b_1\) to minimise the reported \(\ssres\) by hand. It bottoms out when the orange line sits on the grey dashed one.
2. At the least squares fit, read the two sums under the title. Both are zero. Away from it, they are not. These are the normal equations, and we derive them in §7.
3. Raise \(\sigma\). The fitted line moves around more between seeds. That is the sampling distribution of \(\hat\beta_1\).
4. Raise \(n\) with \(\sigma\) fixed and cycle the seed. The line is more stable.

Decomposition of variability

At the least squares fit,

\[ \underbrace{\sum_{i=1}^n (y_i - \bar y)^2}_{\mathrm{SS_{tot}}} \ = \ \underbrace{\sum_{i=1}^n (\hat y_i - \bar y)^2}_{\mathrm{SS_{reg}}} \ + \ \underbrace{\sum_{i=1}^n (y_i - \hat y_i)^2}_{\ssres} \]

and

\[ R^2 \ = \ \frac{\mathrm{SS_{reg}}}{\mathrm{SS_{tot}}} \ = \ 1 - \frac{\ssres}{\mathrm{SS_{tot}}} \]

The cross term vanishes because \(\sum_i e_i = 0\) and \(\sum_i \hat y_i e_i = 0\) at the least squares fit. In §7 this is the statement that the residual vector is orthogonal to the fitted vector.

[60:00, 2 min] Actually do the cross-term cancellation on the board. It is two lines and a good first chalkboard moment for the linear-models half of the lecture.

\(R^2\) is the proportion of \(\mathrm{SS_{tot}}\) accounted for by the fit. In §8 the same decomposition, with more terms, is the ANOVA table.

Estimates and their sampling distributions

\[ \hat\beta_1 \ = \ \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 \]

Under A1 to A3,

\[ \E[\hat\beta_1] = \beta_1, \qquad \mathrm{Var}(\hat\beta_1) = \frac{\sigma^2}{\sum_i (x_i - \bar x)^2} \]

Under A4 as well,

\[ \frac{\hat\beta_1 - \beta_1}{\hat\sigma / \sqrt{\sum_i (x_i - \bar x)^2}} \ \sim \ t_{n-2}, \qquad \hat\sigma^2 = \frac{\ssres}{n-2} \]

[62:00, 2 min] This is and applied to a regression coefficient: parameter, statistic, sampling distribution, standard error, and structurally the same \(t\) statistic as §4, just with \(\hat\beta_1\) in place of a mean difference.

The variance formula is worth a comment for a design course. The denominator is the spread of the \(x_i\), which in a designed experiment is chosen. Spreading the design points out reduces \(\mathrm{Var}(\hat\beta_1)\).

\(n-2\) because two parameters were estimated. This becomes \(n-p\) in §7.

§6. Multiple regression in matrix form

The model

With \(p-1\) predictors,

\[ y_i \ = \ \beta_0 + \beta_1 x_{i1} + \cdots + \beta_{p-1} x_{i,p-1} + \e_i \]

Stacking the \(n\) equations:

\[ \underbrace{\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}}_{\yv,\ n \times 1} \ = \ \underbrace{\begin{pmatrix} 1 & x_{11} & \cdots & x_{1,p-1} \\ 1 & x_{21} & \cdots & x_{2,p-1} \\ \vdots & \vdots & & \vdots \\ 1 & x_{n1} & \cdots & x_{n,p-1} \end{pmatrix}}_{\Xv,\ n \times p} \underbrace{\begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{p-1}\end{pmatrix}}_{\betav,\ p \times 1} \ + \ \underbrace{\begin{pmatrix} \e_1 \\ \e_2 \\ \vdots \\ \e_n \end{pmatrix}}_{\epsilonv,\ n \times 1} \]

\[ \yv \ = \ \Xv\betav + \epsilonv \]

[64:00, 1.5 min] The matrix form is notation, not a new model. Row \(i\) is the equation from the previous slide.

\(\Xv\) is the design matrix. In this course it is the object we choose.

Assumptions in matrix form

\[ \begin{aligned} \text{A1--A3.} \quad & \E[\epsilonv] = \zerov, \quad \mathrm{Var}(\epsilonv) = \sigma^2\Iv_n \\[6pt] \text{A4.} \quad & \epsilonv \sim \mathcal{N}_n(\zerov, \sigma^2\Iv_n) \end{aligned} \]

so that \(\E[\yv] = \Xv\betav\) and \(\mathrm{Var}(\yv) = \sigma^2\Iv_n\).

Two facts used repeatedly, valid for any random vector and fixed \(\boldsymbol{A}, \bv\):

\[ \E[\boldsymbol{A}\yv + \bv] = \boldsymbol{A}\E[\yv] + \bv, \qquad \mathrm{Var}(\boldsymbol{A}\yv + \bv) = \boldsymbol{A}\,\mathrm{Var}(\yv)\,\boldsymbol{A}^\top \]

If \(\yv\) is multivariate normal then \(\boldsymbol{A}\yv + \bv\) is multivariate normal as well.

[65:30, 1 min] \(\sigma^2\Iv_n\) packs constant variance into the diagonal and uncorrelated errors into the off-diagonal.

The two identities require no distributional assumption. You have seen both before, move quickly. Only the last line requires normality. This is why the optimality result in §7 does not need A4.

Encoding a factor in \(\Xv\)

A treatment factor with \(g\) levels is encoded as columns of \(\Xv\). For \(g = 3\), one observation per level:

Cell means

\[ \Xv = \begin{pmatrix} 1&0&0\\ 0&1&0\\ 0&0&1 \end{pmatrix} \]

\(\betav = (\mu_1, \mu_2, \mu_3)^\top\)

Dummy coding

\[ \Xv = \begin{pmatrix} 1&0&0\\ 1&1&0\\ 1&0&1 \end{pmatrix} \]

\(\betav = (\mu_1,\ \mu_2 - \mu_1,\ \mu_3 - \mu_1)^\top\)

Sum to zero

\[ \Xv = \begin{pmatrix} 1&1&0\\ 1&0&1\\ 1&-1&-1 \end{pmatrix} \]

\(\betav = (\bar\mu,\ \mu_1 - \bar\mu,\ \mu_2 - \bar\mu)^\top\)

All three have the same column space, so they give the same fitted values, residuals and \(\ssres\). They differ in what the individual entries of \(\betav\) mean.

\[ \mathcal{C}(\Xv) \ \text{is the design.} \qquad \betav \ \text{is a coordinate system chosen for describing it.} \]

Group means are \((10, 14, 12)\). Write \(\betav\) under each coding. Which entries change if level 3 is relabelled as the baseline, and which do not?

[66:30, 2.5 min] This is the encoding question from §1, now inside \(\Xv\).

Let them work the fragment out before you confirm it: cell means gives \((10, 14, 12)\), dummy coding gives \((10, 4, 2)\), sum to zero gives \((12, -2, 2)\). Relabelling level 3 as the baseline changes the dummy-coded vector to \((12, -2, 2)\) and leaves cell means and sum to zero alone. The quantities that survive every relabelling are the differences between group means, which is exactly the estimability question in §8.

R and Python default to dummy coding; classical ANOVA texts use sum to zero. Reading a coefficient requires knowing which.

The \(g\)-level version with an intercept and \(g\) treatment effects is rank deficient. We come back to that in §8.

The plane of best fit

[69:00, 3.5 min] Points are the data, the orange surface is the plane \(b_0 + b_1x_1 + b_2x_2\) set by the three sliders, and the green segments are the residuals measured against that plane. The plot is a plotly scene, so drag it to rotate rather than reaching for a slider.

The data were generated with \(\betav = (2, 3, -1)\) exactly, so the sliders can be dialled onto the truth.

Things to do:

1. Start at \((0, 0, 0)\): a flat plane at zero, long residuals, large \(\ssres\). Move \(b_1\) up and watch the plane tilt along \(x_1\) and \(\ssres\) fall. Do \(b_2\) next and note it tilts the other way, since \(\beta_2 < 0\).
2. Get \(\ssres\) as low as you can by hand, then switch the radio to OLS. It will be lower than whatever they reached.
3. Now set the sliders to the true \((2, 3, -1)\) and compare. \(\ssres\) is \(61.2\) at the truth and \(59.0\) at OLS. Least squares fits this sample better than the true plane does, and that is not a paradox: least squares minimises \(\ssres\) by construction, so it must win on the observed data. The gap is the noise it has fitted.
4. Rotate the scene. The residual segments stay vertical, in the \(y\) direction. Least squares minimises vertical distance, not perpendicular distance to the plane.

The reported OLS fit is \((1.88, 2.89, -0.91)\) against a truth of \((2, 3, -1)\), which is the sampling variability from §5 showing up in three dimensions.

Compressible: if short on time, do steps 1 and 2 only, and state the step 3 comparison without dialling it in.

§7. Estimating the parameters

Least squares and the normal equations

\[ \hat\betav \ = \ \argmin_{\betav} \ \norm{\yv - \Xv\betav}^2 \]

Differentiating and setting to zero gives the normal equations

\[ \Xv^\top\!\Xv\,\hat\betav \ = \ \Xv^\top\!\yv \]

If \(\Xv\) has full column rank \(p\), then \(\Xv^\top\!\Xv\) is invertible and

\[ \hat\betav \ = \ \pa{\Xv^\top\!\Xv}^{-1}\Xv^\top\!\yv \]

Equivalently, \(\Xv^\top(\yv - \Xv\hat\betav) = \zerov\): the residual vector is orthogonal to every column of \(\Xv\). For simple regression these are the two sums from the widget in §5.

[72:30, 2 min] Derive the normal equations on the board. This is the chalkboard moment of §7, don’t shortcut it by reading the slide.

Then rewrite them as \(\Xv^\top\ev = \zerov\) and connect back to the widget: the first column of \(\Xv\) is the intercept, giving \(\sum_i e_i = 0\); the second is \(x\), giving \(\sum_i x_i e_i = 0\).

Fitted values as a projection

\[ \hat\yv \ = \ \Xv\hat\betav \ = \ \Hv\yv, \qquad \Hv \ = \ \Xv\pa{\Xv^\top\!\Xv}^{-1}\Xv^\top \]

\(\Hv\) is symmetric and idempotent, so it is the orthogonal projection onto the column space \(\mathcal{C}(\Xv)\):

\[ \Hv = \Hv^\top, \qquad \Hv^2 = \Hv, \qquad \Hv\Xv = \Xv \]

and \(\ev = (\Iv - \Hv)\yv\) is the projection onto \(\mathcal{C}(\Xv)^\perp\), with

\[ \mathrm{rank}(\Hv) = \mathrm{tr}(\Hv) = p, \qquad \mathrm{rank}(\Iv - \Hv) = n - p \]

Since \(\hat\yv \perp \ev\), \(\ \norm{\yv}^2 = \norm{\hat\yv}^2 + \norm{\ev}^2\), which is the decomposition from §5.

[74:30, 2 min] Fitting a linear model is projecting \(\yv\) onto the column space of \(\Xv\).

The ranks are where degrees of freedom come from: \(p\) for the fit, \(n-p\) for the residual. In §5 with \(p=2\) this was \(n-2\).

Note that \(\Hv\) depends only on \(\Xv\), not on \(\yv\). Everything about the precision of the fit is fixed once the design is chosen.

The multivariate normal distribution

For \(\xv \in \R^p\), if \(\boldsymbol{\Sigma}\) is positive definite,

\[ \Xv \sim \mathcal{N}_p(\muv, \boldsymbol{\Sigma}) \]

has density

\[ \begin{aligned} f_{\Xv}(\xv) &= \frac{1}{(2\pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}} \\[-0.2em] &\quad {}× \exp\pc{ -\frac{1}{2}(\xv-\muv)^\top \boldsymbol{\Sigma}^{-1}(\xv-\muv) }. \end{aligned} \]

The contours of equal density are ellipsoids centered at \(\muv\); their shape and orientation are determined by \(\boldsymbol{\Sigma}\).

[76:30, 2 min] This is the multivariate version of the bell curve. Point to the quadratic form in the exponent: fixing its value gives an ellipse, and the covariance matrix controls its axes. The plot is the standard bivariate case, \(\mathcal{N}_2(\zerov,\Iv_2)\).

The normalizing constant is less important today than the geometry of \((\xv-\muv)^\top\boldsymbol{\Sigma}^{-1}(\xv-\muv)\).

Standard multivariate normal: facts

Let \(\zv \sim \mathcal{N}_p(\zerov,\Iv_p)\).

Property	Consequence
Components	\(Z_1,\ldots,Z_p \iid \mathcal{N}(0,1)\)
Spherical symmetry	\(\boldsymbol{Q}\zv \stackrel{d}{=} \zv\) for every orthogonal \(\boldsymbol{Q}\)
Linear transformations	\(\boldsymbol{A}\zv+\bv \sim \mathcal{N}_m(\bv,\boldsymbol{A}\boldsymbol{A}^\top)\)
Orthogonal projections	If \(\boldsymbol{P}=\boldsymbol{P}^\top=\boldsymbol{P}^2\), then \(\boldsymbol{P}\zv \indep (\Iv-\boldsymbol{P})\zv\)
Quadratic forms	\(\zv^\top\boldsymbol{P}\zv \sim \chi^2_{r}\) when \(\mathrm{rank}(\boldsymbol{P})=r\)
Centered general normal	If \(\Xv\sim\mathcal{N}_p(\muv,\boldsymbol{\Sigma})\), then \((\Xv-\muv)^\top\boldsymbol{\Sigma}^{-1}(\Xv-\muv)\sim\chi^2_p\)

With \(\zv=\epsilonv/\sigma\) and \(\boldsymbol{P}=\Hv\) or \(\Iv-\Hv\), these facts produce the \(\chi^2\), independence, \(t\), and \(F\) results used in regression.

[78:30, 2 min] This is a reference slide to talk through, not a derivation. For linear transformations, dimensions matter: \(\boldsymbol{A}\) is \(m\times p\). For quadratic forms, symmetric idempotent means orthogonal projection, and rank is the number of independent standard normals being squared.

The last line is the squared Mahalanobis distance. The regression connection at the bottom previews why fitted values and residuals are independent under normality, and why residual sums of squares are chi-squared.

Properties of \(\hat\betav\)

Under A1 to A3,

\[ \E[\hat\betav] = \betav, \qquad \mathrm{Var}(\hat\betav) = \sigma^2\pa{\Xv^\top\!\Xv}^{-1} \]

Least squares is unbiased under A1 to A3 alone, with no normality needed. Without proof: it’s also the best such choice. Among all linear unbiased estimators of \(\cv^\top\betav\), \(\cv^\top\hat\betav\) has the smallest variance, for every \(\cv \in \R^p\).

For \(\sigma^2\): a short calculation with \(\Iv - \Hv\) gives

\[ \E\pc{\norm{\ev}^2} = \sigma^2 (n-p) \qquad\Longrightarrow\qquad \hat\sigma^2 = \frac{\ssres}{n - p} \]

[76:30, 2 min] \(\mathrm{Var}(\hat\betav) = \sigma^2(\Xv^\top\Xv)^{-1}\) has two factors: \(\sigma^2\), reduced by blocking and better measurement, and \((\Xv^\top\Xv)^{-1}\), reduced by choosing the design. Orthogonal designs make \(\Xv^\top\Xv\) diagonal, which is one reason balance is used. That is most of the second half of the course.

The optimality claim is the classical result usually named after Gauss and Markov. The name isn’t worth board time, and neither is the proof (a short quadratic-form argument). What’s worth keeping is that least squares needs no normality to be the best linear choice; A4 buys exact \(t\)/\(F\) distributions on top of that, next slide.

The collinearity demonstration in was \((\Xv^\top\Xv)^{-1}\) becoming large.

\(n-p\) in \(\hat\sigma^2\) is the rank of \(\Iv - \Hv\) from the previous slide.

Distribution theory under normality

Adding A4, \(\epsilonv \sim \mathcal{N}_n(\zerov, \sigma^2\Iv_n)\):

\[ \hat\betav \ \sim \ \mathcal{N}_p\!\pa{\betav, \ \sigma^2(\Xv^\top\!\Xv)^{-1}}, \qquad \frac{\ssres}{\sigma^2} \ \sim \ \chi^2_{n-p}, \qquad \hat\betav \indep \ssres \]

Independence holds because \(\hat\yv = \Hv\yv\) and \(\ev = (\Iv-\Hv)\yv\) are jointly normal with

\[ \mathrm{Cov}\pa{\Hv\yv, (\Iv - \Hv)\yv} \ = \ \sigma^2\Hv(\Iv - \Hv) \ = \ \boldsymbol{0} \]

Combining the normal and the \(\chi^2\) gives, for any \(\cv\),

\[ \frac{\cv^\top\hat\betav - \cv^\top\betav} {\hat\sigma\sqrt{\cv^\top(\Xv^\top\!\Xv)^{-1}\cv}} \ \sim \ t_{n-p} \]

[78:30, 2 min] The covariance calculation is one line: \(\Hv(\Iv - \Hv) = \Hv - \Hv^2 = \boldsymbol{0}\). Put it on the board, it is short.

Two things make it work: joint normality, so that zero covariance implies independence, and \(\mathrm{Var}(\epsilonv) = \sigma^2\Iv\). If the errors are correlated, as in a split plot design, this argument fails and the usual \(t\) and \(F\) statistics are wrong.

The \(t\) statistic at the bottom is the general version of and , and the one reported for every coefficient in software output.

§8. ANOVA and general linear hypotheses

Three regression questions, three tests

Question	Null hypothesis	Statistic	What rejection establishes
Does one coefficient or contrast matter?	\(H_0:\cv^\top\betav=d\)	\(t_{n-p}\)	This particular linear combination differs from \(d\)
Does the whole regression improve on the mean alone?	\(H_0:\beta_1=\cdots=\beta_{p-1}=0\)	overall \(F_{p-1,\,n-p}\)	At least one predictor contributes signal
Do these added terms improve a smaller model?	\(H_0:\Cv\betav=\dv\), \(\mathrm{rank}(\Cv)=q\)	partial \(F_{q,\,n-p}\)	The added terms improve fit beyond noise

The first test is the \(q=1\) case of the third, and then \(F=t^2\). A categorical factor with several indicator columns needs a joint \(F\)-test; its individual \(t\)-tests do not answer whether the factor matters as a whole.

[80:30, 1.5 min] These are the three tests students see in regression software: the coefficient table, the overall model \(F\) near \(R^2\), and an ANOVA comparison between nested models. They are not three unrelated inventions. All are restrictions on \(\betav\), and the final two slides of the section derive the general form.

Emphasize the factor row verbally. Testing three dummy coefficients one at a time is not a replacement for the single three-degree-of-freedom question, “does this factor matter?”

Overall regression ANOVA

The intercept-only model predicts \(\bar y\). The full model adds \(p-1\) slopes:

\[ H_0:\beta_1=\cdots=\beta_{p-1}=0 \qquad\text{versus}\qquad H_a:\text{at least one }\beta_j\neq 0 \]

Source	df	SS	MS
Regression	\(p-1\)	\(\mathrm{SS_{reg}}=\mathrm{SS_{tot}}-\ssres\)	\(\mathrm{SS_{reg}}/(p-1)\)
Residual	\(n-p\)	\(\ssres\)	\(\ssres/(n-p)\)
Total	\(n-1\)	\(\mathrm{SS_{tot}}\)

\[ F\ =\ \frac{\mathrm{MS_{reg}}}{\mathrm{MS_{res}}} \ =\ \frac{\mathrm{SS_{reg}}/(p-1)}{\ssres/(n-p)} \ \sim\ F_{p-1,\,n-p}\quad\text{under }H_0 \]

[82:00, 2 min] ANOVA is comparing a full model to a deliberately weak baseline: predicting every response by the grand mean. The numerator is signal per slope degree of freedom; the denominator is noise per residual degree of freedom. “Large” means large relative to an \(F\) reference distribution, not large by eye.

The overall test is omnibus: rejection says at least one slope is nonzero, but does not say which one. In simple regression there is one slope, so this \(F\) is exactly the square of its \(t\) statistic.

One-way ANOVA: a refresher

Application: plant growth

Thirty plants receive one of three conditions: control, treatment 1, or treatment 2. The response is dried plant weight.

One-way ANOVA asks whether the population means are equal:

\[ H_0:\mu_{\mathrm{ctrl}}=\mu_{\mathrm{trt1}}=\mu_{\mathrm{trt2}} \]

It compares the separation among the group means with the variation of plants inside each group.

PlantGrowth experiment; 10 plants per condition

[84:00, 1.5 min] This is the PlantGrowth example distributed with R: one categorical factor, treatment condition, and one quantitative response, dried weight. The boxplots suggest treatment 2 may produce larger plants, but visual separation alone does not quantify uncertainty.

One-way ANOVA compares variation among the three sample means with variation among plants receiving the same condition. The null concerns population means, not the medians drawn inside the boxes. Independence comes from assigning a separate experimental unit to one treatment condition.

One-way ANOVA is a linear model

\(g\) treatments, \(n_i\) units on treatment \(i\), \(N = \sum_i n_i\).

\[ y_{ij} \ = \ \mu + \tau_i + \e_{ij}, \qquad j = 1, \ldots, n_i \]

This is \(\yv = \Xv\betav + \epsilonv\) with \(\betav = (\mu, \tau_1, \ldots, \tau_g)^\top\) and \(\Xv\) an \(N \times (g+1)\) matrix of indicators.

\[ \Xv \ = \ \begin{pmatrix} \mathbf{1}_{n_1} & \mathbf{1}_{n_1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{1}_{n_2} & \mathbf{0} & \mathbf{1}_{n_2} & \cdots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{1}_{n_g} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{1}_{n_g} \end{pmatrix} \]

One column for the intercept, then one indicator column per treatment.

For \(g = 3\) with \(n_1 = n_2 = n_3 = 2\), writing \(\yv = \Xv\betav + \epsilonv\) out in full:

\[ \begin{pmatrix} y_{11} \\ y_{12} \\ y_{21} \\ y_{22} \\ y_{31} \\ y_{32} \end{pmatrix} \ = \ \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \mu \\ \tau_1 \\ \tau_2 \\ \tau_3 \end{pmatrix} \ + \ \begin{pmatrix} \e_{11} \\ \e_{12} \\ \e_{21} \\ \e_{22} \\ \e_{31} \\ \e_{32} \end{pmatrix} \]

[80:30, 1 min] Read the block form first: column one is the intercept, and each remaining column indicates membership of one treatment. Row \(ij\) picks out \(\mu + \tau_i\).

Then step to the concrete \(6 \times 4\) case and let them look at it. Ask what the first column equals. It is the sum of the other three, which is the point of the next slide, so do not resolve it here.

With \(g = 2\) this is exactly the two-sample setup from §4: \(\tau_1 - \tau_2\) is the \(\mu_A - \mu_B\) we tested by hand.

The one-way design matrix is rank deficient

\(\Xv\) has rank \(g\), not \(g+1\): the intercept column is the sum of the \(g\) treatment columns.

Equivalently, adding \(c\) to \(\mu\) and subtracting \(c\) from every \(\tau_i\) leaves \(\Xv\betav\) unchanged. So \(\Xv^\top\!\Xv\) is singular and \(\hat\betav\) is not unique.

The usual fixes are the constraints from §6: set \(\tau_1 = 0\), or \(\sum_i \tau_i = 0\), or drop \(\mu\) and use cell means.

[81:30, 1 min] Answer the question left hanging on the previous slide, then give the numerical version: \(\mu = 10\), \(\tauv = (2,-1,-1)\) and \(\mu = 12\), \(\tauv = (0,-3,-3)\) give identical fitted values.

The fitted values \(\hat\yv = \Hv\yv\) are unique regardless, since they depend on \(\mathcal{C}(\Xv)\) and not on the parameterization. So \(\ssres\), \(\hat\sigma^2\) and the \(F\) tests are all unaffected. Only the individual parameters are ambiguous, which is what the next slide is about.

The ANOVA table

\[ \R^N \ = \ \underbrace{\mathrm{span}(\mathbf{1})}_{1} \ \oplus \ \underbrace{\mathcal{C}(\Xv) \ominus \mathrm{span}(\mathbf{1})}_{g-1} \ \oplus \ \underbrace{\mathcal{C}(\Xv)^\perp}_{N-g} \]

Source	df	SS	\(\E[\mathrm{MS}]\)
Treatments	\(g-1\)	\(\sum_i n_i(\bar y_{i\cdot} - \bar y_{\cdot\cdot})^2\)	\(\sigma^2 + \dfrac{\sum_i n_i \tau_i^2}{g-1}\)
Residual	\(N-g\)	\(\sum_{i,j}(y_{ij} - \bar y_{i\cdot})^2\)	\(\sigma^2\)
Total	\(N-1\)	\(\sum_{i,j}(y_{ij} - \bar y_{\cdot\cdot})^2\)

Under \(H_0 : \tau_1 = \cdots = \tau_g\), the two mean squares have the same expectation and

\[ F \ = \ \frac{\mathrm{MS_{treatments}}}{\mathrm{MS_{residual}}} \ \sim \ F_{g-1,\,N-g} \]

[84:30, 2 min] The decomposition on the top line gives the df column directly as dimensions, and the SS column as squared lengths of projections, the same \(\norm{\yv}^2\) split as §7, with one more piece.

The two sums of squares are independent by the same covariance argument as §7, applied to orthogonal subspaces. That is what allows the ratio to be \(F\).

The expected mean square column explains why the ratio is built this way. In more complicated designs the expected mean squares change, and they are what tell you which mean square belongs in the denominator.

ANOVA as a right triangle

Between-group scale is the separation among treatment means; within-group scale is the spread of observations around their own treatment mean.

[86:30, 2 min] Read the triangle as the vector decomposition \(\yv - \bar y\mathbf{1} = (\hat\yv - \bar y\mathbf{1}) + \ev\). The horizontal leg is the treatment component and the vertical leg is the residual component. They meet at a right angle because fitted values and residuals are orthogonal.

The between-group slider moves the treatment means farther from or closer to the grand mean. The within-group slider changes how dispersed observations are around their own treatment mean without moving those means.

The lengths are square roots of sums of squares. Pythagoras therefore gives \(\mathrm{SS_{total}} = \mathrm{SS_{treatments}} + \ssres\). Move the two sliders separately: one changes only the treatment leg and the other changes only the residual leg, while the hypotenuse follows both.

General linear hypotheses

Every test in this course has the form

\[ H_0 : \Cv\betav = \dv, \qquad \Cv \in \R^{q \times p} \text{ of rank } q, \ \text{rows estimable} \]

Hypothesis	\(\Cv\)
\(\beta_j = 0\)	one row, \(\cv^\top = \ev_j^\top\)
All treatments equal	\(g-1\) contrast rows
No interaction	the interaction rows
Dose response is linear	the curvature contrasts

Fitting under \(H_0\) means projecting onto the smaller subspace \(\mathcal{C}(\Xv_0) \subset \mathcal{C}(\Xv)\) of dimension \(p - q\), with hat matrix \(\Hv_0\) and residual sum of squares \(\ssres_0\).

[86:30, 1.5 min] Filling in the anatomy one more time: the assumption is A1 to A4, the parameter is \(\betav\), the null is \(\Cv\betav = \dv\), the test statistic is on the next slide, and the rejection region comes from an \(F\) distribution. Both worked examples in §4 and the two-sample test in are \(q=1\) special cases of this row.

A constraint on \(\betav\) corresponds to a smaller column space.

Examples of \(\Cv\betav=\dv\)

For example 1, \(\betav=(\beta_0,\beta_1,\beta_2)^\top\). For examples 2 and 3, \(\betav=(\mu_A,\mu_B,\mu_C)^\top\) contains three treatment means.

Scientific question	\(\Cv\)	\(\dv\)	Equivalent null hypothesis
1. Does \(x_2\) contribute after accounting for \(x_1\)?	\(\begin{pmatrix}0&0&1\end{pmatrix}\)	\(0\)	\(H_0:\beta_2=0\)
2. Do two active treatments average to the control?	\(\begin{pmatrix}-1&\tfrac12&\tfrac12\end{pmatrix}\)	\(0\)	\(H_0:\dfrac{\mu_B+\mu_C}{2}=\mu_A\)
3. Are all three treatment means equal?	\(\begin{pmatrix}1&-1&0\\0&1&-1\end{pmatrix}\)	\(\begin{pmatrix}0\\0\end{pmatrix}\)	\(H_0:\mu_A=\mu_B=\mu_C\)

Each row of \(\Cv\) is one linear restriction. Thus \(q=1\) in examples 1 and 2, while \(q=2\) in example 3.

[88:00, 2 min] Read each example from right to left: start with the scientific question, write it as one or more equations, and then read the coefficients of those equations into the rows of \(\Cv\). The first two examples have \(q=1\); the third has \(q=2\) and therefore asks a genuinely joint question.

In most null hypotheses \(\dv=\zerov\), but it need not be. For example, replacing \(d=0\) by a prespecified clinically meaningful difference would test against that target rather than equality.

Scheffé’s method

For \(H_0:\Cv\betav=\dv\), with \(\mathrm{rank}(\Cv)=q\), use a full-rank parameterization and let

\[ \boldsymbol{V}_C=\Cv(\Xv^\top\Xv)^{-1}\Cv^\top. \]

Then

\[ \Cv\hat\betav \sim \mathcal{N}_q\!\left(\Cv\betav,\sigma^2\boldsymbol{V}_C\right), \qquad \frac{ (\Cv\hat\betav-\dv)^\top \boldsymbol{V}_C^{-1} (\Cv\hat\betav-\dv) }{ q\hat\sigma^2 } \sim F_{q,n-p} \quad\text{under }H_0. \]

Scheffé’s extension. In a one-way ANOVA, the same overall \(F\) distribution simultaneously protects every contrast \(\sum_i a_i\mu_i\), even one chosen after seeing the data. Reject \(H_0:\sum_i a_i\mu_i=0\) when

\[ \frac{ \left(\sum_i a_i\bar y_{i\cdot}\right)^2 }{ \hat\sigma^2\sum_i a_i^2/n_i } > (g-1)F_{g-1,N-g;\,1-\alpha}, \qquad \sum_i a_i=0. \]

[90:00, 2 min] The first display is the sampling distribution for a general linear hypothesis. Scheffé is not a different null distribution; it is a simultaneous-inference rule that uses the larger, overall-ANOVA critical value for every possible contrast.

Compare it with a single preplanned contrast, whose squared \(t\) statistic would be compared with \(F_{1,N-g}\). Scheffé instead compares every such statistic with \((g-1)F_{g-1,N-g}\), paying for the freedom to choose any contrast, including after looking at the data.

For next time

• Be able to write \(\Xv\) for a one-way layout in all three parameterizations of §6, by hand.

• Redo the §4 two-sample \(t\)-test with unequal \(n_A, n_B\) of your choosing, and check that \(S_p^2\) still makes sense as a weighted average.

• Review estimability and contrasts (Christensen 2011; Searle and Gruber 2016).

• Bring the experiment you would like to run.

Two things to carry forward:

\[ \hat\yv = \Hv\yv \qquad\text{and}\qquad \mathrm{Var}(\hat\betav) = \sigma^2\pa{\Xv^\top\!\Xv}^{-1} \]

[90:00, 1 min] Both expressions depend only on \(\Xv\). Next week we start choosing \(\Xv\).

References

Christensen, Ronald. 2011. Plane Answers to Complex Questions: The Theory of Linear Models. 4th ed. Springer. https://doi.org/10.1007/978-1-4419-9816-3.

Searle, Shayle R., and Marvin H. J. Gruber. 2016. Linear Models. 2nd ed. John Wiley & Sons.