Could Bessel's correction make sample variance estimation even more biased?

Question

It is well known that Bessel's correction creates an unbiased estimator of variance. What it basically does is divide by $n-1$ instead of $n$.

Now what I did is that I chose a few number, like $1,2,3,4,5,60$ and calculated it's population variance which is $452.92$. I then took all possible 4-combinations (15 altogether) and calculated their sample variance (dividing by $n-1$) respectively. The average sample variance is $543.5$ which is off by $90.58$.

When I take the population variance (dividing by $n$) of the samples instead I get an average variance of $407.63$ which is off by only $45.29$!

I did several other experiments of the same kind with different numbers, population and sample sizes, all with this strange result that the population variance of the samples is less biased than the so called unbiased sample variance.

How can that be? What am I missing?

EDIT
Because of the illuminating discussion in the comments I posted this follow-up question:
Unbiased estimator of variance for samples *without* replacement

Answer 1

You sampled without replacement (otherwise there would be a lot more than 15 combinations of sample sizes of 4), and, as a consequence, the individual draws are not independent. That messes things up, as one might suspect. To see this, consider drawing 6 times from (1,2,3,4,5,60) without replacement; all your samples are the same, and the sample variance using n in the denominator equals the population variance every time.

For sampling without replacement, the sample mean is still unbiased. It's actually less variable than in the sampling with replacement case. Math helps, but so do examples:

x <- c(1,2,2,3,4,5,5,6,7,8)
xbar <- rep(0,10000)
for (j in 1:length(xbar)) {
  xbar[j] <- mean(sample(x,4))
}

# The population mean and average sample mean

> mean(x)
[1] 4.3
> mean(xbar)
[1] 4.303725

# The variance of the sample means and the true variance of the
# sample mean under sampling with replacement
> var(xbar)
[1] 0.7922241
> var(x)/4
[1] 1.336111

I am being imprecise about the conditions required for the variance of the mean to equal the population variance / the sample size by implicitly equating "sampling with replacement" with the real condition, namely, "the observations are independent and identically distributed", where the latter is clearly stronger than the former, so consider this statement a clarification of that point.

Answer 2

The problem here is that the metric for "population variance" that divides by the finite population size (without using Bessel's correction) is a bad measure of variance to begin with. If you incorporate Bessel's correction into the (finite) population variance for the data vector $\mathbf{x} = (1,2,3,4,5,60)$ you get:

$$s_N^2 = \frac{1}{N-1} \sum_i (x_i-\bar{x})^2 = \frac{2717.5}{5} = 543.5 \neq 452.92.$$

Now, if you sample with replacement from your data then the sample variance will average to this value. In particular, averaging the sample variance over all 15 combinations of 4 distinct data points gives this value. (Try this and you will see.)

The incompatibility you are getting in your question is hardly surprising when you consider that any finite set of data can be considered as a random sample of an infinite sequence, and in this context, the sample variance (incorporating Bartlett's correction) of any amount of data is an unbiased estimator of the superpopulation variance (i.e., the variance parameter). This, it is unsurprising that a sample variance of four data points does not average to the same thing as a "population variance" of six data points, when the latter does even use the same variance formula as the former.

---

As suggested in O'Neill (2014), Bartlett's correction should be incorporated into the "population variance" for a finite population (contrary to the standard approach presented in text-books and introductory statistics courses). It makes no sense to have a metric for the (finite) population variance that uses a different denominator than the sample variance. Not only does this create the strange situation you highlight in your question, but it also then means that the finite population variance is not an unbiased estimator of the variance for a superpopulation into which it is embedded.)

"The reason we choose to incorporate Bartlett’s correction is that it makes sense to consider the finite population variance as an estimator of a larger infinite superpopulation variance in the context of a superpopulation model. With this consideration, Bartlett’s correction ensures that the sample variance and population variance both have the same expected value and therefore function as unbiased estimators of the superpopulation variance parameter." (O'Neill 2014, pp 282-283)