Day 4: One-Way ANOVA and pairwise significant testing

1. Setup

# core
library(dplyr)
library(ggplot2)
library(readr)
library(tidyr)
library(stringr)
# inference helpers
library(car)        # Levene's test
library(emmeans)    # estimated marginal means & Tukey

1.1 Load data

##Import dataset
# Read the combined dataset (ensure the CSV file is in your working directory)
dat <- read.csv("MultiFactor_data.csv")

2. Create the analysis subset (control + single-factor treatments)

In the dataset:

  • a column giving number of factors applied (given as Lv), where 0 = control, 1 = single-factor;
  • a column giving treatment identity for single-factor rows (given as remark).
treatments <- unique(dat$remark)[1:13] # extract control + single-factor treatments labels
single <- dat %>%
  filter(remark %in% treatments)         # control and single-factor only

# Choose responses of interest (edit as needed)
responses <- c("Decom") # soil pH, water-stable aggregates, decomposition rate etc.

summary(select(single, all_of(responses)))
##      Decom       
##  Min.   :0.5882  
##  1st Qu.:0.6425  
##  Median :0.6592  
##  Mean   :0.6608  
##  3rd Qu.:0.6781  
##  Max.   :0.7570  
##  NA's   :2

3. Exploratory plots

3.1 Boxplots by treatment

single_long <- single %>%
  pivot_longer(all_of(responses), names_to = "Response", values_to = "Value")%>%
  mutate(remark = factor(remark, levels = treatments))

ggplot(single_long, aes(remark, Value, fill = remark)) +
  geom_boxplot(outlier.alpha = 0.5) +
  labs(title = "Single-factor treatments vs Control",
       x = "Treatment (CT = Control; others = single GCFs)",
       y = "Observed value") +
  theme_bw() + theme(legend.position = "none",
                     axis.text.x = element_text(angle = 45, hjust = 1))

Tip: Look for visibly different medians/spreads across treatments. Mark imbalanced replication or potential outliers.

4. One-way ANOVA for one response

4.1 Compute One-way ANOVA

We’ll demonstrate for soil pH; repeat the same steps for WSA and Soil Decomposition rate.

# Compute the analysis of variance
res.aov <- aov(PH ~ remark, data = single)

# Summary of the analysis
summary(res.aov)

Interpretation: As the p-value is less than the significance level 0.05, we can conclude that there are significant differences between control and single-factor treatment groups.

4.2 Multiple pairwise-comparison between the means of groups

In one-way ANOVA test, a significant p-value indicates that some of the group means are different, but we don’t know which pairs of groups are different.

It’s possible to perform multiple pairwise-comparison, to determin if the mean difference between specific pairs of group are statistically significant.

Tukey multiple pairwise-comparisons

As the ANOVA test is significant, we can compute Tukey HSD (Tukey Honest Significant Differences, R function: TukeyHSD()) for performing multiple pairwise-comparison between the means of groups.

The function TukeyHD() takes the fitted ANOVA as an argument.

TukeyHSD(res.aov)

We have calculated adjusted p-values for all pairwise-groups. Given our research question, we’ll mainly focus on the comparisons of control treatment and single factor treatments (e.g., D-CT)

Pairwise t-test

The function pairewise.t.test() can be also used to calculate pairwise comparisons between group levels with corrections for multiple testing.

pairwise.t.test(single$PH, single$remark,
                 p.adjust.method = "BH")

The result is a table of p-values for the pairwise comparisons. Here, the p-values have been adjusted by the Benjamini-Hochberg method.

5. Check ANOVA assumptions: test validity

The ANOVA test assumes that, the data are normally distributed and the variance across groups are homogeneous. We can check that with some diagnostic plots.

5.1 Check the homogeneity of variance assumption

The residuals versus fit plot can be used to check the homogeneity of variances.

# 1. Homogeneity of variances
plot(res.aov, 1)

Points 110,99,56 are detected as outliers, which can severely affect normality and homogeneity of variance. It can be useful to remove outliers to meet the test assumptions.

In the plot above, there is no evident relationships between residuals and fitted values (the mean of each groups), which is good. So, we can assume the homogeneity of variances.

5.2 Check the normality assumption

Normality plot of residuals. In the plot below, the quantiles of the residuals are plotted against the quantiles of the normal distribution. A 45-degree reference line is also plotted.

The normal probability plot of residuals is used to check the assumption that the residuals are normally distributed. It should approximately follow a straight line.

# 2. Normality
plot(res.aov, 2)

Given all the points do not fall along this reference line, the normality is probably violated.

Then we conduct the Shapiro-Wilk test on the ANOVA residuals, which finds significant result (indicates the normality is violated).

# Extract the residuals
aov_residuals <- residuals(object = res.aov )
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals )

5.3 Non-parametric alternative to one-way ANOVA test

a non-parametric alternative to one-way ANOVA is Kruskal-Wallis rank sum test, which can be used when ANNOVA assumptions are not met:

kruskal.test(PH ~ remark, data = single)

Genearl linear Model