Linear Regression

WarningDraft

This lecture note is still being edited: the text may contain typos, inaccuracies, and locally underdeveloped passages. If you find something, please let the author know (contact information on the course page).

In this chapter \(\vec{X}\) denotes the observations.

1 Problem statement

In this chapter we consider a fairly special but widely applicable problem in the theory of estimating statistical parameters. The observation in this problem is no longer a sample (a collection of independent random variables), but a vector made up of random variables whose distributions may differ.

1.1 Linear regression without an intercept

ExampleExample 1

There are 2 coins of different denominations, whose masses \(\theta_1\) and \(\theta_2\) are unknown. Using the same scale, the first coin was weighed separately, then the second, and then both coins together. Set up a statistical model for the experiment. Try to estimate \(\theta_1\) and \(\theta_2\).

SolutionSolution

Denote the weighing results \(X_{1}, X_{2}, X_{3}\). If the scale gave error-free measurements, then we would have

\[ \begin{pmatrix} X_1 \\ X_2 \\ X_3 \end{pmatrix} = \begin{pmatrix} \theta _1 \\ \theta _2 \\ \theta _1 + \theta _2 \end{pmatrix} = \theta _1 \cdot \begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix} + \theta _2 \cdot \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix} \]

or, in short,

\[ \vec{X} = \theta _1 \cdot \vec{z}_{\cdot 1} + \theta _2 \cdot \vec{z}_{\cdot 2} = \mathbf{z} \cdot \vec{\theta } \]

where

\[ \vec{X} = \begin{pmatrix} X_1 \\ X_2 \\ X_3 \end{pmatrix}, \quad \mathbf{z} = (\vec{z}_{\cdot 1} \; \vec{z}_{\cdot 2} ) = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix}, \quad \vec{\theta } = \begin{pmatrix} \theta _1 \\ \theta _2 \end{pmatrix} \]

In this case \(\vec{X}\) is non-random, linearly dependent on \(\vec{\theta }\). The 3rd measurement is superfluous; the first 2 would have sufficed to estimate \(\vec{\theta }\).

In practice, any scale admits some error, and it will differ slightly from one experiment to another. The errors in different experiments can be treated as random variables. Denote them \(\vec{\varepsilon } = (\varepsilon_1, \varepsilon_2, \varepsilon_3)^T\). Then the statistical model is:

\[ \vec{X} = \mathbf{z} \cdot \vec{\theta } + \vec{\varepsilon }, \]

where \(\mathbf{z}\vec{\theta }\) is an unknown non-random vector, and \(\vec{\varepsilon }\) is a random vector of errors. Thanks to \(\vec{\varepsilon }\), the observation vector \(\vec{X}\) becomes random.

In general, we are given a bag of coins of \(k\) different types, \(\theta_{1}\) is the mass of a type-1 coin, \(\ldots\), \(\theta_{k}\) is the mass of a type-\(k\) coin. The researcher estimates the masses of the coins. To do this, they weigh combinations of coins \(n\) times: at step \(i \in \left\{ 1,\ldots , n\right\}\) they weigh \(z_{i,1}\) coins of type 1, \(\ldots\), \(z_{i,k}\) coins of type \(k\). If the scale operates without error, then the measurement values will equal linear combinations of \(\theta_{1}, \ldots , \theta_{k}\). But since errors are nevertheless present, the observed quantity obtained at the \(i\)-th weighing \((i \in \{ 1, \ldots , n\} )\) equals

\[ X_i = z_{i, 1} \theta _{1}+\ldots +z_{i, k} \theta _{k}+\varepsilon _{i}, \]

where \(z_{i, j} \in \mathbb {R}\) are known coefficients of the linear combination, \(k\) is the number of coins participating in the \(i\)-th weighing, and \(\varepsilon_{i}\) is the measurement error at the \(i\)-th weighing. The full sample \(\vec{X}\):

\[ \begin{align} \begin{pmatrix} X_1 \\ X_2 \\ \vdots \\ X_n \end{pmatrix} & = \underbrace{\begin{pmatrix} z_{1, 1} \theta _{1}+\ldots +z_{1, k} \theta _{k} \\ z_{2, 1} \theta _{1}+\ldots +z_{2, k} \theta _{k} \\ \vdots \\ z_{n, 1} \theta _{1}+\ldots +z_{n, k} \theta _{k} \end{pmatrix}}_{= \vec{z}_{\cdot 1} \theta _1 + \ldots + \vec{z}_{\cdot k} \theta _k} + \underbrace{\begin{pmatrix} \varepsilon _1 \\ \varepsilon _2 \\ \vdots \\ \varepsilon _n \end{pmatrix}}_{=\vec{\varepsilon }} \end{align} \]

in short:

\[ \vec{X} = \mathbf{z}\vec{\theta } + \vec{\varepsilon }, \]

where

\[ \begin{align} \vec{z}_{\cdot 1} & = \begin{pmatrix} z_{1,1} \\ z_{2,1} \\ \vdots \\ z_{n,1} \end{pmatrix}, \ldots , \vec{z}_{\cdot k} = \begin{pmatrix} z_{1,k} \\ z_{2,k} \\ \vdots \\ z_{n,k} \end{pmatrix}, \\ \mathbf{z} = (\vec{z}_{\cdot 1} \; \vec{z}_{\cdot 2} \; \ldots \; \vec{z}_{\cdot k} ) & = \begin{pmatrix} z_{1,1} & z_{1,2} & \ldots & z_{1,k} \\ z_{2,1} & z_{2,2} & \ldots & z_{2,k} \\ \vdots & \vdots & \vdots & \vdots \\ z_{n,1} & z_{n,2} & \ldots & z_{n,k} \end{pmatrix}, \\ \vec{\theta } & = \begin{pmatrix} \theta _1 \\ \vdots \\ \theta _k \end{pmatrix}, \quad \vec{l} = \mathbf{z}\vec{\theta } \end{align} \]

\(X_i\) is called the dependent variable, \(z_{i,1}, \ldots , z_{i,k}\) are called the independent variables or regressors, and \(\varepsilon_i\) is called the error or residual. Note that \(z_{i,1}, \ldots , z_{i,k}\) are known non-random values, even though they are called independent variables. The only random component in the model is the error \(\varepsilon_i\), and it is through this that \(X_i\) is random.

