Test for normal distribution

There are various tests and the outcome might differ!

The data

A tibble with two variables: normal and non_normal

set.seed(123)
mydata <- tibble(
  normal = rnorm(
    n = 200,
    mean = 50,
    sd = 5
  ),
  non_normal = runif(
    n = 200,
    min = 45,
    max = 55
  )
)

mydata
#> # A tibble: 200 × 2
#>    normal non_normal
#>     <dbl>      <dbl>
#>  1   47.2       54.9
#>  2   48.8       46.4
#>  3   57.8       54.1
#>  4   50.4       50.8
#>  5   50.6       49.0
#>  6   58.6       49.5
#>  7   52.3       52.1
#>  8   43.7       45.8
#>  9   46.6       48.4
#> 10   47.8       51.8
#> # ℹ 190 more rows

mydata |>
  pivot_longer(
    cols = 1:2,
    names_to = "type",
    values_to = "value"
  ) |>
  ggplot() +
  geom_histogram(
    alpha = 0.5,
    aes(
      x = value,
      fill = type,
      # plot probability instead of count on y axis
      y = ..density..
    )
  ) +
  stat_function(
    fun = dnorm,
    args = list(
      mean = 50,
      sd = 5
    ),
    color = "darkorange",
    size = 1
  ) +
  stat_function(
    fun = dnorm,
    args = list(
      mean = mean(mydata$non_normal),
      sd = sd(mydata$non_normal)
    ),
    color = "cyan4",
    size = 1
  ) +
  scale_fill_manual(values = c("cyan4", "darkorange")) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.85, 0.85)
  )

Shapiro-Wilk-Test

\(H_0\): Data does not differ from a normal distribution

shapiro.test(mydata$normal)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  mydata$normal
#> W = 0.99076, p-value = 0.2298

shapiro.test(mydata$non_normal)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  mydata$non_normal
#> W = 0.95114, p-value = 2.435e-06

Visual test with QQ-Plot

Points should match the straight line. Small deviations are okay.

# ggplot(
#   mydata,
#   aes(sample = normal)
# ) +
#   stat_qq() +
#   stat_qq_line()
ggpubr::ggqqplot(mydata$normal)

# ggplot(
#   mydata,
#   aes(sample = non_normal)
# ) +
#   stat_qq() +
#   stat_qq_line()
ggpubr::ggqqplot(mydata$non_normal)

The data

Counts of insects in agricultural units treated with different insecticides.

InsectSprays <- InsectSprays |>
  filter(spray %in% c("A", "B", "E")) |>
  mutate(spray = recode(spray, "E" = "C"))

TreatA <- filter(InsectSprays, spray == "A")$count
TreatB <- filter(InsectSprays, spray == "B")$count
TreatC <- filter(InsectSprays, spray == "C")$count

ggplot(InsectSprays, aes(x = spray, y = count, fill = spray)) +
  geom_boxplot(notch = TRUE) +
  theme(legend.position = "none")

Test for equal variance

First, test for normal distribution!

F-Test

• Normal distribution of groups
• Calculates ratio of variances (if equal, ratio = 1)
• p: How likely is this (or a more extreme) ratio of variances if the variances were truly equal?

Levene test

• Non-normal distribution of groups
• Compare difference between data sets with difference within data sets

Test for equal variance

First, test for normal distribution

shapiro.test(TreatA)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  TreatA
#> W = 0.95757, p-value = 0.7487
shapiro.test(TreatB)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  TreatB
#> W = 0.95031, p-value = 0.6415
shapiro.test(TreatC)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  TreatC
#> W = 0.92128, p-value = 0.2967

Result: All 3 treatments are normally distributed.

F-Test

\(H_0\): Variances do not differ between groups

var.test(TreatA, TreatB)
#> 
#>  F test to compare two variances
#> 
#> data:  TreatA and TreatB
#> F = 1.2209, num df = 11, denom df = 11, p-value = 0.7464
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.3514784 4.2411442
#> sample estimates:
#> ratio of variances 
#>            1.22093

F-Test

\(H_0\): Variances do not differ between groups

var.test(TreatA, TreatC)
#> 
#>  F test to compare two variances
#> 
#> data:  TreatA and TreatC
#> F = 7.4242, num df = 11, denom df = 11, p-value = 0.002435
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>   2.137273 25.789584
#> sample estimates:
#> ratio of variances 
#>           7.424242

