Quadratic Forms and Geometry of Least Squares 1 Basic Geometry

[figure: \underline{y}; Span \underline{1}; Span(X_{1}, X_{2}, …, X_{k}); SpanX = Span(\underline{1}, X_{1}, …, X_{k}); \underline{ŷ} = P \underline{y}; \underline{E} = Q\underline{y} = (I − P)\underline{y}]

Figure 1: Geometry of Linear Regression

We have seen that fitting the linear regression model

\underline{y} = X\underline{β} + \underline{ε}, \underline{ε} ∼ N(\underline{0}, σ^{2}I_{n×n}),

where

X = [ \begin{matrix}| & | & | & & | \\ \underline{1} & \underline{X}_{1} & \underline{X}_{2} & ⋯ & \underline{X}_{k} \\ | & | & | & & | \end{matrix} ] , \underline{β}^{T} = (β_{0}, β_{1}, …, β_{k}),

by least squares was equivalent to decomposing the vector \underline{y} ∈ R^{n} into orthogonal components

\underline{y} = P \underline{y} + Q\underline{y} = ˆ \underline{y} + \underline{E}

where

ˆ\underline{y} = P \underline{y} ∈ ColSpace(P) = ColSpace(X)
\underline{E} = Q\underline{y} ⊥ ColSpace(P)

and

P = X(X^{T} X)^{−1}X^{T}
Q = I_{n×n} − P
ˆ\underline{β} = (X^{T} X)^{−1}X^{T} \underline{y}
ˆ\underline{y} = X ˆ\underline{β}  = P \underline{y}
\underline{E} = \underline{y} − ˆ \underline{y}  = Q\underline{y}

since

P^{2} = P; P^{T} = P  (idempotent, symmetric)
Q^{2} = Q; Q^{T} = Q  (idempotent, symmetric)
P Q = QP = 0  (orthogonal)
P + Q = I  (complete)
and P X = X  (projection onto ColSpace(X))

This gives rise to the picture in Figure 1. It is not hard to show for example that the Pythagorean Theorem holds,

∥\underline{y}∥^{2} = \underline{y}^{T} \underline{y} = [(P +Q)\underline{y}]^{T} [(P +Q)\underline{y}] = (\underline{y}^{T} P +\underline{y}^{T} Q)(P y+Qy) = \underline{y}^{T} P \underline{y}+\underline{y}^{T} Q\underline{y} = ∥ˆ\underline{y}∥^{2}+∥\underline{E}∥^{2}

and it turns out, according to a theorem called “Cochran’s Theorem” which we shall see below, that \underline{y}^{T} P \underline{y} and \underline{y}^{T} Q\underline{y} are statistically independent of one another. Much of the theory of F -statistic testing in linear models is built on elaborations of this Pythagorean Theorem, and on the corresponding statistical indepedence facts about the component sums of squares.

2 Some Results on Quadratic Forms

Let \underline{y} ∈ R^{n}, and suppose A is an n × n matrix.

Definition 1 \underline{y}^{T} A\underline{y} is a quadratic form if

1. A^{T} = A  (symmetric)
2. \underline{y}^{T} A\underline{y} ≥ 0, ∀\underline{y} ∈ R^{n} (positive semidefinite)

In this case, A is the defining matrix of \underline{y}^{T} A\underline{y}.

Definition 2 A is a projection matrix if A^{2} = A (we say A is idempotent).

Note that if A is a projection matrix then A(I_{n×n} − A) = \underline{0}, and every \underline{y} ∈ R^{n} has a unique decomposition

\underline{y} = A\underline{y} + (I − A)\underline{y} = \underline{y}_{A} + \underline{y}_{A^{⊥}}

into a vector \underline{y}_{A} ∈ ColSpace(A) and a vector \underline{y}_{A^{⊥}} ⊥ ColSpace(A). If A is also symmetric, then

trace(A) = rank(A) = dim[ColSpace(A)].

Note that, therefore, trace(A) is the “number of dimensions” that A\underline{y} can wander around in, and \underline{y}^{T} A\underline{y} = (A\underline{y})^{T} (A\underline{y}) = ∥A\underline{y}∥^{2} is the squared length of A\underline{y}.

Facts:

1. If \underline{y} is a random vector with E[\underline{y}] = \underline{μ}, Var(\underline{y}) = V, and if \underline{y}^{T} A\underline{y} is a quadratic form, then

E[\underline{y}^{T} A\underline{y}] = trace AV + \underline{μ}^{T} A\underline{μ}

2. If y_{i} ∼ N(μ_{i}, σ^{2}), independent, and A is symmetric and idempotent, then \underline{y}^{T} A\underline{y}/σ^{2} has a noncentral χ^{2} distribution on k degrees of freedom, where k = rank A. The quantity \underline{μ}^{T} A\underline{μ}/2σ^{2} is called the noncentrality parameter of the distribution. Note: here, V = Var \underline{y} = σ^{2}I.

Theorem 1 (Cochran’s Theorem, Scheffé, 1959, pp 419–421)

Let y_{i} ∼ N(μ_{i}, 1), i = 1 … n, independent, and suppose