Suppose the column vectors \(\vec{z}_{\cdot 1} , \vec{z}_{\cdot 2} , \ldots , \vec{z}_{\cdot k}\) are linearly dependent. Let \(\vec{z}_{\cdot i_1}, \vec{z}_{\cdot i_2}, \ldots , \vec{z}_{\cdot i_{k'}}\) be a maximal linearly independent subset of \(\vec{z}_{\cdot 1} , \ldots , \vec{z}_{\cdot k}\). One can prove that the original model and the model \(\vec{X} = \vec{x}_{i_1}\theta_{\cdot i_1} + \ldots + \vec{z}_{\cdot i_{k'}}\theta_{i_{k'}} + \vec{\varepsilon }\) are equivalent, in the sense that whatever method we choose for estimating the unknown coefficients and the quality of the regression (more on this below), the quality of the model turns out the same in both cases. Going forward, for brevity we will assume that the original collection \(\vec{z}_{\cdot 1} , \ldots , \vec{z}_{\cdot k}\) is already linearly independent by itself. Moreover, we will assume that the regressors are collectively nontrivial, namely that the space \(L := \operatorname {span}(\vec{z}_{\cdot 1} , \ldots , \vec{z}_{\cdot k} )\) contains no vector consisting solely of constants. Regression with a constant regressor is discussed in the next section.

As a rule, \(k\) is much smaller than \(n\), and the matrix \(\mathbf{z}\) ends up extremely stretched vertically. The problem consists of estimating the true values \(\theta_{1}, \ldots , \theta_{k}\) from the observation \(X_{1}, \ldots , X_{n}\).

The observation \(\vec{X} = (X_1, \ldots , X_n)^T\) is a random vector from \(\mathbb {R}^{n}\), where \(\vec{l} = \vec{l}(\vec{\theta }) := \mathbf{z}\vec{\theta }\) is a fixed unknown vector, and \(\vec{\varepsilon }\) is a random vector.

About the vector \(\vec{l} = \mathbf{z}\vec{\theta }\) we know that \(\vec{l} = \mathbf{z}\vec{\theta } \in L\), where

\[ L=\left\langle \vec{z}_{\cdot 1}, \ldots , \vec{z}_{\cdot k}\right\rangle = \operatorname {span}\left(\vec{z}_{\cdot 1}, \ldots , \vec{z}_{\cdot k}\right) \]

is a given linear subspace of \(\mathbb {R}^{n}\) of dimension \(k<n\) (\(\vec{z}_{\cdot 1}, \ldots , \vec{z}_{\cdot k}\) are basis column vectors of the space \(L\)). The parameters of the linear regression model described above are \(\vec{\theta }\) and \(\sigma^{2}\), whose estimation is discussed below.

Note that in this case, strictly speaking, according to our definition \(\vec{X}\) is not a sample of size \(n\), since the components of \(\vec{X}\) may have different distributions. However, \(\vec{X}\) can be viewed as a sample of size \(1\) from an \(n\)-dimensional distribution with parameters \(\vec{\theta }, \sigma^2\). Therefore, even in the absence of independence, one can, for example, speak of the existence of a likelihood function, apply factorization criteria to search for sufficient statistics, and so on. The likelihood function in this case is simply the density of the \(n\)-dimensional vector \(\vec{X}\).

1.2 Linear regression with an intercept

The model

\[ X_i = \theta _0 + z_{i, 1} \theta _{1}+\ldots +z_{i, k} \theta _{k}+\varepsilon _{i}. \]

is called linear regression with an intercept. \(\theta_0\) is called the intercept of the regression. We have

\[ \begin{align} \vec{X} = \begin{pmatrix} X_1 \\ X_2 \\ \vdots \\ X_n \end{pmatrix} & = \begin{pmatrix} 1 & z_{1,1} & \ldots & z_{1,k} \\ 1 & z_{2,1} & \ldots & z_{2,k} \\ \vdots & \vdots & \vdots \\ 1 & z_{n,1} & \ldots & z_{n,k} \end{pmatrix} \cdot \begin{pmatrix} \theta _0 \\ \theta _1 \\ \vdots \\ \theta _k \end{pmatrix} + \begin{pmatrix} \varepsilon _1 \\ \varepsilon _2 \\ \vdots \\ \varepsilon _n \end{pmatrix} = \\ & = \vec{z}_{\cdot 0}\theta _0 + \vec{z}_{\cdot 1} \theta _1 + \ldots + \vec{z}_{\cdot k} \theta _k + \vec{\varepsilon } = \\ & = \mathbf{z}\vec{\theta } + \vec{\varepsilon } \end{align} \]

The vector of ones \(\vec{z}_{\cdot 0} = (1, \ldots , 1)^T\) is conveniently denoted \(\overrightarrow {1}\).

1.3 Conditions on the errors \(\vec{\varepsilon }\)

The following conditions are typically imposed on the error vector \(\vec{\varepsilon }\): it is square-integrable (expectations and variances exist), and moreover

1. \(\mathbb {E}\left[\vec{\varepsilon }\right] = \overrightarrow {0} \in \mathbb {R}^{n}\), i.e. all errors are centered;

2. \(\operatorname {Var}\left[\vec{\varepsilon }\right] = \sigma^2 \mathbf{I}_n\), where \(\operatorname {Var}\left[\vec{\varepsilon }\right]\) is the covariance matrix of the vector \(\vec{\varepsilon }\), and \(\mathbf{I}_{n}\) is the \(n \times n\) identity matrix. That is, all errors have the same variance \(\sigma^2\) and different errors are uncorrelated. This property is also called homoscedasticity, from Greek ὁμός (homos, “same”) + σκεδαστός (skedastos, “scattered”).

In practice these conditions on the errors are far from always satisfied, but often the original model can be reduced to one in which the conditions on the errors do hold. For example, under homoscedasticity the condition of centered errors can be relaxed: suppose \(\mathbb {E}\left[\vec{\varepsilon }\right] = \mu \overrightarrow {1}\), i.e. all errors have the same, but now unknown, nonzero expectation. Then from the model \(X_i = \theta_0 + z_{i, 1} \theta_{1}+\ldots +z_{i, k} \theta_{k}+\varepsilon_{i}\) we can pass to the model

\[ X_i = \theta '_0 + z_{i, 1} \theta _{1}+\ldots +z_{i, k} \theta _{k}+\varepsilon '_{i} \]

where

\[ \theta '_0 := \theta _0 + \mu , \quad \varepsilon _i' := \varepsilon _i - \mu . \]

In this model, the errors \(\varepsilon '\) have become centered. For more on reducing a model to standard form, see the problems below.

2 Ordinary least squares (OLS) in parameter estimation

ExampleExample 2

Under the conditions of example Example 1, try to estimate \(\theta_1\) and \(\theta_2\).

SolutionSolution

For simplicity of exposition, suppose \(\theta_1, \theta_2 \in \mathbb {R}\), i.e. the masses may be negative. Ranging over all possible values \(\theta_1, \theta_2 \in \mathbb {R}\), \(\vec{l}(\vec{\theta }) := \mathbf{z}\vec{\theta } = \vec{z}_{\cdot 1} \theta_1 + \vec{z}_{\cdot 2} \theta_2\) forms a linear subspace \(L:= \operatorname {span}(\vec{z}_{\cdot 1} ,\vec{z}_{\cdot 2} )\) in the space \(\mathbb {R}^{3}\).

Suppose \(\vec{X} \in L\), i.e. there are no measurement errors. It is natural to take \(\widehat{\vec{\theta }}_1, \widehat{\vec{\theta }}_2\) equal to the coefficients of the expansion of \(\vec{X}\) in the basis \(\vec{z}_{\cdot 1} , \vec{z}_{\cdot 2}\).

Suppose a more realistic situation: \(\vec{X} \not\in L\), i.e. measurement errors are present. Denote

\[ \widehat{\vec{X}} := \mathbf{z} \cdot \widehat{\vec{\theta }}, \qquad \widehat{\vec{\varepsilon }} := \vec{X} - \widehat{\vec{X}} = \vec{X} - \mathbf{z}\widehat{\vec{\theta }}. \]

• \(\widehat{\vec{X}}\) is the estimate of the observations, i.e. what the observations would have been had there been no errors and the parameter vector \(\vec{\theta }\) coincided with our estimate \(\widehat{\vec{\theta }}\). Note that \(\widehat{\vec{X}}\) is a linear combination of the columns of \(\mathbf{z}\), i.e. \(\widehat{\vec{X}} \in L\).

• \(\widehat{\vec{\varepsilon }}\) is the vector of error estimates.

Given the conditions on the errors, we want to minimize \(\widehat{\vec{\varepsilon }}\), to find the vector closest to \(\vec{X}\) in the subspace \(L\). \(\widehat{\vec{\varepsilon }}\) is an \(n\)-dimensional quantity, so accordingly we need to choose some norm on \(\mathbb {R}^{n}\) that we will minimize. The most natural and convenient choice is the Euclidean norm

\[ \left\| \vec{x}\right\| _2 = \sqrt{x_1^2 + \ldots + x_n^2} \]

Why is it convenient? On the one hand, the Euclidean norm is easy to differentiate to find extrema. On the other hand, the Euclidean norm is generated by an inner product, so we can work with \(\mathbb {R}^{n}\) as a Euclidean space and use the geometry of Euclidean space. The vector \(\widehat{\vec{X}}\) from the subspace \(L\) closest to \(\vec{X}\) is simply the orthogonal projection of \(\vec{X}\) onto \(L\):

\[ \widehat{\vec{X}} = \operatorname {proj}_{L}(X), \]

and \(\widehat{\vec{\theta }}\) is the coordinates of this projection in the basis \(\vec{z}_{\cdot 1} , \vec{z}_{\cdot 2}\).

This method of estimating the parameter \(\vec{\theta }\) and constructing the model is called the method of ordinary least squares, since we minimize the sum of squares of the errors (the Euclidean norm).

So, the problem of estimating the vector \(\vec{\theta }\) can be solved using the method of ordinary least squares (OLS). The least-squares estimate of the parameters equals

\[ \widehat{\vec{\theta }}_{\text{OLS}} \; := \; \operatorname *{arg\, min}_{\vec{\theta } \in \mathbb {R}^{k}}\left\| \vec{X} -\mathbf{z}\vec{\theta }\right\| ^2 \]

Further,

• \(\widehat{\vec{X}}_{\text{OLS}} =\mathbf{z} \cdot \widehat{\vec{\theta }}_{\text{OLS}} = \vec{X}_{L}\) (the orthogonal projection of \(\vec{X}\) onto \(L\)),

• \(\widehat{\vec{\varepsilon }}_{\text{OLS}} = \vec{X} - \widehat{\vec{X}}_{\text{OLS}} = \vec{X} - \vec{X}_{L} = \vec{X}_{L^\perp }\) (the orthogonal projection onto \(L^{\perp }\)).

Going forward, for brevity, whenever we speak of estimates of the parameters, observations, and errors, we will assume that estimation is performed by OLS and write \(\widehat{\vec{\theta }}, \widehat{\vec{X}}, \widehat{\vec{\varepsilon }}\) instead of \(\widehat{\vec{\theta }}_{\text{OLS}} , \widehat{\vec{X}}_{\text{OLS}} , \widehat{\vec{\varepsilon }}_{\text{OLS}}\), although, generally speaking, there are other ways to estimate \(\vec{\theta }\), model the sample, and model the errors.

2.1 Explicit formula for the OLS estimate

Let’s find an explicit formula for \(\widehat{\vec{\theta }}\). Let’s solve a problem from linear algebra.

ExampleExample 3

Let \(\vec{y} \in \mathbb {R}^{n}\), and let \(\vec{z}_{\cdot 1} , \ldots , \vec{z}_{\cdot k} \in \mathbb {R}^{n}\) be a collection of linearly independent vectors (\(k\leq n\)). Find \(\vec{t} \in \mathbb {R}^{k}\) — the coordinates of the projection of \(\vec{y}\) onto \(L = \operatorname {span}(\vec{z}_{\cdot 1} , \ldots , \vec{z}_{\cdot k} )\) in the basis \(\mathbf{z} = (\vec{z}_{\cdot 1} , \ldots , \vec{z}_{\cdot k} )\):
\[ \operatorname {proj}_L(\vec{y}) = \mathbf{z} \cdot \vec{t} = t_1\vec{z}_{\cdot 1} + \ldots + t_k\vec{z}_{\cdot k} \]

SolutionSolution

We have

\[ \begin{align} \left\| \vec{y} - \mathbf{z}\vec{t}\right\| ^2 & = \left\langle \vec{y} - \mathbf{z}\vec{t}, \vec{y} - \mathbf{z}\vec{t} \right\rangle = \\ & = \left\| \vec{y}\right\| ^2 - \underbrace{ \left\langle \vec{y}, \mathbf{z}\vec{t} \right\rangle }_{=\vec{y}^T \cdot \mathbf{z} \cdot \vec{t}} + \left\| \mathbf{z}\vec{t}\right\| ^2 = \\ & = \left\| \vec{y}\right\| ^2 - \left\langle \mathbf{z}^T\vec{y}, \vec{t} \right\rangle + \left\| \mathbf{z}\vec{t}\right\| ^2 = \\ & = \sum _{i=1}^{n}y_i^2 - 2\sum _{i=1}^{k}\left(\mathbf{z}^T \vec{y}\right)_i \cdot t_i + \sum _{i=1}^{n}\left(\mathbf{z}\vec{t}\right)^2_i = \\ & = \sum _{i=1}^{n}y_i^2 - 2\sum _{i=1}^{k}\left(\mathbf{z}^T \vec{y}\right)_i \cdot t_i + \sum _{i=1}^{n}\left(\sum _{j=1}^{k}z_{i,j}t_j\right)^2 \end{align} \]

Let’s differentiate the resulting expression with respect to \(t_{m}\) for each \(m \in \{ 1, \ldots , k\}\) and then simplify:

\[ \begin{align} \frac{\partial }{\partial t_m} \left(\left\| \vec{y} - \mathbf{z}\vec{t}\right\| ^2\right) & = -2\left(\mathbf{z}^T\vec{y}\right)_m + \sum _{i=1}^{n}2\left(\sum _{j=1}^{k}z_{i,j}t_j\right) \cdot z_{i,m} = \\ & = -2\left(\mathbf{z}^T\vec{y}\right)_m + 2\sum _{i=1}^{n} z_{i,m} \cdot \left(\mathbf{z} \cdot \vec{t}\right)_i = \\ & = -2\left(\mathbf{z}^T\vec{y}\right)_m + 2\left(\mathbf{z}^T \cdot \left(\mathbf{z} \cdot \vec{t}\right)\right)_m \end{align} \]

Setting these derivatives to zero to find the argument at which the minimum is attained:

\[ \begin{cases} \frac{\partial }{\partial t_1}\left(\left\| \vec{y} - \mathbf{z}\vec{t}\right\| ^2\right) = 0 \\ \vdots \\ \frac{\partial }{\partial t_k}\left(\left\| \vec{y} - \mathbf{z}\vec{t}\right\| ^2\right) = 0 \\ \end{cases} \iff \begin{cases} -2\left(\mathbf{z}^T\vec{y}\right)_1 + 2\left(\mathbf{z}^T \cdot \left(\mathbf{z} \cdot \vec{t}\right)\right)_1 = 0 \\ \vdots \\ -2\left(\mathbf{z}^T\vec{y}\right)_k + 2\left(\mathbf{z}^T \cdot \left(\mathbf{z} \cdot \vec{t}\right)\right)_k = 0 \\ \end{cases} \]

In matrix-vector form:

\[ -2\mathbf{z}^T\vec{y} + 2\mathbf{z}^T \cdot \mathbf{z} \cdot \vec{t} = 0 \]

Consequently, \(\vec{t} = \left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T\vec{y}\). We obtain

\[ \operatorname {proj}_L(\vec{y}) = \mathbf{z}\left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T \cdot \vec{y} \]

The matrix \(\mathbf{z}\left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T\) is called the projection matrix onto the subspace \(L\).

If the columns of \(\mathbf{z}\) are orthonormal, then \(\mathbf{z}^T\mathbf{z} = \left(\mathbf{z}^T\mathbf{z}\right)^{-1} = \mathbf{I}_{k}\), and the projection matrix simplifies to \(\mathbf{z}\mathbf{z}^T\).

Check question. Where and how was the linear independence of the column vectors \(\vec{z}_{\cdot 1}, \ldots , \vec{z}_{\cdot k}\) used in the solution?

AnswerAnswer

\(\vec{t} = \left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T\vec{y}\).

Altogether, the OLS estimation formula for linear regression is:

\[ \widehat{\vec{\theta }} = \widehat{\vec{\theta }}(\vec{X}) = \left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T \cdot \vec{X} \]

2.2 Geometry of OLS: regression without an intercept

Geometry of OLS regression without an intercept (Three.js version): \(\overrightarrow X\) is the vector of observations, \(\widehat{\overrightarrow X}=\overrightarrow X_L=\mathbf z\widehat{\overrightarrow\theta}\) is its orthogonal projection onto \(L=\operatorname{span}(\overrightarrow z_{\cdot 1},\ldots,\overrightarrow z_{\cdot k})\), and \(\widehat{\overrightarrow\varepsilon}\) is the error vector. The right angle at vertex \(B\) of triangle \(OBA\) illustrates that \(X-\widehat X\perp L\); the coefficient of determination is \(R^2=\cos^2\varphi=\operatorname{ESS}/\operatorname{TSS}\), where \(\varphi=\angle AOB\) — the larger it is, the closer \(\varphi\) is to zero, and the better the linear model. The labels are genuinely rendered LaTeX (KaTeX), anchored to their line/vector and rotated along it; they stay in place as you rotate the scene with the mouse. The mouse wheel zooms.

Consider the graph @LinearRegr:GraphRegrWoParam. By the Pythagorean theorem we have

\[ \left\| \vec{X}\right\| ^2 = \left\| \widehat{\vec{X}}\right\| ^2 + \left\| \widehat{\vec{\varepsilon }}\right\| ^2 \]

• The total sum of squares of the observations \(\left\| \vec{X}\right\|^2\) is called the total sum of squares, denoted \(\operatorname {TSS}\):

  \[ \operatorname {TSS} := \left\| \vec{X}\right\| ^2 = \sum _{i=1}^n X_i^2 \]

• The sum of squares of the estimated observations \(\left\| \widehat{\vec{X}}\right\|^2\) is called the explained sum of squares, denoted \(\operatorname {ESS}\):

  \[ \operatorname {ESS} := \left\| \widehat{\vec{X}}\right\| ^2 = \sum _{i=1}^n \widehat{X}_i^2 \]

• The sum of squares of the errors \(\left\| \widehat{\vec{\varepsilon }}\right\|^2\) is called the residual sum of squares, denoted \(\operatorname {RSS}\):

  \[ \operatorname {RSS} := \left\| \widehat{\vec{\varepsilon }}\right\| ^2 = \left\| \vec{X} - \widehat{\vec{X}}\right\| ^2 = \sum _{i=1}^n \left(X_i - \widehat{X}_i\right)^2 \]

By the Pythagorean theorem, the total sum of squares splits into the explained sum of squares and the sum of squares of the error estimates:

\[ \operatorname {TSS} = \operatorname {ESS} + \operatorname {RSS} \]

2.2.1 Coefficient of determination

The better the OLS-fitted regression describes the behavior of the dependent variable \(\vec{X}\), the smaller the length of the error vector \(\widehat{\vec{\varepsilon }} = \vec{X} - \vec{X}_{L}\) should be compared to the length of the observation vector \(\vec{X}\). The quantity

\[ R^2 := 1 - \frac{\left\| \widehat{\vec{\varepsilon }}\right\| ^2}{\left\| \vec{X}\right\| ^2} = 1 - \frac{\sum _{i=1}^n \left(X_i - \widehat{X}_i\right)^2}{\sum _{i=1}^n X_i^2} \]

is called the coefficient of determination.

The coefficient of determination takes values between \(0\) and \(1\). The closer it is to \(1\), the better the regression.

Various formulas for the coefficient of determination:

\[ \begin{align} R^2 & = 1 - \frac{\left\| \widehat{\vec{\varepsilon }}\right\| ^2}{\left\| \vec{X}\right\| ^2} = 1 - \frac{\operatorname {RSS}}{\operatorname {TSS}} = \\ & = 1 - \frac{\sum _{i=1}^n \left(X_i - \widehat{X}_i\right)^2}{\sum _{i=1}^n X_i^2} = \\ & = \frac{\operatorname {ESS}}{\operatorname {TSS}} = \frac{\left\| \widehat{\vec{X}}\right\| ^2}{\left\| \vec{X}\right\| ^2} = \frac{\sum _{i=1}^n \widehat{X}_i^2}{\sum _{i=1}^n X_i^2} \end{align} \]

The coefficient of determination has several closely related interpretations. On one hand it is the fraction of the spread in the dependent variable \(\vec{X}\) that our model has managed to “explain”: \(\operatorname {ESS}/\operatorname {TSS}\).

On the other hand, one can say that our fitted model \(\widehat{\vec{X}}\) is being compared to the only model available for explaining \(\vec{X}\) in the absence of regressors — the model \(\vec{X} = \overrightarrow {0} + \vec{\varepsilon }\), i.e. \(\vec{X}\) is simply noise made of independent, identically distributed, centered errors. In this approach, the coefficient of determination is one minus the ratio of the spread in the first model to the spread in the second:

\[ 1 - \frac{\operatorname {RSS}}{\operatorname {TSS}} \]

The closer it is to one, the better our model performs compared to the most trivial model.

A simple model against which more complex models are compared is called a reference model or baseline model (English benchmark model). In this case the baseline model is independent, identically distributed noise. \(\operatorname {TSS}\) is the sum of squares of the errors when using such a model. \(\operatorname {ESS}\) is the sum of squares of the difference between the baseline model and our model.

2.3 Geometry of OLS: regression with an intercept

With an intercept \(\theta_0\) present, the baseline model for comparisons changes. Now, in the absence of regressors, the best model for explaining \(\vec{X}\) is the model of a constant \(\theta_0\) plus noise \(\vec{\varepsilon }\). The estimate of \(\theta_0\) minimizing the sum of squared errors in such a model is the sample mean \(\overline{X}\). The baseline model for modeling \(\vec{X}\):

\[ \left(\overline{X}, \ldots , \overline{X}\right)^T = \overline{X} \cdot \overrightarrow {1}, \]

where \(\overrightarrow {1} = (1, \ldots , 1)^T\) is the vector in \(\mathbb {R}^{n}\) made up of ones.

\(\operatorname {TSS}\), \(\operatorname {ESS}\) are, respectively, the sum of squares of the errors when using the baseline model and the sum of squares of the difference between the baseline model and our model. Their computation now also changes. The computation of \(\operatorname {RSS}\) remains the same.

• Total sum of squares:

  \[ \operatorname {TSS} = \left\| \vec{X} - \overline{X}\cdot \overrightarrow {1}\right\| ^2 = \sum _{i=1}^n (X_i - \overline{X})^2 \]

• Explained sum of squares:

  \[ \operatorname {ESS} = \left\| \widehat{\vec{X}} - \overline{X}\cdot \overrightarrow {1}\right\| ^2 = \sum _{i=1}^n (\widehat{X}_i- \overline{X})^2 \]

• Residual sum of squares:

  \[ \operatorname {RSS} = \left\| \vec{X} - \widehat{\vec{X}}\right\| ^2 = \sum _{i=1}^n \left(X_i - \widehat{X}_i\right)^2 \]

Geometry of OLS regression with an intercept (Three.js version): \(D\) is the point of the baseline model \(\overline X\cdot\overrightarrow 1\), \(A\) is the vector of observations \(\overrightarrow X\), \(B\) is its projection \(\widehat{\overrightarrow X}=\overrightarrow X_L=\mathbf z\widehat{\overrightarrow\theta}\) onto \(L\). The right angle at vertex \(B\) of triangle \(DBA\) illustrates that \(X-\widehat X\perp L\); the coefficient of determination is \(R^2=\cos^2\varphi=\operatorname{ESS}/\operatorname{TSS}\), where \(\varphi=\angle ADB\) — the larger it is, the closer \(\varphi\) is to zero, and the better the linear model. The labels are genuinely rendered LaTeX (KaTeX), anchored to their line/vector and rotated along it; they stay in place as you rotate the scene with the mouse. The mouse wheel zooms.

Consider the graph @LinearRegr:GraphRegrWithParam. By the Pythagorean theorem for triangle \(ABD\) we again have (recall that \(\widehat{\vec{X}} = \vec{X}_{L}\) when the model is fitted by OLS):

\[ \begin{align} \operatorname {TSS} & = \left\| \vec{X} - \overline{X}\cdot \overrightarrow {1}\right\| ^2 = \left|AD\right|^2 = \\ & = \left|BD\right|^2 + \left|AB\right|^2 = \left\| \widehat{\vec{X}} - \overline{X}\cdot \overrightarrow {1}\right\| ^2 + \left\| \vec{X} - \widehat{\vec{X}}\right\| ^2 = \\ & = \operatorname {ESS} + \operatorname {RSS} \end{align} \]

Coefficient of determination:

\[ \begin{align} R^2 & = 1 - \frac{\operatorname {RSS}}{\operatorname {TSS}} = 1 - \frac{\sum _{i=1}^n \left(X_i - \widehat{X}_i\right)^2}{\sum _{i=1}^n (X_i - \overline{X})^2} = \\ & = \frac{\operatorname {ESS}}{\operatorname {TSS}} = \frac{\sum _{i=1}^n \left(\overline{X} - \widehat{X}_i\right)^2}{\sum _{i=1}^n (X_i - \overline{X})^2} \end{align} \]

For the same sample, the coefficient of determination for the model with an intercept can turn out to be either larger or smaller than the coefficient of determination of the model without one, despite the fact that the first model is obviously more flexible.

Linear regression with one regressor \(\vec{z}_{\cdot 1}\) and intercept \(\theta_0\), \(n=6\) (regression parameters denoted \(\theta_0, \theta_1\)). Orange points are the sample points \((z_{1,i}, X_i)\). Green points are the regression points \((z_{1,i}, \widehat{X}_i)\). The total area of the violet squares equals \(\operatorname{TSS}\), and of the red ones — \(\operatorname{RSS}\).

2.4 Properties of the OLS estimates

Recall, \(\widehat{\vec{\theta }} = \left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T\vec{X}\). Denote by \(\widetilde{k}\) the number of columns of the matrix \(\mathbf{z}\):

\[ \widetilde{k} := \begin{cases} k, & \text{ if the regression has no intercept }\theta _0;\\ k+1 , & \text{ if the regression has an intercept }\theta _0. \end{cases} \]

1. Expectation and covariance matrix of \(\widehat{\vec{\theta }}\):

   \[ \mathbb {E}_{\vec{\theta }, \sigma ^2}\left[\widehat{\vec{\theta }}\right] = \vec{\theta }, \qquad \operatorname {Var}_{\vec{\theta }, \sigma ^2}\left[\widehat{\vec{\theta }}\right] = \sigma ^2 \left(\mathbf{z}^T \mathbf{z}\right)^{-1} \]

   Thus, the ordinary least squares estimator is unbiased.

2. Knowing the covariance matrix of the vector \(\widehat{\vec{\theta }}\), we can find the expectation of the (squared) distance from \(\vec{X}\) to the subspace \(L \subset \mathbb {R}^{n}\):

   \[ \mathbb {E}_{\vec{\theta }, \sigma ^2}\left[\left\| \vec{X} - \underbrace{\mathbf{z}\widehat{\vec{\theta }}}_{=\vec{X}_{L} = \widehat{\vec{X}}}\right\| ^2\right] = \mathbb {E}_{\vec{\theta }, \sigma ^2}\left[\left\| \underbrace{\widehat{\vec{\varepsilon }}}_{=\vec{X}_{L^\perp }}\right\| ^2\right] = \mathbb {E}_{\vec{\theta }, \sigma ^2}\left[\operatorname {RSS}\right] = \left(n-\widetilde{k}\right)\sigma ^2 \]

   This last equality gives an unbiased estimator of the parameter \(\sigma^{2}\):

   \[ \widehat{\sigma ^{2}} := \frac{1}{n-\widetilde{k}}\left\| \vec{X}-\mathbf{z} \widehat{\vec{\theta }}\right\| ^{2} = \frac{\operatorname {RSS}}{n-\widetilde{k}}. \]

Sampling distribution of the OLS estimator

The cloud shows estimates \((\hat\theta_0,\hat\theta_1)\) over \(800\) independent repetitions of the experiment \(X_i=\theta_0+\theta_1 z_i+\varepsilon_i\), \(\varepsilon_i\sim\mathcal N(0,\sigma^2)\); the star is the true \(\theta=(2,\,0.5)\); the red ellipse is the theoretical \(95\%\) ellipse of the covariance matrix \(\sigma^2(\vec z^T\vec z)^{-1}\).

\(n=\)20 \(\sigma=\)1.0 Design \(z_i\)UniformClustered near zero

Fraction of points inside the \(95\%\) ellipse: — (theoretically \(\approx 0.95\))

ExampleExample 4

There are 2 pieces of cheese with weights \(\theta_1\) and \(\theta_2\). Using the same scale, the first, then the second, and then both pieces together were weighed. Find the least-squares estimates of \(\theta_1\) and \(\theta_2\), and also an unbiased estimate of the variance of the measurement error. Find the coefficient of determination.

SolutionSolution

Suppose the weighings gave results \(X_{1}, X_{2}, X_{3}\). Set \(\vec{X}=\left(X_{1}, X_{2}, X_{3}\right)^{T}\). Then

\[ \vec{X} =\mathbf{z}\vec{\theta } + \vec{\varepsilon }, \quad \text{where} \quad \mathbf{z} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix}, \quad \vec{\theta } = \begin{pmatrix} \theta _1 \\ \theta _2 \end{pmatrix}, \quad \vec{\varepsilon } = \begin{pmatrix} \varepsilon _1 \\ \varepsilon _2 \\ \varepsilon _3 \end{pmatrix} \]

Let’s first compute \(\mathbf{z}^T\mathbf{z}\), its inverse, and \(\mathbf{z}^T\vec{X}\):

\[ \mathbf{z}^T\mathbf{z} = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}, \qquad \left(\mathbf{z}^T\mathbf{z}\right)^{-1} = \frac13\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}, \qquad \mathbf{z}^T\vec{X} = \begin{pmatrix} X_1+X_3 \\ X_2+X_3 \end{pmatrix}. \]

Then

\[ \begin{align} \widehat{\vec{\theta }} & = \left(\begin{array}{c} \widehat{\vec{\theta }}_1 \\ \widehat{\vec{\theta }}_2 \end{array}\right)=\left(\mathbf{z}^{T} \mathbf{z}\right)^{-1} \mathbf{z}^{T} \vec{X}= \frac13\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}\begin{pmatrix} X_1+X_3 \\ X_2+X_3 \end{pmatrix} = \left(\begin{array}{c} \frac{2}{3} X_{1}-\frac{1}{3} X_{2}+\frac{1}{3} X_{3} \\ -\frac{1}{3} X_{1}+\frac{2}{3} X_{2}+\frac{1}{3} X_{3} \end{array}\right) \\ \widehat{\vec{X}} & = \mathbf{z} \widehat{\vec{\theta }} = \begin{pmatrix} \widehat{\vec{\theta }}_1 \\ \widehat{\vec{\theta }}_2 \\ \widehat{\vec{\theta }}_1 + \widehat{\vec{\theta }}_2 \\ \end{pmatrix} = \begin{pmatrix} \frac{2}{3} X_{1}-\frac{1}{3} X_{2}+\frac{1}{3} X_{3} \\ -\frac{1}{3} X_{1}+\frac{2}{3} X_{2}+\frac{1}{3} X_{3} \\ \frac{1}{3} X_{1}+\frac{1}{3} X_{2}+\frac{2}{3} X_{3} \\ \end{pmatrix} \\ \widehat{\sigma ^{2}} & = \frac{1}{3-2}\operatorname {RSS} = \| \vec{X}-\mathbf{z} \widehat{\vec{\theta }}\| ^{2}=\left(X_{1}-\widehat{\vec{\theta }}_1\right)^{2}+\left(X_{2}-\widehat{\vec{\theta }}_2\right)^{2}+\left(X_{3}-\widehat{\vec{\theta }}_1-\widehat{\vec{\theta }}_2\right)^{2} \end{align} \]

