Skip to contents

Basic Idea

In this post, I will document how I conduct planned contrasts analyses to test whether experimental groups differ from each other. Planned contrasts are similar to t tests, but provide more power when you have several groups:

Statistical power is lower with the standard t test compared than it is with the planned contrast version for two reasons: a) the sample size is smaller with the t test, because only the cases in the two groups are selected; and b) in the planned contrast the error term is smaller than it is with the standard t test because it is based on all the cases (source)

I will first demonstrate a simple example for conducting the analysis, exporting the table of results to Microsoft Word, and producing the figure. After demonstrating it for a single variable, we will then look at a full workflow, i.e., how we can automatize the process for several variables simultaneously.

Getting started

For the minimal example, we will use the iris dataset, which is installed with R by default.

Using nice_contrasts()

Load the rempsyc package:

Note: If you haven’t installed this package yet, you will need to install it via the following command: install.packages("rempsyc").

Let’s test out the function on existing data

table.stats <- nice_contrasts(response = "Sepal.Length",
                              group = "Species",
                              data = iris)
table.stats
##   Dependent Variable             Comparison  df          t            p
## 1       Sepal.Length     setosa - virginica 147 -15.365506 2.214821e-32
## 2       Sepal.Length versicolor - virginica 147  -6.332686 2.765638e-09
## 3       Sepal.Length    setosa - versicolor 147  -9.032819 8.770194e-16
##          dR  CI_lower   CI_upper
## 1 -2.929466 -3.557067 -2.2204369
## 2 -1.205879 -1.690957 -0.6616629
## 3 -1.723587 -2.165855 -1.2561530

This will give us the ‘robust’ version of Cohen’s d (i.e., robust to deviations from parametric assumptions). It will also provide us with bootstrapped 95% confidence intervals of that effect size.

Make it publication-ready using nice_table().

(my_table <- nice_table(table.stats))

Save the table to Word

save_as_docx(my_table, path = "contrasts.docx")

Make a violin plot comparing groups through nice_violin().

(figure <- nice_violin(group = "Species", 
                       response = "Sepal.Length", 
                       data = iris,
                       ytitle = "Length of Sepal", 
                       signif_annotation = c("***", "***", "***"),
                       signif_yposition = c(8.7, 7.3, 8.2),
                       signif_xmin = c("setosa", "setosa", "versicolor"),
                       signif_xmax = c("virginica", "versicolor", "virginica")))

Let’s save a high-resolution/vector version of this figure in PDF (it could be .png or any other format too if necessary).

ggplot2::ggsave('Figure 1.pdf', figure, width = 7, height = 7, 
                unit = 'in', dpi = 300)

Congratulations for making it that far! You’ve done it! From A to Z, you now know how to make things work! If you would like to integrate this with a broader workflow, read on for the next part!

Full Workflow

For the full workflow, we will use the dataset from one of my studies published at the Quarterly Journal of Experimental Psychology and available on the Open Science Framework.

Load the data file from the Open Science Framework

data <- read.csv("https://osf.io/qkmnp//?action=download", header=TRUE)

Specify the order of factor levels for “Group” (otherwise R will alphabetize them)

(data$Group <- factor(data$Group, levels = c("Embodied", "Mental", "Control")))
##  [1] Embodied Embodied Control  Control  Mental   Embodied Control  Control 
##  [9] Control  Embodied Mental   Mental   Control  Control  Mental   Control 
## [17] Mental   Mental   Embodied Embodied Control  Embodied Mental   Control 
## [25] Embodied Embodied Mental   Control  Mental   Embodied Embodied Mental  
## [33] Mental   Mental   Embodied Mental   Embodied Mental   Control  Embodied
## [41] Mental   Mental   Embodied Mental   Mental   Embodied Control  Control 
## [49] Control  Mental   Embodied Control  Embodied Control  Mental   Embodied
## [57] Embodied Embodied Control  Control  Control  Control  Mental   Control 
## [65] Mental   Embodied Embodied Mental   Embodied Control  Control  Mental  
## [73] Mental   Control  Control  Embodied Mental   Control  Embodied Embodied
## [81] Embodied Mental   Control  Control  Embodied Mental   Mental   Mental  
## [89] Embodied Control 
## Levels: Embodied Mental Control

Define our dependent variables

library(dplyr)
(DV <- data %>% select(QCAEPR:IOS) %>% names)
## [1] "QCAEPR" "IRIPT"  "IRIFS"  "IRIEC"  "IRIPD"  "IOS"
table.stats <- nice_contrasts(response = DV,
                              group = "Group",
                              data = data)
