[BioC] question about limma ebayes Fit lod score

[Gordon Smyth]

At 06:14 AM 31/08/2004, magates at u.washington.edu wrote:
>Hello all,
>
>With limma I am confused about the meaning of lods values (ebfit$lods) for an
>eBayes Fit Object.
>Many of the genes listed as significant by classifyTestsF have negative
>ebfit$lods scores, which is what is confusing me.

classifyTestsF() doesn't do any adjustment for multiple testing across 
genes. (This is documented, although briefly.) You are using a p-value of 
0.05 without any adjustment for multiple testing. On the other hand the 
B-statistic adjusts for multiple testing by its definition. In one case you 
are adjusting, in the other you are not, so it is not surprising that that 
they don't agree. If you use classifyTests() it is more usual to stick with 
small p-values like 0.00001, depending on the number of genes.

>  Associated ebfit$p.values of the genes listed by classifyTestsF are all 
> indeed
>smaller than the threshold p.value, as I would expect.
>(A plot of log(p.values) vs lods is pretty darn linear, again as I would
>expect, just the intercept is shifted by about -8 units.)
>
>My dim understanding of Lods would be that negative numbers are not
>significant. Do I misunderstand the meaning of the lods being reported, or 
>has
>something gone awry?

lods is a Bayesian measure and so doesn't have a direct interpretation in 
terms of "significance". The lods gives the log-odds of being 
differentially expressed, and it is for you to decide how stringent you 
want to be. The concept of statistical "significance" belongs to classical 
inference, which is a different philosophical framework.

Gordon

>As always, thanks for your help.
>
>Michael Gates
>
>some details:
>fit <- lmFit(MAb1,design)
>ebfit <- eBayes(fit)
>results <- classifyTestsF(ebfit,p.value=0.05)
>o <- order(ebfit$lods[,3],decreasing=TRUE)
>
># now filter thru a boolean function that -among other things -
># selects a subset of the significant genes for which the
># interaction effect is significant, but a particular main effect
># is not significant (details below)
>s <- myResultsFilter()
>
># and look at the fits, p.values, results, lods, etc
># of the desired subset of signficant genes
>data.frame(
>Name=MAb1$genes$Name[o[s[o]]],
>Expression=ebfit$coefficients[o[s[o]],3],
>P=ebfit$p.value[o[s[o]],3],
>Ftest=results05[o[s[o]],3],
>LOD=ebfit$lods[o[s[o]],3]
>)
>
>myResultsFilter <- function(fitObject=ebfit, resultObject=results,foldDiff=0)
>{
># make sure LB is a number
># and that it's larger than the desired fold difference
>   (!is.na(fitObject$coefficients[,3]) &
>   (abs(fitObject$coefficients[,3]) > foldDiff)) &
>   (
># if L is not significant, select significant LB
># that are not significant for B
>    ((resultObject[,1] %in% c(0)) & (resultObject[,3] %in% c(-1,1)) &
>     (resultObject[,2] %in% c(0))) |
># if L is significant, make sure it and LB go same way,
># and that B is not significant
>    ((resultObject[,1] %in% c(1)) &
>     (resultObject[,3] %in% c(1)) &
>     (resultObject[,2] %in% c(0))) |
>    ((resultObject[,1] %in% c(-1)) &
>     (resultObject[,3] %in% c(-1)) &
>     (resultObject[,2] %in% c(0))
>    )
>   )
>}