Let’s substitute the explicit expressions for \(\widehat{\vec{\theta }}_1, \widehat{\vec{\theta }}_2\) (and for their sum — the third coordinate of \(\widehat{\vec{X}}\), found above) and simplify each residual individually:

\[ \begin{align} X_1 - \widehat{\vec{\theta }}_1 & = X_1 - \left(\tfrac {2}{3} X_{1}-\tfrac {1}{3} X_{2}+\tfrac {1}{3} X_{3}\right) = \tfrac 13 \left(X_1+X_2-X_3\right), \\ X_2 - \widehat{\vec{\theta }}_2 & = X_2 - \left(-\tfrac {1}{3} X_{1}+\tfrac {2}{3} X_{2}+\tfrac {1}{3} X_{3}\right) = \tfrac 13 \left(X_1+X_2-X_3\right), \\ X_3 - \widehat{\vec{\theta }}_1 - \widehat{\vec{\theta }}_2 & = X_3 - \left(\tfrac {1}{3} X_{1}+\tfrac {1}{3} X_{2}+\tfrac {2}{3} X_{3}\right) = -\tfrac 13\left(X_1+X_2-X_3\right). \end{align} \]

All three residuals are equal in absolute value to \(\frac{\left|X_1 + X_2 - X_3\right|}{3}\), so

