[BioC] multtest Paired T-test

[Katherine Pollard]

Ken recently posted a question to the list regarding computation of t- and
F-statistics (and corresponding unadjusted p-values) in multtest. Here are
some suggestions.

Ken's data set consists of the following:

m: a 22277x6 matrix of probeset expressions
m.cl=c(0,0,0,1,1,1): the corresponding class label
m_pt: the 22277x6 matrix m with columns rearranged for paired tests
m_pt.cl=c(0,1,0,1,0,1): the corresponding class label for paired tests

(presumably m_pt=m[,c(1,4,2,5,3,6)])

He is doing two sample tests, but wants to run the F tests (typically
used to compare k>2 groups) for practice.

Ken computes the test statistics for each possible test available using
the function mt.teststat: regular and rank-based (nonpara argument) for
each value of the test argument.

teststat_t_para <- mt.teststat(m,m.cl,test="t",nonpara="n")
teststat_t_nonpara <- mt.teststat(m,m.cl,test="t",nonpara="y")
teststat_te_para <- mt.teststat(m,m.cl,test="t.equalvar",nonpara="n")
teststat_te_nonpara <- mt.teststat(m,m.cl,test="t.equalvar",nonpara="y")
teststat_w_para <- mt.teststat(m,m.cl,test="wilcoxon",nonpara="n")
teststat_w_nonpara <- mt.teststat(m,m.cl,test="wilcoxon",nonpara="y")
teststat_f_para <- mt.teststat(m,m.cl,test="f",nonpara="n")
teststat_f_nonpara <- mt.teststat(m,m.cl,test="f",nonpara="y")
teststat_pt_para <- mt.teststat(m_pt,m_pt.cl,test="pairt",nonpara="n")
teststat_pt_nonpara <- mt.teststat(m_pt,m_pt.cl,test="pairt",nonpara="y")
teststat_bf_para <- mt.teststat(m_pt,m_pt.cl,test="blockf",nonpara="n")
teststat_bf_nonpara <- mt.teststat(m_pt,m_pt.cl,test="blockf",nonpara="y")

These look OK, except that for the "blockf" tests, Ken might want to use m
and m_bf.cl=c(0,1,2,0,1,2), since presumably he is thinking about the
two groups as two blocks (each with 3 observations). As written, with
m_pt, the "blockf" test corresponds to looking at three blocks each with
two observations. Note that in the way Ken has this now, the square of the
paired t-statistic is equal to the block F-statistic (for three blocks of
two observations).

Next, Ken wants to use tabled distributions (with the functions pt,
pwilcox, and pf) to compute marginal, unadjusted p-values for each gene
based on the above test statistics. Assuming two-tailed tests, he
proposes the following (which isn't quite right):

rawp0_teststat_t_para <- 2 * (1 - pt(abs(teststat_t_para), df =5))
rawp0_teststat_t_nonpara <- 2 * (1 - pt(abs(teststat_t_nonpara)), df=5)
rawp0_teststat_te_para <- 2 * (1 - pt(abs(teststat_te_para)), df=5)
rawp0_teststat_te_nonpara <- 2 * (1 - pt(abs(teststat_te_nonpara)), df=5)
rawp0_teststat_w_para <- 2 * (1 - pwilcox(abs(teststat_w_para)), df=5)
rawp0_teststat_w_nonpara <- 2 * (1 - pwilcox(abs(teststat_w_nonpara)), df=5)
rawp0_teststat_f_para <- 2 * (1 - pf(abs(teststat_f_para)), df=5)
rawp0_teststat_f_nonpara <- 2 * (1 - pf(abs(teststat_f_nonpara)), df=5)
rawp0_teststat_pt_para <- 2 * (1 - pt(abs(teststat_pt_para)), df=5)
rawp0_teststat_pt_nonpara <- 2 * (1 - pt(abs(teststat_pt_nonpara)), df=5)
rawp0_teststat_bf_para <- 2 * (1 - pf(abs(teststat_bf_para)), df=5)
rawp0_teststat_bf_nonpara <- 2 * (1 - pf(abs(teststat_bf_nonpara)), df=5)

Here are some suggetions:

1. Typo: For the nonpara ones, the parens are in the wrong place. They
should be the same as for the para ones, eg:

rawp0_teststat_t_nonpara <- 2 * (1 - pt(abs(teststat_t_nonpara), df=5))

2. Degrees of freedom: For t-tests with equal variance (t.equalvar),
df=n1+n2-2=4. For unequal variances, the df formula is more complicated
(Welch - or Satterthwaite - approximation). For the paired t-test, df=2.
For F-tests, you need two dfs (see ? pf).
For pwilcox, you just give sample sizes in each group, not df (see ?
pwilcox).

3. You don't need abs() for F-tests: the statistics are always positive
since they are sums of squares. (It doesn't hurt though).

With the correct rawp0 values in hand, Ken could use the function
mt.rawp2adjp to adjust for multiple comparisons.

Let me add that while this is a good exercise in understanding marginal
p-values, the multtest package is designed to compute these without using
these tabled distributions (and making the corresponding assumptions). The
functions MTP and mt.minP (or mt.maxT) provide p-values (raw and
adjusted) through bootstrap and permutation resampling. These p-values
account for correlations between genes and do not make assumptions about
the parametric forms of the distributions for each test statistic. Hence,
one does not have to guess degrees of freedom or distributions. The imput
is simply the data and class labels.

Best,
Katie