# Variance of the mean of correlated variables

I am totally puzzled by this one. Let's say that my data is $$X \sim MVN(\mu,\Sigma)$$ where the variance-covariance matrix $$\Sigma$$ is of size $$m \times m$$ and is constructed in such a way that all variances are the same and all covariances are the same.

I draw a sample from this multivariate normal distribution and compute the variance of the sampled data. My intuition says that this estimated variance should be equal to the variance of the data minus the covariance, because there is dependence so we should take this into account when computing the variance.

My initial thought was that the variance of the sampled data is equal to the variance of the mean of correlated variables. However, this is not the case when I compute this variance with the well-known formula for this (see e.g., this question and this question and page 228 of this book on meta-analysis) $$$$(\frac{1}{m})^2 (\sum_{i=1}^m \sigma^2_i + \sum_{i\neq j}(r_{ij} \sqrt{\sigma^2_i}\sqrt{\sigma^2_j}))$$$$ where $$r_{ij}$$ is the correlation between the scores and $$\sigma^2_i$$ and $$\sigma_j^2$$ are the variances.

My questions are:

1. Is it correct to estimate the variance of the sampled data of the multivariate normal distribution in the way I am doing it, so variance minus covariance (i.e., $$\sigma^2_i-\sigma_{ij}$$)
2. I am apparently not computing the variance of the mean of correlated variables when I subtract the covariance from the variance. What is the name of the quantity that I am computing?

Below some R code where I show that these two variances are not the same.

Thank you in advance for any help!

library(MASS)

m <- 100 # Number of outcomes
sigma2 <- 1 # Variance
r <- 0.9 # Correlation between outcomes

### Create variance-covariance matrix
cov <- r * sqrt(sigma2*sigma2) # Covariance among outcomes
Sigma <- matrix(cov, nrow = m, ncol = m)
diag(Sigma) <- sigma2

### Generate data from multivariate normal distribution
dat <- mvrnorm(n = 1000, mu = rep(0,m), Sigma = Sigma)

### Compute variance per sample from multivariate normal
vars <- apply(dat, 1, var)

### Mean of computed variances
mean(vars)

### Mean of computed variances seems to be equal to
sigma2 - cov

### Formula for computing the variance of the mean of correlated outcomes
(1/m)^2 * (sum(diag(Sigma)) + r*sqrt(sigma2)*sqrt(sigma2)*(m*m-2))


$$\newcommand{\one}{\mathbf 1}$$If we have $$X\sim\mathcal N(\mu,\Sigma)$$ then the sample variance of a single draw can be computed as $$\frac 1{m-1} X^TSX$$ where $$S = I - \frac 1m \one\one^T$$ is the matrix that projects into the space orthogonal to $$\one$$. This is a quadratic form so we can compute its mean as $$\text{E}(X^TSX) = \text{tr}(S\Sigma) + \mu^T\Sigma\mu.$$ In your experiment we have $$\mu = \mathbf 0$$ so the variance is just proportional to that trace.
If you compare mean(vars) to S <- diag(m) - matrix(1/m,m,m); sum(diag(S %*% Sigma)) / (m-1) (and maybe increase n) you'll see that they agree.
Note that $$\text{tr}(S\Sigma) = \text{tr}(\Sigma - \frac 1m \one\one^T\Sigma)$$ so $$\text{tr}(S\Sigma) = \sum_{i=1}^m\left( \sigma^2_{i} - \frac 1m\sum_{j=1}^m \Sigma_{ij}\right)$$ which is the sum of the diagonal entries minus the row means of $$\Sigma$$. With your experiment this is $$\sum_{i=1}^m \left(1 - \frac{1}{m}\left((m-1)\cdot 0.9 + 1\right)\right) \\ = m - 1 - (m-1)\cdot 0.9 \\ = \frac{m-1}{10}$$ so the actual variance is $$\frac{\frac{m-1}{10}}{m-1} = \frac 1{10}$$ which the simulation confirms.
Another way to do this is to use some properties of the trace to get $$\text{tr}(S\Sigma) = \text{tr}(\Sigma) - \text{tr}\left(\frac 1m \one^T\Sigma\one\right) \\ = \text{tr}(\Sigma) - \frac 1m \one^T\Sigma\one \\ = \sum_{i=1}^m \sigma^2_i - \frac 1m \sum_{ij} \Sigma_{ij}$$