Standard Error and confidence interval for multiple linear regression in matrix notation

Lets say we are trying to fit a normal linear model to our data as follows:

$$y = \beta _0 + x_1\beta_1 + x_2\beta_2 + ... + x_p\beta_p + \epsilon$$.

My question is that how we can derive the standard error and confidence interval for each of the coefficient in matrix notation.

Under the usual Normal Gauss-Markov Model where $$\boldsymbol{Y}=\boldsymbol{X\beta}+\boldsymbol{\epsilon}$$, it is assummed $$\boldsymbol{\epsilon}\sim N\left(0,\boldsymbol{I}\sigma^{2}\right)$$, where $$\boldsymbol{I}$$ is an $$n\times n$$ identity matrix, where $$n$$ is the number of observations in your dataset. This implies $$Var(\boldsymbol{Y})=Var\left(\boldsymbol{X\beta}+\boldsymbol{\epsilon}\right)=Var(\boldsymbol{\epsilon})=\boldsymbol{I}\sigma^{2}$$.Then, under this model:

$$\begin{eqnarray*} Var\left(\hat{\boldsymbol{\beta}}\right) & = & Var\left[\boldsymbol{\left(X^{\prime}X\right)^{-}X^{\prime}Y}\right]\\ & = & Var\left[\boldsymbol{AY}\right]\\ & = & \boldsymbol{A}Var(\boldsymbol{Y})\boldsymbol{A^{\prime}}\\ & = & \boldsymbol{A}\left(\boldsymbol{I}\sigma^{2}\right)\boldsymbol{A}^{\prime}\\ & = & \left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]\left(\boldsymbol{I}\sigma^{2}\right)\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]^{\prime}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\left[\boldsymbol{X}^{\prime}\boldsymbol{X}\right]\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\\ & = & \sigma^{2}\boldsymbol{I}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\\ & = & \sigma^{2}\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-} \end{eqnarray*}$$

Normally we don't know $$\sigma^{2}$$ so it has to estimated from the residuals, $$\boldsymbol{\hat{\epsilon}}$$ to give us: $$Var\left(\boldsymbol{\hat{\beta}}\right)=\hat{\sigma}^{2}\left(X^{\prime}X\right)^{-}$$. Without getting into all the gory details, the estimated $$\sigma^2$$ is given by: $$\hat{\sigma}^{2}=\frac{\hat{\epsilon}^{\prime}\hat{\epsilon}}{n-r}$$ where $$n$$ is the number of rows in $$\boldsymbol{X}$$ and $$r$$ is the rank of $$\boldsymbol{X}$$.

Then, again, without going into all the distribution derivations and other gory details, the $$100(1-\alpha)\%$$ confidence interval for $$\hat{\beta}_{j}$$, $$j=0,1,\ldots,k,$$ ($$k=r-1$$) is given by:

$$\begin{eqnarray*} \hat{\beta}_{j} & \pm & t_{\alpha/2,n-r}\hat{\sigma}\sqrt{g_{jj}} \end{eqnarray*}$$

where $$g_{jj}$$ is the $$j$$th diagonal element of $$\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}$$, and $$t_{\alpha/2,n-r}$$ is the $$\left(\alpha/2\right)100$$ percentile of the $$t$$-distribution, with $$n-r$$ degrees of freedom.

Here is an example of the computations done manually and then checked against confidence intervals computed by R:

> library(MASS)
>
> ##create some data
> data("mtcars")
> y<-mtcars$$mpg > y [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 10.4 10.4 14.7 32.4 30.4 33.9 21.5 15.5 15.2 13.3 19.2 27.3 [27] 26.0 30.4 15.8 19.7 15.0 21.4 > x0<-rep(1, dim(mtcars)[1]) > x1<-mtcars$$disp
> x2<-mtcars$$hp > x3<-mtcars$$drat
> X<-cbind(x0, x1, x2, x3)
>
> ##print first rows of design matrix X
x0  x1  x2   x3
[1,]  1 160 110 3.90
[2,]  1 160 110 3.90
[3,]  1 108  93 3.85
[4,]  1 258 110 3.08
[5,]  1 360 175 3.15
[6,]  1 225 105 2.76
>
> ##calculate beta hats - note:  not efficient, just for teaching
> beta.hat<-ginv(t(X)%*%X)%*%t(X)%*%y
> beta.hat
[,1]
[1,] 19.34429256
[2,] -0.01923223
[3,] -0.03122932
[4,]  2.71497521
> ##get predicted values
> y.hat <- X%*%beta.hat
> residuals <- y-y.hat
> ##find n and rank, r
> r <- qr(X)$$rank > n <- dim(X)[1] > ##compute covariance matrix > var.beta.hat <- as.numeric(t(residuals)%*%residuals/(n - r))*ginv(t(X)%*%X) > ##compute confidence interval > alpha <- .05 > t <-qt(1 - alpha/2, n - r) > conf.int <- round(cbind(beta.hat - t*sqrt(diag(var.beta.hat)), beta.hat + t*sqrt(diag(var.beta.hat))), 5) > conf.int [,1] [,2] [1,] 6.29413 32.39445 [2,] -0.03843 -0.00004 [3,] -0.05857 -0.00389 [4,] -0.33176 5.76171 > > #fit model with standard R lm package > r.fit <-lm(y ~ x1 + x2 + x3) > #show confidence intervals in R > round(confint(r.fit, names(r.fit$$coef), level = 0.95), 5)
2.5 %   97.5 %
(Intercept)  6.29413 32.39445
x1          -0.03843 -0.00004
x2          -0.05857 -0.00389
x3          -0.33176  5.76171


@prony, so we are trying to calculate the standard error, so I'm not sure what you mean by the standard error is defined as $$\sigma/\sqrt{n}.$$ But I think I understand your question. To understand where $$\sqrt{g_{jj}}$$ is coming from, first, notice that the standard error is just the square root of the variance of an estimator -- in this case, $$\hat{\boldsymbol{\beta}}$$. So really, all we are trying to do is obtain the square root of each of the diagonal elements of $$Var(\hat{\boldsymbol{\beta}}) = \hat{\sigma}^2(\boldsymbol{X}^\prime \boldsymbol{X})^-$$. Here, $$\hat{\boldsymbol{\beta}}$$ is a vector. So to obtain the standard error for a single $$\hat{\boldsymbol{\beta}}$$, (what I call $$\hat{\beta_j}$$ and note this is not longer bolded as it's just a single element), we need to pick off the element corresponding to the $$j$$th row and $$j$$th column of the $$\hat{\sigma}^2(\boldsymbol{X}^\prime \boldsymbol{X})^-$$ matrix, and take the square root of it (this amounts to taking the square root of the diagonal entries of the matrix, which are the variances before square-rooting. The off-diagonal entries are the covariances). I have simply defined $$\sqrt{g_{jj}}$$ to be equal to the $$j$$th diagonal element of the square matrix $$(\boldsymbol{X}^\prime \boldsymbol{X})^-$$ so that when we obtain the square root of the diagonal element of $$\hat{\sigma}^2(\boldsymbol{X}^\prime \boldsymbol{X})^-$$ to obtain the standard error, we get $$\sqrt{\hat{\sigma}^2g_{jj}}=\hat{\sigma}\sqrt{g_{jj}}.$$