Test for equal means

t-test

• Normal distribution AND equal variance
• Compares if mean values are within range of standard error of each other
• p: How likely is this (or more extreme) difference between means if the population means were truly equal?

t-test

\(H_0\): The samples do not differ in their mean

t.test(TreatA, TreatB, var.equal = TRUE)
#> 
#>  Two Sample t-test
#> 
#> data:  TreatA and TreatB
#> t = -0.45352, df = 22, p-value = 0.6546
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -4.643994  2.977327
#> sample estimates:
#> mean of x mean of y 
#>  14.50000  15.33333

Welch-Test

\(H_0\): The samples do not differ in their mean

t.test(TreatA, TreatC, var.equal = FALSE)
#> 
#>  Welch Two Sample t-test
#> 
#> data:  TreatA and TreatC
#> t = 7.5798, df = 13.91, p-value = 2.655e-06
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>   7.885546 14.114454
#> sample estimates:
#> mean of x mean of y 
#>      14.5       3.5

Wilcoxon-rank-sum Test

\(H_0\): The samples do not differ in their mean

wilcox.test(TreatA, TreatB)
#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  TreatA and TreatB
#> W = 62, p-value = 0.5812
#> alternative hypothesis: true location shift is not equal to 0

Paired values

Are there pairs of data points?

Example: samples of invertebrates across various rivers before and after sewage plants.

t.test(TreatA, TreatB, var.equal = TRUE, paired = TRUE)
t.test(TreatA, TreatB, var.equal = FALSE, paired = TRUE)
wilcox.test(TreatA, TreatB, paired = TRUE)

Plot test results with ggsignif

The ggsignif package offers a geom_signif() layer that can be added to a ggplot to annotate significance levels

# install.packages("ggsignif")
library(ggsignif)

Plot test results with geom_signif()

ggplot(
  InsectSprays,
  aes(x = spray, y = count)
) +
  geom_boxplot(notch = TRUE) +
  geom_signif(
    comparisons = list(
      c("A", "B"),
      c("B", "C"),
      c("A", "C")
    ),
    map_signif_level = TRUE,
    y_position = c(23, 24, 25)
  )

• By default, a Wilcoxon test is performed

Plot test results with geom_signif()

ggplot(
  InsectSprays,
  aes(x = spray, y = count)
) +
  geom_boxplot(notch = TRUE) +
  geom_signif(
    comparisons = list(
      c("A", "B"),
      c("B", "C"),
      c("A", "C")
    ),
    test = "t.test",
    test.args = list(
      var.equal = TRUE
    ),
    map_signif_level = TRUE,
    y_position = c(23, 24, 25)
  )

• test: run specific test
• test.args: pass additional arguments in a list
• ?geom_signif for more options

Plot mean +- se using stat_summary

Another way to plot the results is to plot mean and standard error of the mean:

ggplot(
  InsectSprays,
  aes(x = spray, y = count)
) +
  stat_summary()

Plot mean +- se using stat_summary

Another way to plot the results is to plot mean and standard error of the mean:

ggplot(
  InsectSprays,
  aes(x = spray, y = count)
) +
  stat_summary(
    fun.data = mean_se,
    geom = "errorbar"
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    color = "#28a87d",
    size = 4
  )

Plot mean +- se using stat_summary

Another way to plot the results is to plot mean and standard error of the mean:

ggplot(
  InsectSprays,
  aes(x = spray, y = count)
) +
  stat_summary(
    fun.data = mean_se,
    geom = "errorbar",
    width = 0.3
  ) +
  stat_summary(
    fun = mean,
    geom = "bar",
    size = 4
  )

Plot mean +- se using stat_summary

Just like before, you can also add a geom_signif to a barplot:

ggplot(
  InsectSprays,
  aes(x = spray, y = count)
) +
  stat_summary(
    fun.data = mean_se,
    geom = "errorbar",
    width = 0.3
  ) +
  stat_summary(
    fun = mean,
    geom = "bar"
  ) +
  ggsignif::geom_signif(
    comparisons = list(
      c("A", "B"),
      c("B", "C"),
      c("A", "C")
    ),
    test = "t.test",
    map_signif_level = TRUE,
    y_position = c(17, 18, 19)
  )