[BioC] limma posterior variance - revisited

[Gordon Smyth]

At 06:43 AM 9/06/2004, Naomi Altman wrote:
>The problem remains, but I have added a few lines of code that were 
>missing in the original posting.
>
>I have just run limma and am getting p-values after eBayes that are 
>smaller than the p-values before, leading to 100% of my genes being 
>declared significant at any value of FDR you care to use.

It seems very surprising to get 100% of genes significant, but nothing in 
the output that you give below suggests that anything is wrong. It all 
seems as it should be. You should tend to get smaller p-values after eBayes 
than before because the degrees of freedom increase, but not uniformly so 
because many of the residual standard deviations also increase.

>The design is a 1-way ANOVA with 6 treatments and 2 reps/treatment (which 
>I know is not great but ...)
>
>I thought that the denominator adjustment would make the posterior 
>sigma^2 > unadjusted MSE,  but this is not the case.

Empirical Bayes methods, like all shrinkage methods, shrink estimators 
towards a common value. This means that some values will go up, and some 
will go down. The help page says that eBayes() "uses an empirical Bayes 
method to      shrink the gene-wise sample variances towards a common 
value". What is happening is that the precisions (the inverse sample 
variances) are being set to their posterior means. You can see the 
complete, pretty simple, formula by following the URL for the reference 
given on the help page.

Gordon

>   Here are the commands I used to fit the model and do the ebayes 
> adjustment.
>
>design=model.matrix(~-1+factor(c(1,1,2,2,3,3,4,4,5,5,6,6)))
>colnames(design)=c("trt1","trt2","trt3","trt4","trt5","trt6")
>
>fitRMA=lmFit(RMAdata,design)
>
>contrast.matrix
>
>      c1 c2 c3 c4 c5
>trt1  1  1 1  1 1
>trt2 -1  1 1  1 1
>trt3  0 -2  1  1 1
>trt4  0  0 -3  1  1
>trt5  0  0 0 -5  1
>trt6  0  0 0  0 -5
>
>fitCont=contrasts.fit(fitRMA,contrast.matrix)
>fitAdj=eBayes(fitCont)
>
>ls.print(lsfit(fitRMA$sigma^2,fitAdj$s2.post))
>Residual Standard Error=0
>R-Square=1
>F-statistic (df=1, 22744)=1.632754e+35
>p-value=0
>
>           Estimate Std.Err      t-value Pr(>|t|)
>Intercept   0.0093       0 1.026963e+17        0
>X              0.5628       0 4.040735e+17        0
>
>mean(fitAdj$s2.post)
>[1] 0.02991697
>
>mean(fitRMA$sigma^2)
>[1] 0.03656270
>
>fitAdj$s2.prior
>[1] 0.02136298
>
>
>Naomi S. Altman                                814-865-3791 (voice)
>Associate Professor
>Bioinformatics Consulting Center
>Dept. of Statistics                              814-863-7114 (fax)
>Penn State University                         814-865-1348 (Statistics)
>University Park, PA 16802-2111