\[ \widehat{\sigma ^2} = 3 \cdot \left(\frac{X_1 + X_2 - X_3}{3}\right)^2 = \frac{\left(X_{1}+X_{2}-X_{3}\right)^2}{3}. \]

Further, we have a regression without an intercept, so by the Pythagorean theorem \(\operatorname {TSS} = \operatorname {ESS} + \operatorname {RSS}\), whence \(\operatorname {ESS} = \operatorname {TSS} - \operatorname {RSS}\); here \(\operatorname {RSS} = (3-2)\widehat{\sigma^2} = \widehat{\sigma^2}\) was already found above. We have

\[ \begin{align} \operatorname {TSS} & = \sum _{i=1}^3 X_i^2 = X_1^2 + X_2^2 + X_3^2, \\ \operatorname {ESS} & = \operatorname {TSS} - \operatorname {RSS} = \left(X_1^2+X_2^2+X_3^2\right) - \frac{\left(X_1+X_2-X_3\right)^2}{3} = \\ & = \left(X_1^2+X_2^2+X_3^2\right) - \frac{X_1^2+X_2^2+X_3^2+2X_1X_2-2X_1X_3-2X_2X_3}{3} = \\ & = \frac{2}{3}X_1^2 + \frac{2}{3}X_2^2 + \frac{2}{3}X_3^2 - \frac{2}{3}X_1X_2 + \frac{2}{3}X_1X_3 + \frac{2}{3}X_2X_3, \\ R^2 & = \frac{\operatorname {ESS}}{\operatorname {TSS}} = \frac{2}{3}\frac{X_1^2 + X_2^2 + X_3^2 - X_1X_2 + X_1X_3 + X_2X_3}{X_1^2 + X_2^2 + X_3^2}. \end{align} \]

