One-way analysis of variance (ANOVA) in R

Anna Shirokanova and Olesya Volchenko

March 2, 2020 + update 2021

This seminar

ANOVA is used:

See more: https://statistics.laerd.com/statistical-guides/one-way-anova-statistical-guide.php

F-test

Source: https://www.youtube.com/watch?v=0Vj2V2qRU10

Example from the World Values Survey data, Russia

Context to the analysis:

For an entertaining overview, see: http://www.worldvaluessurvey.org/WVSContents.jsp?CMSID=Findings

Our research question is: “Do Russians of different levels of education have similar levels of emancipative values?”

Let’s explore the issue. First, get the data:

World Value Survey, Russia subset

library(foreign)
wvs <- read.spss("WV6_Russia.sav", to.data.frame = T, use.value.labels = T)
wvsrus <- subset(wvs, V2 == "Russia")

dim(wvsrus) # there are 2,500 observations of 430 variables
## [1] 2500  430

Variables

wvsrus$RESEMAVAL1 <-  as.numeric(as.character(wvsrus$RESEMAVAL)) # emancipative values is a continuous variable

Emancipative Values

library(ggplot2)
ggplot(wvsrus, aes(x = RESEMAVAL1)) + 
  geom_density() +
  xlim(0, 1)

summary(wvsrus$RESEMAVAL1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.2998  0.3917  0.3931  0.4850  0.9259      22

Education levels in the data:

par(mar = c(5, 26, 2, 2))
barplot(table(wvsrus$V248)/nrow(wvsrus)*100, las = 2, horiz = T, xlim = c(0, 60))

Groups as they are now are very different in numbers, which may affect the F-test.

Recoding education levels (for the purpose of this example)

Let’s merge smaller groups.

wvsrus$educ <- rep(NA, length(wvsrus$V248)) # create a new, empty variable

wvsrus$educ[wvsrus$V248 == "No formal education" | 
            wvsrus$V248 == "Incomplete primary school" |
            wvsrus$V248 == "Complete primary school"|
            wvsrus$V248 == "Incomplete secondary school: university-preparatory type"|
            wvsrus$V248 == "Incomplete secondary school: technical/ vocational type"] <- "primary" 

wvsrus$educ[wvsrus$V248 == "Complete secondary school: university-preparatory type"] <- "sec general"

wvsrus$educ[wvsrus$V248 == "Complete secondary school: technical/ vocational type"] <- "sec vocational"

wvsrus$educ[wvsrus$V248 == "Some university-level education, without degree" | 
            wvsrus$V248 == "University - level education, with degree" ] <- "tertiary" # merge levels
wvsrus$educ <- as.factor(wvsrus$educ)

library(dplyr)
wvsrus <- select(wvsrus, RESEMAVAL1, educ) # select the variables you need to work with
summary(wvsrus)
##    RESEMAVAL1                 educ    
##  Min.   :0.0000   primary       :364  
##  1st Qu.:0.2998   sec general   :367  
##  Median :0.3917   sec vocational:978  
##  Mean   :0.3931   tertiary      :779  
##  3rd Qu.:0.4850   NA's          : 12  
##  Max.   :0.9259                       
##  NA's   :22

Values’ descriptives across new groups

library(psych)
library(magrittr)
library(kableExtra)

describeBy(wvsrus$RESEMAVAL1, wvsrus$educ, mat = TRUE) %>% #create a table
  select(
    Education = group1,
    N = n,
    Mean = mean,
    SD = sd,
    Median = median,
    Min = min,
    Max = max,
    Skew = skew,
    Kurtosis = kurtosis,
    st.error = se
  ) %>%
  kable(
    align = c("lrrrrrrrrr"),
    digits = 2,
    row.names = FALSE,
    caption = "Emancipative Values by Education Level"
  ) %>%
  kable_styling(bootstrap_options = c("bordered", "responsive", "striped"), full_width = FALSE)
Emancipative Values by Education Level
Education N Mean SD Median Min Max Skew Kurtosis st.error
primary 356 0.35 0.13 0.34 0.06 0.74 0.41 -0.08 0.01
sec general 363 0.38 0.13 0.38 0.08 0.76 0.05 -0.59 0.01
sec vocational 973 0.39 0.13 0.39 0.00 0.93 0.07 0.04 0.00
tertiary 774 0.42 0.13 0.43 0.03 0.79 0.03 -0.16 0.00

New groups

par(mar = c(3,20,1,3))
barplot(table(wvsrus$educ)/nrow(wvsrus)*100, 
        horiz = T, 
        xlim = c(0, 60), 
        las = 2,
        main = "Education Levels in the Sample, percent") # we get 4 groups of comparable size, which is good

The groups are of comparable size now, and the continuous variable is normally distributed within the groups.

Comment: If group sizes are unequal, the power of the test is based on the smallest group, see more: https://www.theanalysisfactor.com/when-unequal-sample-sizes-are-and-are-not-a-problem-in-anova/

Boxplot: Emancipative Values by Education Level

ggplot(wvsrus, aes(x = educ, y = RESEMAVAL1)) +
  geom_boxplot() +
  theme_classic() + 
  coord_flip()

From this boxplot, we see that the Y variable is distributed rather normally in across the education groups and that the level of emancipative values is slightly higher among the better educated respondents.

ANOVA test

ANOVA is a parametric test that rests on three assumptions:

If those assumptions are not met by the data, the classic ANOVA test is not reliable.

Homogeneity of variances

library(car)
leveneTest(wvsrus$RESEMAVAL1 ~ wvsrus$educ)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    3  0.4393 0.7249
##       2462

P-value is equal to 0.72, which means that variances are equal, i.e., you can indicate in the ANOVA test that var.equal = TRUE.

F-test

oneway.test(wvsrus$RESEMAVAL1 ~ wvsrus$educ, var.equal = T) 
## 
##  One-way analysis of means
## 
## data:  wvsrus$RESEMAVAL1 and wvsrus$educ
## F = 31.738, num df = 3, denom df = 2462, p-value < 2.2e-16
aov.out <- aov(wvsrus$RESEMAVAL1 ~ wvsrus$educ)
# another function of ANOVA to be used here; implies that var.equal = T
summary(aov.out)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## wvsrus$educ    3   1.58  0.5258   31.74 <2e-16 ***
## Residuals   2462  40.79  0.0166                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 34 observations deleted due to missingness

F(3, 2462) = 31.738, p-value < 0.001 means that the difference in the level of emancipative values across education groups is statistically significant.

Normality of residuals

plot(aov.out, 2)

The data points lie mostly along the diagonal line, which means that the distribution of residuals is close to normal.

layout(matrix(1:4, 2, 2))
plot(aov.out)

In normally distributed residuals, you will see a straight red line in the two upper graphs, and a straight line along the diagonal in the Q-Q plot.

A more formal way to check normality of residuals:

anova.res <- residuals(object = aov.out)
describe(anova.res) # look at skew (<|0.5|) and kurtosis (<|1|)
##    vars    n mean   sd median trimmed  mad  min  max range skew kurtosis se
## X1    1 2466    0 0.13      0       0 0.14 -0.4 0.53  0.93  0.1    -0.12  0
shapiro.test(x = anova.res) # if the p-value is > .05, the distribution is normal. Shapiro-Wilk test is too conservative on larger samples like this one.
## 
##  Shapiro-Wilk normality test
## 
## data:  anova.res
## W = 0.99824, p-value = 0.008871
hist(anova.res)

Here, skew and kurtosis are OK for a normal distribution. However, the Shapiro-Wilk test says that it is not normal. Visual analysis and skew and kurtosis agree here, which is why we can conclude that the normality assumption holds.

Post hoc tests

To inspect which groups contribute to the statistical significance of this test, run post hoc (“after analysis”) tests.

Here, you have two options:

  1. if variances are equal across groups, use Tukey’s post hoc test;

  2. if variances are unequal, use Games-Howell test or pairwise t-test with Bonferroni’s correction.

Tukey’s Honestly Significant Differences

TukeyHSD(aov.out)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = wvsrus$RESEMAVAL1 ~ wvsrus$educ)
## 
## $`wvsrus$educ`
##                                  diff          lwr        upr     p adj
## sec general-primary        0.03380885  0.009126274 0.05849143 0.0024662
## sec vocational-primary     0.04705791  0.026561133 0.06755468 0.0000000
## tertiary-primary           0.07771337  0.056522529 0.09890421 0.0000000
## sec vocational-sec general 0.01324905 -0.007102517 0.03360062 0.3379127
## tertiary-sec general       0.04390452  0.022854091 0.06495494 0.0000005
## tertiary-sec vocational    0.03065546  0.014717825 0.04659310 0.0000048
plot(TukeyHSD(aov.out))

