image from Giphy

1.1 Testing for normality

But hang on!

There are several assumptions behind the classic t-test we haven’t examined properly, namely:

1. independence – sample is independent
2. normality – data for each group is normally distributed
3. homoscedasticity – data across samples have equal variance

We can at least be sure of (1), as we know that senior ICs and Managers are separate populations. However, (2) and (3) are assumptions that we have to validate and address specifically. To test whether our data is normally distributed, we can use the Shapiro-Wilk test of normality, with the function shapiro.test():

pq_data_grouped %>%
  group_by(ManagerIndicator) %>%
  summarise(
    p = shapiro.test(Multitasking_hours)$p.value,
    statistic = shapiro.test(Multitasking_hours)$statistic
  )
## # A tibble: 2 x 3
##   ManagerIndicator      p statistic
##   <fct>             <dbl>     <dbl>
## 1 Manager          0.146      0.936
## 2 Senior IC        0.0722     0.941

As both p-values show up as less than 0.05, the test implies that we should reject the null hypothesis that the data are normally distributed (i.e. not normally distributed). To confirm, you can also perform a visual check for normality using a histogram or a Q-Q plot.

# Multitasking hours - IC
mth_ic <-
  pq_data_grouped %>%
  filter(ManagerIndicator == "Senior IC") %>%
  pull(Multitasking_hours) 

qqnorm(mth_ic, pch = 1, frame = FALSE)
qqline(mth_ic, col = "steelblue", lwd = 2)

< !-- -->

# Multitasking hours - Manager
mth_man <-
  pq_data_grouped %>%
  filter(ManagerIndicator == "Manager") %>%
  pull(Multitasking_hours) 

qqnorm(mth_man, pch = 1, frame = FALSE)
qqline(mth_man, col = "steelblue", lwd = 2)

< !-- -->

In the Q-Q plots, the points broadly adhere to the reference line. Therefore, the graphical approach suggests that the Shapiro-Wilk test may have been slightly over-sensitive. Below is a good thing to bear in mind:³

Statistical tests have the advantage of making an objective judgment of normality but have the disadvantage of sometimes not being sensitive enough at low sample sizes or overly sensitive to large sample sizes. Graphical interpretation has the advantage of allowing good judgment to assess normality in situations when numerical tests might be over or undersensitive.

In other words, the sample sizes may have well played a role in the significant result in our Shapiro-Wilk test.⁴ As our data isn’t conclusively normal – this in turn makes the unpaired t-test less conclusive. When we cannot safely assume normality, we can consider other alternatives such as the non-parametric two-samples Wilcoxon Rank-Sum test. This is covered further down below.

1.2 Testing for equality of variance (homoscedasticity)

Asides from normality, another assumption of the t-test that we hadn’t properly test for prior to running t.test() is to check for equality of variance across the two groups (homoscedasticity). Thankfully, this was not something we had to worry about as we used the Welch’s t-test. Recall that the classic Student’s t-test assumes equality between the two variances, but the Welch’s t-test already takes the difference in variance into account.

If required, however, here is an example on how you can test for homoscedasticity in R, using var.test():

# F test to compare two variances
var.test(
  Multitasking_hours ~ ManagerIndicator,
  data = pq_data_grouped
  )
## 
##  F test to compare two variances
## 
## data:  Multitasking_hours by ManagerIndicator
## F = 4.5726, num df = 22, denom df = 32, p-value = 0.0001085
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   2.146082 10.318237
## sample estimates:
## ratio of variances 
##           4.572575

The var.test() function ran above is an F-test (i.e. uses the F-distribution) used to compare whether the variances of two samples are the same. Under the null hypothesis of the tests, there should be homoscedasticity and as the f-statistic is a ratio of variances, the f-statistic would tend towards 1. The arguments are provided in a similar format to t.test().

It appears that homoscedasticity does not hold: since the p-value is less than 0.05, we should reject the null hypothesis that variances between the manager and IC dataset are equal. The Student’s t-test would not have been appropriate here, and we were correct to have used the Welch’s t-test.

Homoscedasticity can also be examined visually, using a boxplot or a dotplot (using graphics::dotchart() – suitable for small datasets). The code to do so would be as follows. For this example, visual examination is a bit more challenging as the senior IC and Manager groups have starkly different levels of multi-tasking hours.

dotchart(
  x = pq_data_grouped$Multitasking_hours,
  groups = pq_data_grouped$ManagerIndicator
)

< !-- -->

boxplot(
  Multitasking_hours ~ ManagerIndicator,
  data = pq_data_grouped
)

< !-- -->

2.1 Wilcoxon Rank-Sum Test

Previously, we could not safely rely on the unpaired two-sample t-test because we are not fully confident that the data satisfies the normality condition. As an alternative, we can use the Wilcoxon Rank-Sum test (aka Mann Whitney U Test). The Wilcoxon test is described as a non-parametric test, which in statistics typically means that there is no specification on a distribution, or the parameters of a distribution. In this case, the Wilcoxon test does not assume a normal distribution.

