STAT 6160 • Experimental Design • Fall 2026
\[ \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. Notation
§2. The completely randomized design
§3. Randomization
§4. One-way ANOVA
§5. Multiple comparisons
§6. Contrasts
\[y_{ij} = \text{observation of unit $j$ on treatment $i$}, \\ i = 1,\ldots,g, \quad j = 1,\ldots,n_i\]
\(g\) treatments
placebo, treatment 1 and treatment 2)\(n_i\) units on each treatment \(i\)
placebo, \(n_2\) units on treatment 1, and \(n_3\) units on treatment 2)The total number of experimental units is \(\displaystyle N = \sum_{i=1}^{g} n_i\)
Balanced means when:
A dot replaces an index that has been averaged over.
\[ \bar y_{i\b} \ = \ \frac{1}{n_i}\sum_{j=1}^{n_i} y_{ij}, \qquad \bar y_{\b j} \ = \ \frac{1}{g}\sum_{i=1}^{g} y_{ij}, \qquad \bar y_{\b\b} \ = \ \frac{1}{N}\sum_{i=1}^{g}\sum_{j=1}^{n_i} y_{ij} \]
Two treatments, \(n_1\) and \(n_2\) units.
\[ \begin{aligned} y_{ij} &= \mu + \tau_i + \e_{ij}, &\qquad \e_{ij} &\iid \mathcal{N}(0,\sigma^2) \\ &i = 1,2, & j &= 1,\ldots,n_i \end{aligned} \]
\[ t \ = \ \frac{\hat\Delta - \Delta}{s_p\sqrt{1/n_1 + 1/n_2}} \ \sim \ t_{(n_1+n_2-2)} \]
\[ y_i \ = \ \beta_0 + \beta_1 w_i + \e_i, \qquad w_i = \begin{cases} 0 & \text{treatment 1} \\ 1 & \text{treatment 2}\end{cases} \]
\[ \Xv \ = \ \pc{\ \onev \ \ \ \wv \ }, \qquad \hat\beta_1 = \bar y_{2\b} - \bar y_{1\b}, \qquad \pa{(\Xv^\top\Xv)^{-1}}_{22} = \frac{1}{n_1} + \frac{1}{n_2} \]
The two treatment means are
\[ \mu_1 = \mu + \tau_1, \qquad \mu_2 = \mu + \tau_2. \]
For any number \(\eta\), shifting the parameters gives exactly the same two means:
\[ \mu + \tau_i = (\mu - \eta) + (\tau_i + \eta), \qquad i = 1,2. \]
Thus \(\tau_1-\tau_2\) is unchanged, but neither \(\tau_1\) nor \(\tau_2\) has a unique value. The treatment-effects design matrix is
\[ \underbrace{ \begin{pmatrix} \onev_{n_1} & \onev_{n_1} & \zerov_{n_1} \\ \onev_{n_2} & \zerov_{n_2} & \onev_{n_2} \end{pmatrix} }_{\Xv_{\tau}\ \in\ \R^{N\times3},\quad \rank(\Xv_{\tau})=2} \qquad \begin{pmatrix}\mu\\\tau_1\\\tau_2\end{pmatrix}. \]
The shorthand \(\Xv_{\tau} = [\onev\ \zv_1\ \zv_2]\) is convenient. Since
\[ \qquad \onev = \zv_1 + \zv_2, \]
the first column is already the sum of the other two.
With the usual reference-cell dummy encoding, instead drop the second column (\(\zv_1\)):
\[ \Xv_{\mathrm{ref}} = \pc{\ \onev\ \ \zv_2\ }, \qquad \beta_0 = \mu_1, \qquad \beta_1 = \mu_2 - \mu_1. \]
HW: Why?
\[ (\text{A1 }) \e_i \sim \mathcal{N}(0,\sigma^2) \qquad (\text{A2 }) \var(\e_i) \text{ free of } w_i \qquad (\text{A3 }) \e_i \text{ independent} \\ (\text{A4 }) w_i \text{ is unrelated to omitted } u_i \]
\[ \E(y_i) = \beta_0 + \beta_1 w_i + \beta_2 u_i. \]
Leaving out \(u_i\) is harmful when it is systematically associated with treatment. * Complete randomization makes that association average to zero over the randomization; it does not make \(u_i\) disappear from an individual experiment
Two thumbnail designs, A and B, on one account. Ten posting slots \(i = 1,\ldots,10\), one per hour, five slots per design. Response \(y_i\) is engagement in the hour after the post.
S1. Blocked
S2. Alternating
S3. Randomized
\[ \hat\Delta - \Delta \ = \ \frac{1}{5}\sum_{i=1}^{10} v_i\,\e_i \]
\[ v_i \ = \ \begin{cases} +1 & \text{slot } i \text{ gets A}\\ -1 & \text{slot } i \text{ gets B} \end{cases} \]
S1. Blocked
\[ \frac{1}{5}\left(\sum_{i=1}^{5}\e_i - \sum_{i=6}^{10}\e_i\right) \]
S2. Alternating
\[ \frac{1}{5}\sum_{i=1}^{10}(-1)^{i+1}\e_i \]
S3. Randomized
\[ \frac{1}{5}\sum_{i=1}^{10}v_i\e_i \]
Engagement drifts through the day: \(u_i = c\,(i - 5.5)\), absorbed into \(\e_i\).
\[ \E\big(\hat\Delta\big) - \Delta \ = \ \frac{c}{5}\sum_{i=1}^{10} v_i\,(i - 5.5) \]
\[ \sum_i v_i(i-5.5) = -25 \] \[ \text{bias} = -5c \]
\[ \sum_i v_i(i-5.5) = -5 \] \[ \text{bias} = -c \]
\[ \sum_i v_i(i-5.5) = 3 \] \[ \text{bias} = 0.6c \]
Suppose engagement alternates for an unrecorded reason: \(u_i = c\,(-1)^i\). Then
\[ \E\big(\hat\Delta\big) - \Delta \ = \ \frac{c}{5}\sum_{i=1}^{10} v_i(-1)^i \]
| Schedule | \(\sum_i v_i(-1)^i\) | Bias |
|---|---|---|
| S1: blocked | \(0\) | \(0\) |
| S2: alternating | \(-10\) | \(-2c\) |
| S3: randomized | \(-2\) | \(-0.4c\) |
Randomization is not magic for one realized allocation. Its protection is that, over the draw, the bias is centered at zero for any fixed pattern \(u_i\).
\[ \mathrm{corr}(y_i, y_{i'}) \ = \ \begin{cases} 1 & i = i'\\ \rho & |i - i'| = 1\\ 0 & \text{else} \end{cases} \qquad\Longrightarrow\qquad \var\big(\hat\Delta\big) \ = \ \frac{\sigma^2}{25}\pa{10 + 2\rho\sum_{i=1}^{9} v_i v_{i+1}} \]
\[ \textstyle\sum_i v_i v_{i+1} = 7 \] \[ \tfrac{2\sigma^2}{5} + \tfrac{14\rho\sigma^2}{25} \]
\[ \textstyle\sum_i v_i v_{i+1} = -9 \] \[ \tfrac{2\sigma^2}{5} - \tfrac{18\rho\sigma^2}{25} \]
\[ \textstyle\sum_i v_i v_{i+1} = -1 \] \[ \tfrac{2\sigma^2}{5} - \tfrac{2\rho\sigma^2}{25} \]
The reported standard error uses \(\rho = 0\) in all three cases.
\[ \begin{aligned} y_{ij} &= \mu + \tau_i + \e_{ij}, &\qquad \e_{ij} &\iid \mathcal{N}(0,\sigma^2) \\ &i = 1,\ldots,g, & j &= 1,\ldots,n_i \end{aligned} \]
\[ y_{ij} = \mu_i + \e_{ij}, \qquad \mu_i = \mu + \tau_i \]
Let \(\zv_i\) indicate membership in treatment \(i\). The treatment-effects design matrix is
\[ \underbrace{ \begin{pmatrix} \onev_{n_1} & \onev_{n_1} & \zerov_{n_1} & \cdots & \zerov_{n_1} \\ \onev_{n_2} & \zerov_{n_2} & \onev_{n_2} & \cdots & \zerov_{n_2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \onev_{n_g} & \zerov_{n_g} & \zerov_{n_g} & \cdots & \onev_{n_g} \end{pmatrix} }_{\Xv_{\tau}\ \in\ \R^{N\times(g+1)},\quad \rank(\Xv_{\tau})=g} \qquad \begin{pmatrix}\mu\\\tau_1\\\tau_2\\\vdots\\\tau_g\end{pmatrix}. \]
The shorthand \(\Xv_{\tau} = [\onev\ \zv_1\ \cdots\ \zv_g]\) is convenient. Since \(\onev = \zv_1 + \zv_2 + \cdots + \zv_g\), it has \(g+1\) columns but rank only \(g\): one column is redundant. Equivalently, for every \(\eta\),
\[ \mu + \tau_i = (\mu - \eta) + (\tau_i + \eta), \qquad i = 1,\ldots,g. \]
The fitted group means do not change. A constraint simply chooses one coordinate system for the same fitted values:
\[ \sum_i n_i\tau_i=0 \quad \text{(deviations from the weighted grand mean)}, \qquad \tau_1=0 \quad \text{(treatment 1 as reference)}, \]
or use the \(g\) cell means directly and omit \(\mu\).
With reference-cell dummy encoding, drop the second column (\(\zv_1\)):
\[ \Xv_{\mathrm{ref}} = \pc{\ \onev\ \ \zv_2\ \ \cdots\ \ \zv_g\ }, \qquad \gamma_0 = \mu_1, \qquad \gamma_{i-1} = \mu_i - \mu_1, \quad i = 2,\ldots,g. \]
The intercept is the treatment-1 mean; every other coefficient is that treatment’s difference from treatment 1. The fitted group means are unchanged.
Does restricting food intake extend the lives of mice?
349 female mice were randomly assigned to six diets; the response is lifetime in months. (Weindruch et al. 1986; Ramsey and Schafer 2013)
| Diet | What the mice received | \(n_i\) |
|---|---|---|
| A: (N/N85) | Normal before weaning; 85 kcal/week afterward | 57 |
| B: (NP) | Unlimited standard, nonpurified food | 49 |
| C: (N/R50) | Normal before weaning; 50 kcal/week afterward | 71 |
| D: (R/R50) | Restriction before and after weaning; 50 kcal/week afterward | 56 |
| E: (lopro) | As C, with protein reduced as mice aged | 56 |
| F: (N/R40) | Normal before weaning; 40 kcal/week afterward | 60 |
We use A as the reference for comparisons of restricted regimens. B is a separate standard-diet benchmark.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
d = pd.read_csv("data/longevity.csv")
order = ["A", "B", "C", "D", "E", "F"]
colors = ["gray", "#0072B2", "#E69F00",
"#009E73", "#CC79A7", "#D55E00"]
rng = np.random.default_rng(6160)
fig, ax = plt.subplots()
for i, (diet, color) in enumerate(zip(order, colors)):
y = d.loc[d.Diet == diet, "Lifetime"]
x = i + rng.uniform(-0.15, 0.15, len(y))
ax.scatter(x, y, color=color, alpha=0.65)
ax.set(xticks=range(6), xticklabels=order,
ylabel="Lifetime (months)")
fig.tight_layout()Small horizontal offsets separate overlapping points within each diet.
groups = [
d.loc[d.Diet == diet, "Lifetime"]
for diet in order
]
fig, ax = plt.subplots()
boxes = ax.boxplot(
groups, tick_labels=order,
patch_artist=True
)
for box, color in zip(boxes["boxes"], colors):
box.set_facecolor(color)
box.set_alpha(0.65)
ax.set_ylabel("Lifetime (months)")
fig.tight_layout()Which differences look large relative to the spread within a diet?
For mouse \(j\) receiving diet \(i\), \[ y_{ij}=\mu_i+\e_{ij},\qquad \e_{ij}\iid\mathcal N(0,\sigma^2), \quad i=1,\ldots,6,\quad j=1,\ldots,n_i. \]
\[ H_0:\mu_1=\cdots=\mu_6 \qquad\text{versus}\qquad H_A:\text{at least one mean differs}. \] This tests equality of all six means, including the B group.
Order the groups as A, B, C, D, E, F. With A as reference, the design matrix in block form is
\[ \Xv = \underbrace{ \begin{pmatrix} \onev_{57} & \zerov_{57} & \zerov_{57} & \zerov_{57} & \zerov_{57} & \zerov_{57} \\ \onev_{49} & \onev_{49} & \zerov_{49} & \zerov_{49} & \zerov_{49} & \zerov_{49} \\ \onev_{71} & \zerov_{71} & \onev_{71} & \zerov_{71} & \zerov_{71} & \zerov_{71} \\ \onev_{56} & \zerov_{56} & \zerov_{56} & \onev_{56} & \zerov_{56} & \zerov_{56} \\ \onev_{56} & \zerov_{56} & \zerov_{56} & \zerov_{56} & \onev_{56} & \zerov_{56} \\ \onev_{60} & \zerov_{60} & \zerov_{60} & \zerov_{60} & \zerov_{60} & \onev_{60} \end{pmatrix} }_{349\times6,\quad \rank(\Xv)=6}. \]
Here \(\onev_n\) and \(\zerov_n\) are \(n\times1\) vectors of ones and zeros.
The intercept is \(\mu_{\mathrm{A}}\); the other coefficients are \(\mu_i-\mu_{\mathrm{A}}\). The fitted value for every mouse is its diet’s sample mean.
\[ \underbrace{y_{ij}-\bar y_{\b\b}}_{\text{total deviation}} = \underbrace{\bar y_{i\b}-\bar y_{\b\b}}_{\text{between diets}} + \underbrace{y_{ij}-\bar y_{i\b}}_{\text{within diet}}. \]
\[ \underbrace{\sum_{i,j}(y_{ij}-\bar y_{\b\b})^2}_{\sstot} = \underbrace{\sum_i n_i(\bar y_{i\b}-\bar y_{\b\b})^2}_{\sstrt} + \underbrace{\sum_{i,j}(y_{ij}-\bar y_{i\b})^2}_{\ssres}. \]
Why does the cross term vanish? \[ \sum_{i,j}(\bar y_{i\b}-\bar y_{\b\b})(y_{ij}-\bar y_{i\b}) =\sum_i(\bar y_{i\b}-\bar y_{\b\b}) \underbrace{\sum_j(y_{ij}-\bar y_{i\b})}_{0}=0. \]
\[ \mathbb R^{349} = \underbrace{\mathrm{span}(\onev)}_{1} \oplus \underbrace{\bigl(\mathcal C(\Xv)\cap\onev^\perp\bigr)}_{5} \oplus \underbrace{\mathcal C(\Xv)^\perp}_{343}. \]
The centered response has \(348=5+343\) degrees of freedom. Orthogonal components give the sums-of-squares identity by Pythagoras.
Under \(H_0\) and the normal, common-variance model, \[ \frac{\sstrt}{\sigma^2}\sim\chi^2_{5},\qquad \frac{\ssres}{\sigma^2}\sim\chi^2_{343}, \quad\text{independently}. \]
Therefore both mean squares estimate the same \(\sigma^2\) under \(H_0\): \[ F=\frac{\sstrt/5}{\ssres/343}\sim F_{5,343}. \]
When diet means differ, the expected treatment mean square increases: \[ \E(\mathrm{MS_{trt}})=\sigma^2+ \frac{\sum_i n_i(\mu_i-\mu_w)^2}{g-1}, \qquad \mu_w=\frac{\sum_i n_i\mu_i}{N}. \] Large values of \(F\) provide evidence against equality of the diet means.
| Source | df | Sum of squares | Mean square |
|---|---|---|---|
| Diet | \(g-1=5\) | \(\sstrt\) | \(\sstrt/5\) |
| Residual | \(N-g=343\) | \(\ssres\) | \(\ssres/343\) |
| Total | \(N-1=348\) | \(\sstot\) |
\[ F=\frac{\mathrm{MS_{trt}}}{\mathrm{MS_{res}}}, \qquad \hat\sigma=\sqrt{\mathrm{MS_{res}}}, \qquad R^2=\frac{\sstrt}{\sstot}. \]
Unequal group sizes change the weights, not this decomposition.
import statsmodels.formula.api as smf
# First category becomes the reference.
d["Diet"] = pd.Categorical(
d.Diet, categories=order
)
fit = smf.ols(
"Lifetime ~ C(Diet)", data=d
).fit()
coef = pd.DataFrame({
"Estimate": fit.params,
"SE": fit.bse
})
coef.index = ["A mean"] + [
diet + " - A" for diet in order[1:]
]
coef.round(2)| Estimate | SE | |
|---|---|---|
| A mean | 32.69 | 0.88 |
| B - A | -5.29 | 1.30 |
| C - A | 9.61 | 1.19 |
| D - A | 10.19 | 1.26 |
| E - A | 6.99 | 1.26 |
| F - A | 12.43 | 1.24 |
The intercept is a mean in months; every other estimate is a difference in months. These standard errors pool within-diet variation from all six groups.
| df | sum_sq | mean_sq | F | PR(>F) | |
|---|---|---|---|---|---|
| C(Diet) | 5.0 | 12733.942 | 2546.788 | 57.104 | 0.0 |
| Residual | 343.0 | 15297.415 | 44.599 | NaN | NaN |
\(F_{5,343}=57.10\), \(p=4.11e-43\); \(\hat\sigma=6.68\) months.
We have evidence that the six mean lifetimes are not all equal. Comparisons with A require tests of the individual mean differences.
Does reducing the ration from 85 to 50 kcal/week after weaning change mean lifetime?
Groups C and A both received normal feeding before weaning: \[ \theta=\mu_{\mathrm{C}}-\mu_{\mathrm{A}}. \]
In our six-diet order, the coefficients are \[ \cv=(-1,\ 0,\ 1,\ 0,\ 0,\ 0). \] The unrestricted B group is not part of this particular comparison.
\[ \theta=\sum_i c_i\mu_i,\qquad \sum_i c_i=0. \]
| Question | Coefficients in our diet order |
|---|---|
| C versus A: 50 versus 85 kcal after weaning | \((-1,0,1,0,0,0)\) |
| D versus C: earlier versus later restriction at 50 kcal | \((0,0,-1,1,0,0)\) |
| F versus C: 40 versus 50 kcal after weaning | \((0,0,-1,0,0,1)\) |
| E versus C: reduced protein | \((0,0,-1,0,1,0)\) |
Because the coefficients sum to zero, \(\sum_i c_i\mu_i=\sum_i c_i\tau_i\): the answer is unchanged by the identifiability constraint.
\[ \hat\theta=\sum_i c_i\bar y_{i\b},\qquad \mathrm{se}(\hat\theta)= \sqrt{\mathrm{MS_{res}}\sum_i\frac{c_i^2}{n_i}}. \]
Under \(H_0:\theta=\theta_0\), \[ t=\frac{\hat\theta-\theta_0}{\mathrm{se}(\hat\theta)} \sim t_{343}. \]
A pointwise 95% confidence interval is \[ \hat\theta\pm t_{343}^{0.975}\mathrm{se}(\hat\theta). \]
For C versus A (50 versus 85 kcal), only \(1/71+1/57\) enters the variance multiplier.
| c0 | |
|---|---|
| coef | 9.606 |
| std err | 1.188 |
| t | 8.088 |
| P>|t| | 0.000 |
| Conf. Int. Low | 7.270 |
| Conf. Int. Upp. | 11.942 |
Estimated change: 9.61 months; pointwise 95% CI: (7.27, 11.94) months.
Why does the coefficient vector in the Python call differ from \(\cv\)?
How does the average of groups C–F (the four restricted regimens) compare with A?
\[ \cv=(-1,0,\tfrac14,\tfrac14,\tfrac14,\tfrac14). \]
This gives each regimen equal weight, despite unequal sample sizes. Pooling all restricted mice would give a different, sample-size-weighted estimand.
For a null contrast \(\theta=0\), \[ \mathrm{SS}(\cv)=\frac{\hat\theta^2}{\sum_i c_i^2/n_i}, \qquad F_{\text{one contrast}}=\frac{\mathrm{SS}(\cv)}{\mathrm{MS_{res}}}=t^2. \]
Two estimated contrasts are uncorrelated precisely when \[ \sum_i\frac{c_i c_i'}{n_i}=0. \]
C versus A and F versus A share the same reference mean. Their covariance is \(\sigma^2/57>0\): their sums of squares do not simply add.
For \(H_0:\theta=\sum_i c_i\mu_i=0\), with \(\sum_i c_i=0\), write \[ t_{\cv}=\frac{\hat\theta} {\sqrt{\mathrm{MS_{res}}\sum_i c_i^2/n_i}}. \]
Scheffé’s method rejects \(H_0\) when \[ |t_{\cv}|>\sqrt{(g-1)F_{g-1,N-g}^{1-\alpha}}, \qquad\text{equivalently}\qquad \frac{t_{\cv}^2}{g-1}>F_{g-1,N-g}^{1-\alpha}. \]
This controls the probability of any false rejection across all contrasts at level \(\alpha\), under the normal, common-variance model. It therefore also covers contrasts selected after examining the data.
Simultaneous \(100(1-\alpha)\%\) intervals for all contrasts are \[ \hat\theta\ \pm\ \underbrace{\sqrt{(g-1)F_{g-1,N-g}^{1-\alpha}}}_{\text{Scheffé multiplier}} \sqrt{\mathrm{MS_{res}}\sum_i\frac{c_i^2}{n_i}}. \]
For the six diets, the multiplier is \(\sqrt{5F_{5,343}^{0.95}}\).
These intervals are wider than pointwise \(t\) intervals because they cover all contrasts simultaneously. For comparisons only with A, Dunnett’s method gives narrower intervals.
Compare the average of groups C–F (the restricted diets) with A: \[ \theta=\frac{\mu_{\mathrm{C}}+\mu_{\mathrm{D}} +\mu_{\mathrm{E}}+\mu_{\mathrm{F}}}{4} -\mu_{\mathrm{A}}. \]
# Weights on the reference-coded parameters.
result = fit.t_test([0, 0, .25, .25, .25, .25])
effect = result.effect.item()
se = result.sd.item()
df = fit.df_resid
k = len(order) - 1
critical = np.sqrt(k * stats.f.ppf(.95, k, df))
p_scheffe = stats.f.sf((effect / se)**2 / k, k, df)
pd.Series({
"Estimate (months)": effect,
"SE (months)": se,
"Scheffe lower": effect - critical * se,
"Scheffe upper": effect + critical * se,
"Adjusted p": p_scheffe
}).to_frame("Value")| Value | |
|---|---|
| Estimate (months) | 9.805092e+00 |
| SE (months) | 9.837196e-01 |
| Scheffe lower | 6.512712e+00 |
| Scheffe upper | 1.309747e+01 |
| Adjusted p | 2.122759e-17 |
We test groups B–F against A. Together, these tests form a family of comparisons.
\[ \mathrm{FWER}=\Pr(\text{at least one false rejection in the family}). \]
Rejecting the omnibus null does not control the familywise error rate of unadjusted follow-up tests. Dunnett’s procedure controls this rate without a preliminary omnibus test.
| Comparisons included | Procedure | Scope of simultaneous inference |
|---|---|---|
| Groups B–F versus A | Dunnett | All five reference comparisons |
| All 15 pairwise diet differences | Tukey–Kramer | All pairs, allowing unequal \(n_i\) |
| \(K\) contrasts specified before examining the data | Bonferroni | All \(K\) contrasts |
| Any contrast, including one selected after seeing data | Scheffé | The entire contrast family |
Bonferroni uses \(t_{343}^{1-\alpha/(2K)}\) for two-sided intervals; Scheffé uses \(\sqrt{5F_{5,343}^{1-\alpha}}\).
from scipy import stats
control = groups[0] # A
dunnett = stats.dunnett(
*groups[1:], control=control,
alternative="two-sided",
rng=np.random.default_rng(6160)
)
ci = dunnett.confidence_interval()
comparison = pd.DataFrame({
"Difference": [
y.mean() - control.mean()
for y in groups[1:]
],
"Lower": ci.low,
"Upper": ci.high,
"Adjusted p": dunnett.pvalue
}, index=order[1:])
comparison.round(3)| Difference | Lower | Upper | Adjusted p | |
|---|---|---|---|---|
| B | -5.289 | -8.572 | -2.007 | 0.0 |
| C | 9.606 | 6.609 | 12.602 | 0.0 |
| D | 10.194 | 7.024 | 13.365 | 0.0 |
| E | 6.994 | 3.824 | 10.165 | 0.0 |
| F | 12.425 | 9.309 | 15.542 | 0.0 |
Intervals have joint 95% coverage under the common-variance normal model. Positive differences mean longer average life than A.
Which diets have longer mean lifetimes than A, and by how many months? An interval containing zero is not evidence that the two diets are equivalent.
These estimates apply to the mouse population and feeding conditions studied.