To make all labels visible in the plot, arrange margins and change the direction of labels

par(mar = c(5, 15, 3, 1))
Tukey <- TukeyHSD(aov.out)
plot(Tukey, las = 2)

We can see that all but one differences between pairs of groups are statistically significant. The difference between “secondary gen” and “secondary vocational” is not statistically significant.

Bonferroni post hoc

pairwise.t.test(wvsrus$RESEMAVAL1, wvsrus$educ, 
                adjust = "bonferroni", pool.sd = T) # for unequal variances, write 'pool.sd = FALSE'
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  wvsrus$RESEMAVAL1 and wvsrus$educ 
## 
##                primary sec general sec vocational
## sec general    0.00087 -           -             
## sec vocational 2.0e-08 0.09434     -             
## tertiary       < 2e-16 3.6e-07     2.4e-06       
## 
## P value adjustment method: holm

We see that all but one pairs have a statistically significant difference between the means. If you have trouble reading such p-values, use the options parameter:

options(scipen = 999)
pairwise.t.test(wvsrus$RESEMAVAL1, wvsrus$educ, 
                adjust = "bonferroni", pool.sd = TRUE) 
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  wvsrus$RESEMAVAL1 and wvsrus$educ 
## 
##                primary              sec general sec vocational
## sec general    0.00087              -           -             
## sec vocational 0.00000002           0.09434     -             
## tertiary       < 0.0000000000000002 0.00000036  0.00000244    
## 
## P value adjustment method: holm

