CODEWITHSUNDEEP
r
Ping Tang
Ping Tang

Reputation: 425

How to conduct pairwise comparison in R like that in SPSS with "multcomp" package

I have three factors: word, type and register. In SPSS, it is very easy to conduct a pairwise comparison (or simple comparison) in SPSS, the syntax is:

/EMMEANS=TABLES(word*register*type) COMPARE(type) ADJ (BONFERRONI)

And it will give me a result like this:

enter image description here

But how can I achieve this in R with Multcomp package? Say if I have a anova result like this: 

library (multcomp)
anova <- aov (data, min~word*register*type)
summary (anova)

How can I do SPSS-like pairwise comparisons?

Upvotes: 0

Views: 2174

Answers (1)

Hack-R
Hack-R

Reputation: 23200

hsb2<-read.table("http://www.ats.ucla.edu/stat/data/hsb2.csv", sep=",", header=T)
attach(hsb2) # attach is a horrible practice, but just roll with it for this example
pairwise.t.test(write, ses, p.adj = "bonf")

That gives you the Bonferonni pairwise comparison that you see in SPSS.

This may help further and in general UCLA provides some good resources that relate commands in SAS, SPSS, Stata, Mplus and R:

http://statistics.ats.ucla.edu/stat/r/faq/posthoc.htm

http://www.ats.ucla.edu/stat/dae/

So, that's for pairwise comparisons. You can also use p.adjust with multiple comparisons (multi-way). See this manual page "Adjust P-values for Multiple Comparisons".

The example they give showing the different adjustments which can be made, including Bonferroni, is: 

require(graphics)

set.seed(123)
x <- rnorm(50, mean = c(rep(0, 25), rep(3, 25)))
p <- 2*pnorm(sort(-abs(x)))

round(p, 3)
round(p.adjust(p), 3)
round(p.adjust(p, "BH"), 3)

## or all of them at once (dropping the "fdr" alias):
p.adjust.M <- p.adjust.methods[p.adjust.methods != "fdr"]
p.adj    <- sapply(p.adjust.M, function(meth) p.adjust(p, meth))
p.adj.60 <- sapply(p.adjust.M, function(meth) p.adjust(p, meth, n = 60))
stopifnot(identical(p.adj[,"none"], p), p.adj <= p.adj.60)
round(p.adj, 3)
## or a bit nicer:
noquote(apply(p.adj, 2, format.pval, digits = 3))


## and a graphic:
matplot(p, p.adj, ylab="p.adjust(p, meth)", type = "l", asp = 1, lty = 1:6,
        main = "P-value adjustments")
legend(0.7, 0.6, p.adjust.M, col = 1:6, lty = 1:6)

## Can work with NA's:
pN <- p; iN <- c(46, 47); pN[iN] <- NA
pN.a <- sapply(p.adjust.M, function(meth) p.adjust(pN, meth))
## The smallest 20 P-values all affected by the NA's :
round((pN.a / p.adj)[1:20, ] , 4)

This page also gives some nice pairwise comparison examples with ANOVA where they have mental, physical, and medical factors and adjust the p-values in various ways.

This is another excellent resource by Bretz.

Finally, here's a SO question on N-way ANOVA in R where they show how to do 3-way and up to 100-way.

Upvotes: 1

Related Questions