Ordinary least squares estimation explained in "Optimal Design of Experiments" by Goos/Jones, questions about some equationsEq I am currently trying to become familiar with design of experiments with the book “Optimal design of experiments” by Goos and Jones. In chapter 2, they discuss the use of a Plackett–Burman type design (6 factors, 12 experiments), and in chapter 2.3 (“Peek into the black box”), they introduce the math. I am currently struggling with chapter 2.3.4 (“Ordinary least squares estimates”), more specifically with some of the equations.
Quote:

The ordinary least squares estimator of the vector of unknown model >coefficients $\beta $ is
$$\hat{\beta }=({X}'X)^{-1}{X}'Y.$$

Ok, I understand that, assuming X' is the transposed of X. Then they continue

The variance-covariance matrix of this estimator is
$$var(\hat{\beta })=\sigma _{\varepsilon } ^{2}({X}'X)^{-1}.$$

Here, I am getting confused. I thought that $\beta$ (the "real" values) and $\hat{\beta}$ (the calculated values) are vectors, and that the elements $\beta _i$ are scalars. However, the equation and the wording suggest that $\hat{\beta}$ is a matrix, and that the elements $\hat{\beta}_i$ are vectors. Moreover, how is it possible to calculate anything related to $\beta$ without the regressand vector Y. What am I missing here?
Further down in the text, they write

Note that the variance-covariance matrix of the estimator is directly
proportional to the error variance, which is unknown. We can estimate
the error variance using the mean squared error: $$\hat{\sigma}
> _{\varepsilon } ^{2} = \frac{1}{n-p}(Y- X \hat{\beta} )' (Y- X \hat{\beta} ).$$

I can understand that equation. At least, it makes sense that the difference between Y and $X\hat{\beta }$ (or $\hat{Y}$) is proportional to $\hat{\sigma} _{\varepsilon }$.
But I would really like to understand the other equations as well.
Regards, Soltub
 A: I think your main point of confusion might just be a misunderstanding of notation -- you are indeed correct that $\hat{\beta}$ should be a vector. However, $Var(\hat{\beta})$ is a matrix called (as you have noted above) the variance-covariance (or just covariance) matrix. More specifically, for parameter vector $\hat{\beta} \in \mathbb{R}^n$, $Var(\hat{\beta}) \in \mathbb{R}^{n \times n}$.
Why is it a matrix instead of a vector?
Well, there's nothing really stopping you from defining some quantity of $\hat{\beta}$ that measures element-wise variance of each element of the vector, but it's a more informative measure of how a random vector varies, since variances in one element may correlate to variances in another element.
In particular, the OLS estimator has sampling distribution $\mathcal{N}(\beta, \sigma_\varepsilon^2 I)$ -- a multivariate Gaussian. At risk of some hand-waviness -- you can't parameterize a multivariate Gaussian with only the element-wise variances of the vector elements; you need to take into account the dependencies between elements of the vector. To explore this further, you can try playing with the covariance matrix of a 2-D multivariate Gaussian, and observe how the shape changes.
Derivation
There's a lot of possible derivations for the variance of the OLS estimator; here's one I like. I'll try to keep the same notation you used.
We have our usual setup: $$Y = X\beta + \varepsilon, \;\varepsilon\sim \mathcal{N}(0, \sigma_\varepsilon^2 I)$$
By definition:
$$Var(\hat{\beta}) = \mathbb{E}[(\hat{\beta} - \beta)(\hat{\beta} - \beta)'].$$
Let's find $\hat{\beta} - \beta$ as an intermediate step. Recall that
$$\hat{\beta} = (X'X)^{-1}X' Y,$$
so, substituting, we have
$$\hat{\beta} = (X'X)^{-1}X' (X\beta + \varepsilon)$$
$$ = \beta + (X'X)^{-1}X'\varepsilon$$
so that $$\hat{\beta} - \beta = (X'X)^{-1}X'\varepsilon.$$
Then, substituting into the original expression for variance:
$$Var(\hat{\beta}) = \mathbb{E}[(X'X)^{-1}X'\varepsilon\varepsilon'X(X'X)^{-1}]$$
$$=(X'X)^{-1}X'\mathbb{E}[\varepsilon\varepsilon']X(X'X)^{-1}$$.
Note that the $\mathbb{E}[\varepsilon\varepsilon']$ term is simply the covariance of random vector $\varepsilon$, which is given to us as $\sigma_\varepsilon^2 I$. We can substitute, multiply out the $I$, and rearrange constants:
$$=(X'X)^{-1}X'\sigma_\varepsilon^2 I X(X'X)^{-1}$$
$$=\sigma_\varepsilon^2 \cdot (X'X)^{-1}X'X(X'X)^{-1}$$
$$=\sigma_\varepsilon^2 (X'X)^{-1}$$
as needed. $\blacksquare$
I skipped a step where I sneakily removed the transpose from $(X'X)^{-1}$; since $X'X$ is symmetric, $(X'X)^{-1}$ is as well, which is proven here.
