[BioC] Questions about multtest
michael watson (IAH-C)
michael.watson at bbsrc.ac.uk
Tue Feb 24 15:23:06 MET 2004
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)))
However, rawp0:
> rawp0[1]
[1] 0
So according to this procedure, my raw (unadjusted) p-value is zero, whereas according to t.test(), my p-value is 0.0006061.
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
Does this mean that my adjusted p-value from the Benjamini and Hochberg procedure for my first gene is zero? Haven't I lost some information along the way when 0.0006061 was converted to just plain old zero? How do I get accurate raw p-values from the multtest package?
Or am I doing something wrong?
Thanks
Mick
More information about the Bioconductor
mailing list