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.
 A: 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