∑_{i=1}^{n} y_{i}^{2} = \underline{y}^{T} A_{1}\underline{y} + \underline{y}^{T} A_{2}\underline{y} + … + \underline{y}^{T} A_{s}\underline{y}

where each \underline{y}^{T} A_{j}\underline{y} is a quadratic form, with rank A_{j} = n_{j}. \underline{Then:}

(a) Each A_{j}^{2} = A_{j} (\underline{idempotent}, so rank A_{j} = trace A_{j}); and A_{i}A_{j} = 0 for all i ≠ j (A_{i} and A_{j} are orthogonal).

(b) The quadratic forms \underline{y}^{T} A_{j}\underline{y} are distributed as independent, noncentral χ^{2}’s on n_{j} degrees of freedom.

\underline{if and only if}n_{1} + n_{2} + … + n_{s} = n.

Corollary 1 If y_{i} ∼ N(μ_{i}, σ^{2}), the theorem still holds but now \underline{y}^{T} A_{j}\underline{y}/σ^{2} will be independent noncentral χ_{n_{j}}^{2}’s.

Cochran’s Theorem or its corollary tell us immediately, for example, that since

∑_{i}y_{i}^{2} = ∥\underline{y}∥^{2} = \underline{y}^{T} P \underline{y} + \underline{y}^{T} Q\underline{y}

and since trace(P) + trace(Q) = trace(P + Q) = trace(I) = n, then \underline{y}^{T} P \underline{y} and \underline{y}^{T} Q\underline{y} are (proportional to) independent χ^{2}’s on trace(P) and trace(Q) degrees of freedom.

3 Overall ANOVA for regression with an intercept

Let’s continue now with the model

\underline{y} = X\underline{β} + \underline{ε}, \underline{ε} ∼ N(\underline{0}, σ^{2}I_{n×n}),  (1)

where

X = [ \begin{matrix}| & | & | & & | \\ \underline{1} & \underline{X}_{1} & \underline{X}_{2} & ⋯ & \underline{X}_{k} \\ | & | & | & & | \end{matrix} ] , \underline{β}^{T} = (β_{0}, β_{1}, …, β_{k}), k′ = k + 1,

with projection matrix P and Q = I − P.

We wish to compare this model with the simpler model

\underline{y} = \underline{1}β_{0} + \underline{ε}  (2)

of regression onto just the intercept β_{0}. You can verify that the projection matrix for this regression is J/n where J is the n × n matrix of all 1’s. It is easy to see that

∑_{i=1}^{n} y_{i}^{2} =  ∥\underline{y}∥^{2}  =  \underline{y}^{T} I \underline{y}
= \underline{y}^{T} (I − P)\underline{y} + \underline{y}^{T} (P − J/n)\underline{y} + \underline{y}^{T} J/n \underline{y}
= ∥Q\underline{y}∥^{2} + ∥(P − J/n)\underline{y}∥^{2} + ∥(J/n)\underline{y}∥^{2}
=  RSS  +  RegSS  + n(ȳ)^{2}

and we can readily verify^{1} that the matrices A_{1} = (I − P), A_{2} = (P − J/n), and A_{3} = J/n, are symmetric (A_{i}^{T} = A), idempotent (A_{i}^{2} = A_{i}) and orthogonal (A_{i}A_{j} = 0), and trace(A_{1}) + trace(A_{2}) + trace(A_{3}) = n. So, by Cochran’s Theorem and its corollary, the quadratic forms

\underline{y}^{T} (J/n)\underline{y}/σ^{2}, \underline{y}^{T} (P − J/n)\underline{y}/σ^{2}, \underline{y}^{T} Q\underline{y}/σ^{2}

are independent noncentral χ^{2}’s. Moreover you can easily show that

T SS = ∑_{i=1}^{n}(y_{i} − ȳ)^{2} = ∑ y_{i}^{2} − n(ȳ)^{2} = \underline{y}^{T} (I − J/n)\underline{y}

so the RegSS and RSS terms are in fact independent χ^{2}’s in the overall ANOVA table

Source  df  SS
Regression k′ − 1 RegSS = \underline{y}^{T} (P − J/n)\underline{y}
Residual n − k′ RSS = \underline{y}^{T} Q\underline{y}
Total  n − 1 TSS = \underline{y}^{T} (I − J/n)\underline{y}

The geometry for this ANOVA table is given in Figure 2.

The overall F -statistic for the regression is

F = \frac{(T SS − RSS)([(n − 1) − (n − k′)])}{RSS/(n − k′)} = \frac{RegSS/(k′ − 1)}{RSS/(n − k′)} = \frac{RegSS/k}{RSS/(n − (k + 1))},

and tests H_{0} : eq. (2) holds, vs. H_{1} : equation (1) holds. We note that

1Use the facts that P and J/n are the projection matrices for the regressions (1) and (2) and the fact that \underline{1} is a column of X, so that (I − P)(J/n) = 0.

[figure: \underline{y}; Span \underline{1}; Span(X_{1}, X_{2}, …, X_{k}); Span(\underline{1}, X_{1}, …, X_{k}); \underline{ŷ}; \underline{ȳ}; ∥(P − J/n)\underline{y}∥^{2} = RegSS = T SS − RSS; ∥Q\underline{y}∥^{2} = ∑(y_{i} − ŷ_{i})^{2} = RSS; ∥(I − J/n)\underline{y}∥^{2} = ∑(y_{i} − ȳ)^{2} = T SS]