The results are the same as in the TukeyHSD test.

Non-parametric equivalent of ANOVA on ranks

kruskal.test(RESEMAVAL1 ~ educ, data = wvsrus)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  RESEMAVAL1 by educ
## Kruskal-Wallis chi-squared = 92.577, df = 3, p-value <
## 0.00000000000000022

With KW chi-square(3) = 92.6, p-value is < .001, which means that the mean ranks of the education groups are not the same. This results confirms what we saw earlier in the ANOVA test.

If Kruskal-Wallis test is significant, you can run non-parametric post hoc tests, e.g. Games-Howell or Dunn’s test:

! Games-Howell post hoc test is the test of choice for parametric ANOVA with unequal variances.

library(rstatix)
games_howell_test(RESEMAVAL1 ~ educ, data = wvsrus)
## # A tibble: 6 x 8
##   .y.     group1     group2    estimate conf.low conf.high    p.adj p.adj.signif
## * <chr>   <chr>      <chr>        <dbl>    <dbl>     <dbl>    <dbl> <chr>       
## 1 RESEMA… primary    sec gene…   0.0338  0.00905    0.0586  3.00e-3 **          
## 2 RESEMA… primary    sec voca…   0.0471  0.0264     0.0677  4.49e-8 ****        
## 3 RESEMA… primary    tertiary    0.0777  0.0561     0.0994  0.      ****        
## 4 RESEMA… sec gener… sec voca…   0.0132 -0.00680    0.0333  3.23e-1 ns          
## 5 RESEMA… sec gener… tertiary    0.0439  0.0229     0.0650  6.28e-7 ****        
## 6 RESEMA… sec vocat… tertiary    0.0307  0.0146     0.0467  5.68e-6 ****
library(DescTools) 
DunnTest(RESEMAVAL1 ~ educ, data = wvsrus)
## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                            mean.rank.diff                 pval    
## sec general-primary              196.7600              0.00042 ***
## sec vocational-primary           266.8025         0.0000000073 ***
## tertiary-primary                 426.8090 < 0.0000000000000002 ***
## sec vocational-sec general        70.0425              0.10971    
## tertiary-sec general             230.0490         0.0000015181 ***
## tertiary-sec vocational          160.0065         0.0000092217 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

