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

  • $\begingroup$ If you sort(sample(100)) you will see it is identical to 1:100 and hence their variances are identical. Can't help you with the first bit of your post - I would have thought correlated quantities have lower variance (e.g. intra-cluster correlations) but then I don't know what LOOCV is. $\endgroup$ – Matthew Jul 12 '16 at 21:23
  • $\begingroup$ I know they are the same and just the ordering is changed. However 1:100 are correlated numbers but sample(100) are not. $\endgroup$ – Farid Cheraghi Jul 12 '16 at 21:27
  • $\begingroup$ The vector 1:100 is no more correlated than the vector sample(100). They may have been generated differently but, apart from the ordering, are identical. Certainly the calculation of variance does not take ordering into account. There are examples on-line of how to simulate correlated data which are probably what you need. $\endgroup$ – Matthew Jul 12 '16 at 21:38
  • $\begingroup$ Try acf(x) and acf(y) and see yourself! $\endgroup$ – Farid Cheraghi Jul 12 '16 at 21:47
  • $\begingroup$ Ah, I didn't think about autocorrelation. Still, the reason the variances are equal is because the ordering is irrelevant to the var function, just as it would be if you calculated the variance by had. Also, and possibly more helpfully, see link $\endgroup$ – Matthew Jul 12 '16 at 22:04

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.
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.

  • 3
    $\begingroup$ I couldn't find the quote in ESLII (10th printing) using a text search in the pdf, but it's word-for-word in ISLR (so that's the reference I inserted in the question just now) $\endgroup$ – Glen_b -Reinstate Monica Jul 13 '16 at 1:02
  • $\begingroup$ @whuber, thank you very much. I used to create two dimensional correlated quantities using mvtnorm package. Sigma <- outer(1:n, 1:n, function(i,j) rho^abs(i-j)) xy = mvtnorm::rmvnorm(5000,mean=c(0,0),sigma = Sigma) How can I understand the way Q %*% Z leads to correlated quantities? books, references? $\endgroup$ – Farid Cheraghi Jul 13 '16 at 8:59
  • 2
    $\begingroup$ The simplest way is a calculation: $$\operatorname{Cov}(QZ)=Q\operatorname{Cov}(Z)Q^\prime=QQ^\prime=V\sqrt{D} \sqrt{D}V^\prime=VDV^\prime=S.$$ I chose this method, instead of mvtnorm, because it enabled me to use exactly the same values of $Z$ for both simulated datasets, forcibly demonstrating the result was not an artifact of the random number generation. (Thank you, @glen_b, for correcting the reference.) $\endgroup$ – whuber Jul 13 '16 at 13:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.