R Tutorial Series

By John M Quick

The R Tutorial Series provides a collection of user-friendly tutorials to people who want to learn how to use R for statistical analysis.

My Statistical Analysis with R book is available from Packt Publishing and Amazon.

R Tutorial Series: ANOVA Pairwise Comparison Methods

When we have a statistically significant effect in ANOVA and an independent variable of more than two levels, we typically want to make follow-up comparisons. There are numerous methods for making pairwise comparisons and this tutorial will demonstrate how to execute several different techniques in R.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 30 participants who are divided into three stress reduction treatment groups (mental, physical, and medical). The values are represented on a scale that ranges from 1 to 5. This dataset can be conceptualized as a comparison between three stress treatment programs, one using mental methods, one using physical training, and one using medication. The values represent how effective the treatment programs were at reducing participant's stress levels, with higher numbers indicating higher effectiveness.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.
  1. > #read the dataset into an R variable using the read.csv(file) function
  2. > dataPairwiseComparisons <- read.csv("dataset_ANOVA_OneWayComparisons.csv")
  3. > #display the data
  4. > dataPairwiseComparisons

The first ten rows of our dataset

Omnibus ANOVA

For the purposes of this tutorial, we will assume that the omnibus ANOVA has already been conducted and that the main effect for treatment was statistically significant. For details on this process, see the One-Way ANOVA with Pairwise Comparisons tutorial, which uses the same dataset.

Means

Let's also look at the means of our treatment groups. Here, we will use the tapply() function, along with the following arguments, to generate a table of means.
  • X: the data
  • INDEX: a list() of factor variables
  • FUN: the function to be applied
  1. > #use tapply(X, INDEX, FUN) to generate a table displaying each treatment group mean
  2. > tapply(X = dataPairwiseComparisons$StressReduction, INDEX = list(dataPairwiseComparisons$Treatment), FUN = mean)

The treatment group means

Pairwise Comparisons

We will cover five major techniques for controlling Type I error when making pairwise comparisons. These methods are no adjustment, Bonferroni's adjustment, Holm's adjustment, Fisher's LSD, and Tukey's HSD. All of these techniques will be demonstrated on our sample dataset, although the decision as to which to use in a given situation is left up to the reader.

pairwise.t.test()

Our first three methods will make use of the pairwise.t.test() function, which has the following major arguments.
  • x: the dependent variable
  • g: the independent variable
  • p.adj: the p-value adjustment method used to control for the family-wise Type I error rate across the comparisons; one of "none", "bonferroni", "holm", "hochberg", "hommel", "BH", or "BY"

No Adjustment

Using p.adj = "none" in the pairwise.t.test() function makes no correction for the Type I error rate across the pairwise tests. This technique can be useful for employing methods that are not already built into R functions, such as the Shaffer/Modified Shaffer, which use different alpha level divisors based on the number of levels composing the independent variable. The console results will contain no adjustment, but the researcher can manually consider the statistical significance of the p-values under his or her desired alpha level.
  1. > #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the treatment group means
  2. > #no adjustment
  3. > pairwise.t.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, p.adj = "none")

Pairwise comparisons of treatment group means with no adjustment

With no adjustment, the mental-medical and physical-medical comparisons are statistically significant, whereas the mental-physical comparison is not. This suggests that both the mental and physical treatments are superior to the medical treatment, but that there is insufficient statistical support to distinguish between the mental and physical treatments.

Bonferroni Adjustment

The Bonferroni adjustment simply divides the Type I error rate (.05) by the number of tests (in this case, three). Hence, this method is often considered overly conservative. The Bonferroni adjustment can be made using p.adj = "bonferroni" in the pairwise.t.test() function.
  1. > #Bonferroni adjustment
  2. > pairwise.t.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, p.adj = "bonferroni")

Pairwise comparisons of treatment group means using Bonferroni adjustment

Using the Bonferroni adjustment, only the mental-medical comparison is statistically significant. This suggests that the mental treatment is superior to the medical treatment, but that there is insufficient statistical support to distinguish between the mental and physical treatments and the physical and medical treatments. Notice that these results are more conservative than with no adjustment.

Holm Adjustment

