Suppose that we have samples from $m$ different $p$-dimensional normal multivariate distributions, where they share a common covariance matrix $\Sigma$ but the mean vectors may be different for each population. That is, the distribution of the populations are $\mathrm{N}_p(\mu_i, \Sigma)$ , for $i = 1,2, \ldots, m.$ Each sample is made up of $n_i$ observations, depending on the group.
So, $\mathbf{y}_{ij}$ is observation $j$ (for $j = 1,\ldots,n_i$) coming from group $i$ (for $i = 1,\ldots,m$). Notice that $\mathbf{y}_{ij}$ is a $p-$dimensional vector that comes from the multivariate normal distribution $\mathrm{N}_p(\mu_i, \Sigma)$.
An estimator for the covariance matrix for all of the data (ignoring the potential difference between groups) is $$\mathbf{S}=\frac{1}{N-1} \sum_{i=1}^m \sum_{j=1}^{n_i}\left(\mathbf{y}_{i j}-\bar{\mathbf{y}}\right)\left(\mathbf{y}_{i j}-\bar{\mathbf{y}}\right)^{\prime} = \dfrac{SS_T}{N-1},$$ where $N= \sum_{i=1}^m n_i$ is the total number of observations and $SS_T$ refers to the total variation (the sum of squares). Now, the total variation $SS_T$ can be broken down into the within-group and between-group variation components:
\begin{align} \tag{1} SS_T= & \sum_{i=1}^m \sum_{j=1}^{n_i}\left(\mathbf{y}_{i j}-\bar{\mathbf{y}}_i\right)\left(\mathbf{y}_{i j}-\bar{\mathbf{y}}_i\right)^{\prime}+ \sum_{i=1}^m n_i\left(\bar{\mathbf{y}}_i-\bar{\mathbf{y}}\right)\left(\bar{\mathbf{y}}_i-\bar{\mathbf{y}}\right)^{\prime} =SS_W + SS_B, \end{align} where $SS_W$ corresponds to the within-group variation and $SS_B$ to the between-batch variation. From here, a weighted average matrix can be used for the within-group covariance matrix (especially useful if the mean vectors of the distributions are very different): $$\mathbf{S}_W=\frac{\sum_{i=1}^m\left(n_i-1\right) \mathbf{S}_i}{N-m}=\frac{SS_W}{N-m},\tag{2}$$ where $\mathbf{S}_i$ is the covariance matrix estimate from group $i$. This makes sense to me since it's just a weighted average of individual covariance matrices, so it will only model the within-group variation. Similarly, a between-group covariance matrix can be defined as $$\mathbf{S}_B = \dfrac{SS_B}{m-1} = \frac{\sum_{i=1}^m n_i\left(\bar{\mathbf{y}}_i-\bar{\mathbf{y}}\right)\left(\bar{\mathbf{y}}_i-\bar{\mathbf{y}}\right)^{\prime}}{m-1}. \tag{3}$$
The formulas make sense to me, at least mathematically. The total variation can be broken down into its between and within-group components as in Equation (1). Intuitively, using the within-group covariance estimator $S_W$ makes sense, since we are averaging individual covariance matrices, and giving a larger weight to groups with more observations. This would be a good estimate for the matrix $\Sigma$.
However, I don't understand the intuition behind Equation (3), because of the term $n_i$. If the term were not there, the expression would look like this: $$\mathbf{S}_B^* = \frac{\sum_{i=1}^m \left(\bar{\mathbf{y}}_i-\bar{\mathbf{y}}\right)\left(\bar{\mathbf{y}}_i-\bar{\mathbf{y}}\right)^{\prime}}{m-1}. $$ and I could explain it as the covariance matrix of the group means. Thus, the diagonal of $\mathbf{S}_B^*$ would correspond to the between-group variances for each of the $p$ variables or characteristics being modelled. Adding $n_i$ could be understood as a weighting mechanism (where groups with more observations should be considered more, I understand). But then, what is the intepretation of the resulting matrix $\mathbf{S}_B$? What does each entry tell us? Notice that $\mathbf{S}$ and $\mathbf{S}_W$ have $N$ in their denominator, while $\mathbf{S}_B$ does not. This means that, the bigger the groups $n_i$, the bigger the entries of $\mathbf{S}_B$. If a process has a specific between-group variation, why would making groups bigger affect those values?
