One-way ANCOVA

The following code chunk will produce the same ANCOVA output with the rstatix package. One slight difference here though is that we are going to save the output to an anova_test() object, then print the output with the get_anova_table() function. This approach will come in handy when generating the plot of means. In the anova_test() formula argument we enter the covariate first. This is primarily done to pull out the correct F ratio for plotting. In our example, we are setting the sums of squares type = 3 which means we will get the same output whether we use formula Post ~ Pre + Condition or formula = Post ~ Condition + Pre. However, if using Type II sums of squares enter the covariate first.

library(rstatix)
library(ggpubr)

# Evaluated ANCOVA test w rstatix
aocv.model <- anova_test(formula = Post ~  Pre + Condition, 
                         data = chapter_9_table_7, 
                         type = 3,
                         effect.size = "pes")

# Print ANOVA table
get_anova_table(aocv.model)

## ANOVA Table (type III tests)
## 
##      Effect DFn DFd      F     p p<.05   pes
## 1       Pre   1  26 10.772 0.003     * 0.293
## 2 Condition   2  26  3.732 0.038     * 0.223

Post hoc tests

To conduct the post hoc tests, we will use the emmeans_test() function specifying the formula Post ~ Condition, setting covariate = Pre, and the Bonferroni correction with p.adjust.method = "bonferroni". In the following code chunk, we will also save the output to an object called pwc. The pwc object can be used with the get_emmeans() function to disply the estimated marginal means for each group after accounting for the covariate. Finally, we can print pwc to the console with the print() to display the results of the post hoc tests.

# Pairwise comparisons
pwc <- chapter_9_table_7 %>% 
  emmeans_test(Post ~ Condition, 
               covariate = Pre,
               p.adjust.method = "bonferroni")

# Print estimated marginal means
get_emmeans(pwc) %>% kableExtra::kable()

# Print post hoc tests
#print(pwc)
pwc %>% kableExtra::kable()

Effect sizes

The process to calculate the effect sizes of the post hoc comparisons is a little more drawn out in this case. First, we will need the MBESS package. Next we will need three pieces of information, group means, standard deviations, and sample sizes. We will want the estimated marginal means for this situation because these means have been adjusted to take into consideration the covariate. We can also calculate the standard deviations by multiplying the standard error of each group by the square root of the group sample size. Finally, we can enter those values into the smd() function which stands for the standardized mean difference. A nice feature of the smd() is that we can calculate a standardized mean difference as biased or unbiased. In practice, it is recommended to set Unbiased = TRUE with small samples.

library(MBESS)

# Get sample sizes. For these data, each group had 10 so n <- c(10, 10,10) would have worked,
# but the following code here will generalize to other datasets
n <- chapter_9_table_7 %>% 
  group_by(Condition) %>% 
  get_summary_stats() %>% 
  filter(variable == "Post") %>% 
  select(n)

# Need emmeans from pwc
emm <- get_emmeans(pwc) %>% select(emmean)

# Get standard deviations
std <- get_emmeans(pwc)$se * sqrt(n)

# Group 1 vs group 2
smd(Mean.1 = emm[1,], 
    Mean.2 = emm[2,], 
    s.1 = std[1,], 
    s.2 = std[2,], 
    n.1 = n[1,], 
    n.2 = n[2,],
    Unbiased = FALSE)

##       emmean
## 1 -0.8240768

# Group 1 vs group 3
smd(Mean.1 = emm[1,], 
    Mean.2 = emm[3,], 
    s.1 = std[1,], 
    s.2 = std[3,], 
    n.1 = n[1,], 
    n.2 = n[3,],
    Unbiased = FALSE)

##      emmean
## 1 -1.193841

# Group 2 vs group 3
smd(Mean.1 = emm[2,], 
    Mean.2 = emm[3,], 
    s.1 = std[2,], 
    s.2 = std[3,], 
    n.1 = n[2,], 
    n.2 = n[3,],
    Unbiased = FALSE)

##       emmean
## 1 -0.3694917

Plot of means

While the jmv approach is nice because we can get quite a bit of output (ANCOVA test, post hoc tests, and plots) from one function, the benefit of using rstatix is that we gain some additional features when it comes to plotting. The plot produced for this guide displays the estimated marginal means, confidence intervals, a significance marker, the results of the omnibust test, and includes a caption displaying the tests used to compare means and the multiple correction procedure.

# Add the x and y coordinates for plotting significance markers
pwc <- pwc %>% add_xy_position(x = "Condition", 
                               fun = "mean_se")

# Generate an error plot with significance markers, label, and caption
ggerrorplot(get_emmeans(pwc), 
            x = "Condition", 
            y = "emmean", 
            ylab = "Mean Depression Scores") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) + 
  stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) +
  labs(subtitle = get_test_label(aocv.model, detailed = TRUE), caption = get_pwc_label(pwc))

Figure 2: Plot of post-treatment depression score means by Condition.