The Holm adjustment sequentially compares the lowest p-value with a Type I error rate that is reduced for each consecutive test. In our case, this means that our first p-value is tested at the .05/3 level (.017), second at the .05/2 level (.025), and third at the .05/1 level (.05). This method is generally considered superior to the Bonferroni adjustment and can be employed using p.adj = "holm" in the pairwise.t.test() function.
  1. > #Holm adjustment
  2. > pairwise.t.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, p.adj = "holm")

Pairwise comparisons of treatment group means using Holm adjustment

Using the Holm procedure, our results are practically (but not mathematically) identical to using no adjustment.

LSD Method

The Fisher Least Significant Difference (LSD) method essentially does not correct for the Type I error rate for multiple comparisons and is generally not recommended relative to other options. However, should the need arise to employ this method, one should seek out the LSD.test() function in the agricolae package, which has the following major arguments.
  • y: the dependent variable
  • trt: the independent variable
  • DFerror: the degrees of freedom error
  • MSerror: the mean squared error
Note that the DFerror and MSerror can be found in the omnibus ANOVA table.
  1. > #load the agricolae package (install first, if necessary)
  2. > library(agricolae)
  3. #LSD method
  4. #use LSD.test(y, trt, DFerror, MSerror) to test the pairwise comparisons between the treatment group means
  5. > LSD.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, 30.5, 1.13)

Pairwise comparisons of treatment group means using LSD method

Using the LSD method, our results are practically (but not mathematically) identical to using no adjustment or the Holm procedure.

HSD Method

The Tukey Honest Significant Difference (HSD) method controls for the Type I error rate across multiple comparisons and is generally considered an acceptable technique. This method can be executed using the TukeyHSD(x) function, where x is a linear model object created using the aov(formula, data) function. Note that in this application, the aov(formula, data) function is identical to the lm(formula, data) that we are already familiar with from linear regression.
  1. > #HSD method
  2. > #use TukeyHSD(x), in tandem with aov(formula, data), to test the pairwise comparisons between the treatment group means
  3. TukeyHSD(aov(StressReduction ~ Treatment, dataPairwiseComparisons))

Pairwise comparisons of treatment group means using HSD method

Using the HSD method, our results are practically (but not mathematically) identical to using the Bonferroni, Holm, or LSD methods.

Complete Pairwise Comparisons Example

To see a complete example of how various pairwise comparison techniques can be applied in R, please download the ANOVA pairwise comparisons example (.txt) file.
Posted by John 10 comments
Labels: , , , , , , , , ,

R Tutorial Series: Two-Way ANOVA with Pairwise Comparisons

By extending our one-way ANOVA procedure, we can test the pairwise comparisons between the levels of several independent variables. This tutorial will demonstrate how to conduct pairwise comparisons in a two-way ANOVA.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 27 participants who are divided into three stress reduction treatment groups (mental, physical, and medical) and three age groups (young, mid, and old). The stress reduction values are represented on a scale that ranges from 0 to 10. This dataset can be conceptualized as a comparison between three stress treatment programs, one using mental methods, one using physical training, and one using medication, across three age groups. The stress reduction values represent how effective the treatment programs were at reducing participant's stress levels, with higher numbers indicating higher effectiveness. Note that the numbers in this dataset are not very realistic and are simply used to make this example possible.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.
  1. > #read the dataset into an R variable using the read.csv(file) function
  2. > dataTwoWayComparisons <- read.csv("dataset_ANOVA_TwoWayComparisons.csv")
  3. > #display the data
  4. > dataTwoWayComparisons

The first ten rows of our dataset.

Omnibus Test

Let's run a general omnibus test to assess the main effects and interactions present in the dataset.
  1. > #use anova(object) to test the omnibus hypothesis
  2. > #Are main effects or interaction effects present in the independent variables?
  3. > anova(lm(StressReduction ~ Treatment * Age, dataTwoWayComparisons))

The omnibus ANOVA test

Pairwise Comparisons

Since the omnibus test was significant for both variables and no interaction effect was present, we can proceed to testing the main effect pairwise comparisons. To accomplish this, we will apply our pairwise.t.test() function to each of our independent variables. For more details on the pairwise.t.test() function, see the One-Way ANOVA with Pairwise Comparisons tutorial.
  1. > #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the treatment group means
  2. > #What significant differences are present amongst the treatment means?
  3. > pairwise.t.test(dataTwoWayComparisons$StressReduction, dataTwoWayComparisons$Treatment, p.adj = "none")
  4. > #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the age group means
  5. > #What significant differences are present amongst the age group means?
  6. > pairwise.t.test(dataTwoWayComparisons$StressReduction, dataTwoWayComparisons$Age, p.adj = "none")

