Case Study 3
This case study uses the data set from Laerd Statistics for ANOVA.
cs3 <- list()
# Theses are the same, so just use one.
# cs3$kw_dat <- read.spss("./input/kruskal-wallis-h-test.sav", to.data.frame = TRUE)
# cs3$anova_dat <- read.spss("./input/one-way-anova.sav", to.data.frame = TRUE)
cs3$dat <- read.spss("./input/kruskal-wallis-h-test.sav", to.data.frame = TRUE)A study tests whether physically active individuals are better able to cope with workplace stress. The study categorizes \(n\) = 31 participants by physical activity level (“Sedentary”, “Low”, “Moderate”, and “High”) and measures their ability to cope with workplace-related stress (CWWS) as the average score of a series of Likert items on a questionnaire (higher scores indicating a greater CWWS ability).8
cs3$dat %>%
group_by(group) %>%
summarize(
.groups = "drop",
mean_coping_stress = mean(coping_stress),
cl_025 = mean_coping_stress + qnorm(.025) * sd(coping_stress) / sqrt(n()),
cl_975 = mean_coping_stress + qnorm(.975) * sd(coping_stress) / sqrt(n())
) %>%
ggplot(aes(x = group, y = mean_coping_stress)) +
geom_point(shape = 21, fill = "gray80", color = "black", size = 3) +
geom_errorbar(aes(ymin = cl_025, ymax = cl_975, width = 0.1)) +
labs(title = "Distribution of CWWS by Physical Activity Level Group",
x = NULL, y = "Score",
caption = "Error Bars: 95% CI")
(cs3$gt <- cs3$dat %>%
tbl_summary(
by = group,
label = list(coping_stress = "CWWR"),
type = coping_stress ~ "continuous2",
statistic = coping_stress ~ c("{median} ({p25}, {p75})",
"{mean}, {sd}")
) %>%
add_n())| Characteristic | N | Sedentary, N = 7 | Low, N = 9 | Moderate, N = 8 | High, N = 7 |
|---|---|---|---|---|---|
| CWWR | 31 | ||||
| Median (IQR) | 4.12 (3.55, 4.77) | 5.50 (4.45, 7.85) | 7.10 (6.03, 8.11) | 7.47 (6.97, 8.19) | |
| Mean, SD | 4.15, 0.77 | 5.88, 1.69 | 7.12, 1.57 | 7.51, 1.24 |
inline_text(cs3$gt, variable = coping_stress, column = "Sedentary", pattern = "{mean}, {sd}")## [1] "NA, NA"inline_text(cs3$gt, variable = coping_stress, level = "Mean, SD", column = "Sedentary")## [1] "4.15, 0.77"cs3$gt$table_body## # A tibble: 3 × 10
## variable row_type label n var_t…¹ var_l…² stat_1 stat_2 stat_3 stat_4
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 coping_stress label CWWR 31 contin… CWWR <NA> <NA> <NA> <NA>
## 2 coping_stress level Medi… <NA> contin… CWWR 4.12 … 5.50 … 7.10 … 7.47 …
## 3 coping_stress level Mean… <NA> contin… CWWR 4.15,… 5.88,… 7.12,… 7.51,…
## # … with abbreviated variable names ¹var_type, ²var_labelCWWS score increased from the sedentary (4.15, 0.77), to low (5.88, 1.69), to moderate (7.12, 1.57) to high (7.51, 1.24) physical activity groups, in that order.
5.6.1 Assumptions
Recall that the one-way ANOVA test is valid under three assumptions. One, there are no significant outliers that influence the group mean. Two, the dependent variable is at least approximately (ANOVA is robust to this assumption) normally distributed for each group if the sample size is small (for large sample sizes the Central Limit Theorem shows normality is unnecessary). Three, the dependent variable should have equal variances across groups. ANOVA is only sensitive to this condition if the group sample sizes are not similar.
Kruskal-Wallis has no assumptions per se, but the interpretation of its results depend on the distribution of the dependent variable. If the distributions are similar, then the test results tell you whether the medians differ. Otherwise, the test results tell you whether the distributions differ.
Use a boxplot to assess outliers for ANOVA and the data distribution (if you revert to Kruskal-Wallis). Values greater than 1.5 IQR from the hinges (values beyond the whiskers) are outliers. Outliers might occur from data entry errors or measurement errors, so investigate and fix or throw them out. If the outlier is a genuinely extreme, you still have a couple options before reverting to Kruskal-Wallis. You can transform the dependent variable, but don’t do this unless the data is also non-normal. Transforming the variable also has the downside of making interpretation more difficult. You can also leave the outlier(s) in if it doesn’t affect the conclusion. There are no outliers here.
cs3$dat %>%
ggplot(aes(x = group, y = coping_stress)) +
geom_boxplot(outlier.color = "goldenrod", outlier.size = 2) +
theme_minimal() +
labs(title = "Boxplot of CWWR vs Group",
y = "Score", x = "Group")
There is no accepted practice for determining whether distributions are similar. The boxplot reveals a wider range of values for “Low” group, but this is close enough to conclude the distributions are similar.
You can assume the populations are normally distributed if \(n_j >= 30\). Otherwise, try the Q-Q plot, or skewness and kurtosis values, or histograms. If you still don’t feel confident about normality, run the Shapiro-Wilk test of normality or the Kolmogorov-Smirnov test. Definitely do not use Shapiro-Wilk for \(n_j >= 30\) because it is too sensitive. The Normal Q-Q plot below looks good for all groups except perhaps the “Low” group. The Shapiro-Wilk test confirms this, with all p-values over .05.
cs3$dat %>%
ggplot(aes(sample = coping_stress)) +
stat_qq() +
stat_qq_line(col = "goldenrod") +
facet_wrap(facets = vars(group)) +
theme_minimal() +
labs(title = "Q-Q Plot", x = "Theoretical", y = "Sample")
with(cs3$dat, by(coping_stress, group, shapiro.test)) %>%
map(tidy) %>%
bind_rows(.id = "group")## # A tibble: 4 × 4
## group statistic p.value method
## <chr> <dbl> <dbl> <chr>
## 1 Sedentary 0.928 0.538 Shapiro-Wilk normality test
## 2 Low 0.841 0.0589 Shapiro-Wilk normality test
## 3 Moderate 0.976 0.940 Shapiro-Wilk normality test
## 4 High 0.944 0.671 Shapiro-Wilk normality testHad the data failed the normality test, you could probably carry on anyway since the test is fairly robust to deviations from normality, particularly if the sample sizes are nearly equal. You can also try transforming the dependent variable. Transformations will generally only work when the distribution of scores in all groups are the same shape. Otherwise, revert to the Kruskal-Wallis H test.
ANOVA’s equality of sample variances condition is less critical when sample sizes are similar among the groups (as they are here). A rule of thumb is that no group’s standard deviation should be more than double that of any other. In this case, “Moderate” and “Low” are more than double “Sedentary”.
## # A tibble: 4 × 4
## group n sd multiple
## <fct> <int> <dbl> <dbl>
## 1 Sedentary 7 0.771 1
## 2 High 7 1.24 1.61
## 3 Moderate 8 1.57 2.04
## 4 Low 9 1.69 2.19There are two common tests, Bartlett and Levene9. Levene is less sensitive to departures from normality. Neither test rejects the null hypothesis of equality of variance here.
(cs3$levene <- car::leveneTest(coping_stress ~ group, data = cs3$dat, center = "mean"))
## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 3 2.129 0.1199
## 27
(cs3$bartlet <- bartlett.test(coping_stress ~ group, data = cs3$dat))
##
## Bartlett test of homogeneity of variances
##
## data: coping_stress by group
## Bartlett's K-squared = 3.7489, df = 3, p-value = 0.2899Heterogeneity of variances is a common problem in ANOVA. The Box-Cox procedure can help find a good transformation to remove heterogeneity. MASS::boxcox() calculates a profile of log-likelihoods for a power transformation of the dependent variable \(Y^\lambda\).
| \(\lambda\) | \(Y^\lambda\) | Transformation |
|---|---|---|
| 2 | \(Y^2\) | Square |
| 1 | \(Y^1\) | (no transformation) |
| .5 | \(Y^{.5}\) | Square Root |
| 0 | \(\ln(Y)\) | Log |
| -.5 | \(Y^{-.5}\) | Inverse Square Root |
| -1 | \(Y^{-1}\) | Inverse |
The Box-Cox procedure does not recommend any particular transformation of the data in this case.
MASS::boxcox(aov(coping_stress ~ group, data = cs3$dat), plotit = TRUE)
Had the data failed the homogeneity assumption, you can use a modified version of ANOVA called Welch’s ANOVA with Games-Howell post hoc test, or revert to the nonparametric Kruskal-Wallis test.
5.6.2 ANOVA
If the dependent variable conforms to the three ANOVA assumptions of no outliers, normality, and homogeneity, then you can run a one-way ANOVA with aov(). If the dependent variable only violates the homegeneity assumption, you can run Welch’s ANOVA with oneway.test(..., var.equal = FALSE)
cs3$aov <- aov(coping_stress ~ group, data = cs3$dat)
(cs3$anova <- anova(cs3$aov))## Analysis of Variance Table
##
## Response: coping_stress
## Df Sum Sq Mean Sq F value Pr(>F)
## group 3 49.033 16.3443 8.316 0.0004454 ***
## Residuals 27 53.066 1.9654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(cs3$welch <- oneway.test(coping_stress ~ group, data = cs3$dat, var.equal = FALSE))##
## One-way analysis of means (not assuming equal variances)
##
## data: coping_stress and group
## F = 14.821, num df = 3.000, denom df = 14.574, p-value = 0.0001058The ability to cope with workplace-related stress (CWWS score) was statistically significantly different for different levels of physical activity group, F(3, 27) = 8.3, p = 0.0004.
tibble(
f_stat = seq(0, 10, .01),
d_val = df(f_stat, 3, 14.574),
p_f = pf(f_stat, cs3$anova$Df[1], cs3$anova$Df[2], lower.tail = FALSE),
region = if_else(p_f < .05, "reject", "accept")
) %>%
ggplot(aes(x = f_stat, y = d_val)) +
geom_area(aes(fill = region), show.legend = FALSE) +
geom_line() +
geom_vline(xintercept = cs3$anova$"F value"[1], linetype = 2, color = "firebrick") +
scale_fill_manual(values = c(reject = "firebrick", accept = "white")) +
labs(
title = glue::glue("F({paste(cs3$anova$Df, collapse = ', ')}) = {comma(cs3$anova$'F value'[1], .1)}, p = {comma(cs3$anova$'Pr(>F)'[1], .0001)}"),
x = "F", y = "P(F)"
)
The F test does not indicate which populations cause the rejection of \(H_0\). Conduct a Tukey post hoc test if you have no specific hypothesis about two groups differing or want to see all group differences.10 If you want to compare two groups or set of groups, then use a custom contrast. Tukey is valid for balanced designs. If you have different sample sizes per group, use the Tukey-Kramer post hoc test. If the dependent variable failed the homogeneity of variances assumption, you would run the Games-Howell post hoc test instead.
(cs3$tukey <- TukeyHSD(cs3$aov))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = coping_stress ~ group, data = cs3$dat)
##
## $group
## diff lwr upr p adj
## Low-Sedentary 1.7276175 -0.2057757 3.661011 0.0923527
## Moderate-Sedentary 2.9715262 0.9859704 4.957082 0.0018413
## High-Sedentary 3.3540854 1.3034122 5.404759 0.0006806
## Moderate-Low 1.2439086 -0.6202750 3.108092 0.2835038
## High-Low 1.6264679 -0.3069254 3.559861 0.1226045
## High-Moderate 0.3825593 -1.6029965 2.368115 0.9517285
(cs3$games_howell <- rstatix::games_howell_test(cs3$dat, coping_stress ~ group))
## # A tibble: 6 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.…¹
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 coping_stress Sedentary Low 1.73 -0.163 3.62 0.077 ns
## 2 coping_stress Sedentary Moderate 2.97 1.07 4.88 0.003 **
## 3 coping_stress Sedentary High 3.35 1.66 5.05 0.000584 ***
## 4 coping_stress Low Moderate 1.24 -1.04 3.53 0.423 ns
## 5 coping_stress Low High 1.63 -0.508 3.76 0.167 ns
## 6 coping_stress Moderate High 0.383 -1.76 2.52 0.951 ns
## # … with abbreviated variable name ¹p.adj.signifTukey post hoc analysis revealed that the increase from sedentary to moderate (2.97, 95% CI (0.99 to 4.96)) was statistically significant (p = 0.002), as well as the increase from sedentary to high (3.35, 95% CI (1.30 to 5.40)) was statistically significant (p = 0.001), but no other group differences were statistically significant.
Games-Howell post hoc analysis revealed that the increase from sedentary to moderate (2.97, 95% CI (1.07 to 4.88)) was statistically significant (p =0.003), as well as the increase from sedentary to high (3.35, 95% CI (1.66 to 5.05, p = 0.001).
If you have specific hypotheses about the differences between the groups of your independent variable, e.g., whether the mean CWWS differs between the “Low” and “Sedentary” groups, \(H_0: \sum_i^K{c_i u_i} = 0\) where \(c_i = (1, -1, 0, 0)\) or “Sedentary” and average of all others, \(c_i = (1, -1/3, -1/3, -1/3)\), set up a contrast using the multcomp package.
cs3$glht_1 <- multcomp::glht(cs3$aov, linfct = multcomp::mcp(group = c(-1, 1, 0, 0)))
summary(cs3$glht_1)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = coping_stress ~ group, data = cs3$dat)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 1.7276 0.7065 2.445 0.0213 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
cs3$glht_2 <- multcomp::glht(cs3$aov, linfct = multcomp::mcp(group = c(-1, 1/3, 1/3, 1/3)))
summary(cs3$glht_2)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = coping_stress ~ group, data = cs3$dat)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 2.6844 0.6029 4.452 0.000133 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)There are three groups, so you need to adjust the p-value and 95% CI for them.
confint(cs3$glht_1, level = 1-.05/3)
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = coping_stress ~ group, data = cs3$dat)
##
## Quantile = 2.5525
## 98.3333333333333% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 == 0 1.7276 -0.0757 3.5309
confint(cs3$glht_2, level = 1-.05/3)
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = coping_stress ~ group, data = cs3$dat)
##
## Quantile = 2.5525
## 98.3333333333333% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 == 0 2.6844 1.1454 4.2234There was no statistically significant increase in CWWS score from the sedentary group (4.15, 0.77) to the group performing a low level of physical activity (5.88, 1.69), a mean increase of 1.7 (95% CI, -0.08, 3.53), p = 0.021.
CWWS score was statistically significantly higher in the non-sedentary groups (M = 6.8) compared to the sedentary group (4.15, 0.77), a mean increase of 1.7 (95% CI, 1.15, 4.22), p = 0.000.
You may also want to report the \(\omega^2\) effect size,
\[\omega^2 = \frac{SSR - df_R \cdot MSE}{MSE + SST}\]
where SSR is the between groups sum of squares, 49.0, MSE is the within groups mean square, 1.97, and SST is the total sum of squares, 102.1.
(cs3$anova_stats <- sjstats::anova_stats(cs3$aov))## term | df | sumsq | meansq | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## -----------------------------------------------------------------------------------------------------------------------------------------
## group | 3 | 49.033 | 16.344 | 8.316 | < .001 | 0.480 | 0.480 | 0.415 | 0.415 | 0.423 | 0.961 | 0.993
## Residuals | 27 | 53.066 | 1.965 | | | | | | | | |\(\omega^2\) estimates the population effect size. It \(\omega^2\) ranges from -1 to +1. Here, \(\omega^2\) is 0.415.
comma(cs3$anova$`Sum Sq`[1], .1)## [1] "49.0"Alternatively, the partial eta squared statistic, \(\eta^2\), measures the effect size in the sample. Here \(\eta^2\) is 0.48.
Now you can report your results.
A one-way ANOVA was conducted to determine if the ability to cope with workplace-related stress (CWWS score) was different for groups with different physical activity levels. Participants were classified into four groups: sedentary (n = 7), low (n = 9), moderate (n = 8) and high levels of physical activity (n = 7). There were no outliers, as assessed by boxplot; data was normally distributed for each group, as assessed by Shapiro-Wilk test (p > .05); and there was homogeneity of variances, as assessed by Levene’s test of homogeneity of variances (p = 0.120). Data is presented as mean, standard deviation. CWWS score was statistically significantly different between different physical activity groups, F(3, 27) = 8.3, p = 0.0004, \(\omega^2\) = 0.415`. CWWS score increased from the sedentary (4.15, 0.77), to low (5.88, 1.69), to moderate (7.12, 1.57) to high (7.51, 1.24) physical activity groups, in that order. Tukey post hoc analysis revealed that the mean increase from sedentary to moderate (2.97, 95% CI (0.99 to 4.96)) was statistically significant (p = 0.002), as well as the increase from sedentary to high (3.35, 95% CI (1.30 to 5.40), p = 0.001), but no other group differences were statistically significant.
Had the dependent variable failed the homogeneity of variances assumption, you would report the results from Welch’s ANOVA,
The ability to cope with workplace-related stress (CWWS score) was statistically significantly different for different levels of physical activity group, Welch’s F(3, 14.6) = 14.8, p < .0005.
and the Games-Howell post hoc test,
Games-Howell post hoc analysis revealed that the increase from sedentary to moderate (2.97, 95% CI (1.07 to 4.88)) was statistically significant (p =0.003), as well as the increase from sedentary to high (3.35, 95% CI (1.66 to 5.05, p = 0.001).
BTW, it is worth noting the relationship with linear regression. The regression model intercept is the overall mean and the coefficient estimators indirectly indicate the group means. The analysis of variance table in a regression model shows how much of the overall variance is explained by those coefficient estimators. It’s the same thing.
5.6.3 Kruskal-Wallis Test
Run a Kruskal-Wallis H test with kruskal.test().
(cs3$kruskal <- kruskal.test(coping_stress ~ group, data = cs3$dat))##
## Kruskal-Wallis rank sum test
##
## data: coping_stress by group
## Kruskal-Wallis chi-squared = 14.468, df = 3, p-value = 0.002332The dependent variable has similarly shaped distributions for all groups of the independent variable, so you can conclude the median CWWS scores were statistically significantly different between groups, \(\chi^2\)(3) = 14.5, p = 0.0023. Otherwise you would conclude the distributions differ. You rejected the null hypothesis, so continue on with a post hoc test to determine which medians (similar distributions) or mean ranks (dissimilar distributions) differ with the Dunn procedure using a Bonferroni correction for multiple comparisons.
(cs3$dunn <- FSA::dunnTest(coping_stress ~ group, data = cs3$dat, method = "bonferroni"))## Registered S3 methods overwritten by 'FSA':
## method from
## confint.boot car
## hist.boot car## Dunn (1964) Kruskal-Wallis multiple comparison## p-values adjusted with the Bonferroni method.## Comparison Z P.unadj P.adj
## 1 High - Low 1.6801430 0.0929294869 0.557576921
## 2 High - Moderate 0.2163067 0.8287486760 1.000000000
## 3 Low - Moderate -1.5121301 0.1305007724 0.783004635
## 4 High - Sedentary 3.3216144 0.0008949829 0.005369897
## 5 Low - Sedentary 1.8429610 0.0653346998 0.392008199
## 6 Moderate - Sedentary 3.2142419 0.0013078945 0.007847367P.adj equals P.unadj multiplied by the number of comparisons (6). You could report the adjusted p or the unadjusted p with a note that you accepted statistical significance at the p < .05 / 6 = 0.0083 level.
Now you can report your results.
A Kruskal-Wallis test was conducted to determine if there were differences in CWWS scores between groups that differed in their level of physical activity: the “sedentary” (n = NA), “low” (n = NA), “moderate” (n = NA) and “high” (n = NA) physical activity level groups. Distributions of CWWS scores were similar for all groups, as assessed by visual inspection of a boxplot. Median CWWS scores were statistically significantly different between the different levels of physical activity group, \(\chi^2\)(3) = 14.5, p = 0.0023. Subsequently, pairwise comparisons were performed using Dunn’s (1964) procedure with a Bonferroni correction for multiple comparisons. Adjusted p-values are presented. This post hoc analysis revealed statistically significant differences in CWWS scores between the sedentary (Mdn = NA) and moderate (Mdn = NA) (p = 0.0078) and sedentary and high (Mdn = NA) (p = 0.0054) physical activity groups, but not between the low physical activity group (Mdn = NA) or any other group combination.
Had the distributions been different, you would report “CWWS scores” instead of “Median CWWS scores” and report the mean ranks instead of Mdn. Unfortunately, you cannot retrieve those ranks from the test object, so you would have to calculate them yourself.