Another difference between the Wilcoxon Rank-Sum test and the unpaired t-test is that the former tests whether two populations have the same shape via comparing medians, whereas the latter parametric test compares means between two independent groups.

This is run using wilcox.test()

wilcox.test(
  Multitasking_hours ~ ManagerIndicator,
  data = pq_data_grouped,
  paired = FALSE
)
## 
##  Wilcoxon rank sum exact test
## 
## data:  Multitasking_hours by ManagerIndicator
## W = 752, p-value = 2.842e-14
## alternative hypothesis: true location shift is not equal to 0

The p-value of the test is less than the significance level (alpha = 0.05), which allows us to conclude that Managers’ median multitasking hours is significantly different from the ICs’.

Note that the Wilcoxon Rank-Sum test is different from the similarly named Wilcoxon Signed-Rank test, which is the equivalent alternative for the paired t-test. To perform the Wilcoxon Signed-Rank test instead, you can simply specify the argument to be paired = TRUE. Similar to the decision of whether to use the paired or the unpaired t-test, you should ensure that the one-sample condition applies if you use the Wilcoxon Signed-Rank test.

2.2 Kruskal-Wallis test

So far, we have only been looking at tests which compare exactly two populations. If we are looking for a test that works with comparisons across three or more populations, we can consider the Kruskal-Wallis test.

Let us create a new data frame that is grouped at the PersonId level, but filtering out fewer values in LevelDesignation:

pq_data_grouped_2 <-
  pq_data %>%
  filter(LevelDesignation %in% c(
    "Support",
    "Senior IC",
    "Junior IC",
    "Manager",
    "Director"
  )) %>%
  mutate(ManagerIndicator = factor(LevelDesignation)) %>%
  group_by(PersonId, ManagerIndicator) %>%
  summarise(Multitasking_hours = mean(Multitasking_hours), .groups = "drop")
  
glimpse(pq_data_grouped_2)
## Rows: 198
## Columns: 3
## $ PersonId           <chr> "0049ef24-ec83-356d-89f7-46b67364e677", "00f6d464-b~
## $ ManagerIndicator   <fct> Support, Senior IC, Manager, Support, Support, Supp~
## $ Multitasking_hours <dbl> 0.3812649, 0.2813373, 0.5980080, 0.2918829, 0.42288~

We can then run the Kruskal-Wallis test:

kruskal.test(
  Multitasking_hours ~ ManagerIndicator,
  data = pq_data_grouped_2
)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Multitasking_hours by ManagerIndicator
## Kruskal-Wallis chi-squared = 91.061, df = 4, p-value < 2.2e-16

Based on the Kruskal-Wallis test, we reject the null hypothesis and we conclude that at least one value in LevelDesignation is different in terms of their weekly hours spent multitasking. The most obvious downside to this method is that it does not tell us which groups are different from which, so this may need to be followed up with multiple pairwise-comparison tests (also known as post-hoc tests).

3.1 ANOVA

What if we want to run the t-test across more than two groups?

Analysis of Variance (ANOVA) is an alternative method that generalises the t-test beyond two groups, so it is used to compare three or more groups.

There are several versions of ANOVA. The simple version is the one-way ANOVA, but there is also two-way ANOVA which is used to estimate how the mean of a quantitative variable changes according to the levels of two categorical variables (e.g. rain/no-rain and weekend/weekday with respect to ice cream sales). In this example we will focus on one-way ANOVA.

There are three assumptions in ANOVA, and this may look familiar:

• The data are independent.
• The responses for each factor level have a normal population distribution.
• These distributions have the same variance.

These assumptions are the same as those required for the classic t-test above, and it is recommended that you check for variance and normality prior to ANOVA.

ANOVA calculates the ratio of the between-group variance and the within-group variance (quantified using sum of squares), and then compares this with a threshold from the Fisher distribution (typically based on a significance level). The key function is aov():

res_aov <-
  aov(
    Multitasking_hours ~ ManagerIndicator,
    data = pq_data_grouped_2
  )

summary(res_aov)
##                   Df Sum Sq Mean Sq F value Pr(>F)    
## ManagerIndicator   4  40.55   10.14   504.6 <2e-16 ***
## Residuals        193   3.88    0.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The interpretation is as follows:⁵

• Df: degrees of freedom for…

  • the outcome variable, i.e. the number of levels in the variable minus 1
  • the residuals, i.e. the total number of observations minus one and minus the number of levels in the outcome variables

• Sum Sq: sum of squares, i.e. the total variation between the group means and the overall mean

• Mean Sq: mean of the sum of squares, calculated by dividing the sum of squares by the degrees of freedom for each parameter

• F value: test statistic from the F test. This is the mean square of each independent variable divided by the mean square of the residuals. The larger the F value, the more likely it is that the variation caused by the outcome variable is real and not due to chance.

• Pr(>F): p-value of the F-statistic. This shows how likely it is that the F-value calculated from the test would have occurred if the null hypothesis of no difference among group means were true.

