Two-way ANOVA

To perform a two-way ANOVA with the rstatix package, the Feedback and Drug variables first need to be converted to factor so the numerical values get treated as a categorical variable. Next, the anova_test() function can be used to conduct the test and display the results. However, it may be advantageous to first fit an aov() model that can later be used with the emmeans_test() function.

library(tidyverse)
library(rstatix)
library(ggpubr)

# Convert Drug and Feedback to factor
chapter_7_table_5 <- chapter_7_table_5 %>%
  mutate_at(c("Drug", "Feedback"), as.factor)

# Create aov() model
model <- aov(formula = Score ~ Feedback * Drug,
    data = chapter_7_table_5)

# anova_test on  model, with partial eta squared and type III ss
anova_test(model, effect.size = "pes", type = 3)

## ANOVA Table (type III tests)
## 
##          Effect DFn DFd      F        p p<.05   pes
## 1      Feedback   1  24  6.934 0.015000     * 0.224
## 2          Drug   2  24 10.979 0.000411     * 0.478
## 3 Feedback:Drug   2  24  2.504 0.103000       0.173

# Alternative form of using anova_test()
# anova_test(Score ~ Feedback * Drug,
#            data = chapter_7_table_5,
#            effect.size = "pes",
#            type = 3)

The results from the omnibus test reveal that there is no significant interaction between Drug and Feedback. On the other hand, there are significant main effects for Drug and for Feedback. Because the interaction is not significant, we may proceed to perform tests of marginal means. Had the interaction been significant, we could have opted to perform tests of simple effects of Drug within Feedback, or Feedback within drug. However, these are primarily suggestions and the approach to analyzing the data should be guided by the research question.

Tests of marginal means

The rstatix package includes a function, emmeans_test(), that can perform tests of estimated marginal means. To perform these tests, we will use the aov model that we created in the previous step and we will conduct two separate tests, one for Feedback and one for Drug. The correction method in this example is set to "none", but this can be easily changed according to your situation. The available correction methods can be found by typing help(anova_test) in the Console.

We will also use the get_summary_stats() function and illustrate how the tidyverse package meshes with rstatix to produce descriptive statistics to aid in the interpretation of the output. First we will start with the chapter_7_table_5 data and pipe it to the group_by() function which will create subsets of data according to combinations of Drug and Feedback, then that output is fed into the get_summary_stats() which will calculate several descriptive statistics for each group. Finally, to display only a part of this output, the select() function will display data by specifying column names.

To plot the data we will rely on the ggpubr package. The ggpubr and rstatix packages are developed by the same individual and are aimed at simplifying the syntax for conducting statistical tests and generating plots in R. For our purposes we will use the ggerrorplot() function to plot means and confidence intervals.

Feedback

# Test of marginal means for Feedback
pwc <- chapter_7_table_5 %>% 
  emmeans_test(Score ~ Feedback, 
               p.adjust.method = "none", 
               model = model)

## NOTE: Results may be misleading due to involvement in interactions

# Print the estimated marginal means of Feedback
get_emmeans(pwc) %>% kableExtra::kable(., "html")

# Modify pwc to include x and y positions for plotting signficance markers
pwc <- pwc %>% add_xy_position(x = "Feedback", fun = "mean_ci")

# Generate an error plot uses ggpubr AND ggplot functions to layer the plot
ggerrorplot(get_emmeans(pwc), 
       x = "Feedback", 
       y = "emmean",
       ylab = "Mean BP Scores") +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
                  width = 0.1) +
    stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE)

Plot cell means

#Produce an error plot with cell means and ci instead of estimated marginal means
#The same code chunk can be used to get cell means of Drug by changing the "x = " argument
ggerrorplot(chapter_7_table_5,
       x = "Feedback",
       y = "Score",
       add = "mean",
       desc_stat = "mean_ci",
       error.plot = "errorbar",
       width = .1)

#Print cell mean, standard deviation, and standard error for each level of Feedback
chapter_7_table_5 %>%
  group_by(Feedback) %>%
  get_summary_stats(Score) %>%
  select(Feedback, mean, sd, se)

Drug

# Test of marginal means for Drug
pwc <- chapter_7_table_5 %>% 
  emmeans_test(Score ~ Drug, 
               p.adjust.method = "none", 
               model = model)

## NOTE: Results may be misleading due to involvement in interactions

# Print the results of the pairwise comparisons of Feedback
pwc %>% kableExtra::kable(., "html")

# Print the estimated marginal means of Feedback
get_emmeans(pwc) %>% kableExtra::kable(., "html")

# Modify pwc to include x and y positions for plotting signficance markers
pwc <- pwc %>% add_xy_position(x = "Feedback", fun = "mean_ci")

# Generate an error plot uses ggpubr AND ggplot functions to layer the plot
ggerrorplot(get_emmeans(pwc), 
       x = "Drug", 
       y = "emmean",
       ylab = "Mean BP Scores") +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
                  width = 0.1) +
    stat_pvalue_manual(pwc, hide.ns = FALSE, tip.length = FALSE)

Figure 3: Means and confidence intervals for Drug collapsed across Feedback. *, p<0.05; **, p<0.01; ***, p<0.001; ns, not significant

# Generate an interaction plot
ggline(chapter_7_table_5, 
       x = "Feedback", 
       y = "Score",
       color = "Drug",
       add = "mean_se",
       palette = "jama",
       position = position_dodge(.2))

The two tests of marginal means will produce a couple of messages to remind us that the results could be misleading because of interactions. However, as we saw in the omnibus test, there was no significant interaction and can disregard the messages. The results of the emmeans_test() on Feedback indicate a significant effect and suggests that the participants undergoing biofeedback had lower bloodpressure scores than those without biofeedback. When examining the output of the results of the emmeans_test() on Drug, we see that the mean of Drug 1 is significantly lower than that of Drug 2 and Drug 3. However, the mean of Drug 2 is not significantly different than that of Drug 3.