ProblemProblem 1

In a quadrilateral \(A B C D\), independent measurements of the angles, all of equal precision,
\[ A B D, D B C, A B C, B C D, C D B, B D A, C D A, D A B \] in degrees gave the results
\[ 50.78, 30.25,78.29,99.57,50.42,40.59,88.87,89.86 \] respectively. Assuming that the measurement errors are independently distributed, have the same variance \(\sigma^2\) and zero expectation, find the OLS estimates of the angles \(\theta_{1}=A B D, \theta_{2}=D B C, \theta_{3}=C D B, \theta_{4}=B D A\) and of the unknown variance \(\sigma^{2}\). Find the coefficient of determination.

ProblemProblem 2

Let
\[ X_{i}=\theta _0+i \theta _1+\varepsilon _{0}+\ldots +\varepsilon _{i}, \] \(i=0,1, \ldots , n\), where \(\theta_0, \theta_1\) are unknown parameters, and the random variables \(\varepsilon_{0}, \ldots , \varepsilon_{n}\) are independently distributed, have the same variance \(\sigma^2\) and zero expectation. Reduce the problem to a linear model and find the least-squares estimates of \(\theta_0\) and \(\theta_1\), as well as an unbiased estimate of \(\sigma^{2}\). Find the coefficient of determination \(R^2\).

ProblemProblem 3

1. Prove that

\[ \mathbb {E}_{\vec{\theta }, \sigma ^2}\left[\widehat{\vec{\theta }}\right] = \vec{\theta }, \qquad \operatorname {Var}_{\vec{\theta }, \sigma ^2}\left[\widehat{\vec{\theta }}\right] = \sigma ^2 \left(\mathbf{z}^T \mathbf{z}\right)^{-1} \]

2. Prove that

\[ \mathbb {E}_{\vec{\theta }, \sigma ^2}\left[\underbrace{\left\| \vec{X} - \mathbf{z}\widehat{\vec{\theta }}\right\| ^2}_{=\operatorname {RSS}}\right] = (n-\widetilde{k})\sigma ^2 \qquad \]

and

\[ \widehat{\sigma ^{2}} = \frac{1}{n-\widetilde{k}} \left\| \vec{X} - \mathbf{z}\widehat{\vec{\theta }}\right\| ^2 \]

is an unbiased estimate of \(\sigma^2\).

HintHint

If \(\vec{\xi }\) is a random vector of dimension \(n\), then its covariance matrix can be computed as:
\[ \operatorname {Var}\left[\vec{\xi }\right] = \mathbb {E}\left[\vec{\xi } \cdot \vec{\xi }^T\right] - \mathbb {E}\left[\vec{\xi }\right] \cdot \left(\mathbb {E}\left[\vec{\xi }\right]\right)^T. \] Moreover, if \(A \in \mathbb {R}^{m \times n}\) is a matrix of dimension \(m \times n\), then \(A\vec{\xi }\) is a new random vector of dimension \(m\), whose covariance matrix and the covariance matrix of the original vector are related as follows:
\[ \operatorname {Var}\left[A \vec{\xi }\right] = A \cdot \operatorname {Var}\left[\vec{\xi }\right] \cdot A^T \]

3 Gaussian linear regression

We will devote the second half of this chapter exclusively to a special case of the linear regression model of exceptional practical importance, namely the Gaussian linear model. This is the linear regression model \(\vec{X}=\mathbf{z}\vec{\theta }+\vec{\varepsilon }\), in which the errors are distributed according to a multivariate Gaussian law:

\[ \vec{\varepsilon } \sim \mathscr {N}\left(\overrightarrow {0}, \sigma ^2 \mathbf{I}_n\right) \]

Recall that \(\mathbf{z}\vec{\theta }\) is a non-random \(n\)-dimensional vector, so

\[ \vec{X} \sim \mathscr {N}\left(\mathbf{z}\vec{\theta }, \; \sigma ^{2} \mathbf{I}_{n}\right) \]

As already mentioned, \(\vec{X}\) can be viewed as a sample of size \(1\) from an \(n\)-dimensional distribution, in this case Gaussian. It has density with respect to Lebesgue measure on \(\mathbb {R}^{n}\):

\[ \begin{align} f_{\vec{\theta }, \sigma ^2}(\vec{x}) & = \frac{1}{\sqrt{(2\pi \sigma ^2)^n}} \mathrm{exp}\left(-\frac{1}{2\sigma ^2}\left(\vec{x} - \mathbf{z}\vec{\theta }\right)^T \left(\vec{x} - \mathbf{z}\vec{\theta }\right)\right) = \\ & = \frac{1}{\sqrt{(2\pi \sigma ^2)^n}} \mathrm{exp}\left(-\frac{1}{2\sigma ^2} \left\| \vec{x} - \mathbf{z}\vec{\theta }\right\| ^2\right), \qquad \vec{x} \in \mathbb {R}^{n} \end{align} \]

3.1 Complete, sufficient statistic. Optimal estimators

Theorem 1 In the Gaussian linear model, the statistic \(\left(\vec{X}_{L}, \left\| \vec{X}_{L^\perp }\right\|^2\right)\) is complete and sufficient.

Corollary 1 In the Gaussian linear model, the statistic \(\left(\widehat{\vec{\theta }}, \widehat{\sigma^2}\right)\) is complete and sufficient.

ProblemProblem 4

Prove theorem Theorem 1. Prove corollary Corollary 1.

From the standard regression model it is known that the statistics

\[ \begin{align} \widehat{\vec{\theta }} & = \operatorname *{arg\, min}_{\vec{t} \in \mathbb {R}^{k}} \left\| \vec{X} - \mathbf{z}\vec{t}\right\| = \left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T \vec{X}, \\ \widehat{\sigma ^2} & = \frac{1}{n-\widetilde{k}} \left\| \vec{X} - \mathbf{z}\widehat{\vec{\theta }}\right\| ^2 = \frac{1}{n-\widetilde{k}} \left\| \vec{X} - \vec{X}_{L}\right\| ^2 = \\ & = \frac{1}{n-\widetilde{k}}\left\| \vec{X}_{L^\perp }\right\| ^2 = \frac{1}{n -\tilde{k}}\operatorname {RSS} \end{align} \]