Given that the p-value is smaller than 0.05, we reject the null hypothesis, so we reject the hypothesis that all means are equal. Therefore, we can conclude that at least one value in LevelDesignation is different in terms of their weekly hours spent multitasking.

Antoine Soetewey’s blog recommends the use of the report package, which can help you make sense of the results more easily:

library(report)

report(res_aov)
## The ANOVA (formula: Multitasking_hours ~ ManagerIndicator) suggests that:
## 
##   - The main effect of ManagerIndicator is statistically significant and large
## (F(4, 193) = 504.61, p < .001; Eta2 = 0.91, 95% CI [0.90, 1.00])
## 
## Effect sizes were labelled following Field's (2013) recommendations.

The same drawback that applies to the Kruskall-Wallis test also applies to ANOVA, in that doesn’t actually tell you which exact group is different from which; it only tells you whether any group differs significantly from the group mean. This ANOVA test is hence sometimes also referred to as an ‘omnibus’ test.

3.2 Next steps after ANOVA

A pairwise t-test (note: pairwise, not paired!) is likely required to provide more information, and it is recommended that you review the p-value adjustment methods when doing so.⁶ Type I errors are more likely when running t-tests pairwise across many variables, and therefore correction is necessary. Here is an example of how you might run a pairwise t-test:

pairwise.t.test(
  x = pq_data_grouped_2$Multitasking_hours,
  g = pq_data_grouped_2$ManagerIndicator,
  paired = FALSE,
  p.adjust.method = "bonferroni"
)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pq_data_grouped_2$Multitasking_hours and pq_data_grouped_2$ManagerIndicator 
## 
##           Director Junior IC Manager Senior IC
## Junior IC <2e-16   -         -       -        
## Manager   <2e-16   <2e-16    -       -        
## Senior IC <2e-16   1         <2e-16  -        
## Support   <2e-16   1         <2e-16  1        
## 
## P value adjustment method: bonferroni

It may not be surprising that a pairwise method also exists as a follow-up for the Kruskall-Wallis test – which is the pairwise Wilcoxon test! This can be run using pairwise.wilcox.test(). The API for the pairwise.wilcox.test() is very similar to pairwise.t.test() where you can change the p-value adjustment method using the argument p.adjust.method:

pairwise.wilcox.test(
  x = pq_data_grouped_2$Multitasking_hours,
  g = pq_data_grouped_2$ManagerIndicator,
  paired = FALSE,
  p.adjust.method = "bonferroni"
)
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  pq_data_grouped_2$Multitasking_hours and pq_data_grouped_2$ManagerIndicator 
## 
##           Director Junior IC Manager Senior IC
## Junior IC 5.3e-09  -         -       -        
## Manager   3.3e-09  1.3e-09   -       -        
## Senior IC 5.9e-11  1         2.8e-13 -        
## Support   1.3e-08  1         8.6e-13 1        
## 
## P value adjustment method: bonferroni

4.1 Should I use a t-test or ANOVA for comparing exactly two groups?

One question worth discussing is the scenario at (1). Suppose that normality is observed in both groups, does it make a difference whether I use the t-test or ANOVA if I am comparing exactly two groups?

The textbook recommendation is that whenever one is comparing exactly two groups one should use the t-test, and ANOVA whenever there are more than two groups being compared. What can get confusing here is that there is the classic Student’s t-test and the Welch’s t-test.

When ANOVA is used to compare two groups, the results will be equivalent to a classic (Student’s) t-test with equal variances.⁷ However, if we are talking about the Welch’s t-test instead, it may be preferable over ANOVA because the Welch’s t-test takes into account heteroscedasticity. When there is heteroscedasticity, ANOVA (as well as Kruskall-Wallis) would become unstable and produce Type I errors, such as:

• conservative estimates for large sample sizes
• inflated estimates for small sample size⁸

To further complicate matters, there is also a method called Welch’s ANOVA which is like classic ANOVA but handles unequal variances better. This can be done in R using oneway.test(), but there is some debate around best practice that is beyond the scope of this post. ⁹ It would be prudent to run the Welch versions of the tests whenever we suspect the data to be heteroscedastic.

The recurring themes here are: (1) to check for heteroscedasticity and normality, and (2) to run multiple tests to acquire a more comprehensive view.

4.2 t-tests, ANOVA, and linear regression – are they completely different?

The common assumptions shared by the three methods may have gave it away, but the t-test, ANOVA, and linear regression are actually related in the sense that one is a special case of another.

The t-test is considered a special case of ANOVA, since the classic Student’s t-test is the same as ANOVA in comparing two groups when variances are equal. When the t-test statistic is squared, you get the corresponding f-statistic in the ANOVA.¹⁰

On the other hand, an ANOVA model is the same as a regression with a dummy variable. In fact, the aov() function in R is a wrapper around the linear regression function lm(). Steve Midway’s Analysis in R has a chapter which compares the outputs when running ANOVA using lm() versus aov().

All of these procedures are subsumed under the General Linear Model and share the same assumptions.

Flowchart for inferential statistics from Grosofsky (2009)