Figure 2: Geometry of the overall ANOVA table in linear regression

• RegSS/σ^{2} = \underline{y}^{T} (P − J/n)\underline{y}/σ^{2} ∼ noncentral χ_{k−1}^{2} :

df = rank (P − J/n) = trace (P − J/n) = trace P − trace (J/n)
= trace [X(X^{T} X)^{−1}X^{T} ] − 1 = trace (X^{T} X(X^{T} X)^{−1})−1
= trace I_{k′×k′} − 1 = k′ − 1 = k
E [\underline{y}^{T} (P − J/n)\underline{y}] = trace [(P − J/n)σ^{2}I] + (X\underline{β})^{T} (P − J/n)(X\underline{β})
= σ^{2}(k′ − 1) + \underline{β}^{T} X^{T} (P − J/n)X\underline{β}
= σ^{2}(k′ − 1) + \underline{β}^{T} X^{T} (I − J/n)X\underline{β}
since P X = X.

• RSS/σ^{2} = \underline{y}^{T} Q\underline{y} ∼ noncentral χ_{n−k′}^{2} :

df = rank Q = trace Q = trace (I − P) = trace I − trace P
= n − k′ = n − (k + 1)
E [\underline{y}^{T} Q\underline{y}] = trace [Q(σ^{2}I)] + (X\underline{β})^{T} Q(X\underline{β})
= σ^{2}(n − k′) + \underline{β}^{T} X^{T} QX\underline{β}
= σ^{2}(n − k′) = σ^{2}(n − (k + 1))

(Note: \underline{β}^{T} X^{T} QX\underline{β} = 0 since QX = 0.)

After some algebra we see that the ratio of the expected values in the F statistic is

\frac{E [RegSS/(k′ − 1)]}{E [RSS/(n − k′)]} = 1 + [\underline{β}^{T} X^{T} (I − J/n)X\underline{β}] /(k′ − 1)σ^{2} ≥ 1

The thing that prevents this from being exactly 1 is \underline{β}^{T} X^{T} (I − J/n)X\underline{β}. If you multiply out X^{T} (I − J/n)X you will see^{2} that

\underline{β}^{T} X^{T} (I − J/n)X\underline{β} = (β_{0}, β_{1} ⋯ β_{k}) [ \begin{matrix}0 & 0 & … & 0 \\ 0 & stuff & … & stuff \\ ⋮ & ⋮ & & ⋮ \\ 0 & stuff & … & stuff \end{matrix} ] [ β_{0} ; β_{1} ; ⋮ ; β_{k} ]

= (β_{1} ⋯ β_{k}) [ \begin{matrix}stuff & … & stuff \\ ⋮ & & ⋮ \\ stuff & … & stuff \end{matrix} ] [ β_{1} ; ⋮ ; β_{k} ]

which will equal zero when H_{0} : β_{1} = β_{2} = ⋯ = β_{k} = 0 is true.

Distribution of F = \frac{RegSS/(k′−1)}{RSS/(n−k′)} under H_{0}:

According to Cochran’s theorem, F is a ratio of independent noncentral χ^{2}’s divided by their respective degrees of freedom. If they were central χ^{2}’s (noncentrality parameter = 0), F would follow the F distribution we always use.

But they are central χ^{2}’s:

• Noncentrality param. of RSS is \frac{1}{2σ^{2}} (X\underline{β})^{T} Q(X\underline{β}) = \underline{β}^{T} X^{T} QX\underline{β} = 0 since QX = 0

• Noncentrality param. of RegSS is \frac{1}{2σ^{2}} \underline{β}^{T} X^{T} (I − J/n)X\underline{β} = 0 under H_{0} (See above!).

Distribution of F under H_{1}:

The numerator χ^{2}, RegSS/σ^{2}, now has noncentrality parameter

\frac{1}{2σ^{2}} \underline{β}^{T} X^{T} (I − J/n)X\underline{β} ≥ 0

The larger this parameter, the easier H_{1} is to detect.

2Hint: Apply the formula ∑_{i=1}^{n} x_{i}^{2} − n(x̄)^{2} = \underline{x}^{T} (I − J/n)\underline{x} to situations where \underline{x}^{T} and \underline{x} represent a column of X, and recall that the first column of X is all 1’s.

4 Using sums of squares to test nested models

Model: \underline{y} = X\underline{β} + \underline{ε} , \underline{ε} ∼ N(\underline{0}, σ^{2}I) with some small changes of notation: Let k′ = k_{1} + k_{2}, and

X = [ \begin{matrix}| & | & & | \\ \underline{X}_{1} & \underline{X}_{2} & ⋯ & \underline{X}_{k′} \\ | & | & & | \end{matrix} ] = [\underbrace{\begin{matrix}| & | & & | \\ \underline{X}_{1} & \underline{X}_{2} & ⋯ & \underline{X}_{k_{1}} \\ | & | & & | \end{matrix}} \underbrace{\begin{matrix}| & | & & | \\ \underline{X}_{k_{1}+1} & \underline{X}_{k_{1}+2} & ⋯ & \underline{X}_{k_{1}+k_{2}} \\ | & | & & | \end{matrix}}]
= [ X_{1} X_{2} ]