• Thanks for the answer. Actually, I was mostly interested in the gory details since I wanted to get the derivation. I found the derivation for estimation $\hat{\sigma}$. However, i could not understand why standar error of the estimator of a coefficient is $\hat{\sigma}\sqrt(g_{jj})$. – prony Feb 5 '19 at 16:49
• No problem. Do you mean you are confused how we get $\hat{\sigma}\sqrt{g_{jj}}$ as the standard error of a coefficient or did you mean to say that you are wanting to understand how it is that we determine $\frac{\hat{\beta}_{j}}{\hat{\sigma}\sqrt{g_{jj}}}$ is distributed as a central $t$ distribution? – StatsStudent Feb 5 '19 at 20:58
• Sorry for the ambiguity of my question. I meant the first one. Because, to the best of my knowledge (forgive my ignorance), the standard error is defined as $\frac{\sigma}{\sqrt{n}}$. I understand here we use estimator for the sigma but I don't understand where the $\sqrt{g_{jj}}$ is coming from. Btw, the reason why $\frac{\hat{\beta}}{\hat{\sigma}\sqrt{g_{jj}}}$ follows t distribution is that $\hat{\beta}$ is normally distributed and $\hat{\sigma}$ follow chi-square distribution right? – prony Feb 6 '19 at 4:08
• @prony, I just added "Update 1" at the end of the question to address your additional question about the derivation of the standard error. Is that clear now? – StatsStudent Feb 6 '19 at 5:51
• Also, @prony, the distribution is $t$ because (1) $\hat{\boldsymbol{\beta}}\sim N_r[\boldsymbol{\beta},\sigma^2(\boldsymbol{X}^\prime \boldsymbol{X})^-]$, (2) $(n-r)\hat{\sigma}^2/\sigma^2\sim \chi^2_{n-r}$ and (3) $\hat{\boldsymbol{\beta}}$ and $\hat{\sigma}^2$ are independent. – StatsStudent Feb 6 '19 at 6:28

Estimate the covariance matrix of the coefficients via

$$\widehat{\bf{\Sigma}} = \hat{\sigma}^2(X'X)^{-1}$$

The diagonal elements are the variances

$$\operatorname{diag}(\widehat{\Sigma)} = \left[ \hat{\sigma}^2_0, \dots, \hat{\sigma}_p^2 \right]$$

Here,

$$\hat{\sigma}^2_j = \hat{\sigma}^2 (X'X)^{-1}_{jj}$$

Though it isn't typical to write the square root of a vector, let $$\sqrt{ \operatorname{diag}(\widehat{\Sigma)}}$$ be the diagonal of $$\widehat{\Sigma}$$ where each component has been square rooted.

I suppose the confidence intervals are then

$$\beta\pm t_{1-\alpha/2, n-q-1}\sqrt{ \operatorname{diag}(\widehat{\Sigma)}}$$

Let's try it in R

#model
set.seed(0)
x = rnorm(9)
y = 2*x+1 +rnorm(9)
model = lm(y~x)

sigma = vcov(model)
tr.sigma = diag(sigma)

betas = coef(model)

#right
betas + qt(0.975,model$df.residual)*sqrt(tr.sigma) Intercept) x 1.859659 3.059125 #left betas - qt(0.975,model$df.residual)*sqrt(tr.sigma)

(Intercept)           x
0.01696051  1.17202056

confint(model)
2.5 %   97.5 %
(Intercept) 0.01696051 1.859659
x           1.17202056 3.059125


EDIT: Faraway examines this very question on pg. 38 of Linear Models With R. Link here

• Can you explain why the trace of the square root of the matrix equals to the standard error? – prony Feb 5 '19 at 16:51
• How is the trace a vector as you have it defined? The trace is a scalar since it's the sum of diagonal elements. – StatsStudent Feb 6 '19 at 0:16
• @StatsStudent Its been a long time since I've done linear algebra. I think I just meant $\operatorname{diag}$ – Demetri Pananos Feb 6 '19 at 2:27
• @StatsStudent I've made some edits, which should now agree with the R code I've presented. – Demetri Pananos Feb 6 '19 at 2:29
• That's looking better, but while your ultimate confidence interval formula is correct, the RHS of your penultimate equation is not: the $diag(\boldsymbol{\Sigma}) \ne \left[ \sigma^2_0, \dots, \sigma_p^2 \right]$. How did your $(X'X)^{-1}$ term vanish? – StatsStudent Feb 6 '19 at 5:00