<p>I assume that you are talking about the implementation of these methods in the <a href="https://bioconductor.org/packages/limma/" rel="nofollow noreferrer">limma package</a>. Otherwise this answer does not apply. </p>

<p>I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship:</p>

<pre><code>library("limma")
set.seed(123)

# Create some genes and samples
nGenes &lt;- 40
nSamples &lt;- 10
expr &lt;- matrix(rnorm(nGenes*nSamples), ncol = nSamples, nrow = nGenes)
genes &lt;- combn(LETTERS, 2, FUN = paste, collapse = "")
rownames(expr) &lt;- genes[seq_len(nGenes)]
colnames(expr) &lt;- letters[seq_len(nSamples)]

# Create some other genes with slightly higher values 
expr2 &lt;- expr+0.5
rownames(expr2) &lt;- genes[seq(41, to = 80)]
expr &lt;- rbind(expr, expr2)

# Some genes are more expressed in the last 5 samples
expr[41:80, 5:10] &lt;- expr[41:80, 5:10] + 2

# Heatmap of the simulated genes and samples
heatmap(expr, scale = "row", main = "Heatmap")
</code></pre>

<p><a href="https://i.sstatic.net/9doOR.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/9doOR.png" alt="Heatmap of the genes and samples"></a></p>

<pre><code># Create a design matrix and 
design &lt;- matrix(c(rep(1, 10), rep(0, 5), rep(1, 5)), 
                 ncol = 2, nrow = 10)
fit &lt;- lmFit(expr, design)
fit2 &lt;- eBayes(fit)
tt &lt;- topTable(fit2, coef = 2, number = Inf, sort.by = "none")

randomGeneSet &lt;- sample(1:80, 5, replace = FALSE)
# Calculate the enrichment with geneSetTest
geneSetTest(1:40, tt$t)
## [1] 1
geneSetTest(41:80, tt$t)
## [1] 9.043842e-13
geneSetTest(randomGeneSet, tt$t)
## [1] 0.9730717

# Explore the interGeneCorrelation and the  variance inflation factor
interGeneCorrelation(expr[1:40, ], design)
## $vif
## [1] 0.2055766
## 
## $correlation
## [1] -0.02036983
## 
interGeneCorrelation(expr[41:80, ], design)
## $vif
## [1] 12.79611
## 
## $correlation
## [1] 0.3024643
## 
interGeneCorrelation(expr[randomGeneSet, ], design)
## $vif
## [1] 0.3526011
## 
## $correlation
## [1] -0.1618497

