6
$\begingroup$

I'm performing a relatively simple one-way ANOVA as part of a class exercise. I'd like to get pairwise uncorrected p-values out of R, so I can do a sequential FDR test in another package (I realize there are FDR packages in R as well). I've set up my ANOVA as shown below, and it works fine, produces results, but I can't seem to figure out how to get the raw, uncorrected p-values back. Ultimately, I'd like to do the pairwise tests for both FDR and sequential Bonferroni in R, but this seems like a first step.

It looks like pairwiseCI or multcomp might get me where I'm trying to go, but I'm having a difficult time figuring out which will do what I'm looking for.

R> head(d10)
  time  breed
1 27.4 type.A
2 18.3 type.A
3 24.3 type.B
4 19.6 type.B
5 21.6 type.C
6 30.3 type.D

a10 <- aov(time~breed,data=d10)

# reports the overall significance, but nothing on the pairs  
summary(a10)

# reports corrected values only
TukeyHSD(a10)
$\endgroup$
7
$\begingroup$

For the multcomp package, see the help page for glht; you want to use the "Tukey" option; this does not actually use the Tukey correction, it just sets up all pairwise comparisons. In the example section there's an example that does exactly what you want.

This calculates the estimates and se's for each comparison but doesn't do p-values; for that you need summary.glht; on the help page for that, notice in particular the test parameter, that allows you to set the actual test that is run. For doing multiply-adjusted p-values, you use the adjusted function for this parameter, and to not multiply-adjust, you use test=adjusted("none") (which isn't mention specifically in the help page, though it does say it takes any of the methods in p.adjust, which is where you'd find none.)

You can also compute the estimates and se's by hand using matrix multiplication, and then get the p-values however you want; this is what the glht function is doing behind the scenes. To get the matrices you need to start you'd use coef and vcov.

I didn't put complete code as you say it's for a class project (thanks for being honest, by the way!) and the policy here is to provide helpful hints but not solutions.

$\endgroup$
  • $\begingroup$ thanks for the detailed answer! Re: homework; this is a graduate course and this work isn't being turned in, I know how to perform the same in the class-used JMP but would like to stick to R as much as possible. $\endgroup$ – scw Jun 2 '11 at 18:42
  • $\begingroup$ Thanks for the update about the homework status. If you need more specific help, let us know. $\endgroup$ – Aaron - Reinstate Monica Jun 2 '11 at 19:04
7
$\begingroup$

You can use pairwise.t.test() with one of the available options for multiple comparison correction in the p.adjust.method= argument; see help(p.adjust) for more information on the available option for single-step and step-down methods (e.g., BH for FDR or bonf for Bonferroni). Of note, you can directly give p.adjust() a vector of raw p-values and it will give you the corrected p-values.

So, I would suggest to run something like

pairwise.t.test(time, breed, p.adjust.method ="none")  # uncorrected p-value
pairwise.t.test(time, breed, p.adjust.method ="bonf")  # Bonferroni p-value

The first command gives you t-test-based p-values without controlling for FWER or FDR. You can then use whatever command you like to get corrected p-values.

$\endgroup$
  • 2
    $\begingroup$ By the name, one might think it computes each pairwise t test separately (I did, anyway), but it actually uses the pooled SD by default, so the result is the same as using the anova result. $\endgroup$ – Aaron - Reinstate Monica Jun 2 '11 at 15:12
  • 2
    $\begingroup$ @Aaron You're right, it uses pooled SD by default, but it can be disabled. We could also construct every contrast we want with simple algebra. I like your response btw (+1). $\endgroup$ – chl Jun 2 '11 at 15:38
2
$\begingroup$
library(multcomp)
df = mtcars
df$am = as.factor(df$am)
m1 <- aov(mpg ~ am, data= df)
ht = glht(m1, linfct = mcp(am = "Tukey"))
summary(ht, test = adjusted("none"))
# Linear Hypotheses:
#            Estimate Std. Error t value Pr(>|t|)
# 1 - 0 == 0    7.245      1.764   4.106 0.000285 ***
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.