Why is the expectation of a random vector still a vector? In my previous post on derivation of covariance between y and random effect, for the following linear model:

Frank's anwer proved cov(y, u) = ZD as below:

Frank's proof involved E[u]. As far as I could understand, the random effect variable u in linear mixed model is a vector of 1*q, where q is the number of covariates in the model. My question is why the E[u] is not a single number?
What I do not understand:
If I have 3 random covariates in the model, then u MUST BE a vector of random variable of size 3*1. An example of u is as follows:

then what is E[u]? I think u is only a column of numbers, so E[u] should be the average of c(1,2,3), which is 2. However, frank says the expectation of a random vector is still a vector. I know he is right but I do not see why. I want to see the answer E[u] for my example so that I could really see what is E[u], including its dimension and how each number in E[u] is derived.
 A: It seems to me that you're overthinking things :). For any vector or matrix of any configuration, its expectation is simply the expectation of each of its elements. If we have:
$$
\mathbb{E}[\mathbf{A}] = 
\begin{pmatrix}
x_1 & y_1 \\
x_2 & y_2
\end{pmatrix}
$$
Then to calculate its expectation:
$$
\mathbb{E}[\mathbf{A}] = 
\mathbb{E}\bigg[\begin{pmatrix}
x_1 & y_1 \\
x_2 & y_2
\end{pmatrix}\bigg] = 
\begin{pmatrix}
\mathbb{E}[x_1] & \mathbb{E}[y_1] \\
\mathbb{E}[x_2] & \mathbb{E}[y_2]
\end{pmatrix}
$$
So now we see why that the expectation of any vector is going to some other vector of the same shape, because $\mathbb{E}[\mathbf{x}]$ has one element $\mathbb{E}[x_i]$ for each element $x_i$ of $\mathbf{x}$ in the same place.
For your specific question, you ask about $1\times 3$ vectors, but just so you know, the vector you have posted:
$$
\mathbf{u} = \begin{pmatrix}
1 \\
2 \\
3
\end{pmatrix}
$$
would be denoted in the common notation as a matrix of size $3\times 1$, or would be called a column vector of dimension 3.
This is a constant matrix, so its not random, so if we took its expectation we would get the same thing back $\mathbb{E}[\mathbf u]=\mathbf{u}$. But if we had a vector $\mathbf{v}$ with random variables $x_1,x_2,x_3$ in it:
$$
\mathbb{E}[\mathbf{v}] = 
\mathbf{E}\bigg[\begin{pmatrix}
x_1\\x_2\\x_3\end{pmatrix}\bigg] = 
\begin{pmatrix}
\mathbb{E}[x_1]\\\mathbb{E}[x_2]\\\mathbb{E}[x_3]\end{pmatrix}
$$
we would get an estimate for the population mean of $\mathbf{v}$ (the population mean and sample mean are vectors of the same shape, otherwise the sample mean would make a pretty darn bad estimator!).
Oh, right, also, so far we've been talking about the population quantity: the expectation. You're question title mentions moments but you specifically deal with "sample moments". The sample mean vector is simply the vector of means for each variable:
$$\bar{\mathbf{x}} = \begin{pmatrix}
\bar{x_1} \\ \bar{x_2} \\ \bar{x_3} \end{pmatrix}
$$
You mention covariances too, I couldn't really tell, do you have a question about that or was it just to give background?
A: Consider $X \in \mathbb{R}^m$, i.e., $X = \left(\begin{array}{c} X_1 \\ \vdots \\ X_m \end{array}\right)$. Then, we can see that
\begin{align*}
    \mathbb{E}(X) & = \int_{\mathbb{R}^{m}} X f_{X}(X) d X\\
    & = \int_{\mathbb{R}^{m}} \left(\begin{array}{c}
       x_1\\
       \vdots \\
       x_m
    \end{array}\right) f_{X}\left(x_{1}, \ldots, x_{m}\right) d x_{1} \ldots d x_{m}\\
    & = \left(\begin{array}{c}
      \int_{\mathbb{R}^{m}} x_1 f_{X}\left(x_{1}, \ldots, x_{m}\right) d x_{1} \ldots d x_{m}\\
      \vdots\\
      \int_{\mathbb{R}^{m}} x_m f_{X}\left(x_{1}, \ldots, x_{m}\right) d x_{1} \ldots d x_{m}
    \end{array}\right)\\
    & = \left(\begin{array}{c}
      \int_{\mathbb{R}} x_1 f_{X_1}\left(x_{1}\right) d x_{1}\\
      \vdots\\
      \int_{\mathbb{R}} x_m f_{X_m}\left(x_{m}\right) d x_{m}
    \end{array}\right)\\
    & = \left(\begin{array}{c}
       \mathbb{E}(X_1)\\
       \vdots\\
       \mathbb{E}(X_m)
    \end{array}\right)
    \end{align*}

Notice that $\int_{\mathbb{R}^{m}} x_1 f_{X}\left(x_{1}, \ldots, x_{m}\right) d x_{1} \ldots d x_{m} = \int_{\mathbb{R}} \dots \int_{\mathbb{R}} x_1 f_{X}\left(x_{1}, \ldots, x_{m}\right) d x_{1} \ldots d x_{m} = \int_{\mathbb{R}} x_1 \underbrace{\left(\int_{\mathbb{R}} \dots \int_{\mathbb{R}}f_{X}\left(x_{1}, \ldots, x_{m}\right) d x_{2} \ldots d x_{m}\right)}_{f_{X_1}\left(x_{1}\right)}d x_{1}$