# Calculate the enrichment with camera
camera(expr, 1:40, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1     40 -0.02036983      Down 3.0445e-06
camera(expr, 41:80, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction      PValue
## set1     40   0.3024643        Up 0.002402933
camera(expr, randomGeneSet, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1      5  -0.1618497      Down 0.04272392
</code></pre>

<p>In this simulation we created three gene sets were the p-value of both camera and geneSetTest where significant, one where only camera was significant and another one where both were below the 0.05 threshold of significance. </p>

<p>As you can see we can calculate the variance inflation factor for each gene set, which also yields the intercorrelation for each gene set we can summarize it into a table:</p>

<pre><code>               camera    VIF    correlation  geneSetTest
GS1 (1:40)     3.0e-06   0.203     -0.020     1.0
GS2 (41:80)    0.024     12.8       0.30      9.0e-13
randomGeneSet  0.043     0.353     -0.16      0.97
</code></pre>

<p>As GS1 and GS2 are the complete opposite we see that they have quite different inter correlation, which leads to a different VIF. We can see that higher VIF or small doesn't mean that the gene set has lower p-value for camera.</p>

<blockquote>
  <p>How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? </p>
</blockquote>

<p>I am unsure how to interpret that comparison. It would be better to explore this  with real data (for example from ALL library), and contact the authors, either via the <a href="https://support.bioconductor.org" rel="nofollow noreferrer">support site in Bioconductor</a> or through personal communications. </p>

<blockquote>
  <p>Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another?</p>
</blockquote>

<p>As you might have noticed I used <code>inter.gene.cor = NA</code> in order to estimate the correlations and avoid guessing the correlation between the genes but still the correction of intercorrelations has lead to higher p-values. So no, a wrong estimation of the correlation of the query set doesn't explain why the p-value is lower. </p>

<p>Note: I only used the unranked method of the limma 3.34.1, I didn't use the preranked camera (<code>cameraPR</code>) because it doesn't have the same results as in <code>camera</code> and it doesn't allow to estimate the intercorrelation.</p>

CC BY-SA 3.0
2017-12-04T20:31:51.523Z
Noticed an error in your code - in the call to topTable, you need to include 'sort.by = "none"' or else it reorders by significance, making the 1:40 and 41:80 indices for geneSetTest no longer applicable.The geneSetTest p-values for the three gene sets now come out to 1, 9.04e-13, and 0.97. It seems that the geneSetTest p-value is arguably more "correct" here, since only GS2 truly shows differential expression between samples 1-5 and 6-10. CAMERA says they're all significant!
Nuclear Hoagie
CC BY-SA 3.0
2017-12-04T21:13:15.640Z
@NuclearWang Nice catch, that's why I usually use `ids2indeces`! I am not sure if this is the real case, I'll try to expand the answer using the ALL dataset in a couple of days. Hope this helps! If you find a better answer, please let me know, I am interested in those methods :)
llrs
CC BY-SA 4.0
2018-07-30T18:09:02.900Z
Why is CAMERA giving a significant p-value for the first comparison?
flies
CC BY-SA 4.0
2018-07-30T18:31:24.210Z
@flies There is only one comparison. Do you refer to the enrichment in the GS2?
llrs
CC BY-SA 4.0
2018-07-31T15:40:39.767Z
@Llopis I refer to GS1. AFAICT, rows 1 through 40 are just a uniform normal distribution across samples, so the test should be negative.
flies
CC BY-SA 4.0
2018-07-31T15:51:05.213Z
@flies Yes, you are right. The [GSEABenchmarkeR](http://bioconductor.org/packages/release/bioc/html/GSEABenchmarkeR.html) aims to find a consensus. My understanding is that the low correlation and the VIF makes strange things. But it might be an effect of a small dataset...
llrs
CC BY-SA 4.0
2018-07-31T17:45:20.673Z
@Llopis thanks for your response on this old question! It's really odd and rather disturbing that this is what you found. I will say that I recently had a case where I used `cameraPR` where the `ids2indices` was scrambled with respect to the expression input (I used t-statistics output by limma), and the significance was all negative. Once I fixed the bug that scrambled the indices, I got good results that made biological sense. So camera will give negative results in some cases - at least with fixed `inter.gene.cor`!
flies
CC BY-SA 3.0
2017-12-02T15:46:24.577Z
2017-12-04T21:08:23.367Z
2017-12-04T21:08:23.367Z
llrs
<p>I am performing a gene set enrichment analysis to determine if particular gene sets are coherently up- or down-regulated. I have seen several statistics for computing a p-value of GSEA-style enrichment. I'm particularly interested in the differences between a mean-rank gene set test, like the one implemented in the geneSetTest R function, and the correlation-corrected version of the test implemented with the CAMERA method.</p>

<p>As I understand it, CAMERA corrects for correlation within a gene set, because a highly correlated query set violates the underlying assumptions and results in a bias toward more significant p-values. So, if I find a significant mean-rank gene set test p-value, but an insignificant CAMERA p-value, I conclude that my query set is correlated, hence the correction factor and reduction in significance.</p>

<p>How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another? I feel like there's a clear connection between CAMERA's variance inflation factor and a reduction in significance, but I struggle to see how the correction factor can increase the significance over the mean-rank test.</p>

<gene><r><differential-expression><pathway><gsea>
Comparison of gene set enrichment statistics
# Comparison of gene set enrichment statistics

(Asked by: username_0 on 2017-12-01)

I am performing a gene set enrichment analysis to determine if particular gene sets are coherently up- or down-regulated. I have seen several statistics for computing a p-value of GSEA-style enrichment. I'm particularly interested in the differences between a mean-rank gene set test, like the one implemented in the geneSetTest R function, and the correlation-corrected version of the test implemented with the CAMERA method.

As I understand it, CAMERA corrects for correlation within a gene set, because a highly correlated query set violates the underlying assumptions and results in a bias toward more significant p-values. So, if I find a significant mean-rank gene set test p-value, but an insignificant CAMERA p-value, I conclude that my query set is correlated, hence the correction factor and reduction in significance.

How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another? I feel like there's a clear connection between CAMERA's variance inflation factor and a reduction in significance, but I struggle to see how the correction factor can increase the significance over the mean-rank test.

**Tags:** gene, r, differential-expression, pathway, gsea

---

## 1 answers

### Answer by username_1 on 2017-12-02. Score: 4

I assume that you are talking about the implementation of these methods in the limma package. Otherwise this answer does not apply. 
I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship:
library("limma")
set.seed(123)

# Create some genes and samples
nGenes <- 40
nSamples <- 10
expr <- matrix(rnorm(nGenes*nSamples), ncol = nSamples, nrow = nGenes)
genes <- combn(LETTERS, 2, FUN = paste, collapse = "")
rownames(expr) <- genes[seq_len(nGenes)]
colnames(expr) <- letters[seq_len(nSamples)]

# Create some other genes with slightly higher values 
expr2 <- expr+0.5
rownames(expr2) <- genes[seq(41, to = 80)]
expr <- rbind(expr, expr2)

# Some genes are more expressed in the last 5 samples
expr[41:80, 5:10] <- expr[41:80, 5:10] + 2

# Heatmap of the simulated genes and samples
heatmap(expr, scale = "row", main = "Heatmap")


# Create a design matrix and 
design <- matrix(c(rep(1, 10), rep(0, 5), rep(1, 5)), 
                 ncol = 2, nrow = 10)
fit <- lmFit(expr, design)
fit2 <- eBayes(fit)
tt <- topTable(fit2, coef = 2, number = Inf, sort.by = "none")

randomGeneSet <- sample(1:80, 5, replace = FALSE)
# Calculate the enrichment with geneSetTest
geneSetTest(1:40, tt$t)
## [1] 1
geneSetTest(41:80, tt$t)
## [1] 9.043842e-13
geneSetTest(randomGeneSet, tt$t)
## [1] 0.9730717

# Explore the interGeneCorrelation and the  variance inflation factor
interGeneCorrelation(expr[1:40, ], design)
## $vif
## [1] 0.2055766
## 
## $correlation
## [1] -0.02036983
## 
interGeneCorrelation(expr[41:80, ], design)
## $vif
## [1] 12.79611
## 
## $correlation
## [1] 0.3024643
## 
interGeneCorrelation(expr[randomGeneSet, ], design)
## $vif
## [1] 0.3526011
## 
## $correlation
## [1] -0.1618497

# Calculate the enrichment with camera
camera(expr, 1:40, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1     40 -0.02036983      Down 3.0445e-06
camera(expr, 41:80, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction      PValue
## set1     40   0.3024643        Up 0.002402933
camera(expr, randomGeneSet, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1      5  -0.1618497      Down 0.04272392

In this simulation we created three gene sets were the p-value of both camera and geneSetTest where significant, one where only camera was significant and another one where both were below the 0.05 threshold of significance. 
As you can see we can calculate the variance inflation factor for each gene set, which also yields the intercorrelation for each gene set we can summarize it into a table:
               camera    VIF    correlation  geneSetTest
GS1 (1:40)     3.0e-06   0.203     -0.020     1.0
GS2 (41:80)    0.024     12.8       0.30      9.0e-13
randomGeneSet  0.043     0.353     -0.16      0.97

As GS1 and GS2 are the complete opposite we see that they have quite different inter correlation, which leads to a different VIF. We can see that higher VIF or small doesn't mean that the gene set has lower p-value for camera.

How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? 

I am unsure how to interpret that comparison. It would be better to explore this  with real data (for example from ALL library), and contact the authors, either via the support site in Bioconductor or through personal communications. 

Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another?

As you might have noticed I used inter.gene.cor = NA in order to estimate the correlations and avoid guessing the correlation between the genes but still the correction of intercorrelations has lead to higher p-values. So no, a wrong estimation of the correlation of the query set doesn't explain why the p-value is lower. 
Note: I only used the unranked method of the limma 3.34.1, I didn't use the preranked camera (cameraPR) because it doesn't have the same results as in camera and it doesn't allow to estimate the intercorrelation.

> **Comments:**
> **username_0**: Noticed an error in your code - in the call to topTable, you need to include 'sort.by = "none"' or else it reorders by significance, making the 1:40 and 41:80 indices for geneSetTest no longer applicable.The geneSetTest p-values for the three gene sets now come out to 1, 9.04e-13, and 0.97. It seems that the geneSetTest p-value is arguably more "correct" here, since only GS2 truly shows differential expression between samples 1-5 and 6-10. CAMERA says they're all significant!  
> **username_1**: @NuclearWang Nice catch, that's why I usually use `ids2indeces`! I am not sure if this is the real case, I'll try to expand the answer using the ALL dataset in a couple of days. Hope this helps! If you find a better answer, please let me know, I am interested in those methods :)  
> **username_2**: Why is CAMERA giving a significant p-value for the first comparison?  
> **username_1**: @flies There is only one comparison. Do you refer to the enrichment in the GS2?  
> **username_2**: @Llopis I refer to GS1. AFAICT, rows 1 through 40 are just a uniform normal distribution across samples, so the test should be negative.  
> **username_1**: @flies Yes, you are right. The [GSEABenchmarkeR](http://bioconductor.org/packages/release/bioc/html/GSEABenchmarkeR.html) aims to find a consensus. My understanding is that the low correlation and the VIF makes strange things. But it might be an effect of a small dataset...  
> **username_2**: @Llopis thanks for your response on this old question! It's really odd and rather disturbing that this is what you found. I will say that I recently had a case where I used `cameraPR` where the `ids2indices` was scrambled with respect to the expression input (I used t-statistics output by limma), and the significance was all negative. Once I fixed the bug that scrambled the indices, I got good results that made biological sense. So camera will give negative results in some cases - at least with fixed `inter.gene.cor`!  

---