These results show that all but one pair of groups have statistically significant differences in their medians.

See: https://stats.stackexchange.com/questions/126686/how-to-read-the-results-of-dunns-test

This also confirms what we saw above for the ANOVA post hoc tests.

Effect size

Now that you have established the statistical significance, you may think, what is the substantive significance of group membership in determining the level of emancipative values.

To answer this question, look at the share of variance explained by education groups in all the variance of emancipative values.

This is called the ‘effect size’ in ANOVA, or omega-squared. The larger it is, the more important your finding.

Omega-squared varies from 0 to 1. Any effect between 0.01 and 0.06 is considered small, 0.06 and 0.14 medium, and above 0.14 large.

library(sjstats)
anova_stats(aov.out) # replace the 'aov.out' with the name of your ANOVA resulting object
##                    term   df  sumsq meansq statistic p.value etasq
## wvsrus$educ wvsrus$educ    3  1.578  0.526    31.738       0 0.037
## ...2          Residuals 2462 40.791  0.017        NA      NA    NA
##             partial.etasq omegasq partial.omegasq epsilonsq cohens.f power
## wvsrus$educ         0.037   0.036           0.036     0.036    0.197     1
## ...2                   NA      NA              NA        NA       NA    NA
library(effectsize)
omega_squared(aov.out)
## Parameter   | Omega2 (partial) |    9e+01% CI
## ---------------------------------------------
## wvsrus$educ |             0.04 | [0.02, 0.05]

The omega-squared is 0.036 - this is a small effect (between 0.01 and 0.06). It means that even though there is a statistically significant difference in the level of emancipative values across education groups and this difference is statistically significant across all but one pairs of groups, in practical terms this effect is modest.

Extra: How is the effect size for ANOVA calculated?

\[ \omega^2 = \frac{(SS_{model} - (df_{model})*(MS_{error}))}{ SS_{model} + MS_{error} + SS_{error}}\]

See: http://www.statisticshowto.com/omega-squared/

What about the effect size for a non-parametric test?

Consider the epsilon-squared: \[ \epsilon^2 = \frac{{\chi^2} - k + 1}{n - k}\] where the chi-squared is the statistic from the Kruskal-Wallis test, k is the number of groups, and n is the total number of observations.

The effect size was already provided in the anova_stats() output, but there is also a separate function for it if you need it:

library(rcompanion)
epsilonSquared(x = wvsrus$RESEMAVAL1, g = wvsrus$educ)
## epsilon.squared 
##           0.037
library(rstatix)
kruskal_effsize(data = wvsrus, RESEMAVAL1 ~ educ)
## # A tibble: 1 x 5
##   .y.            n effsize method  magnitude
## * <chr>      <int>   <dbl> <chr>   <ord>    
## 1 RESEMAVAL1  2500  0.0359 eta2[H] small

“The epsilon-squared shows the percentage of variance in the dependent variable explained by the independent variable. The interpretation values commonly in published literature are: 0.01- < 0.06 (small effect), 0.06 - < 0.14 (moderate effect) and >= 0.14 (large effect).”

Read more: https://rpkgs.datanovia.com/rstatix/reference/kruskal_effsize.html

Here, the epsilon-squared equals 0.036, which means a small effect.

Sources to read and think over:

To summarise the key points