variance of the mean of correlated and uncorrelated data

Question

I read this paragraph in James et al, Introduction to Statistical Learning, p183-184 [1]:

Since the mean of many highly correlated quantities has higher variance than does the mean of many quantities that are not as highly correlated, the test error estimate resulting from LOOCV tends to have higher variance than does the test error estimate resulting from k-fold CV.

Can you give me a numerical example e.g. in R to check the validity of this claim?

I tried to check it using the following code:

x = 1:100 #highly correlated data
y = sample(100) #same data without correlation
var(x) == var(y) # TRUE

What is wrong with this code?

• LOOCV stands for "leave one out cross-validation"

[1]: James, G., Witten, D., Hastie, T., Tibshirani, R. (2013),
An Introduction to Statistical Learning with Applications in R,
Springer Texts in Statistics, Springer Science+Business Media, New York

Answer 1

The variance computed in the code views each array as if it were one sample of 100 separate values. Because both the array and its permuted version contain the same 100 values, they have the same variance.

The right way to simulate the situation in the quotation requires repetition. Generate a sample of values. Compute its mean. (This plays the role of "test error estimate.") Repeat many times. Collect all these means and look at how much they vary. This is the "variance" referred to in the quotation.

We may anticipate what will happen:

• When the elements of each sample in this process are positively correlated, when one value is high the others tend to be high, too. Their mean will then be high. When one value is low the others tend to be low, too. Their mean will then be low. Thus, the means tend either to be high or low.

• When elements of each sample are not correlated, the amount by which some elements are high is often balanced (or "canceled out") by other low elements. Overall the mean tends to be very close to the average of the population from which the samples are drawn--and rarely much greater or much less than that.

R makes it easy to put this into action. The main trick is to generate correlated samples. One way is to use standard Normal variables: linear combinations of them can be used to induce any amount of correlation you might like.

Here, for instance, are the results of this repeated experiment when it was conducted 5,000 times using samples of size $n=2$. In one case the samples were obtained from a standard Normal distribution. In the other they were obtained in a similar way--both with zero means and unit variances--but the distribution they were drawn from had a correlation coefficient of $90\%$.

The top row shows the frequency distributions of all 5,000 means. The bottom row shows the scatterplots generated by all 5,000 pairs of data. From the difference in spreads of the histograms, it is clear the set of means from the uncorrelated samples is less scattered than the set of means from the correlated samples, exemplifying the "canceling out" argument.

The difference in the amount of spread becomes more pronounced with higher correlation and with larger sample sizes. The R code allows you to specify these as rho and n, respectively, so you can experiment. Like the code in the question, its aim is to produce arrays x (from the uncorrelated samples) and y (from the correlated samples) for further comparison.

n <- 2
rho <- 0.9
n.sim <- 5e3
#
# Create a data structure for making correlated variables.
#
Sigma <- outer(1:n, 1:n, function(i,j) rho^abs(i-j))
S <- svd(Sigma)
Q <- S$v %*% diag(sqrt(S$d))
#
# Generate two sets of sample means, one uncorrelated (x) and the other correlated (y).
#
Z <- matrix(rnorm(n*n.sim), nrow=n)
x <- colMeans(Z)
y <- colMeans(Q %*% Z)
#
# Display the histograms of both.
#
par(mfrow=c(2,2))
h.y <- hist(y, breaks=50, plot=FALSE)
h.x <- hist(x, breaks=h.y$breaks, plot=FALSE)
ylim <- c(0, max(h.x$density))
hist(x, main="Uncorrelated", freq=FALSE, breaks=h.y$breaks, ylim=ylim)
hist(y, main="Correlated", freq=FALSE, breaks=h.y$breaks, ylim=ylim)
#
# Show scatterplots of the first two elements of the samples.
#
plot(t(Z)[, 1:2], pch=19, col="#00000010", xlab="x.1", ylab="x.2", asp=1)
plot(t(Q%*%Z)[, 1:2], pch=19, col="#00000010", xlab="x.1", ylab="x.2", asp=1)

Now when you compute the variances of the arrays of means x and y, their values will differ:

> var(x)
[1] 0.5035174
> var(y)
[1] 0.9590535

Theory tells us these variances will be close to $(1+1)/2^2 = 0.5$ and $(1 + 2\times 0.9 + 1)/2^2 = 0.95$. They differ from the theoretical values only because just 5,000 repetitions were done. With more repetitions, the variances of x and y will tend closer to their theoretical values.