are unbiased estimators of \(\vec{\theta }, \sigma^2\). They are functions of the complete, sufficient statistic, and hence, by the Lehmann-Scheffé theorem, \(\widehat{\vec{\theta }}, \widehat{\sigma^2}\) are optimal estimators of \(\vec{\theta }, \sigma^2\) in the Gaussian regression model.

ProblemProblem 5

Two loads of masses \(a\) and \(b\) are weighed as follows: the first load is weighed \(n_{1}\) times (all measurement errors having distribution \(\left.\mathscr {N}\left(0, \sigma^{2}\right)\right)\), the second load is weighed \(n_{2}\) times on the same scale, then the first and second loads together are weighed \(n_{3}\) times on a different scale, all of whose measurement errors have distribution \(\mathscr {N}\left(0, 3 \sigma^{2}\right)\). Reduce the problem to a linear model and find the least-squares estimates of \(a\) and \(b\), as well as the optimal estimate of \(\sigma^{2}\).

ProblemProblem 6

Let the vector \(\vec{X}=\left(X_{1}, \ldots , X_{n}\right)^T\) have distribution \(\mathscr {N}\left(a \cdot \vec{b}, \sigma^{2} \vec{\Sigma }\right)\), where \(\vec{b}\) is a known vector of dimension \(n\), \(\vec{\Sigma }\) is a known positive-definite matrix, and \(a\) and \(\sigma^{2}\) are unknown parameters. Reduce the problem to a Gaussian linear model and find the least-squares estimates of the parameters \(a\) and \(\sigma^{2}\).

3.2 Confidence intervals for the parameters

We have constructed point unbiased estimates of the parameters of the linear model, which in the Gaussian case turned out to be optimal. Let us now construct confidence intervals for these parameters.

Consider a motivating example.

ExampleExample 5

Let \(X_1, \ldots , X_n\) be a sample from \(\mathscr {N}\left(\mu , \sigma^2\right)\), where

1. \(\mu\) is unknown, \(\sigma^2\) is known;

2. \(\mu\) is known, \(\sigma^2\) is unknown;

3. \(\mu\) is unknown, \(\sigma^2\) is unknown. Construct confidence intervals of confidence level \(\alpha\) for the unknown parameters.

SolutionSolution

This model is the simplest example of Gaussian linear regression. Namely, let

\[ \vec{X} = \mathbf{z}\vec{\theta } + \vec{\varepsilon } \]

where we set

\[ \vec{X} := \begin{pmatrix} X_1 \\ \vdots \\ X_n \end{pmatrix}, \quad \mathbf{z} := \begin{pmatrix} 1 \\ \vdots \\ 1 \end{pmatrix}, \quad \vec{\theta } := (\mu ), \quad \vec{\varepsilon } := \begin{pmatrix} X_1 - \mu \\ \vdots \\ X_n - \mu \end{pmatrix} \]

Why is this Gaussian regression? Most importantly, we need to check the conditions on the errors. From the assumption we know that the components of \(\vec{X}\) are independent and Gaussian. So \(\vec{X}\) is a Gaussian vector, and

\[ \vec{X} \sim \mathscr {N}\left(\mu \cdot \vec{1}, \sigma ^2 \cdot \mathbf{I}_{n}\right) \]

Hence

\[ \vec{\varepsilon } \sim \mathscr {N}\left(\vec{0}, \sigma ^2 \cdot \mathbf{I}_{n}\right) \]

The conditions on the errors (independence, centeredness, identical Gaussian distribution) are satisfied.

Case c), where both \(\vec{\theta } = (\mu )\) and \(\sigma^2\) are unknown, corresponds fully to the regression model.

Let’s first note that \(\sum_{i=1}^n X_i \sim \mathscr {N}\left(n\mu , n\sigma^2\right)\), so the sample mean has distribution

\[ \overline{X} = \frac{1}{n}\sum _{i=1}^n X_i \sim \mathscr {N}\left(\mu , \sigma ^2/n\right) \]

If we center and normalize (divide by the standard deviation), we get the standard normal distribution:

\[ \frac{\overline{X} - \mu }{\sigma /\sqrt{n}} \sim \mathscr {N}\left(0, 1\right) \]

Recall that to construct a confidence interval, we need to find a pivotal statistic: its distribution must not depend on the unknown parameters.

Enum-item(1)

\(\mu\) is unknown, \(\sigma^2\) is known, so we can freely use \(\sigma\). Our pivotal statistic will be exactly

\[ \frac{\overline{X} - \mu }{\sigma /\sqrt{n}} \]

Enum-item(2)

\(\mu\) is known, \(\sigma^2\) is unknown. We need to estimate the variance. A good estimate for the variance \(\sigma^2\) can be given by the sample variance computed using the true expectation:

\[ \widehat{\sigma ^2} := \frac{1}{n}\sum _{i=1}^n \left(X_i - \mathbb {E}\left[X_i\right]\right)^2 = \frac{1}{n}\sum _{i=1}^n \left(X_i - \mu \right)^2 \]

To construct a confidence interval, we need to find the distribution of \(\widehat{\sigma^2}\). Note that \(X_i - \mu \sim \mathscr {N}\left(0, \sigma^2\right)\), so \(\frac{X_i - \mu }{\sigma } \sim \mathscr {N}\left(0, 1\right)\). Consequently

\[ n\frac{\widehat{\sigma ^2}}{\sigma ^2} = \sum _{i=1}^n \left(\frac{X_i - \mu }{\sigma }\right)^2 \]

is a sum of squares of \(n\) independent standard normal random variables.

Definition 1 The distribution of a random variable equal to the sum of squares of \(n\) independent standard normal random variables is called the chi-squared distribution with \(n\) degrees of freedom. Notation: \(\chi^2_n\). One can prove that the chi-squared distribution is a special case of the gamma distribution: \(\chi^2_n = \Gamma \left(\frac{n}{2}, \frac{1}{2}\right)\).

The distribution of the statistic \(n\frac{\widehat{\sigma^2}}{\sigma^2} \sim \chi^{2}_{n}\) does not depend on the parameter, so it is a pivotal statistic. Let \(Q_{\chi^2_n}(\alpha )\) be the quantile of level \(\alpha\). Then

\[ \mathbb {P}\left(0 \quad \leq \quad n\frac{\widehat{\sigma ^2}}{\sigma ^2} \quad \leq \quad Q_{\chi ^2_n}(\alpha )\right) = \alpha \]

Let’s solve the inequality inside the probability with respect to \(\sigma^2\). We get

\[ \mathbb {P}\left(n \widehat{\sigma ^2} \cdot \frac{1}{Q_{\chi ^2_n}(\alpha )} \quad \leq \quad \sigma ^2 \quad < \quad +\infty \right) = \alpha \]

Enum-item(3)

\(\mu\) is unknown, \(\sigma^2\) is unknown. If \(\sigma^2\) were known, constructing a confidence interval for \(\mu\) would present no difficulty (see the part above). In our case \(\sigma^2\) is unknown, so let’s try to estimate it. A good estimate for the variance \(\sigma^2\) is given by the sample variance \(S^2\):

\[ \widehat{\sigma ^2} := S^2 = \frac{1}{n}\sum _{i=1}^n \left(X_i - \overline{X}\right)^2 \]

Now we need to find the distribution of \(\widehat{\sigma^2}\). Let \(L = \operatorname {span}(\vec{z}_{\cdot 1} ) = \operatorname {span}(\vec{1})\) be the one-dimensional linear subspace of \(\mathbb {R}^{n}\) generated by the vector of ones. As we established above, the projection of \(\vec{X}\) onto \(L\) is the vector made up of sample means:

\[ \vec{X}_{L} = \begin{pmatrix} \overline{X} \\ \vdots \\ \overline{X} \end{pmatrix} = \overline{X} \cdot \vec{1} \]

Then

\[ n\widehat{\sigma ^2} = \sum _{i=1}^n \left(X_i - \overline{X}\right)^2 = \left\| \vec{X} - \overline{X} \cdot \vec{1}\right\| ^2 = \left\| \vec{X}_{L^{\perp }}\right\| ^2 \]

\(\widehat{\sigma^2}\) is a function of the vector \(\vec{X}_{L^{\perp }}\), while \(\overline{X}\) is a function of \(\vec{X}_{L}\).

Let \(\mathbf{z}_{L_2} \in \mathbb {R}^{n \times (n-1)}\) be the matrix whose columns form an orthonormal basis of the subspace \(L_2 := L^{\perp }\) (then \(\mathbf{z}_{L_2}\mathbf{z}_{L_2}^T\) is the projection matrix onto \(L_2\), see the remark above about orthonormal columns). We have

\[ \begin{align} n\widehat{\sigma ^2} & = \left\| \vec{X}_{L^{\perp }}\right\| ^2 = \left\langle \vec{X}_{L^{\perp }}, \vec{X}_{L^{\perp }} \right\rangle = \\ & = \left\langle \mathbf{z}_{L_2}\mathbf{z}_{L_2}^T \vec{X}, \mathbf{z}_{L_2}\mathbf{z}_{L_2}^T \vec{X} \right\rangle = \\ & = \left(\mathbf{z}_{L_2}\mathbf{z}_{L_2}^T \vec{X}\right)^T \cdot \mathbf{z}_{L_2}\mathbf{z}_{L_2}^T \vec{X} = \\ & = \vec{X}^T \mathbf{z}_{L_2} \underbrace{\mathbf{z}_{L_2}^T \cdot \mathbf{z}_{L_2}}_{=\mathbf{I}_{n-1}}\mathbf{z}_{L_2}^T \vec{X} = \\ & = \left\langle \mathbf{z}_{L_2}^T\vec{X}, \mathbf{z}_{L_2}^T\vec{X} \right\rangle = \\ & = \left\| \mathbf{z}_{L_2}^T\vec{X}\right\| ^2 \end{align} \]

