Contrast analysis

Source: vignettes/estimate_contrasts.Rmd

Warning: we will go full Bayesian in this vignette. If you’re not familiar with the Bayesian framework, we recommend starting with this gentle introduction.

Testing pairwise differences

In the previous tutorial, we computed marginal at the 3 different Species levels from the iris dataset. However, one might also want to statistically test the differences between each levels, which can be achieved through contrast analysis. Although the procedure is much more powerful, its aim is somehow analog to the post hoc analysis (pretty much consisting of pairwise t-tests) heavily used in psychological science to palliate the uselessness of ANOVAs.

Let’s do that based on the simple model from the previous tutorial:

library(ggplot2)
library(see)
library(rstanarm)
library(modelbased)


model <- stan_glm(Sepal.Width ~ Species, data = iris)
means <- estimate_means(model)
ggplot(iris, aes(x = Species, y = Sepal.Width, fill = Species)) +
  geom_violin() +
  geom_jitter2(width = 0.05, alpha = 0.5) +
  
  geom_line(data = means, aes(y = Mean, group = 1), size = 1) +
  geom_pointrange(data = means, aes(y = Mean, ymin = CI_low, ymax = CI_high), size = 1, color = "white") +
  theme_modern()

Contrast analysis can be achieved through the estimate_contrasts function:

> Level1     |     Level2 | Difference |         95% CI |     pd | % in ROPE | Difference (std.)
> ----------------------------------------------------------------------------------------------
> setosa     | versicolor |       0.66 | [ 0.52,  0.79] |   100% |        0% |              1.50
> setosa     |  virginica |       0.45 | [ 0.32,  0.59] |   100% |        0% |              1.04
> versicolor |  virginica |      -0.20 | [-0.34, -0.08] | 99.88% |     6.33% |             -0.46

As we can see here, all pairwise differences can be considered as significant.

Complex model

Again, as contrast analysis is based on marginal means, it can be applied to more complex models:

model <- stan_glm(Sepal.Width ~ Species * Petal.Width, data = iris)
estimate_contrasts(model)
> Level1     |     Level2 | Difference |        95% CI |     pd | % in ROPE | Difference (std.)
> ---------------------------------------------------------------------------------------------
> setosa     | versicolor |       1.52 | [ 0.80, 2.23] |   100% |        0% |              3.48
> setosa     |  virginica |       1.70 | [ 0.98, 2.50] |   100% |        0% |              3.89
> versicolor |  virginica |       0.17 | [-0.11, 0.43] | 88.20% |    27.30% |              0.39

For instance, if we add Petal.Width in the model, we can see that the difference between versicolor and virginica becomes not significant (and even changes sign).

Note that we can plot simple contrast analysis through lighthouse plots with the help of the see package:

These represent the estimated means and their CI range (in black), while the gray areas show the CI range of the difference (as compared to the point estimate).

Changes in difference

Interestingly, we can also see how these differences are modulated by another continuous variable. Based on the model above (including the interaction with Petal.Width), we will compute the contrasts at 100 equally-spaced points of Petal.Width, that we will then visualise.

contrasts <- estimate_contrasts(model, modulate = "Petal.Width", length = 100)

# Create a variable with the two levels concatenated 
contrasts$Contrast <- paste(contrasts$Level1, "-", contrasts$Level2) 

# Visualise the changes in the differences
ggplot(contrasts, aes(x = Petal.Width, y = Difference)) +
  geom_ribbon(aes(fill = Contrast, ymin = CI_low, ymax = CI_high), alpha = 0.2) +
  geom_line(aes(colour = Contrast), size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_modern() +
  ylab("Difference")

As we can see, the difference between versicolor and virginica increases as Petal.Width increases. In conclusion, contrast analysis is a powerful tool to interpret and understand statistical models.

References