Questions about multtest
2
0
Entering edit mode
@michael-watson-iah-c-378
Last seen 9.6 years ago
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
multtest multtest • 1.2k views
ADD COMMENT
0
Entering edit mode
@sean-davis-490
Last seen 12 weeks ago
United States
> 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 --
ADD COMMENT
0
Entering edit mode
@oosting-j-path-412
Last seen 9.6 years ago
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. >
ADD COMMENT
0
Entering edit mode
On Tue, 24 Feb 2004, Oosting, J. (PATH) wrote: <snip> > 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 it should be pointed out that with just 3 reps the t-test is not a very powerful test (low sensitivity). you might consider alternatives, such as the one performed by SAM (in the siggenes library) or one of the many implemented in limma. > > > 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. > > > > _______________________________________________ > Bioconductor mailing list > Bioconductor@stat.math.ethz.ch > https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor >
ADD REPLY

Login before adding your answer.

Traffic: 706 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6