table.stats
##    Dependent Variable         Comparison df           t            p
## 1              QCAEPR Embodied - Control 87  2.19599050 3.075313e-02
## 2              QCAEPR   Mental - Control 87  0.23740438 8.129013e-01
## 3              QCAEPR  Embodied - Mental 87  1.95858612 5.336471e-02
## 4               IRIPT Embodied - Control 86  0.12249596 9.027921e-01
## 5               IRIPT   Mental - Control 86  0.01932946 9.846231e-01
## 6               IRIPT  Embodied - Mental 86  0.10212396 9.188960e-01
## 7               IRIFS Embodied - Control 86  0.15139155 8.800214e-01
## 8               IRIFS   Mental - Control 86 -1.38629687 1.692400e-01
## 9               IRIFS  Embodied - Mental 86  1.53639996 1.281114e-01
## 10              IRIEC Embodied - Control 86  2.30118373 2.379970e-02
## 11              IRIEC   Mental - Control 86  1.44763512 1.513548e-01
## 12              IRIEC  Embodied - Mental 86  0.83396371 4.066119e-01
## 13              IRIPD Embodied - Control 86  2.13034255 3.599933e-02
## 14              IRIPD   Mental - Control 86  0.85841114 3.930517e-01
## 15              IRIPD  Embodied - Mental 86  1.25380051 2.133121e-01
## 16                IOS Embodied - Control 87  4.70778837 9.385982e-06
## 17                IOS   Mental - Control 87  2.12312024 3.658408e-02
## 18                IOS  Embodied - Mental 87  2.58466812 1.141344e-02
##             dR    CI_lower  CI_upper
## 1   0.52518512 -0.06478346 1.1231293
## 2   0.03643629 -0.52454848 0.6181501
## 3   0.48874884 -0.01956104 1.0168335
## 4  -0.06831227 -0.69830600 0.4812051
## 5  -0.17448179 -0.71742146 0.2895616
## 6   0.10616953 -0.43394443 0.7091485
## 7   0.09489742 -0.42082445 0.6489127
## 8  -0.35461669 -0.86858844 0.1598174
## 9   0.44951412 -0.12606486 1.1265277
## 10  0.58908525  0.01124882 1.2468039
## 11  0.42094553 -0.11215480 0.9577911
## 12  0.16813972 -0.37170209 0.7579083
## 13  0.47999336 -0.06202268 1.0883890
## 14  0.18562672 -0.28074287 0.6598243
## 15  0.29436664 -0.26676217 0.9764130
## 16  1.08839542  0.49749276 1.8495629
## 17  0.49918888 -0.03121046 1.0386874
## 18  0.58920655  0.01388086 1.2115056

Add more meaningful names for measures

table.stats[1] <- rep(c("Peripheral Responsivity (QCAE)",
                        "Perspective-Taking (IRI)", 
                        "Fantasy (IRI)", "Empathic Concern (IRI)", 
                        "Personal Distress (IRI)",
                        "Inclusion of Other in the Self (IOS)"), each = 3)

Make it publication-ready using nice_table().

Note that with option highlight = TRUE, it will automatically highlight significant results.

(my_table <- nice_table(table.stats, highlight = TRUE))

Save the table to Word

save_as_docx(my_table, path = "contrasts.docx")

Make violin plots comparing groups through nice_violin.

We will have to make the four figures for which comparisons were significant and then combine them in one plot.

This function will throw an error if your dataset contains missing data. We will omit those rows missing data for now.

Data <- na.omit(data)

Make figure panel 1

(EC <- nice_violin(group = "Group",
                   response = "IRIEC", 
                   data = Data,
                   colours = c("#00BA38", "#619CFF", "#F8766D"), 
                   ytitle = "Empathic Concern (from IRI)", 
                   signif_annotation = "*", 
                   signif_yposition = 5.2, 
                   signif_xmin = 1, 
                   signif_xmax = 3))

Make figure panel 2

(PD <- nice_violin(group = "Group", 
                   response = "IRIPD", 
                   data = Data,
                   colours = c("#00BA38", "#619CFF", "#F8766D"), 
                   ytitle = "Personal Distress (from IRI)", 
                   signif_annotation = "*", 
                   signif_yposition = 5, 
                   signif_xmin = 1,
                   signif_xmax = 3))

Make figure panel 3

(PR <- nice_violin(group = "Group", 
                   response = "QCAEPR", 
                   data = Data,
                   colours = c("#00BA38", "#619CFF", "#F8766D"), 
                   ytitle = "Peripheral Responsivity (from QCAE)", 
                   signif_annotation="*",
                   signif_yposition=4.2,
                   signif_xmin=1,
                   signif_xmax=3))

Make figure panel 4

(IOS <- nice_violin(group = "Group", 
                    response = "IOS", 
                    data = Data,
                    colours = c("#00BA38", "#619CFF", "#F8766D"), 
                    ytitle = "Self-Other Merging (from IOS)", 
                    signif_annotation=c("***", "*", "*"),
                    signif_yposition=c(8.25, 7.5, 6.75),
                    signif_xmin=c(1,1,2),
                    signif_xmax=c(3,2,3)))

It’s now time to combine our four plots into one figure! Yeah!

library(ggpubr)
(figure <- ggarrange(EC, PD, PR, IOS,
                     labels = "AUTO",
                     ncol = 2, nrow = 2))

Let’s save a high-resolution/vector version of this figure in PDF (it could be .png or any other format too if necessary).

ggplot2::ggsave('Figure 1.pdf', figure, width = 14, height = 14, 
                unit = 'in', dpi = 300)

Thanks for checking in

Make sure to check out this page again if you use the code after a time or if you encounter errors, as I periodically update or improve the code. Feel free to contact me for comments, questions, or requests to improve this function at https://github.com/rempsyc/rempsyc/issues. See all tutorials here: https://remi-theriault.com/tutorials.