A permutation test works provided you follow it up by a test for an unknown number of high outliers.
Here's an implementation, some illustrations, and a discussion.
Implementation
The test statistic itself consists of the column variances:
stats <- apply(data, 2, var)
Compute a permutation statistic by independently scrambling each row, obtaining the variances of all data columns, and sorting them:
var.prmttn <- function(x, center=FALSE) {
if (center) x <- apply(x, 2, function(y) y - mean(y))
sort(apply(apply(x, 1, sample), 1, var))
}
Repeat this a bunch of times to obtain the null distribution of the permutation statistic:
sim <- replicate(1000, var.prmttn(data))
This is our reference for assessing the actual statistic. (For the use of the center parameter, see below.)
Illustrations
Let's synthesize some data in which all variables have randomly chosen values but have the same expected variance. prob specifies the expected frequencies of the values in values and n.obs stipulates the number of rows:
set.seed(17)
names <- c(116, 111, 11, 109, 108, 107, 105, 104, 55, 52, 58, 1)
vars <- sapply(names, function(i) paste("T", i, sep=""))
values <- c(0, 1/2, 1, -1)
prob <- c(.1, .4, .3, .1)
n.obs <- 60
data <- data.frame(sapply(vars,
function(s) sample(values, n.obs, prob=prob, replace=TRUE)))
row.names(data) <- sapply(runif(n.obs), function(x) sprintf("cg%08d", floor(10^8*x)))
The null distribution in this case is the set of distributions of the sorted variances, which are effectively visualized with boxplots:
extent <- c(min(min(stats), min(sim)), max(max(stats), max(sim)))
extent <- 1.1 * (extent - mean(extent)) + mean(extent)
do.call(boxplot, list(lapply(1:dim(sim)[1], function(i) sim[i,]),
ylim=extent, main="Permutation distribution",
ylab="Variance", xaxt="n"))
Over these we may superimpose the actual vector of test statistics, with appropriate labels to identify them:
points(sort(stats), pch=19, col="Red", cex=1.25)
text(rep(min(extent), length(stats)), labels=vars[order(stats)])
Let's do this again but increase the expected variances of two of the columns (by giving them greater probabilities of exhibiting the extreme values $\pm 1$):
prob.alt <- c(.1, .1, .4, .4)
data$T55 <- sample(values, n.obs, prob=prob.alt, replace=TRUE)
data$T11 <- sample(values, n.obs, prob=prob.alt, replace=TRUE)
We see two kinds of effects. Sure enough, the two different variables are detected by having two large values. At the same time, the values of all the rest of the variances are depressed compared to the null distribution.
Discussion
It may be difficult to obtain a useful p-value for such a result, given how changes in variance in one variable will affect both the null distribution and the apparent extremeness of variances in all the other variables. This is especially the case when the number of high-variance variables is not known in advance.
This situation is ideal for applying a multiple outlier detection procedure to the variances of the variables. One that searches for unusually large gaps among the largest values ought to perform well. The NIST handbook provides some guidance and links to further reading.
Finally, if the means of the variables differ substantially, their variation can make the null distribution less useful as a reference. Instead, permute their residuals in the test rather than the values themselves by setting center=TRUE in var.prmttn (or, more efficiently, by centering the variables individually before running the test).
