2
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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
> head(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

Update 1. to Assist With Your Additional Comment/Question About the Standard Error

@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}}.$

$\endgroup$
  • $\begingroup$ 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})$. $\endgroup$ – prony Feb 5 '19 at 16:49
  • $\begingroup$ 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? $\endgroup$ – StatsStudent Feb 5 '19 at 20:58
  • $\begingroup$ 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? $\endgroup$ – prony Feb 6 '19 at 4:08
  • $\begingroup$ @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? $\endgroup$ – StatsStudent Feb 6 '19 at 5:51
  • $\begingroup$ 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. $\endgroup$ – StatsStudent Feb 6 '19 at 6:28
2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Can you explain why the trace of the square root of the matrix equals to the standard error? $\endgroup$ – prony Feb 5 '19 at 16:51
  • $\begingroup$ How is the trace a vector as you have it defined? The trace is a scalar since it's the sum of diagonal elements. $\endgroup$ – StatsStudent Feb 6 '19 at 0:16
  • $\begingroup$ @StatsStudent Its been a long time since I've done linear algebra. I think I just meant $\operatorname{diag}$ $\endgroup$ – Demetri Pananos Feb 6 '19 at 2:27
  • $\begingroup$ @StatsStudent I've made some edits, which should now agree with the R code I've presented. $\endgroup$ – Demetri Pananos Feb 6 '19 at 2:29
  • $\begingroup$ 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? $\endgroup$ – StatsStudent Feb 6 '19 at 5:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.