and

\underline{β^{T}} = (β_{1} ⋯ β_{k′})  = (\underbrace{β_{1}, β_{2}, ⋯, β_{k_{1}},} \underbrace{β_{k_{1}+1}, β_{k_{1}+2}, ⋯, β_{k_{1}+k_{2}}})
= ( \underline{β}_{1}^{T} , \underline{β}_{2}^{T} )

so that

\underline{y} = X\underline{β} + \underline{ε} = X_{1}\underline{β}_{1} + X_{2}\underline{β}_{2} + \underline{ε}

(We could treat “regression with an intercept” again by making the column \underline{X}_{1} just a column of 1’s.) We want to compare the two models

SMALL: \underline{y} = X_{1}\underline{β}_{1} + \underline{ε}  Hat matrix P_{1} = X_{1}(X_{1}^{T} X_{1})^{−1}X_{1}^{T}
LARGE: \underline{y} = X_{1}\underline{β}_{1} + X_{2}\underline{β}_{2} + ε Hat matrix P_{12} = X(X^{T} X)^{−1}X^{T}
RSS(\underline{β}_{1}, \underline{β}_{2}) = RSS_{1}2 = ∑_{i=1}^{n}(y_{i} − ŷ_{i(12)})^{2}
= (\underline{y} − \underline{ŷ}_{(12)})^{T} (\underline{y} − \underline{ŷ}_{(12)})
= [(I − P_{12})\underline{y}]^{T} [(I − P_{12})\underline{y}]
= \underline{y}^{T} (I − P_{12}) \underline{y} = \underline{y}^{T} Q_{12}\underline{y}  (Q_{12} = I − P_{12})

Similarly

RSS(\underline{β}_{1}) = RSS_{1} = ∑_{i=1}^{n}(y_{i} − y_{i(1)})^{2} = \underline{y}^{T} Q_{1}\underline{y}  (Q_{1} = I − P_{1})

We want a test for H_{0} : \underline{β}_{2} = \underline{0} vs. H_{1} : \underline{β_{2}} free. [See Figure 3.] This is a test for

\begin{aligned}H_{0} : \underline{y} & = X_{1}\underline{β}_{1} + \underline{ε} \\ vs. H_{1} : \underline{y} & = X_{1}\underline{β}_{1} + X_{2}\underline{β}_{2} + \underline{ε} \end{aligned}

Clearly we can write, thinking of the smaller model,

∑ y_{i}^{2} = ∑ ŷ_{i(1)}^{2} + (∑ y_{i}^{2} − ∑ ŷ_{i(1)}^{2})
= \underline{y}^{T} P_{1}\underline{y} + \underline{y}^{T} I\underline{y} − \underline{y}^{T} P_{1}\underline{y}
= \underline{y}^{T} P_{1}\underline{y} + \underline{y}^{T} (I − P_{1})\underline{y}
= \underline{y}^{T} P_{1}\underline{y} + \underline{y}^{T} Q_{1}\underline{y}
= \underline{y}^{T} P_{1}\underline{y} + \underline{y}^{T} (Q_{1} − Q_{12})\underline{y} + \underline{y}^{T} Q_{12}\underline{y}
= \underline{y}^{T} P_{1}\underline{y} + (RSS_{1} − RSS_{12}) + RSS_{12}
= \underline{y}^{T} P_{1}\underline{y} + SS(\underline{β}_{2}|\underline{β}_{1}) + RSS_{12}

[figure: \underline{y}; SpanX_{1}; SpanX_{2}; SpanX = Span[X_{1}, X_{2}]; \underline{ŷ}_{(12)}; \underline{ŷ}_{(1)}; (length)^{2} = RSS_{1} − RSS_{12} = SS(\underline{β}_{2}|\underline{β}_{1}); (length)^{2} = RSS_{12}; (length)^{2} = RSS_{1}]

Figure 3: Geometry of the F -test comparing the nested models X = [X_{1}, X_{2}] and X_{1}

where the third-to-last inequality follows from Q_{1} = (Q_{1} − Q_{12}) + Q_{12} and the fact that (check!) Q_{12}(Q_{1} − Q_{12}) = 0.

Now, if SS(\underline{β}_{2}|\underline{β}_{1}) = \underline{y}^{T} (Q_{1} − Q_{12})\underline{y}:

(1) is a quadratic form (positive-semidefinite)?

(2) has noncentrality parameter = 0 under H_{0}?

then according to Cochran’s theorem, the F -statistic

F = \frac{SS(\underline{β}_{2}|β_{1})/ trace [Q_{1} − Q_{12}]}{RSS(\underline{β}_{1}, \underline{β}_{2})/ trace Q_{12}}

can be used to test H_{0} : \underline{β}_{2} = 0 vs. H_{1} : (H_{0})^{c}. In order to verify (1) and (2) for SS(\underline{β}_{2}|\underline{β}_{1}) = \underline{y}^{T} (Q_{1} − Q_{12})\underline{y} we will reparametrize the larger models to make Q_{12} easier to calculate. The idea is similar to the idea that ˆ \underline{β}_{2} can be obtained by regressing the residuals Q_{1}\underline{y} onto the residuals Q_{1}X_{2} that we saw empirically earlier in the course.

