[Bioperl-l] Re: [BioC] Questions about multtest

[Sean Davis]

> OK, I using the multtest package to analyse my data, following the
> instructions in multtest.pdf.
> 
> I run:
> 
>> t <- mt.teststat(data[,6:12], c(0,0,0,1,1,1,1), test="t")
> 
> which calculates the t statistic for my data.  The t statistic for my first
> gene comes up as:
> 
>> t[1]
> [1] 40.60158
> 
> Presumably, this is equivalent to me running t.test:
> 
>> t.test(data[1,9:12], data[1,6:8], var.equal=FALSE, alternative="two.sided")
> 
>       Welch Two Sample t-test
> 
> data:  data[1, 9:12] and data[1, 6:8]
> t = 40.6016, df = 2, p-value = 0.0006061
> alternative hypothesis: true difference in means is not equal to 0
> 95 percent confidence interval:
> 1.713804 2.120092
> sample estimates:
>   mean of x     mean of y
> -1.596190e-15 -1.916948e+00
> 
> So I want to know how I can get p-values for the t statistics I have just
> calculated using mt.teststat.
> 
> This is where I get confused - multtest.pdf says I should "compute raw nominal
> two-sided p-values for the 3,051 test statistics using the standard Gaussian
> distribution":
> 
>> rawp0 <- 2 * (1 - pnorm(abs(t)))

Keep in mind what this is doing:  computing a p-value based on the normal
approximation of the t-distribution with infinite degrees of freedom.  In
your case, this approximation does not hold (probably) because of the
smaller than infinite (or, in practice 50 or so) number of degrees of
freedom.  

The above is asking what the probability of seeing a z-score of 40, which is
nearly equivalent to 0 (and, to the number of significant digits here, IS
0).  

What you probably want is:
rawp0 <- 2 * (1 - pt(abs(t),df=2))

Like so:

> 2*(1-pt(40.60158,df=2))
[1] 0.000606065

Which agrees with your t-test value.  Then, you can soldier on.

> Soldiering on, I want to calculate adjusted p-values accoridng to Benjamini
> and Hochberg:
> 
>> res <- mt.rawp2adjp(rawp0, "BH")
>> adjp <- res$adjp[order(res$index), ]
>> adjp[1]
> [1] 0



Sean
-- 
Sean Davis, M.D., Ph.D.

Clinical Fellow
National Institutes of Health
National Cancer Institute
National Human Genome Research Institute

Clinical Fellow, Johns Hopkins
Department of Pediatric Oncology
--