Worked Mean-Response Regions from Design Matrices

From Design Points to a Joint Mean-Response Region

A fitted regression model can estimate a mean response at any specified design point. When several means are of interest together, the calculations must account for how their estimates vary jointly. The design matrix supplies both the fitted values and the covariance factors needed for an exact confidence region.

We use the normal linear model with a fixed design matrix \(X\) of full column rank \(p\), \(n\) observations, and \(n>p\). The error variance is unknown, and \(s^2\) denotes the usual residual mean square. Each row of a matrix \(A\) represents one specified design point; the corresponding vector of mean responses is \(A\beta\). The rows of \(A\) must be linearly independent for the full-dimensional region below, so the number \(q\) of targets satisfies \(q\leq p\).

Definition: For specified design-point rows collected in a full-row-rank matrix \(A\), the target vector is \(A\beta\), and its fitted estimate is \(A\hat{\beta}\). Its covariance-factor matrix is \(H=A(X'X)^{-1}A'\). Under the normal linear model, the covariance matrix of \(A\hat{\beta}\), conditional on \(X\), is \(\sigma^2 H\).

The exact joint confidence region uses the \(F\) distribution with \(q\) numerator degrees of freedom and \(n-p\) denominator degrees of freedom. Write \(F_{q,n-p}(1-\alpha)\) for its \(1-\alpha\) quantile.

Theorem: If the errors are independent \(N(0,\sigma^2)\), \(X\) has full column rank, \(n>p\), and \(A\) has full row rank \(q\), then an exact \(1-\alpha\) confidence region for \(A\beta\), conditional on \(X\), is

\left\{\mu: (A\hat{\beta}-\mu)'H^{-1}(A\hat{\beta}-\mu) \leq q s^2 F_{q,n-p}(1-\alpha) \right\}.

The left side measures distance from the fitted target vector while accounting for correlation among its components. In particular, it is not generally correct to treat the several mean-response estimates as independent.

Why the Region Has Exact Coverage

Here is the argument behind the theorem. Let \(M_X=I-X(X'X)^{-1}X'\), the residual projection matrix, and let \(\nu=n-p\). Under the stated model, \(A\hat{\beta}-A\beta=A(X'X)^{-1}X'\varepsilon\) is multivariate normal with mean zero and covariance \(\sigma^2H\). Consequently, \((A\hat{\beta}-A\beta)'H^{-1}(A\hat{\beta}-A\beta)/\sigma^2\) has a chi-squared distribution with \(q\) degrees of freedom.

The residual sum of squares is \(\varepsilon'M_X\varepsilon\), and \(M_X\) is an idempotent projection of rank \(\nu\). Thus \(\nu s^2/\sigma^2\) has a chi-squared distribution with \(\nu\) degrees of freedom. The fitted-error projection and residual projection are orthogonal. Since \(\varepsilon\) is normal, the two projected vectors are independent; hence the two chi-squared quantities just described are independent. Dividing their scaled versions gives

\frac{(A\hat{\beta}-A\beta)'H^{-1}(A\hat{\beta}-A\beta)} {q s^2} \sim F_{q,\nu}.

Therefore the probability that this statistic is no larger than \(F_{q,\nu}(1-\alpha)\) is \(1-\alpha\). Substituting \(\mu=A\beta\) in the region gives exactly that event, proving the claimed conditional coverage. \(\square\)

Projecting the Region to Get Intervals

A joint region can be summarized by its coordinate projections: the range of each mean response among all vectors inside the region. If \(H_{jj}\) is the \(j\)th diagonal entry of \(H\), the projection for coordinate \(j\) is centered at its fitted mean and has half-width \(\sqrt{q s^2F_{q,n-p}(1-\alpha)H_{jj}}\).

Result: The projection of the joint region onto its \(j\)th coordinate is

\left[ (A\hat{\beta})_j-\sqrt{K H_{jj}}, \ (A\hat{\beta})_j+\sqrt{K H_{jj}} \right], \qquad K=q s^2F_{q,n-p}(1-\alpha).

To verify this, set \(d=A\hat{\beta}-\mu\). For any point in the region, \(d'H^{-1}d\leq K\). Write \(d_j=e_j'd=(H^{1/2}e_j)'(H^{-1/2}d)\), where \(e_j\) is the \(j\)th coordinate vector. The Cauchy–Schwarz inequality gives

|d_j|^2 \leq (e_j'He_j)(d'H^{-1}d) \leq H_{jj}K.

Equality is attainable: take \(d\) proportional to \(He_j\), choosing the proportionality constant so that \(d'H^{-1}d=K\). Thus the bound is the exact projection endpoint, not merely an upper bound. \(\square\)

These projected intervals come from one joint region, so they collectively contain the true target vector whenever that vector lies in the region. They are not the same construction as separate pointwise intervals using a \(t\) critical value. A projection interval can be wider because the joint region is calibrated to cover all \(q\) targets together.

Worked Calculations

Worked Example: Two Mean Responses at Extrapolation Points

Suppose six observations follow the model with intercept and slope, with predictor values \(-1,-1,0,0,1,1\). The design matrix is

X= \begin{pmatrix} 1&-1\\1&-1\\1&0\\1&0\\1&1\\1&1 \end{pmatrix}.

Take responses \(1+\sqrt{3},1-\sqrt{3},2,2,3,3\). Their fitted line is \(\hat{\beta}=(2,1)'\): the residual vector is \((\sqrt{3},-\sqrt{3},0,0,0,0)'\), whose components sum to zero and whose inner product with the predictor column is zero. Its squared length is \(6\), so \(s^2=6/(6-2)=1.5\).

We want the means at \(x=-2\) and \(x=2\), so \(A\) has rows \((1,-2)\) and \((1,2)\). Since \(X'X\) is diagonal with entries \(6\) and \(4\),

H=A(X'X)^{-1}A' = \begin{pmatrix} 7/6&-5/6\\ -5/6&7/6 \end{pmatrix}.

The fitted means are \(0\) and \(4\). For a 95% region, \(q=2\), \(\nu=4\), and \(F_{2,4}(0.95)\approx6.9443\). Thus \(K=2(1.5)(6.9443)\approx20.8329\). The region is

\left\{\mu: (\hat{\mu}-\mu)' \begin{pmatrix}7/4&5/4\\5/4&7/4\end{pmatrix} (\hat{\mu}-\mu)\leq20.8329 \right\}, \qquad \hat{\mu}=\begin{pmatrix}0\\4\end{pmatrix}.

Each coordinate has \(H_{jj}=7/6\), so its projected half-width is \(\sqrt{20.8329(7/6)}=\sqrt{24.30505}\approx4.9300\). The projections are approximately \([-4.930,4.930]\) at \(x=-2\) and \([-0.930,8.930]\) at \(x=2\). These are the coordinate ranges of the same joint ellipse.

Worked Example: Computing a Region From Observed Responses

Consider predictor values \(-2,-1,0,1,2\) and responses \(1,2,2,4,6\), with an intercept-and-slope model. Here \(X'X=\begin{pmatrix}5&0\\0&10\end{pmatrix}\). The response mean is \(3\), and the sum of predictor times response is \(12\), giving \(\hat{\beta}=(3,1.2)'\).

The fitted values are \(0.6,1.8,3,4.2,5.4\). The residuals are \(0.4,0.2,-1,-0.2,0.6\), and their squared sum is \(0.16+0.04+1+0.04+0.36=1.60\). Therefore \(s^2=1.60/3=8/15\).

For the means at \(x=-1\) and \(x=1\), let \(A\) have rows \((1,-1)\) and \((1,1)\). Direct multiplication gives

H= \begin{pmatrix} 0.3&0.1\\ 0.1&0.3 \end{pmatrix}, \qquad A\hat{\beta}= \begin{pmatrix}1.8\\4.2\end{pmatrix}.

At 95% confidence, \(F_{2,3}(0.95)\approx9.5521\), so \(K=2(8/15)(9.5521)\approx10.1889\). The coordinate half-width is \(\sqrt{10.1889(0.3)}\approx1.7483\). The projections of the joint region are therefore approximately \([0.0517,3.5483]\) at \(x=-1\) and \([2.4517,5.9483]\) at \(x=1\). Their overlap is not an error: the two means are correlated, and each interval is a projection of the full region.

Worked Example: A Quadratic Model and Two Design Points

For predictor values \(-2,-1,0,1,2\), use the quadratic design columns \(1,x,x^2\), and responses \(3,-2.5,8,-0.5,7\). The fitted polynomial coefficients are \(\hat{\beta}=(2,1,0.5)'\), because the fitted values are \(2,1.5,2,3.5,6\) and the residual vector \((1,-4,6,-4,1)'\) is orthogonal to each design column. Its squared length is \(70\), so with \(n-p=2\), \(s^2=35\).

The design cross-product matrix and its inverse are

X'X= \begin{pmatrix}5&0&10\\0&10&0\\10&0&34\end{pmatrix}, \qquad (X'X)^{-1}= \begin{pmatrix} 34/70&0&-10/70\\ 0&1/10&0\\ -10/70&0&5/70 \end{pmatrix}.

At \(x=0\) and \(x=1\), the design rows are \((1,0,0)\) and \((1,1,1)\). Multiplication gives

H= \begin{pmatrix} 17/35&12/35\\ 12/35&13/35 \end{pmatrix}, \qquad A\hat{\beta}= \begin{pmatrix}2\\3.5\end{pmatrix}.

For 95% confidence, \(F_{2,2}(0.95)=19\), so \(K=2(35)(19)=1330\). The half-widths are \(\sqrt{1330(17/35)}=\sqrt{646}\approx25.4165\) and \(\sqrt{1330(13/35)}=\sqrt{494}\approx22.2261\). Thus the projected intervals are approximately \([-23.4165,27.4165]\) at \(x=0\) and \([-18.7261,25.7261]\) at \(x=1\). Their large widths reflect the small residual degrees of freedom and the estimated variance; the arithmetic follows directly from the exact region.

Common Mistakes and Practical Checks

Using the wrong degrees of freedom. The numerator degrees of freedom are the number \(q\) of jointly estimated targets, not the number \(p\) of regression coefficients. The denominator degrees of freedom are \(n-p\).
Dropping the factor \(q\). The region’s right side is \(q s^2F_{q,n-p}(1-\alpha)\). Omitting \(q\) changes both the region and every projected margin.
Using only the diagonal of \(H\) to describe the joint region. The diagonal gives coordinate variances and projection widths, but the off-diagonal entries determine the region’s orientation and dependence structure.
Calling the projections pointwise intervals. They are coordinate projections of a joint region. Their critical value and simultaneous interpretation differ from separate pointwise intervals.
Rounding too early. Keep \(H_{jj}\), \(s^2\), and the critical value unrounded through the product under the square root, then round the final margin.

In practice, first construct \(A\) from the design rows of interest, then calculate \(A\hat{\beta}\) and \(H=A(X'X)^{-1}A'\). Check that \(H\) is symmetric positive definite, compute \(K\), and only then form the region or its projections. This sequence makes it easier to detect a transposed design row, an incorrect degrees-of-freedom choice, or a variance-factor arithmetic error.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

In the joint mean-response region, what determines the numerator and denominator degrees of freedom?
Why do the off-diagonal entries of \(H\) matter even when only the projected intervals are being reported?
For a region with right side \(K\), what is the projection half-width for the \(j\)th mean response?
Why can coordinate projections from a joint region be wider than separate pointwise intervals?
What condition on \(A\) ensures that \(H^{-1}\) exists in the stated region formula?