Note that in the last 2 lines the inner product is already taken in \(\mathbb {R}^{n-1}\).

We need to find the distribution of \(\mathbf{z}_{L_2}^T\vec{X}\). This is a linear transformation of a Gaussian vector, so it is also a Gaussian vector, with expectation:

\[ \mathbb {E}\left[\mathbf{z}_{L_2}^T\vec{X}\right] = \mathbf{z}_{L_2}^T\mathbb {E}\left[\vec{X}\right] = \mathbf{z}_{L_2}^T \cdot \mu \cdot \underbrace{\vec{1}}_{\in \mathbb {R}^{n}} = \vec{0} \in \mathbb {R}^{n-1} \]

Why did we get \(\vec{0}\)? Because the columns of \(\mathbf{z}_{L_2}\) are orthogonal to \(L_1 = \operatorname {span}(\vec{1})\).

Covariance matrix:

\[ \begin{align} \operatorname {Var}\left[\mathbf{z}_{L_2}^T\vec{X}\right] & = \mathbf{z}_{L_2}^T \cdot \operatorname {Var}\left[\vec{X}\right] \cdot \mathbf{z}_{L_2} = \mathbf{z}_{L_2}^T \cdot \sigma ^2 \cdot \mathbf{I}_{n} \cdot \mathbf{z}_{L_2} = \sigma ^2 \cdot \mathbf{z}_{L_2}^T \cdot \mathbf{z}_{L_2} = \\ & = \sigma ^2 \cdot \mathbf{I}_{n-1} \end{align} \]

In total,

\[ \mathbf{z}_{L_2}^T\vec{X} \sim \mathscr {N}\left(\vec{0}, \sigma ^2 \cdot \mathbf{I}_{n-1}\right) \]

Hence the components of \(\frac{1}{\sigma }\mathbf{z}_{L_2}^T\vec{X}\) are independent standard normal random variables, \(n-1\) of them in total. So

\[ n\cdot \frac{\widehat{\sigma ^2}}{\sigma ^2} = \left\| \frac{1}{\sigma }\mathbf{z}_{L_2}^T \vec{X}\right\| ^2 \sim \chi ^2_{n-1} \]

We have found the distribution of \(\widehat{\sigma^2}\)!

Next we need to find a pivotal statistic for the confidence interval for \(\mu\). In the part above we used \(\frac{\overline{X} - \mu }{\sqrt{\sigma^2/n}}\). In this case \(\sigma^2\) is unknown. Let’s replace \(\sigma^2\) with the estimate \(\widehat{\sigma^2}\):

\[ \frac{\overline{X} - \mu }{\sqrt{\widehat{\sigma ^2}/n}} \]

We need to find the distribution of this expression. Let’s transform it:

\[ \frac{\overline{X} - \mu }{\sqrt{\widehat{\sigma ^2}/n}} = \frac{(\overline{X} - \mu )/\sqrt{\sigma ^2 / n}}{\sqrt{\left(n \widehat{\sigma ^2}/\sigma ^2\right)/n}} = \frac{\xi }{\sqrt{\eta / n}} \]

where

\[ \xi := \frac{\overline{X} - \mu }{\sqrt{\sigma ^2/n}}, \qquad \eta := n\frac{\widehat{\sigma ^2}}{\sigma ^2} \]

Definition 2 Let \(\xi , \eta\) be independent random variables, distributed according to the laws \(\mathcal{N}(0,1), \chi_{n}^{2}\) respectively. Then the random variable
\[ \frac{\xi }{\sqrt{\eta / n}} \] has Student’s distribution with \(n\) degrees of freedom (denoted \(T_n\)).

In our case \(\xi \sim \mathscr {N}\left(0, 1\right)\) and \(\eta \sim \chi^2_{n-1}\) are independent (\(\xi\) is a function of \(\vec{X}_L\), \(\eta\) is a function of \(\vec{X}_{L^\perp }\); independence is justified by Cochran’s theorem below), but in the denominator we have \(\sqrt{\eta /n}\), not \(\sqrt{\eta /(n-1)}\). Let’s correct the normalization:

\[ \frac{\overline{X} - \mu }{\sqrt{\widehat{\sigma ^2}/(n-1)}} = \frac{\xi }{\sqrt{\eta /(n-1)}} \sim T_{n-1} \]

(equivalently, \(\frac{\overline{X} - \mu }{\sqrt{s^2/n}} \sim T_{n-1}\), where \(s^2 = \frac{n}{n-1}\widehat{\sigma^2}\) is the unbiased sample variance). This is a pivotal statistic, and the confidence interval for \(\mu\) at confidence level \(\alpha\) is obtained in the standard way:

\[ \mathbb {P}\left(\mu \in \left(\overline{X} - \sqrt{\tfrac {s^2}{n}}\, Q_{T_{n-1}}\left(\tfrac {1+\alpha }{2}\right), \; \; \overline{X} - \sqrt{\tfrac {s^2}{n}}\, Q_{T_{n-1}}\left(\tfrac {1-\alpha }{2}\right)\right)\right) = \alpha . \]

The interval for \(\sigma^2\) with unknown \(\mu\) is constructed as in part b), except now \(n\widehat{\sigma^2}/\sigma^2 \sim \chi^2_{n-1}\), i.e. with \(n-1\) degrees of freedom.

Densities of \(\chi^2_k\) and \(T_k\)

Left: densities of \(\chi^2_k\) for \(k=1,2,3,5,10\). Right: densities of \(T_k\) for \(k=1,3,10\) and the limiting \(\mathcal N(0,1)\) (dashed). The black/teal curve on each panel is the current \(k\) from the slider, recomputed continuously.

\(k=\)5.0 Log scale for the tails (\(T_k\))

Deviation of the current \(T_k\) from \(\mathcal N(0,1)\): — (max. \(|T_k(x)-\mathcal N(0,1)(x)|\) over the grid \(x\in[-4.6,4.6]\))

In general, the following theorem holds.

Theorem 2 (Cochran’s theorem on the orthogonal decomposition of a Gaussian vector) Let \(X_1, X_2, \ldots , X_n\) be independent standard normal random variables, i.e.

\[ \vec{X} \sim \mathscr {N}\left(\vec{0}, \mathbf{I}_n\right) \]

where \(\vec{X} = (X_1, \ldots , X_n)^T\). Let \(L_{1} \oplus \ldots \oplus L_{r}\) be an arbitrary decomposition of \(\mathbb {R}^{n}\) into a direct sum of orthogonal subspaces. For each \(j \in \{ 1, \ldots , r\}\) set \(\vec{X}_{L_j} := \operatorname {proj}_{L_{j}} \vec{X}\). Then

• \(\vec{X}_{L_1}, \ldots , \vec{X}_{L_r}\) are jointly independent random vectors;

• \(\left\| \vec{X}_{L_j}\right\|^2 \sim \chi^2_{\dim (L_j)}\) for each \(j \in \{ 1, \ldots , r\}\).

Let

\[ \vec{X} \sim \mathscr {N}\left(\vec{a}, \sigma ^2 \mathbf{I}_n\right), \]

and let \(L_{1} \oplus \ldots \oplus L_{r}\) be a decomposition of \(\mathbb {R}^{n}\) into a direct sum of orthogonal subspaces. For each \(j \in \{ 1, \ldots , r\}\) set \(\vec{X}_{L_j} := \operatorname {proj}_{L_{j}} \vec{X}\). Then \(\vec{X}_{L_1}, \ldots , \vec{X}_{L_r}\) are jointly independent, and moreover

\[ \frac{1}{\sigma ^2}\left\| \vec{X}_{L_j} - \mathbb {E}\left[\vec{X}_{L_j}\right]\right\| ^2 \sim \chi ^2_{\dim (L_j)} \]

ProblemProblem 7

Let \(\vec{X} = \left(X_{1}, \ldots , X_{n}\right)^T\) be a sample from the normal distribution with parameters \(\left(a, \sigma^{2}\right)\). Prove that the statistics \(\bar{X}\) and
\[ S^{2}=\frac{1}{n} \sum _{i=1}^{n}\left(X_{i}- \overline{X}\right)^{2} \] are independent, and compute the distribution of the statistic \(n S^{2}\).

In the general case \(\vec{X} \sim \mathscr {N}\left(\mathbf{z}\vec{\theta }, \sigma^2 \mathbf{I}_n\right)\), we need to construct confidence intervals for the components of \(\vec{\theta }\) and for \(\sigma^2\). Let

\[ L = \begin{cases} \operatorname {span}\left(\vec{z}_{\cdot 1}, \vec{z}_{\cdot 2}, \ldots , \vec{z}_{\cdot k}\right), & \text{if the model has no intercept } \theta _0 \\ \operatorname {span}\left(\overrightarrow {1}, \vec{z}_{\cdot 1}, \vec{z}_{\cdot 2}, \ldots , \vec{z}_{\cdot k}\right), & \text{if the model has an intercept } \theta _0 \end{cases} \]

Then the vectors \(\vec{X}_{L}, \vec{X}_{L^\perp }\) are independent Gaussian random vectors.

• \(\widehat{\vec{\theta }}\) is linearly dependent on \(\vec{X}_{L}\),

