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

> 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 --

Michael, I was also struggling with this a few days ago mt.maxT(data[,6:12], c(0,0,0,1,1,1,1), test="t") from multtest gives the raw P values in the 3rd column of the result In my case that gave some concern because there were almost no low values so I basically computed all p-values from the t.test also mytt<-function(x,cl){ t.test(x[cl==0],x[cl==1])$p.val } pvals<-apply(data[,6:12],1,mytt,cl= c(0,0,0,1,1,1,1)) Well, values turned out to be the same when doing this so I had to conclude there was no difference between my groups there Keep in mind that the mt.maxT returns a sorted array with lowest p-values on top Jan > > 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. >