**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**.

In the **previous tutorial**, we computed marginal means 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 analogous to the ** post hoc** analysis (pretty much consisting of pairwise

Let’s carry out contrast analysis on the simple model from the previous tutorial:

```
library(ggplot2)
library(see)
library(rstanarm)
library(modelbased)
# you can ignore the `refresh` argument, it just prevents the function from printing
# a few details to the console
model <- stan_glm(Sepal.Width ~ Species, data = iris, refresh = 0)
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:

`estimate_contrasts(model)`

```
> Level1 | Level2 | Difference | 95% CI | pd | % in ROPE | Std_Difference
> -------------------------------------------------------------------------------------------
> 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
```

If you can recall your Bayesian apprenticeship, we can conclude that all pairwise differences are statistically significant.

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 | Std_Difference
> ------------------------------------------------------------------------------------------
> 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:

```
library(see)
plot(estimate_contrasts(model), estimate_means(model)) +
theme_modern()
```

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).

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.

Contrast analysis can be a powerful tool to interpret and understand statistical models.