• \(\widehat{\sigma^2}\) is a function of \(\vec{X}_{L^\perp }\)

Hence \(\widehat{\vec{\theta }}, \widehat{\sigma^2}\) are independent. We have already computed above the expectation and covariance matrix of \(\widehat{\vec{\theta }} = \left(\mathbf{z}^T\mathbf{z}\right)^{-1}\mathbf{z}^T \vec{X}\): \(\mathbb {E}\left[\widehat{\vec{\theta }}\right] = \vec{\theta }\), \(\operatorname {Var}\left[\widehat{\vec{\theta }}\right] = \sigma^2 (\mathbf{z}^T \mathbf{z})^{-1}\); in the Gaussian case this means

\[ \widehat{\vec{\theta }} \sim \mathscr {N}\left(\vec{\theta }, \sigma ^2 (\mathbf{z}^T \mathbf{z})^{-1}\right). \]

Further,

\[ \begin{align} \vec{X}_{L^\perp } & = \vec{X} - \mathbf{z}\widehat{\vec{\theta }}, \qquad \mathbb {E}\left[\vec{X}_{L^\perp }\right] = \mathbb {E}\left[\vec{X} - \mathbf{z} \widehat{\vec{\theta }}\right] = \mathbf{z}\vec{\theta } - \mathbf{z}\vec{\theta } = \vec{0}, \\ \frac{1}{\sigma ^2} \left\| \vec{X}_{L^\perp } - \mathbb {E}\left[\vec{X}_{L^\perp }\right]\right\| ^2 & = \frac{1}{\sigma ^2 } \left\| \vec{X}_{L^\perp }\right\| ^2 = \frac{n-\widetilde{k}}{\sigma ^2}\widehat{\sigma ^2} \sim \chi ^2_{n-\widetilde{k}}. \end{align} \]

Let \(\left(\mathbf{z}^T\mathbf{z}\right)^{-1} =: (a_{i,j})\). Then

\[ \frac{\widehat{\vec{\theta }} - \vec{\theta }}{\sqrt{\sigma ^2}} \sim \mathscr {N}\left(\overrightarrow {0}, (\mathbf{z}^T \mathbf{z})^{-1}\right), \qquad \frac{\widehat{\vec{\theta }}_i - \theta _i}{\sqrt{a_{i,i}\sigma ^2}} \sim \mathscr {N}\left(0, 1\right), \; \forall i \in \left\{ 1,\ldots ,\widetilde{k}\right\} \]

In total, \(\frac{\widehat{\vec{\theta }}_i - \theta_i}{\sqrt{a_{i,i}\sigma^2}}, \frac{n-\widetilde{k}}{\sigma^2}\widehat{\sigma^2}\) are 2 independent random variables, where the first is distributed as \(\mathscr {N}\left(0, 1\right)\) and the second as \(\chi^2_{n-\widetilde{k}}\). Then, by the definition of Student’s distribution, we have

\[ \begin{align} \frac{\frac{\widehat{\vec{\theta }}_i - \theta _i}{\sqrt{a_{i,i}\sigma ^2}}}{\sqrt{\frac{n-\widetilde{k}}{\sigma ^2}\widehat{\sigma ^2}} / \sqrt{n-\widetilde{k}}} = \frac{\widehat{\vec{\theta }}_i - \theta _i}{\sqrt{a_{i,i}\widehat{\sigma ^2}}} & \sim T_{n-\widetilde{k}} \end{align} \]

ExampleExample 6

In the Gaussian linear model, find confidence intervals for the parameters \(\theta_{i}, i \in \{ 1, \ldots , k\}\), and \(\sigma^{2}\) at confidence level \(\gamma\).

SolutionSolution

• Denote by \(a_{i, i}, i \in \{ 1, \ldots , \widetilde{k}\}\), the diagonal element of the matrix \(\left(\mathbf{z}^{T} \mathbf{z}\right)^{-1}\) in the \(i\)-th position. Then, by the theorem on the orthogonal decomposition, the random variable

  \[ \frac{\widehat{\vec{\theta }}_{i}-\theta _{i}}{\sqrt{a_{i, i} \widehat{\sigma ^{2}}}} \]

  has distribution \(T_{n-\widetilde{k}}\). Let \(Q_{T_{n-\tilde{k}}}(\cdot )\) be the quantile function of the distribution \(T_{n-\tilde{k}}\). Then

  \[ \mathbb {P}_{\vec{\theta }, \sigma ^{2}}\left(Q_{T_{n-\tilde{k}}}\left(\frac{1-\gamma }{2}\right)<\frac{\widehat{\vec{\theta }}_{i}-\theta _{i}}{\sqrt{a_{i, i} \widehat{\sigma ^{2}}}} < Q_{T_{n-\tilde{k}}}\left(\frac{1+\gamma }{2}\right) \right) =\gamma , \]

  from which we get that the confidence interval at confidence level \(\gamma\) for \(\theta_{i}\) equals

  \[ \left(\widehat{\vec{\theta }}_{i}-\sqrt{a_{i, i} \widehat{\sigma ^{2}}}Q_{T_{n-\tilde{k}}}\left(\frac{1+\gamma }{2}\right), \widehat{\vec{\theta }}_{i}-\sqrt{a_{i, i} \widehat{\sigma ^{2}}} Q_{T_{n-\tilde{k}}}\left(\frac{1-\gamma }{2}\right)\right) \]

• Further, let \(Q_{\chi_{n-\tilde{k}}^2}(\cdot )\) be the quantile function of the distribution \(\chi_{n-\tilde{k}}^2\). Then

  \[ \mathbb {P}_{\vec{\theta }, \sigma ^{2}}\left(\frac{(n-\widetilde{k}) \widehat{\sigma ^{2}}}{\sigma ^{2}}>Q_{\chi _{n-\tilde{k}}^2}(1 - \gamma )\right) = \gamma \]

  The confidence interval at confidence level \(\gamma\) for \(\sigma^{2}\) equals

  \[ \left(0, \frac{(n-\widetilde{k}) \widehat{\sigma ^{2}}}{Q_{\chi _{n-\tilde{k}}^2}(1 - \gamma )}\right). \]

ProblemProblem 8

Let \(X_{1}, \ldots , X_{n}\) be a sample from \(\mathscr {N}\left(a, \sigma^{2}\right)\) (both parameters unknown). Construct exact confidence intervals for each of the parameters \(a, \sigma^{2}\).

ProblemProblem 9

Let \(X_{i}, i \in \{ 1, \ldots , n\}\) be independent random variables, distributed according to the normal law with parameters \(\left(a+b i, \sigma^{2}\right), \; a,b \in \mathbb {R}, \; \sigma^2 \in \mathbb {R}_+\). Construct exact confidence intervals for the parameters \(a, b, \sigma^{2}\).

4 Regression with random regressors

The researcher does not always have the ability to control the regressors \(\mathbf{z}\). In general, they too can be treated as random:

\[ \vec{X} = \mathbf{Z}\vec{\theta } + \vec{\varepsilon }, \]

where \(\mathbf{Z}\) is a random element with values in the space \(\mathbb {R}^{n \times \widetilde{k}}\). The assumptions on the errors here are the same as in Section 1.3, but all expectations become conditional on \(\mathbf{Z}\):

1. \(\mathbb {E}\left[ \vec{\varepsilon } \mid \mathbf{Z} \right] = \overrightarrow {0} \in \mathbb {R}^{n}\);

2. \(\operatorname {Var}\left[ \vec{\varepsilon } \mid \mathbf{Z} \right] = \mathbb {E}\left[ \vec{\varepsilon }\vec{\varepsilon }^T \mid \mathbf{Z} \right] - \mathbb {E}\left[ \vec{\varepsilon } \mid \mathbf{Z} \right] \cdot \mathbb {E}\left[ \vec{\varepsilon }^T \mid \mathbf{Z} \right] = \sigma^2 \mathbf{I}_n\).

All the results obtained above remain valid if we replace expectations and distributions with conditional expectations and conditional distributions given \(\mathbf{Z}\). For example, in the case of Gaussian regression we now have \(\operatorname {Law}\left(\widehat{\vec{\theta }} \mid \mathbf{Z}\right) = \mathscr {N}\left(\vec{\theta }, \sigma^2 \left(\mathbf{Z}^T\mathbf{Z}\right)^{-1}\right)\).

The difference is, for example, that it is no longer so simple to obtain the unconditional distribution \(\operatorname {Law}\left(\widehat{\vec{\theta }}\right)\), since the distribution of \(\mathbf{Z}\) may be arbitrary. However, if the conditional distribution does not depend on \(\mathbf{Z}\), as, for example, in these cases:

\[ \operatorname {Law}\left(\frac{n-\widetilde{k}}{\sigma ^2}\widehat{\sigma ^2} \mid \mathbf{Z}\right)[\bigg] = \chi ^2_{n-\widetilde{k}}, \qquad \operatorname {Law}\left(\frac{\widehat{\vec{\theta }}_i - \theta _i}{\sqrt{a_{i,i}\widehat{\sigma ^2}}} \mid \mathbf{Z}\right)[\bigg] = T_{n-\widetilde{k}}, \]

then we can pass from the conditional distribution to the unconditional one:

\[ \operatorname {Law}\left(\frac{n-\widetilde{k}}{\sigma ^2}\widehat{\sigma ^2} \right) = \chi ^2_{n-\widetilde{k}}, \qquad \operatorname {Law}\left(\frac{\widehat{\vec{\theta }}_i - \theta _i}{\sqrt{a_{i,i}\widehat{\sigma ^2}}}\right) = T_{n-\widetilde{k}} \]