Define

P_{1}
X_{2·1} = Q_{1}X_{2} = X_{2} − X_{1} \underbrace{(X_{1}^{T} X_{1})^{−1} X_{1}^{T} X_{2}}
A

Clearly X_{2·1} is the residuals of (each column of) X_{2} after regressing on X_{1}.

• X_{2·1} is the part of X_{2} unexplained by regression on X_{1}.

• X_{2·1} is orthogonal to X_{1} : X_{1}^{T} X_{2·1} = 0. With X_{2·1} and A = (X_{1}^{T} X_{1})^{−1}X_{1}^{T} X_{2} we can write

\underline{y} = X_{1}\underline{β}_{1} + X_{2}\underline{β}_{2} + \underline{ε} = X_{1}\underline{β}_{1} + (X_{1}A + X_{2} − X_{1}A)\underline{β}_{2} + \underline{ε}

= X_{1}(\underbrace{\underline{β}_{1} + A\underline{β}_{2}}) + (\underbrace{X_{2} − X_{1}A})\underline{β}_{2} + \underline{ε} = X_{1} \underline{α}_{1} + X_{2·1} \underline{β}_{2} + \underline{ε}

So the larger model is identical to

\underline{y} = X_{1}\underline{α}_{1} + X_{2·1}\underline{β}_{2} + \underline{ε}

with model matrix X̃ = [X_{1} X_{2·1}] instead of X = [X_{1} X_{2}]. Therefore

\underline{y}^{T} Q_{12}\underline{y} = RSS (\underline{β}_{1}, \underline{β}_{2}) = RSS (\underline{α}_{1}, \underline{β}_{2}) = \underline{y}^{T} [I − X̃ (X̃^{T} X̃)^{−1} X̃^{T} ] \underline{y}

= \underline{y}^{T} \{I − [X_{1} X_{2·1}] [ \begin{matrix}X_{1}^{T} X_{1} & X_{1}^{T} X_{2·1} \\ X_{2·1}^{T}X_{1} & X_{2·1}^{T}X_{2·1} \end{matrix} ]^{−1} [ \begin{matrix}X_{1}^{T} \\ X_{2·1}^{T} \end{matrix} ]\} \underline{y}

= \underline{y}^{T} \{I − [X_{1} X_{2·1}] [ \begin{matrix}(X_{1}^{T} X_{1})^{−1} & 0 \\ 0 & (X_{2·1}^{T}X_{2·1})^{−1} \end{matrix} ] [ \begin{matrix}X_{1}^{T} \\ X_{2·1}^{T} \end{matrix} ]\} \underline{y}

= \underline{y}^{T} (I − P_{1} − P_{2·1}) \underline{y}

Hence

\begin{aligned}\underline{y}^{T} Q_{12}\underline{y} & = \underline{y}^{T} (I − P_{1} − P_{2·1}) \underline{y} \\ & = \underline{y}^{T} Q_{1}\underline{y} − \underline{y}^{T} P_{2·1}\underline{y} \end{aligned}

Therefore

SS (\underline{β}_{2}|\underline{β}_{1}) = \underline{y}^{T} (Q_{1} − Q_{12}) \underline{y} = \underline{y}^{T} P_{2·1}\underline{y}

is a quadratic form with defining matrix P_{2·1} = X_{2·1} (X_{2·1}^{T}X_{2·1})^{−1} X_{2·1}^{T}. (It is positive definite because

\underline{y}_{1}^{T} P_{2·1}\underline{y} = ∑ ŷ_{i(2·1)}^{2} ≥ 0

from the model \underline{y} = X_{2·1} ˜\underline{β} + \underline{ε}). Now, since P_{2·1} projects onto X_{2·1} which is orthogonal to X_{1}, it follows that P_{2·1}X_{1} = 0. Therefore,

E[SS(\underline{β}_{2}|\underline{β}_{1})] = E [\underline{y}^{T} P_{2·1}\underline{y}]
= trace (P_{2·1}σ^{2}I) + ([X_{1} X_{2}] \underline{β})^{T} P_{2·1} ([X_{1} X_{2}] \underline{β})
= σ^{2} trace [X_{2·1}^{T}X_{2·1} (X_{2·1}^{T}X_{2·1})^{−1}] + [\underline{β}_{1}^{T} \underline{β}_{2}^{T} ] [ \begin{matrix}X_{1}^{T} \\ X_{2}^{T} \end{matrix} ] P_{2·1} [X_{1} X_{2}] [ \begin{matrix}\underline{β} & 1 \\ \underline{β} & 2 \end{matrix} ]
= σ^{2}k_{2} + \underline{β}_{2}^{T} X_{2}^{T} P_{2·1}X_{2}\underline{β}_{2}

This last term controls the noncentrality parameter, and equals 0 under H_{0} : \underline{β}_{2} = \underline{0}. Once again, E [RSS (\underline{β}_{1}, \underline{β}_{2})] = E [\underline{y}^{T} Q_{12}\underline{y}] = σ^{2} (n − k_{1} − k_{2}), (verify this!) so

