Welch's one-way ANOVA

In rstatix, the anova_test() function can analyse several types of between and within subjects ANOVA designs. However, to conduct a Welch’s ANOVA , the welch_anova_test() function is required. The syntax to conduct the actual test is rather simple. In the following code chunk, I’ve chosen to explicitly declare formula = Scores ~ Group and data = C3E9. However, one could just as easily specify the dataframe and formula as long as data is the first entered variable. Similarly, data and formula do not need to be explicitly specified for the games_howell_test() function, but hare specified in the example for consistency.

Check normality

To produce the Shapiro-Wilk’s test of normality as in the jmv output, we will need to create an analysis of variance (aov) object with the base R aov() function. The aov() function is actually a wrapper for the lm() which highlights the relationship between linear regression and ANOVA. The rstatix package does contain its own shapiro_test() function, but it will not test normality of the residuals, only the actual values. In the code chunk below, the aov object is piped to the residuals() function, which is then piped to the shapiro.test() function. This will conduct Shapiro-Wilk’s test on the residuals of the aov object and not the values of the dependent variable. Our Shaprio-Wilk’s test result, W = 0.81 is significant with a p-values less than 0.05 which is taken as evidence that our data violate the assumption of normality.

In addition to the Shapiro-Wilk test, we can visualize the residuals of our data and plot them against the expected residuals of a normal distribution. If our data are normally distributed, we would expect the individual data points to hover near the diagonal line. As we can see in the plot, we have quite a few data points fall far away from the line which is additional evidence that our data are not normally distributed. When normality is violated, the Welch’s ANOVA is one option to analyse the data.

# Convert group to factor
C3E9$Group <- as.factor(C3E9$Group)

# Base R Shapiro-Wilk test on residuals of the aov object
aov(Scores ~ Group, data = C3E9) %>% 
  residuals() %>% 
  shapiro.test()

## 
## 	Shapiro-Wilk normality test
## 
## data:  .
## W = 0.81079, p-value = 0.01246

# Create a plot of standardised residuals, indexed at position 2 of plot(aov(x))
plot(aov(formula = Scores ~ Group, data = C3E9), 2)

Figure 1: Q-Q Plot

Check homogeneity of variance

Finally, for the homogeneity of variance test, the rstatix levene_test() function specifying a formula and a data frame does the job. The result is a non-significant p-value for the test statistic. This indicates that the data meet the assumption of homogeneity of variance.

As in the examination of normality above, base R also can also plot the Residuals vs Fitted values to examine homogeneity of variance. The plot maintains a straight red line which is what would be expected for data that meet the homogeneity of variance assumption. When the assumption of homogeneity of variance is violated one can explore the use of a robust ANOVA.

library(rstatix)

# Residuals vs Fitted values plot
plot(aov(formula = Scores ~ Group, data = C3E9), 1)

Figure 2: Residuals vs Fitted Values

# rstatix Levene's test for homogeneity of variance
levene_test(formula = Scores ~ Group,
            data = C3E9)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3     8  6.93e-33     1

Welch’s ANOVA

Considering that the data failed to meet the assumption of normality, but met the assumption of homogeneity of variance, we can proceed to conduct a Welch’s ANOVA.

# Conduct Welch's ANOVA
welch_anova_test(formula = Scores ~ Group,
                 data = C3E9) 

## # A tibble: 1 × 7
##   .y.        n statistic   DFn   DFd     p method     
## * <chr>  <int>     <dbl> <dbl> <dbl> <dbl> <chr>      
## 1 Scores    12      7.69     3  4.44 0.032 Welch ANOVA

# Alternatively
# welch_anova_test(C3E9, Scores ~ Group)

Post hoc tests

# Games-Howell post-hoc tests
games_howell_test(formula = Scores ~ Group,
                  data = C3E9, 
                  conf.level = 0.95, 
                  detailed = FALSE)