Pairwise comparisons of treatment group means


Pairwise comparisons of age group means

Note that the desired p-adjustment method will vary by researcher, study, etc. Here, we will assume an alpha level of .05 for all tests, effectively making no adjustment for the family-wise Type I error rate.
These results indicate that there are are no statistically significant pairwise differences between the treatment groups and that all of the comparisons between age groups are statistically significant. The age group means are 8 for young, 5 for mid, and 2 for old. Consequently, we are inclined to conclude that, regardless of treatment, young patients are going to be most responsive, followed by middle aged patients, followed by older ones. However, there is insufficient support to differentiate between the effectiveness of the treatment methods themselves.

Complete Two-Way ANOVA with Pairwise Comparisons Example

To see a complete example of how two-way ANOVA pairwise comparisons can be conducted in R, please download the two-way ANOVA comparisons example (.txt) file.
Posted by John 3 comments
Labels: , , , , , ,

R Tutorial Series: One-Way ANOVA with Pairwise Comparisons

When we have more than two groups in a one-way ANOVA, we typically want to statistically assess the differences between each group. Whereas a one-way omnibus ANOVA assesses whether a significant difference exists at all amongst the groups, pairwise comparisons can be used to determine which group differences are statistically significant.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 30 participants who are divided into three stress reduction treatment groups (mental, physical, and medical). The values are represented on a scale that ranges from 1 to 5. This dataset can be conceptualized as a comparison between three stress treatment programs, one using mental methods, one using physical training, and one using medication. The values represent how effective the treatment programs were at reducing participant's stress levels, with higher numbers indicating higher effectiveness.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.
  1. > #read the dataset into an R variable using the read.csv(file) function
  2. > dataOneWayComparisons <- read.csv("dataset_ANOVA_OneWayComparisons.csv")
  3. > #display the data
  4. > dataOneWayComparisons

The first ten rows of our dataset.

Omnibus Test

One way to begin an ANOVA is to run a general omnibus test. The advantage to starting here is that if the omnibus test comes up insignificant, you can stop your analysis and deem all pairwise comparisons insignificant. If the omnibus test is significant, you should continue with pairwise comparisons.
  1. > #use anova(object) to test the omnibus hypothesis
  2. > #Is there a significant difference amongst the treatment means?
  3. > anova(lm(StressReduction ~ Treatment, dataOneWayComparisons))
The omnibus ANOVA test

Pairwise Comparisons

Since the omnibus test was significant, we are safe to continue with our pairwise comparisons. To make pairwise comparisons between the treatment groups, we will use the pairwise.t.test() function, which has the following major arguments.
  • x: the dependent variable
  • g: the independent variable
  • p.adj: the p-value adjustment method used to control for the family-wise Type I error rate across the comparisons; one of "none", "bonferroni", "holm", "hochberg", "hommel", "BH", or "BY"
The pairwise.t.test() function can be implemented as follows.
  1. > #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the treatment group means
  2. > #What significant differences are present amongst the treatment means?
  3. > pairwise.t.test(dataOneWayComparisons$StressReduction, dataOneWayComparisons$Treatment, p.adj = "none")

Pairwise comparisons of treatment group means

Note that the desired p-adjustment method will vary by researcher, study, etc. Here, we will assume an alpha level of .05 for all tests, effectively making no adjustment for the family-wise Type I error rate.
These results indicate that there is a statistically significant difference between the mental and medical (p = .004) and physical and medical (p = 0.045), but not the mental and physical (p = 0.302) treatments. The treatment means are 3.5 for mental, 3 for physical, and 2 for medical. Subsequently, we are inclined to conclude based on this study that the mental and physical treatments lead to greater stress reduction than the medical method, but that there is insufficient statistical support to determine that either the mental or physical treatment method is superior.

Complete One-Way ANOVA with Pairwise Comparisons Example

To see a complete example of how one-way ANOVA pairwise comparisons can be conducted in R, please download the one-way ANOVA comparisons example (.txt) file.

References

R Documentation (n.d.). Adjust P-values for Multiple Comparisons. Retrieved January 16, 2011 from http://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html
Posted by John 2 comments
Labels: , , , , , ,
Subscribe to: Posts (Atom)