\frac{E [SS (\underline{β}_{2}|\underline{β}_{1}) /k_{2}]}{E [RSS (\underline{β}_{1}, \underline{β}_{2}) / (n − k_{1} − k_{2})]} = 1 + \underline{β}_{2}^{T} X_{2}^{T} P_{2·1}X_{2}\underline{β}_{2}

equals 1 when H_{0} : \underline{β}_{2} = \underline{0} is true; and ≥ 1 when H_{0} is false. This is suggestive of the behavior of the F statistic

F = \frac{SS (\underline{β}_{2}|\underline{β}_{1}) /k_{2}}{RSS (\underline{β}_{1}, \underline{β}_{2}) / (n − k_{1} − k_{2})}

5 Testing a General Linear Hypothesis

Let us return to the general linear regression model

\underline{y}_{n×1} = X_{n×k′} \underline{β}_{k′×1} + ε_{n×1}

and now consider inference for a set of general linear combinations of the parameters

\underline{ψ}_{q×1} = L_{q×k′}β_{k′×1}

Some examples of L’s that might be of interest include • L selects the last k_{2} parameters (as in the previous section):

L_{k_{2}×(k_{1}+k_{2})} = \underbrace{[ \begin{matrix}0 & ⋯ & 0 \\ 0 & ⋯ & 0 \\ ⋮ & ... & ⋮ \\ 0 & ⋯ & 0 \end{matrix}}_{k_{1} columns} \underbrace{\begin{matrix}1 & 0 & ⋯ & 0 \\ 0 & 1 & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & 1 \end{matrix} ]}_{k_{2} columns}; L [ β_{1} ; ⋮ ; β_{k_{1}+k_{2}} ] = [ β_{k_{1}+1} ; ⋮ ; β_{k_{1}+k_{2}} ]

• L codes the contrasts to test the hypotheses:

\begin{matrix}H_{0} : μ_{1} = (μ_{2} + μ_{3})/2 vs. H_{1} : & μ_{1} ≠ (μ_{2} + μ_{3})/2 \\ H_{0} : μ_{2} = μ_{3} vs. H_{1} : & μ_{2} ≠ μ_{3} \end{matrix}

in a one-way ANOVA model with three cells, i.e.

L_{2×3} = [ \begin{matrix}1 & −1/2 & −1/2 \\ 0 & 1 & −1 \end{matrix} ]; L [ \begin{matrix}μ_{1} \\ μ_{2} \\ μ_{3} \end{matrix} ] = [ \begin{matrix}μ_{1} − (μ_{2} + μ_{3})/2 \\ μ_{2} − μ_{3} \end{matrix} ]

(note that L is the transpose of the “contrast matrix” that we have dealt with previously). Now, since \underline{ψ}_{q×1} = L_{q×k′} \underline{β}_{k′×1}, it follows that

ˆ\underline{ψ} = Lˆ\underline{β} = L(X^{T} X)^{−1}X^{T} y = Ay ∼ N(ψ, Σ̸ψ̂)

(either by least-squares, or more simply from a maximum likelihood argument), where Σ̸ψ̂ = VarAy = AVaryA^{T} = AA^{T} σ_{ϵ}^{2} Now let

B = Σ̸ψ̂/σ_{ϵ}^{2} = AA^{T} = L(X^{T} X)^{−1}X^{T} X(X^{T} X)^{−1}L^{T} = L(X^{T} X)^{−1}L^{T}

It then follows (see for example Scheffé pp. 1–30) that

(ˆ\underline{ψ} − \underline{ψ})^{T} B^{−1}(ˆ\underline{ψ} − \underline{ψ}) ∼ σ_{ϵ}^{2}χ_{q}^{2}, independent of RSS

Intuition:

Chi-squared: Let R be a “square root” of B^{−1}; then R(ˆ\underline{ψ}−\underline{ψ}) uncorrelates (ˆ\underline{ψ}−\underline{ψ}) and the result is a vector of q independent N(0, 1)′s; the quadratic form (ˆ\underline{ψ}−\underline{ψ})^{T} B^{−1}(ˆ\underline{ψ}−\underline{ψ}) = ∥R(ˆ\underline{ψ}−\underline{ψ})∥^{2} is just the sum of squares of these q N(0, 1)’s;

Indep. of RSS: Like our earlier arguments based on Cochran’s Theorem.

Thus, when the true value of \underline{ψ} = \underline{ψ}_{0},

F = \frac{(ˆ\underline{ψ} − \underline{ψ}_{0})^{T} B^{−1}(ˆ\underline{ψ} − \underline{ψ}_{0})/q}{RSS/(n − k′)} = \frac{(ˆ\underline{ψ} − \underline{ψ}_{0})^{T} [L(X^{T} X)^{−1}L^{T}]^{−1}(ˆ\underline{ψ} − \underline{ψ}_{0})/q}{S_{E}^{2}} ∼ F_{q,n−k′}

This can immediately be “inverted” to provide direct tests of hypotheses:

• Under H_{0} : \underline{ψ} = \underline{ψ}_{0}, the F statistic above is F_{q,n−k′};

• Under H_{1} : \underline{ψ} free, the F statistic tends to be too large.

