Multiple comparisons problem

Yanwu Guo

Multiple testing (comparisons) problem occurs when one considers a set of statistical inferences simultaneously. The more inferences are made, the more likely erroneous inferences are to occur, for example, when we perform high through-put data analysis. See more https://en.wikipedia.org/wiki/Multiple_comparisons_problem.

Let’s see an example: we select 100 genes which has the same expression (following the same normal distribution). Each genes are measured 10 times for each group. Then, we perform t.test.

set.seed(1)
pvalues <- c()
for(i in 1:100){ # 100 genes
  pvalues[i] <- t.test(rnorm(10), rnorm(10))$p.value # Each genes are measured 10 times for each "group"
}

plot(pvalues, ylab="p values")
abline(h=0.05, col="red")
text(10, 0.1, "p values= 0.05", col = "red")

At the end we got 100 p values and let’s look at the plot. The red line is the generally used cutoff (p-value equals to 0.05). There are 4 points crossed the line, which is accuturlly what we expected by using the 0.05 cut off.

Multiple testing correction

There are several methods developed to prevent the problem. The first fuction is p.adjust from The R Base Package.

Bonferroni correction is the simplest and most conservative approach. For example, in the example above, with 100 tests and p = 0.05, you’d only reject a null hypothesis if the p-value is less than 0.05/100 = 0.0005. Or the p-values from the t.test are mutipled by 100, then we still use the 0.05 cutoff.

## 
?p.adjust
p.adjust.methods # options for different adjust methods
## [1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"        
## [6] "BY"         "fdr"        "none"
df<- data.frame(1:100)
for (i in 1:length(p.adjust.methods)) {
  df[,i] <- p.adjust(pvalues, method = p.adjust.methods[i] )
}
names(df) <- p.adjust.methods
head(df)
##   holm  hochberg    hommel bonferroni        BH BY       fdr      none
## 1    1 0.9988074 0.9988074          1 0.9697029  1 0.9697029 0.7840382
## 2    1 0.9988074 0.9988074          1 0.9621211  1 0.9621211 0.5287546
## 3    1 0.9988074 0.9988074          1 0.9988074  1 0.9988074 0.9803540
## 4    1 0.9988074 0.9988074          1 0.9272705  1 0.9272705 0.1471865
## 5    1 0.9988074 0.9988074          1 0.9988074  1 0.9988074 0.9717914
## 6    1 0.9988074 0.9988074          1 0.9272705  1 0.9272705 0.2099613
plot(pvalues, ylab="p values")
abline(h=0.05, col="red")
text(10, 0.1, "p values= 0.05", col = "red")
n <- length(p.adjust.methods)
colors <- rainbow(n, s = 1, v = 1, start = 0, end = max(1, n - 1)/n, alpha = 1)
for (i in 1:n) {
  points(df[,i], col= colors[i], pch=20)
}

After the adjustment, there is no p-values less than 0.05.

The second function we are going to use is qvalue from Bioconductor. the qvalue funciton returns a qvalue object.

A list of object type “qvalue” containing:

call Function call.

pi0 An estimate of the proportion of null p-values.

qvalues A vector of the estimated q-values (the main quantity of interest).

pvalues A vector of the original p-values.

lfdr A vector of the estimated local FDR values.

significant If fdr.level is specified, and indicator of whether the q-value fell below fdr.level (taking all such q-values to be significant controls FDR at level fdr.level).

pi0.lambda An estimate of the proportion of null p-values at each lambda value (see vignette).

lambda A vector of the lambda values utilized to obtain pi0.lambda.

library(qvalue)
qobj <- qvalue(pvalues) 
class(qobj)
## [1] "qvalue"
str(qobj)
## List of 8
##  $ call      : language qvalue(p = pvalues)
##  $ pi0       : num 1
##  $ qvalues   : num [1:100] 0.97 0.962 0.999 0.927 0.999 ...
##  $ pvalues   : num [1:100] 0.784 0.529 0.98 0.147 0.972 ...
##  $ lfdr      : num [1:100] 1 1 1 1 1 1 1 1 1 1 ...
##  $ pi0.lambda: num [1:19] 1.01 1.02 1.04 1.02 1.01 ...
##  $ lambda    : num [1:19] 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 ...
##  $ pi0.smooth: num [1:19] 1.03 1.02 1.02 1.02 1.02 ...
##  - attr(*, "class")= chr "qvalue"
head(qobj$qvalues)
## [1] 0.9697029 0.9621211 0.9988074 0.9272705 0.9988074 0.9272705