STAT 306 - Module 09
\[ \mathbf{Y} = \begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{bmatrix} \;, \qquad \mathbf{X} = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{bmatrix} \;, \\ \pmb{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{bmatrix} \;, \qquad \pmb{\varepsilon} = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix} \;, \]
with \(\pmb{\varepsilon}\sim N_n(\mathbf{0}, \sigma^2\mathbf{I}_n)\), i.e., the errors \(\varepsilon_i\) are i.i.d. \(N(0,\sigma^2)\).
\[ \mathbb{E}[\mathbf{Y}] = \begin{bmatrix} \mathbb{E}[Y_1] \\ \mathbb{E}[Y_2] \\ \vdots \\ \mathbb{E}[Y_n] \end{bmatrix}\; \]
\[ \mathrm{Cov}(\mathbf{Y}) = \begin{bmatrix} \mathrm{Var}(Y_1) & \mathrm{Cov}(Y_1,Y_2) & \cdots & \mathrm{Cov}(Y_1,Y_n) \\ \mathrm{Cov}(Y_2,Y_1) & \mathrm{Var}(Y_2) & \cdots & \mathrm{Cov}(Y_2,Y_n) \\ \vdots & \vdots & \ddots & \vdots \\ \mathrm{Cov}(Y_n,Y_1) & \mathrm{Cov}(Y_n,Y_2) & \cdots & \mathrm{Var}(Y_n) \end{bmatrix} \;. \]
A covariance matrix is always symmetric, i.e., \(\mathbf{\Sigma}=\mathbf{\Sigma}^\top\).
A covariance matrix is always positive semi-definite, i.e., for any vector \(\mathbf{a}\), \(\mathbf{a}^\top \mathbf{\Sigma} \mathbf{a} \geq 0\).
If the components of a random vector are independent, then its covariance matrix is diagonal.
If \(\mathbf{A}\) and \(\mathbf{B}\) are constant matrices and \(\mathbf{X}\), \(\mathbf{Y}\) are random vectors, then \[ \mathrm{Cov}(\mathbf{A}\mathbf{X}, \mathbf{B}\mathbf{Y}) = \mathbf{A}\text{Cov}(\mathbf{X}, \mathbf{Y})\mathbf{B}^\top \]
In particular, taking \(\mathbf{X}=\mathbf{Y}\) and \(\mathbf{B}=\mathbf{A}\), \[ \mathrm{Cov}(\mathbf{A}\mathbf{Y}) = \mathbf{A}\mathbf{\Sigma}_\mathbf{Y}\mathbf{A}^\top \]
If \(\mathbf{a}\) is a constant vector of the same dimension as \(\mathbf{Y}\), then \[ \mathrm{Cov}(\mathbf{a}+\mathbf{Y}) = \mathrm{Cov}(\mathbf{Y}) = \mathbf{\Sigma}_{\mathbf{Y}} \;. \]
\[ \mathbb{E}[\pmb{\varepsilon}] = \begin{bmatrix} \mathbb{E}[\varepsilon_1] \\ \mathbb{E}[\varepsilon_2] \\ \vdots \\ \mathbb{E}[\varepsilon_n] \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} = \mathbf{0}_n \; \]
\[ \text{Cov}(\pmb{\varepsilon}) = \mathbf{\Sigma}_{\pmb{\varepsilon}} = \begin{bmatrix} \mathrm{Var}(\varepsilon_1) & \mathrm{Cov}(\varepsilon_1,\varepsilon_2) & \cdots & \mathrm{Cov}(\varepsilon_1,\varepsilon_n) \\ \mathrm{Cov}(\varepsilon_2,\varepsilon_1) & \mathrm{Var}(\varepsilon_2) & \cdots & \mathrm{Cov}(\varepsilon_2,\varepsilon_n) \\ \vdots & \vdots & \ddots & \vdots \\ \mathrm{Cov}(\varepsilon_n,\varepsilon_1) & \mathrm{Cov}(\varepsilon_n,\varepsilon_2) & \cdots & \mathrm{Var}(\varepsilon_n) \end{bmatrix} = \sigma^2\mathbf{I}_n \;. \]
To carry out any inference on the parameters \(\pmb{\beta}\), we need to know the sampling distribution of the least squares estimator \(\widehat{\pmb{\beta}}\).
What is the mean vector of the least squares estimator \(\widehat{\pmb{\beta}}\) in a model with \(p\) covariates?
What is the mean vector of the least squares estimator \(\widehat{\pmb{\beta}}\) in a model with \(p\) covariates?
\[ \begin{align} \mathbf{E}[\widehat{\pmb{\beta}}] &= \mathbf{E}\left[(\mathbf{X^\top X})^{-1}\mathbf{X^\top}\mathbf{Y}\right] \\ &= (\mathbf{X^\top X})^{-1}\mathbf{X^\top}\mathbf{E}\left[\mathbf{Y}\right] \\ &= (\mathbf{X^\top X})^{-1}\mathbf{X^\top}(\mathbf{X}\pmb{\beta}) \\ &= (\mathbf{X^\top X})^{-1}\mathbf{X^\top}\mathbf{X}\pmb{\beta} \\ &= \pmb{\beta} \;. \end{align} \]
Which of the following is the covariance matrix of the least squares estimator \(\widehat{\pmb{\beta}}\) in a model with \(p\) covariates?
Which of the following is the covariance matrix of the least squares estimator \(\widehat{\pmb{\beta}}\) in a model with \(p\) covariates?
\[ \begin{align} \text{Cov}(\widehat{\pmb{\beta}}) &= \text{Cov}\left(\left(\mathbf{X^\top X}\right)^{-1}\mathbf{X^\top Y}\right)\\ &=\left(\mathbf{X^\top X}\right)^{-1}\mathbf{X^\top} \mathbf{Cov}(\mathbf{Y})\left[\left(\mathbf{X^\top X}\right)^{-1}\mathbf{X^\top} \right]^\top\\ &= \left(\mathbf{X^\top X}\right)^{-1}\mathbf{X^\top} \sigma^2\mathbf{I}_n \mathbf{X}\left(\mathbf{X^\top X}\right)^{-1}\\ &= \sigma^2\left(\mathbf{X^\top X}\right)^{-1}\mathbf{X^\top}\mathbf{X}\left(\mathbf{X^\top X}\right)^{-1}\\ &= \sigma^2\left(\mathbf{X^\top X}\right)^{-1} \;. \end{align} \]
\[ \text{Cov}(\widehat{\pmb{\beta}}) = \begin{bmatrix} \mathrm{Var}(\widehat{\beta}_0) & \mathrm{Cov}(\widehat{\beta}_0,\widehat{\beta}_1) & \cdots & \cdots & \mathrm{Cov}(\widehat{\beta}_0,\widehat{\beta}_p) \\ \mathrm{Cov}(\widehat{\beta}_1,\widehat{\beta}_0) & \mathrm{Var}(\widehat{\beta}_1) & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \mathrm{Var}(\widehat{\beta}_{p-1}) & \mathrm{Cov}(\widehat{\beta}_{p-1},\widehat{\beta}_p) \\ \mathrm{Cov}(\widehat{\beta}_p,\widehat{\beta}_0) & \cdots & \cdots & \mathrm{Cov}(\widehat{\beta}_p,\widehat{\beta}_{p-1}) & \mathrm{Var}(\widehat{\beta}_p) \end{bmatrix} \;. \]
\(\widehat{\pmb{\beta}} = \left(\mathbf{X^\top X}\right)^{-1}\mathbf{X^\top}\mathbf{Y}\) is a linear transformation of \(\mathbf{Y}\sim N_n(\mathbf{X}\pmb{\beta}, \sigma^2\mathbf{I}_n)\), so it is also (multivariate) normal.
Therefore, \[ \widehat{\pmb{\beta}} \sim N_{p+1}\left(\pmb{\beta},\; \sigma^2\left(\mathbf{X^\top X}\right)^{-1}\right) \;. \]
In particular, \(\widehat{\beta}_j \sim N\left(\beta_j,\; \sigma^2\left[\left(\mathbf{X^\top X}\right)^{-1}\right]_{jj}\right)\), \(j=0,1,\dotsc,p\).
We have \(\text{Cov}(\widehat{\pmb{\beta}})=\mathbf{\Sigma}_{\widehat{\pmb{\beta}}} = \sigma^2(\mathbf{X^\top X})^{-1}\), which depends on the unknown parameter \(\sigma^2\).
To estimate this covariance matrix, we estimate \(\sigma^2\) with
\[ \hat\sigma^2 = \frac{RSS}{n-(p+1)} = \frac{\mathbf{e^\top e}}{n-(p+1)}\;. \]
\(\hat\sigma^2 = \frac{RSS}{n-(p+1)}\) is frequently referred to as residual mean square and commonly denoted as \(MSE\).
Hence, \(\widehat{\text{Cov}}(\widehat{\pmb{\beta}}) = \widehat{\mathbf{\Sigma}}_{\widehat{\pmb{\beta}}} = \widehat{\sigma}^2(\mathbf{X^\top X})^{-1}\)
Why divide by \(n-(p+1)\)? Same idea as the \(n-2\) in SLR: the residuals must satisfy \(p+1\) constraints (\(\mathbf{X^\top e} = \mathbf{0}\), the normal equations), so only \(n-(p+1)\) of them are “free”.
In Module 10, we will prove that \(\hat\sigma^2\) is unbiased for \(\sigma^2\) and see why replacing \(\sigma^2\) by \(\hat\sigma^2\) leads to a \(t\) distribution.
Suppose we fit a model with \(p=3\) covariates and obtain the estimated covariance matrix of \(\widehat{\pmb{\beta}}\) below:
\[ \widehat{\mathbf{\Sigma}}_{\widehat{\pmb{\beta}}} = \begin{bmatrix} 4 & 2 & 2 & 0 \\ 2 & 5 & 1 & -2 \\ 2 & 1 & 6 & 1 \\ 0 & -2 & 1 & 6 \end{bmatrix} \]
Based on this matrix, what is \(\widehat{\text{Cov}}(\widehat{\beta}_1,\widehat{\beta}_3)\)?
Now, we can create confidence intervals and hypothesis tests for individual parameters \(\beta_j\), \(j=0,1,\dotsc,p\), based on \[ \frac{\widehat{\beta}_j - \beta_j}{\sqrt{\left(\widehat{\mathbf{\Sigma}}_{\widehat{\pmb{\beta}}}\right)_{jj}}} \sim t_{n-(p+1)} \;, \]
Therefore, a \(100(1-\alpha)\%\) confidence interval for \(\beta_j\) is
\[ \left[\widehat{\beta}_j - t_{n-(p+1),\alpha/2} \cdot \sqrt{\left(\widehat{\mathbf{\Sigma}}_{\widehat{\pmb{\beta}}}\right)_{jj}}, \widehat{\beta}_j + t_{n-(p+1),\alpha/2} \cdot \sqrt{\left(\widehat{\mathbf{\Sigma}}_{\widehat{\pmb{\beta}}}\right)_{jj}}\right] \]
Back to penguins! Let’s fit a linear model with response body_mass_g and covariates flipper_length_mm, bill_length_mm, and bill_depth_mm.
I have already computed the least squares estimates for you \(\widehat{\pmb{\beta}}\)
Manually compute the following quantities:
Tip
res: recall the residual vector in matrix form is \(\mathbf{e} = \mathbf{y} - \mathbf{X}\widehat{\pmb{\beta}}\).Sigma_beta: recall \(\widehat{\mathbf{\Sigma}}_{\widehat{\pmb{\beta}}} = \widehat{\sigma}^2(\mathbf{X^\top X})^{-1}\), and \(\widehat{\sigma}^2 = \dfrac{\mathbf{e^\top e}}{n-(p+1)}\). You’ll need nrow(X) and ncol(X) to get \(n\) and \(p+1\).CI_bill_depth: a \(100(1-\alpha)\%\) CI for \(\beta_j\) is \(\widehat{\beta}_j \pm t_{n-(p+1),\alpha/2}\sqrt{\left(\widehat{\mathbf{\Sigma}}_{\widehat{\pmb{\beta}}}\right)_{jj}}\). Check the order of the columns in X to figure out which entry of beta_hat (and which diagonal entry of Sigma_beta) corresponds to bill_depth_mm.Solution:
#################
# Design Matrix #
#################
X <- cbind(1, flipper_length_mm, bill_length_mm, bill_depth_mm)
############
# Response #
############
y <- body_mass_g
###########################
# Least Squares Estimates #
###########################
beta_hat <- solve(t(X)%*%X,t(X)%*%y)
#############
# Residuals #
#############
res <- y - X %*% beta_hat
########################################
# Sample Covariance Matrix of beta-hat #
########################################
Sigma_beta <- solve(t(X)%*%X)*sum(res**2/(nrow(X)-ncol(X)))
#############################
# CI for beta_bill_depth_mm #
#############################
(CI_bill_depth <- beta_hat[4] + c(-1,1)*qt(0.975, nrow(X)-ncol(X))*sqrt(Sigma_beta[4,4]) )
We fit the model
body_mass_g ~ flipper_length_mm + bill_length_mm + bill_depth_mm
We obtain \(\text{CI}\left(\beta_{\text{bill_depth_mm}}, 0.95\right) = (-9.36,\, 45.03)\).
If we test \[H_0: \beta_{\text{bill_depth_mm}} = 0\quad \text{ vs } \quad H_1: \beta_{\text{bill_depth_mm}} \neq 0\] at \(\alpha=0.05\), would we reject or not reject \(H_0\)?
We fit the model
body_mass_g ~ flipper_length_mm + bill_length_mm + bill_depth_mm
We obtain \(\text{CI}\left(\beta_{\text{bill_depth_mm}}, 0.95\right) = (-9.36,\, 45.03)\).
If we test \[H_0: \beta_{\text{bill_depth_mm}} = 0\quad \text{ vs } \quad H_1: \beta_{\text{bill_depth_mm}} \neq 0\] at \(\alpha=0.05\), would we reject or not reject \(H_0\)?
Let \(\mathbf{x}_{new}=\begin{bmatrix} 1 & x_{new,1} & x_{new,2} & \cdots & x_{new,p} \end{bmatrix}^\top\).
The expected response at \(\mathbf{x}_{new}\) is \[ \mathbb{E}[Y_{new}] = \mathbf{x}_{new}^\top \pmb{\beta} \;. \]
An estimator of \(\mathbb{E}[Y_{new}]\) is \[ \widehat{\mathbb{E}}[Y_{new}] = \mathbf{x}_{new}^\top \widehat{\pmb{\beta}} \;. \]
\[ \begin{align} \mathbb{E}[\widehat{\mathbb{E}}[Y_{new}]] &= \mathbb{E}[\mathbf{x}_{new}^\top \widehat{\pmb{\beta}}] \\ &= \mathbf{x}_{new}^\top \mathbb{E}[\widehat{\pmb{\beta}}] \\ &= \mathbf{x}_{new}^\top \pmb{\beta} \\ &= \mathbb{E}[Y_{new}] \;. \end{align} \]
\[ \begin{align} \text{Var}[\mathbf{x}_{new}^\top \widehat{\pmb{\beta}}] &= \mathbf{x}_{new}^\top \text{Cov}[\widehat{\pmb{\beta}}] \mathbf{x}_{new} \\ &= \mathbf{x}_{new}^\top \sigma^2(\mathbf{X^\top X})^{-1} \mathbf{x}_{new} \\ &= \sigma^2 \mathbf{x}_{new}^\top (\mathbf{X^\top X})^{-1} \mathbf{x}_{new} \;. \end{align} \]
\[ CI(\mathbb{E}[Y_{new}], 1-\alpha) = \mathbf{x}_{new}^\top \widehat{\pmb{\beta}} \pm t_{n-(p+1),\alpha/2} \times \sqrt{\widehat{\sigma}^2 \mathbf{x}_{new}^\top (\mathbf{X^\top X})^{-1} \mathbf{x}_{new}} \]
Tip
x_new is already built for you: its entries follow the same order as the columns of X (intercept first, then the covariates in the order used to build X).beta_hat, X, res) is already available from the setup chunk.nrow(X) - ncol(X).Solution:
x_new <- matrix(c(1, 200, 40, 18), ncol = 1)
df <- nrow(X) - ncol(X)
(CI_expected_body_mass <- c(
t(x_new) %*% beta_hat - qt(0.975, df)*sqrt(t(x_new) %*% solve(t(X) %*% X) %*% x_new * sum(res**2)/(df)),
t(x_new) %*% beta_hat + qt(0.975, df)*sqrt(t(x_new) %*% solve(t(X) %*% X) %*% x_new * sum(res**2)/(df))
))Now derive the prediction interval for a new observation \(Y_{new}\) at \(\mathbf{x}_{new}\).
The prediction error is \[ \widehat{Y}_{new}-Y_{new} = \mathbf{x}_{new}^\top \widehat{\pmb{\beta}} - Y_{new} \;. \]
What is \(\text{Var}\left[\widehat{Y}_{new}-Y_{new}\right]\), the variance of the prediction error for a new observation at \(\mathbf{x}_{new}\)?
What is \(\text{Var}\left[\widehat{Y}_{new}-Y_{new}\right]\), the variance of the prediction error for a new observation at \(\mathbf{x}_{new}\)?
The prediction error is \[ \widehat{Y}_{new}-Y_{new} = \mathbf{x}_{new}^\top \widehat{\pmb{\beta}} - Y_{new} \;. \]
\(Y_{new}\) is a new observation, so it is independent of \(\widehat{\pmb{\beta}}\) (which only depends on the training data). Hence \(\text{Cov}\left(\mathbf{x}_{new}^\top \widehat{\pmb{\beta}}, Y_{new}\right) = 0\).
The variance of the prediction error is
\[ \begin{align} \text{Var}[\widehat{Y}_{new}-Y_{new}] &= \text{Var}[\mathbf{x}_{new}^\top \widehat{\pmb{\beta}} - Y_{new}] \\ &= \text{Var}\left(\mathbf{x}_{new}^\top \widehat{\pmb{\beta}}\right) + Var[Y_{new}] \qquad (\text{by independence}) \\ &= \mathbf{x}_{new}^\top \text{Cov}\left(\widehat{\pmb{\beta}}\right)\mathbf{x}_{new} + Var[Y_{new}] \\ &= \sigma^2 \mathbf{x}_{new}^\top (\mathbf{X^\top X})^{-1} \mathbf{x}_{new} + Var[Y_{new}] \\ &= \sigma^2 \mathbf{x}_{new}^\top (\mathbf{X^\top X})^{-1} \mathbf{x}_{new} + \sigma^2 \;. \end{align} \]
\[ PI(Y_{new}|\mathbf{x}_{new}, 1-\alpha) = \mathbf{x}_{new}^\top \widehat{\pmb{\beta}} \pm t_{n-(p+1),\alpha/2} \times \sqrt{\widehat{\sigma}^2 \left(1 + \mathbf{x}_{new}^\top (\mathbf{X^\top X})^{-1} \mathbf{x}_{new}\right)} \]
Tip
x_new is built the same way as in Activity 4: a column vector matching the column order of X (intercept first).Solution:
x_new <- matrix(c(1, 200, 40, 18), ncol = 1)
df <- nrow(X) - ncol(X)
(prediction_interval_body_mass <- c(
t(x_new) %*% beta_hat - qt(0.975, df)*sqrt((t(x_new) %*% solve(t(X) %*% X) %*% x_new + 1)* sum(res**2)/(df)),
t(x_new) %*% beta_hat + qt(0.975, df)*sqrt((t(x_new) %*% solve(t(X) %*% X) %*% x_new + 1)* sum(res**2)/(df))
))© 2025 Rodolfo Lourenzutti – Material Licensed under CC By-SA 4.0