One can show that when L selects k_{2} parameters and we test for \underline{ψ} = \underline{0}, we get exactly the same test that we derived in connection with least-squares geometry in the last section. In this case [L(X^{T} X)^{−1}L^{T}]^{−1} has the effect of “selecting out” the rows and columns of (X^{T} X)^{−1} corresponding to the parameters β_{k_{1}+1}, …, β_{k_{1}+k_{2}} and inverting the resulting submatrix.

6 Confidence Regions

The same basic fact, that

F = \frac{(ˆ\underline{ψ} − \underline{ψ}_{0})^{T} [L(X^{T} X)^{−1}L^{T}]^{−1}(ˆ\underline{ψ} − \underline{ψ}_{0})/q}{S_{E}^{2}} ∼ F_{q,n−k′}

when \underline{ψ} = \underline{ψ}_{0}, can be used to generate confidence regions for parameters or combinations of parameters. The idea of a confidence region is the same as that of a confidence interval, except that it is a region in q-space for all q parameters \underline{ψ} = (ψ_{1}, …, ψ_{q}).

Definition 3 A 100×(1 − α)% confidence ellipsoid for the parameters \underline{ψ} = (ψ_{1}, …, ψ_{q}) is given by

\{\underline{ψ} : (ˆ\underline{ψ} − \underline{ψ})^{T} [L(X^{T} X)^{−1}L^{T}]^{−1}(ˆ\underline{ψ} − \underline{ψ}) ≤ qS_{E}^{2}F_{α,q,n−k′}\}

Once again, L may be designed so that it selects a particular group of parameters β_{j} of interest, or a particular group of linear combinations of β_{j}’s, etc.

Some examples will help to fix ideas here. The prestige data contains various variables that were thought to be predictive of occupational prestige in Canada, taken largely from 1971 Canadian Census archives. The variables are fairly self-explanatory

> names(prestige)
[1] "education"     "income"        "pct.female"   "occ.prestige" "occ.code"
[6] "occ.type"
> dim(prestige)
[1] 102   6

Here we review a couple of the fits that we did when we first encountered this data set.

> summary(fit.ed <- lm(occ.prestige ~ education, data = prestige))
Call: lm(formula = occ.prestige ~ education, data = prestige)
Residuals:
     Min      1Q Median     3Q   Max
 -26.04 -6.523 0.6611 6.743 18.16
Coefficients:
                 Value Std. Error  t value Pr(>|t|)
(Intercept) -10.7320   3.6771     -2.9186   0.0043
  education   5.3609   0.3320     16.1478   0.0000
Residual standard error: 9.103 on 100 degrees of freedom
Multiple R-Squared: 0.7228
F-statistic: 260.8 on 1 and 100 degrees of freedom, the p-value is 0
Correlation of Coefficients:
           (Intercept)
education -0.9695

We see that neither the slope nor the intercept are likely to be zero, according to their own t-tests; we could also see this with 95% t-intervals:

> qt((1-.05/2),100)
[1] 1.983972
> -10.732+c(-1,1)*1.983972*3.6771     # for beta_0
[1] -18.027263  -3.436737
> 5.3609  +c(-1,1)*1.983972*0.332     # for beta_1
[1] 4.702221 6.019579

on the other hand, the parameter estimates are very highly correlated, and these intervals are not corrected for doing two inferences at once.

We can fix the multiple inference problem with the Bonferroni correction, replacing α with α/(number of tests),

> -10.732+c(-1,1)*qt((1-(.05/2)/2),100)*3.6771     # for beta_0
[1] -19.099801  -2.364199
> 5.3609  +c(-1,1)*qt((1-(.05/2)/2),100)*0.332     # for beta_1
[1] 4.605383 6.116417

but we still don’t learn how the parameters covary. Instead we calculate the confidence ellipse

> X _model.matrix(fit.ed)
> dim(X)
[1] 102   2
> X[1:5,]
                       (Intercept) education
 GOV_ADMINISTRATORS             1      13.11
   GENERAL_MANAGERS             1      12.26
         ACCOUNTANTS             1      12.77
PURCHASING_OFFICERS             1      11.42
             CHEMISTS             1      14.62
> SE2 _sum(residuals(fit.ed)^2)/100
> sqrt(SE2) # check against earlier output
[1] 9.103291
> M _t(X)%*%X/(2*SE2)
> b0 _seq(-20,-2,length=50)
> b1 _seq(3,8,length=50)
> bhat0 _coefficients(fit.ed)[1]
> bhat1 _coefficients(fit.ed)[2]
> qform _function(psi, psi0, M)
+ {
+ #
+ # qform
+ #
+ # computes the quadratic form (psi - psi0)' %*% M %*% (Psi - Psi0)
+ # on every row of the matrix psi, using the single row in psi0
+ #
+ #
+        qfn <- function(x, M)
+        {
+                 x %*% M %*% x
+        }
+        n <- dim(psi)[1]
+        delta <- psi - matrix(rep(psi0, n), byrow = T, ncol = 2)
+        result <- apply(delta, 1, qfn, M)
+        return(result)
+ }
> qform(cbind(0,0),c(bhat0,bhat1),M)

Figure 4: 95% Confidence Ellipse for β_{0} and β_{1} from fit.ed

