Variance of an average of random variables

Question

This seems like it should be a pretty common problem. I have four estimates of fishing effort, each with its own variance. For subsequent calculations, I want the mean of the four estimates, and a variance that captures variability of the numbers being averaged, and also propagates the measurement error.

In general terms, I have a series of estimates $(X_1, X_2, \dots X_n)$, each with a variance $(\sigma^2_1, \sigma^2_2, \dots \sigma^2_n)$. I know how to calculate the overall mean, $\bar{X}$ (sum of the $X$s over n). I am unhappy with the formulas I have found for the variance of $\bar{X}$.

It seems that that variance can come in two forms.

First, there is the 'standard' variance that measures the spread of the observed values around the mean. i.e.,

$$ \text{Var}_1 = \frac{\sum((X_i - \bar{X})^2)}{(n-1)}$$

But this variance ignores the fact that each of the $X$ values was measured with error.

Second, I could propagate the error, and calculate the variance as:

$$\text{Var}_2 = \frac{1}{n^2} (\sigma^2_1 + \sigma^2_2 + ... + \sigma^2_n)$$

But this variance ignores the fact that each of the X values differed from each other.

I have also seen:

$$\text{variance} = \frac{\sigma^2}{n}$$

but that appears to be for the situation in which the original $X$s all have the same variance, which doesn't apply here. If anything, I am looking for a 'final variance' that is bigger than any of the original $\sigma^2_i$ values, and also maybe bigger than the thing I called '$\text{Var}_1$', above.

How can I calculate a variance that captures variability of the numbers being averaged, and also propagates the measurement error?

Answer 1

If the $X_i$ are random variables with a variance $\sigma_i^2$, then the variance of $X=\sum_i X_i$ their sum is $\sigma_X^2$ is given by:

$\sigma_X^2=\sum_i \sigma_i^2 + 2 \sum_i \sum_{j<i} cov(X_i, X_j)$.

So if draw a random sample $x_i$ from these distributions, then $x=\sum_i x_i$ will be random (when we draw another sample we will have another value for the sum). The variance of these random outcomes (the sums) will be given by the formula supra.

If I divide a random variable $X$ with variance $\sigma_X^2$ by $n$ then the variance of this variable will be $\frac{\sigma_X^2}{n^2}$.

Consequently, if you draw a random sample $x_i$ from the distributions of $X_i$, then $\bar{x}=\frac{\sum_i x_i}{n}$ is random (with another sample I will have another average), if we draw many samples and each time compute the average, then we find the distribution of this random variable (the average $\bar{X}=\frac{\sum_i X_i}{n}$) and it has a variance equal to:

$\sigma_{\bar{X}}^2=\frac{\sum_i \sigma_i^2 + 2 \sum_i \sum_{j<i} cov(X_i, X_j)}{n^2}$

If the $X_i$ are independent then the covariances are zero and the formula simplifies to.

$\sigma_{\bar{X}}^2=\frac{\sum_i \sigma_i^2 }{n^2}$

If the $X_i$ are all independent and have the same variance $\sigma$ then this becomes:

$\sigma_{\bar{X}}^2=\frac{\sum_i \sigma ^2}{n^2}=\frac{n \sigma^2}{n^2}=\frac{\sigma^2}{n}$

Note: the capital letters represent random variables, the small letters are draws from a random variable and so one particular outcome.

Below some sample code to simulate this in the case of two independent variables

#draw hughe sample from x1 and from x2 
set.seed(1)
#random variable with variance 25
x1<-rnorm(n=1000000, mean=2, sd=5)
#random variable with variance 81
x2<-rnorm(n=1000000, mean=0, sd=9)

# compute the sum
x<-x1+x2
# compute the variance of the sum , should be close to 25+81
sd(x)^2

# compute the average
x.bar<-x/2
# compute the variance of the average , should be close to (25+81)/(2^2)
sd(x.bar)^2

Answer 2

If your estimates are independent and the likelihood on which each parameter estimate is based is reasonably close to a normal likelihood, then this seems like a meta-analysis with estimates $\hat{\theta}_i$ with standard errors $\sigma_i$, $i=1,\ldots,I$ and one reasonable overall estimate (the inverse variance fixed effects estimate) is $$\hat{\theta} = \frac{ \sum_{i=1}^I \frac{\hat{\theta}_i}{\sigma_i^2}}{\sum_{i=1}^I \sigma_i^{-2}} $$ and $$\text{SE}(\hat{\theta}) = \sqrt{ \frac{ 1 }{\sum_{i=1}^I \sigma_i^{-2}}}.$$ This "fixed effects" estimate assumes that the underlying parameter $\theta$ that you are trying to estimate is one fixed value and does not vary across the different estimates you have.

Alternatives include random effects estimates that assume that the true underlying value $\theta_i$ for each experiment you have varies a bit in a way that you cannot explain by observed covariates (if you have observed covariates that explain these differences, you should ideally use those and model the differences rather than "dumping" them into extra variability). That might be what you are looking for.

Additionally, note that if your estimates are not independent, then you need to take the covariance into account, as pointed out in another response.

Answer 3

I'd like to add these details to the answer by f coppens. Basically, my take-home message is that the amount by which the original numbers vary does not affect confidence in the mean of the original numbers.

#Here are the four estimates that I would like to average:
m1=4.8; sd1=0.438178045
m2=3.6; sd2=0.328633533
m3=2.4; sd3=0.219089022
m4=1.2; sd4=0.109544511
n=4

In 10,000 fake simulations of four flights, I get 10,000 estimates of 'average number of boats'

fake.values1 = array (0); fake.means1= array (0)
for (ix in 1:10000) {
  fake.values1[1:4] = c(rnorm(n=1, m1, sd1), rnorm(n=1, m2, sd2), rnorm(n=1, m3, sd3), rnorm(n=1, m4, sd4))
  fake.means1[ix] = mean(fake.values1)
}

# The variance of my 10,000 boat estimates is
sd(fake.means1)^2
# which is darn near the Var2 value from my original data, 0.0225

Next, I've add 1 to m1 and subtracted 1 from m4, thus increasing the overall range (and variance) in the values, but keeping the mean constant.

new.m1 = m1+1
new.m4 = m4-1

fake.values2 = array (0); fake.means2= array (0)
for (ix in 1:10000) {
  fake.values2[1:4] = c(rnorm(n=1, new.m1, sd1), rnorm(n=1, m2, sd2), rnorm(n=1, m3, sd3), rnorm(n=1, new.m4, sd4))
  fake.means2[ix] = mean(fake.values2)
}
# The variance of this set of 10,000 boat estimates is
sd(fake.means2)^2
# which, again is near the Var2 value, 0.0225

So what appears to be the answer to my question is 'Use your Var2 equation'. The Var1 equation captures variability in the estimates. Var2 captures variability in the mean of the estimates, and is not dependent on the spread of the original values being averaged.