[1] 1480.219
> qlev _qform(expand.grid(b0,b1),c(bhat0,bhat1),M) # wait a while...
> contour(b0,b1,matrix(qlev,nrow=length(b0),ncol=length(b1)),
+ levels=qf(.95,2,100))
> abline(h=c(4.702221, 6.019579),v=c(-18.027263, -3.436737))

The resulting confidence ellipse is shown in Figure 4. Even though the individual confidence intervals are fairly wide, the confidence ellipse shows that β̂_{1} can vary only a little, for example, at each fixed value of β̂_{0}. Horizontal and vertical lines indicating the original uncorrected t-based confidence intervals for β_{0} and β_{1} are also overlaid on the plot; one can see that the “box” formed by these two intervals both misses some of the area covered by the ellipse—reflecting the fact that the intervals are uncorrected for multiple inferences—and covers much more area than the ellipse covers—reflecting the inefficiency of using individual intervals.

If one did want individual intervals for the β’s that were guaranteed to cover the β’s with 95% confidence, a natural choice would be the “shadow” of the ellipse on each axis. This generates a larger box (not shown in Figure 4); the two intervals would be guaranteed to contain β_{0} and β_{1} in at least (and in fact more than) 95% of repeated samples from the “model” that generated the data.

Scheffé’s method. An extension of the “shadow” idea is Scheffé’s method of generating simultaneous confidence intervals for all possible linear combinations of the original parameters \underline{β}.

Let ψ = ∑_{j=1}^{k′} w_{j}β_{j} = \underline{w}^{T} \underline{β} be any linear combination of the original k′ model parameters; then as usual

ψ̂ = \underline{w}^{T} ˆ \underline{β} = \underline{w}^{T} (X^{T} X)^{−1}X^{T} \underline{y}

and

Var(ψ̂) = \underline{w}^{T} Var(ˆ\underline{β})\underline{w} = σ_{ϵ}^{2}\underline{w}^{T} (X^{T} X)^{−1}\underline{w}

and hence

SEψ̂ = S_{E}\sqrt{\underline{w}^{T} (X^{T} X)^{−1}\underline{w}}

A simple version^{3} and of Scheffé’s method replaces the t cutoff in an uncorrected confidence interval

ψ̂ ± t_{1−α/2,n−k′} · SEψ̂

with a scaled F cutoff:

ψ̂ ± \sqrt{k′F (1 − α, k′, n − k′)} · SEψ̂

The F cutoff k′F (1 − α, k′, n − k′) comes from calculating a confidence ellipsoid for the original set of k′ parameters \underline{β}; the essential idea of Scheffé’s method is that once you have controlled these, you automatically control any linear combination of them. For a proof that these are really simultaneous 100(1 − α)% intervals, see Scheffé, 1959, section 3.5.

Scheffé’s intervals are very useful for

• Data snooping, or any situation in which you are evaluating a lot of linear contrasts or other linear combinations of the original parameters, and you want to control the lengths of your intervals; and

• Graphing simultanous confidence envelopes around regression lines or surfaces.

We will illustrate the latter with the prestige example. We wish to graph

• The data

• The fitted regression line ŷ = β̂_{0} + xβ̂_{1}

• A simultaneous 95% confidence envelope computed at every fitted value

We recall that ˆ \underline{y} = X ˆ\underline{β} = P \underline{y} = H\underline{y} and Var(y_{i}) = h_{ii}σ_{ϵ}^{2}, so that SEŷ_{i} = S_{E}\sqrt{h_{ii}}. Also n = 102 and k′ = 2 for the simple regression of occ.prestige on education.

> attach(prestige)
> plot(education,occ.prestige)
> abline(fit.ed)
> SEy _sqrt(lm.hat(fit.ed)*sum(residuals(fit.ed)^2)/100)
> l _fitted(fit.ed) - sqrt(2*qf(.95,2,100))*SEy
> h _fitted(fit.ed) + sqrt(2*qf(.95,2,100))*SEy
> ord _order(education)
> lines(education[ord],l[ord],lty=3)
> lines(education[ord],h[ord],lty=3)
> detach()

3If you know that you will only be considering ψ’s that are linear combinations of a subset of q of the β’s then you may replace (X^{T} X)^{−1} with L(X^{T} X)^{−1}L^{T} where L selects the subset, and use the cutoff qF (1 − α, q, n − k′)

instead.

Figure 5: 95% Envelope for Fitted Regression Line from from fit.ed

The resulting graph is displayed in Figure 5. Comparing the t cutoff with the F cutoff in Scheffé’s method,

> sqrt(2*qf(.95,2,100))
[1] 2.484873
> qt(1-.05/2,100)
[1] 1.983972

we see that there is a substantial difference between the corrected and the uncorrected envelopes for the regression line here.

Remark: There are other methods to generate confidence envelopes that may do better. For example, at each x_{i} value, we could calculate and plot the values

ũ_{i} = max{β_{0} + x_{i}β_{1} : (β_{0}, β_{1}) ∈ the 95% confidence ellipsoid in Fig 4 }
ℓ̃_{i} = min{β_{0} + x_{i}β_{1} : (β_{0}, β_{1}) ∈ the 95% confidence ellipsoid in Fig 4 }

along with the ŷ_{i}’s.