# ==============================================================================
# Programmer: [RINKESHKUMAR AHIR]
# Title: Advance statistical programming - Data Analysis
# Purpose: Analysis of air pollutant data for Port Talbot, including data cleaning,
# visualization, correlation, t-test, and ANOVA.
# Date: November 2025
# ==============================================================================
# ==============================================================================
# Q.1 Create a data frame PortTalbotAQ, ensuring that the headers are imported. You should then rename the headers to shorter but still understandable descriptors. Your script should assume that the data file is stored in the same folder as the script file. Display the resulting data frame.
# Setup and Data Loading
# ==============================================================================
# Install needed libraries if not already installed
if (!requireNamespace("tidyverse", quietly = TRUE)) install.packages("tidyverse")
if (!requireNamespace("lubridate", quietly = TRUE)) install.packages("lubridate")
# Load libraries
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
# ------------------------------------------------------------------------------
# Locate the CSV in the same folder as this script.
# ------------------------------------------------------------------------------
args <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", args[grep("^--file=", args)][1])
if (is.na(script_path) || script_path == "") {
script_dir <- getwd()
} else {
script_dir <- dirname(normalizePath(script_path))
}
csv_path <- file.path(script_dir, "Port_Talbot_Hourly_AQ.csv")
# For check if the file is present where expected
if (!file.exists(csv_path)) {
stop(paste0(
"CSV not found at: ", csv_path, "\n",
))
}
# ------------------------------------------------------------------------------
# Load CSV into a data frame. "No data" is converted to NA
# ------------------------------------------------------------------------------
PortTalbotAQ <- read.csv(
csv_path,
header = TRUE,
stringsAsFactors = FALSE,
na.strings = "No data"
)
# ------------------------------------------------------------------------------
# Create a single DateTime column by combining Date and Time
# ------------------------------------------------------------------------------
PortTalbotAQ <- PortTalbotAQ %>%
mutate(
DateTime = as.POSIXct(
paste(Date, Time),
format = "%d-%m-%Y %H:%M:%S",
tz = "GMT"
)
)
# ------------------------------------------------------------------------------
# Select and rename columns to shorter, meaningful names
# ------------------------------------------------------------------------------
PortTalbotAQ <- PortTalbotAQ %>%
select(
Date,
Time,
DateTime,
Ozone,
Nitric_Oxide = `Nitric.oxide`,
Nitrogen_Dioxide = `Nitrogen.dioxide`,
Nitrogen_Oxides = `Nitrogen.oxides.as.nitrogen.dioxide`,
Sulphur_Dioxide = `Sulphur.dioxide`,
Carbon_Monoxide = `Carbon.monoxide`,
PM10 = `PM10.particulate.matter..Hourly.measured.`,
PM2_5 = `PM2.5.particulate.matter..Hourly.measured.`,
Wind_Direction = `Modelled.Wind.Direction`,
Wind_Speed = `Modelled.Wind.Speed`,
Temperature = `Modelled.Temperature`
)
# ------------------------------------------------------------------------------
# Convert measurement columns to numeric
# ------------------------------------------------------------------------------
numeric_cols <- c(
"Ozone", "Nitric_Oxide", "Nitrogen_Dioxide", "Nitrogen_Oxides",
"Sulphur_Dioxide", "Carbon_Monoxide", "PM10", "PM2_5",
"Wind_Direction", "Wind_Speed", "Temperature"
)
PortTalbotAQ[numeric_cols] <- lapply(PortTalbotAQ[numeric_cols], as.numeric)
# ------------------------------------------------------------------------------
# Display first 6 rows of the data frame
# ------------------------------------------------------------------------------
cat("PortTalbotAQ data frame first 6 rows:\n")
## PortTalbotAQ data frame first 6 rows:
print(head(PortTalbotAQ))
## Date Time DateTime Ozone Nitric_Oxide
## 1 01-01-2020 01:00:00 2020-01-01 01:00:00 32.39686 0.65250
## 2 01-01-2020 02:00:00 2020-01-01 02:00:00 27.07500 NA
## 3 01-01-2020 03:00:00 2020-01-01 03:00:00 34.72518 0.46059
## 4 01-01-2020 04:00:00 2020-01-01 04:00:00 39.21551 0.40301
## 5 01-01-2020 05:00:00 2020-01-01 05:00:00 36.12217 0.20151
## 6 01-01-2020 06:00:00 2020-01-01 06:00:00 33.37808 0.40301
## Nitrogen_Dioxide Nitrogen_Oxides Sulphur_Dioxide Carbon_Monoxide PM10 PM2_5
## 1 12.26557 13.26605 0.44178 0.070977 25.1 21.1
## 2 NA NA 0.70684 0.086749 30.6 25.0
## 3 12.30813 13.01435 0.92773 0.106465 27.6 25.9
## 4 7.24397 7.86191 0.66267 0.073934 23.6 19.8
## 5 7.33085 7.63982 0.39760 0.065062 22.9 20.2
## 6 8.88742 9.50536 0.46387 0.065062 20.9 15.9
## Wind_Direction Wind_Speed Temperature
## 1 96.7 4.6 3.0
## 2 101.2 5.3 3.1
## 3 98.9 5.1 3.1
## 4 103.5 4.8 3.2
## 5 97.4 4.4 3.4
## 6 90.0 4.3 3.6
# I used the tidy and usable data frame named PortTalbotAQ. I got the air-quality data as a CSV file and it was in the same folder with my script and ensured that the headers were properly added. Some no data I also transformed into meaningful missing values. I subsequently gave the names of the initial long column names shorter and clearer names to make the variables easier to comprehend and analyze. Then i transformed all the pollutant and weather related data into the numeric form and presented the initial rows of the cleaned data frame. This is for preconditions further analysis and visualisation of the dataset.
# ==============================================================================
# Q.2.Create a new data frame invalidData that shows the percentage of readings where no measurement has been recorded for each pollutant. Display this data frame.
# Analysis the Missing Data
# ==============================================================================
# Pollutant columns
pollutant_cols <- c(
"Ozone", "Nitric_Oxide", "Nitrogen_Dioxide", "Nitrogen_Oxides",
"Sulphur_Dioxide", "Carbon_Monoxide", "PM10", "PM2_5"
)
# Percentage of missing value per pollutant
invalidData <- PortTalbotAQ %>%
summarise(across(all_of(pollutant_cols), ~ mean(is.na(.)) * 100)) %>%
pivot_longer(everything(), names_to = "Pollutant", values_to = "Percentage_Missing") %>%
mutate(Percentage_Missing = round(Percentage_Missing, 2)) %>%
arrange(desc(Percentage_Missing))
cat("\n2.2 invalidData percentage of missing values per pollutant:\n")
##
## 2.2 invalidData percentage of missing values per pollutant:
print(invalidData)
## # A tibble: 8 × 2
## Pollutant Percentage_Missing
## <chr> <dbl>
## 1 PM2_5 6.86
## 2 PM10 6.54
## 3 Carbon_Monoxide 4.74
## 4 Nitrogen_Dioxide 3.98
## 5 Nitrogen_Oxides 3.98
## 6 Nitric_Oxide 3.67
## 7 Sulphur_Dioxide 2.4
## 8 Ozone 1.98
# I analysed the data to determine the degree of missing data of each pollutant. I takk all the pollutant variables and used the percentage of the readings in which no valid measurement was found. Then I added a new data frame, which is named invalidData, that summarises all these percentages in a clear format. This help me to easily identify the most missing data of pollutants and evaluate how good and reliable the dataset is in general.
# ==============================================================================
# Q.3. Create a function displayPollutantGraph() that will take the pollutant name, start month, start year, end month, and end year, and will produce a graph of the variation in pollutant levels measured at noon on the 1st of the month. Display this using the data for nitrogen oxides as nitrogen dioxide between July 2015 and November 2024 inclusive.
# Time Series Plotting for Single Pollutant
# ==============================================================================
# load packages
library(tidyverse)
library(lubridate)
library(rlang)
##
## Attaching package: 'rlang'
## The following objects are masked from 'package:purrr':
##
## flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
# Reading at 12:00:00 on the 1st of each month
monthly_noon_data <- PortTalbotAQ %>%
mutate(
Day = day(DateTime),
Hour = hour(DateTime)
) %>%
filter(Day == 1, Hour == 12)
# plot a pollutant over a date range. start month, start year, end month, end year numerics
displayPollutantGraph <- function(pollutant_name, start_month, start_year, end_month, end_year) {
# Validate the pollutant exists
if (!pollutant_name %in% names(monthly_noon_data)) {
stop(sprintf("Column '%s' not found. Available columns include: %s",
pollutant_name, paste(names(monthly_noon_data), collapse = ", ")))
}
# Zero pad months for parsing
sm <- sprintf("%02d", as.integer(start_month))
em <- sprintf("%02d", as.integer(end_month))
# Build inclusive range for start, first day of month after end
start_date <- as.POSIXct(paste0("01-", sm, "-", start_year, " 12:00:00"),
format = "%d-%m-%Y %H:%M:%S", tz = "GMT")
end_date <- as.POSIXct(paste0("01-", em, "-", end_year, " 12:00:00"),
format = "%d-%m-%Y %H:%M:%S", tz = "GMT")
end_date_inclusive <- end_date %m+% months(1) # safe month addition
# Filter for DateTime and selected pollutant
plot_data <- monthly_noon_data %>%
filter(DateTime >= start_date, DateTime < end_date_inclusive) %>%
select(DateTime, !!sym(pollutant_name)) %>%
drop_na()
if (nrow(plot_data) == 0) {
message("No data found for the requested pollutant and date range at noon on the 1st.")
return(invisible(NULL))
}
# Identify max and min rows for easyily understand.
max_row <- plot_data %>% slice_max(!!sym(pollutant_name), n = 1, with_ties = FALSE)
min_row <- plot_data %>% slice_min(!!sym(pollutant_name), n = 1, with_ties = FALSE)
# Labels
subtitle_txt <- sprintf("From %s %d to %s %d (noon on the 1st)",
month.abb[as.integer(sm)], start_year,
month.abb[as.integer(em)], end_year)
title_txt <- paste("Variation in", gsub("_", " ", pollutant_name))
ylab_txt <- paste(gsub("_", " ", pollutant_name), "Concentration")
# Plot
p <- ggplot(plot_data, aes(x = DateTime, y = !!sym(pollutant_name))) +
geom_line(linewidth = 1.1, colour = "steelblue") +
geom_point(size = 2.3, colour = "black") +
# Red color for max
geom_point(data = max_row, aes(x = DateTime, y = !!sym(pollutant_name)),
colour = "red", size = 4) +
geom_text(data = max_row,
aes(label = round(!!sym(pollutant_name), 1)),
colour = "red", vjust = -1, fontface = "bold") +
# Green color for min
geom_point(data = min_row, aes(x = DateTime, y = !!sym(pollutant_name)),
colour = "green4", size = 4) +
geom_text(data = min_row,
aes(label = round(!!sym(pollutant_name), 1)),
colour = "green4", vjust = 1.8, fontface = "bold") +
labs(title = title_txt, subtitle = subtitle_txt, x = "Date", y = ylab_txt) +
# Show dates on the x-axis
scale_x_datetime(date_labels = "%d-%b-%Y", date_breaks = "6 months") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1)
)
print(p)
}
# Nitrogen oxides as nitrogen dioxide, Jul 2015 – Nov 2024
cat("\nProducing Nitrogen Oxides Graph (Jul 2015 – Nov 2024):\n")
##
## Producing Nitrogen Oxides Graph (Jul 2015 – Nov 2024):
displayPollutantGraph("Nitrogen_Oxides", 7, 2015, 11, 2024)
# I creat a function that would generate a time series graph of any pollutant in the given date range. My initial screening of the dataset was to only serve the measurements that were made at 12:00 noon on the 1 st of every month as this would provide a regular monthly picture. Then I created a function that receives the name of the pollutant and the beginning and ending date, filters the data this way and draws the results. I made the graph better by including a coloured trend line, black markers of each reading, special markers that indicate the highest value in red and the lowest one in green with its respective values. I also made the x-axis format display entire dates to make the graph easy to read. Then I plotted nitrogen oxides in July to November 2015 through this functional.
# ==============================================================================
# Q.4. Create a function comparePollutants() that will take a vector of pollutants and will produce a graph of the variation in each pollutant level measured at noon on the 1st of the month for all ten years on the same graph. Display this using an appropriate vector of 4 pollutants.
# Time Series Plotting for four different Pollutant
# ==============================================================================
library(tidyverse)
library(lubridate)
library(RColorBrewer)
# For ensure monthly noon data exists
if (!exists("monthly_noon_data")) {
monthly_noon_data <- PortTalbotAQ %>%
mutate(Day = day(DateTime), Hour = hour(DateTime)) %>%
filter(Day == 1, Hour == 12)
}
# Choose 4 pollutants which i want to compare
pollutants <- c("Ozone", "Nitric_Oxide", "Nitrogen_Dioxide", "Sulphur_Dioxide")
# For color palette
pal <- setNames(brewer.pal(length(pollutants), "Set1"), pollutants)
clean_time_theme <- function() {
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 0, vjust = 0.5)
)
}
# Combined plot for standardized
combined_df <- monthly_noon_data %>%
select(DateTime, all_of(pollutants)) %>%
pivot_longer(-DateTime, names_to = "Pollutant", values_to = "Concentration") %>%
drop_na() %>%
group_by(Pollutant) %>%
mutate(Value = as.numeric(scale(Concentration))) %>%
ungroup() %>%
mutate(Label = gsub("_", " ", Pollutant))
p_combined <- ggplot(combined_df, aes(DateTime, Value, colour = Pollutant)) +
geom_line(linewidth = 1) +
labs(
title = "4 Diffrent Pollutants Levels at Noon on the 1st",
subtitle = sprintf("From %s to %s",
format(min(combined_df$DateTime), "%d-%b-%Y"),
format(max(combined_df$DateTime), "%d-%b-%Y")),
x = "Date", y = "Standardized Concentration", colour = "Pollutant"
) +
scale_color_manual(values = pal, labels = gsub("_", " ", pollutants)) +
scale_x_datetime(date_breaks = "1 year", date_labels = "%Y") +
clean_time_theme()
print(p_combined)
# Single pollutant plots
make_single_plot <- function(df, colname, colour_hex) {
nice <- gsub("_", " ", colname)
dat <- df %>% select(DateTime, !!sym(colname)) %>%
rename(Value = !!sym(colname)) %>%
drop_na()
ggplot(dat, aes(DateTime, Value)) +
geom_line(linewidth = 1, colour = colour_hex) +
labs(
title = paste0(nice, " — Noon on the 1st"),
subtitle = sprintf("From %s to %s",
format(min(dat$DateTime), "%d-%b-%Y"),
format(max(dat$DateTime), "%d-%b-%Y")),
x = "Date", y = paste(nice, "Concentration")
) +
scale_x_datetime(date_breaks = "1 year", date_labels = "%Y") +
clean_time_theme()
}
for (pname in pollutants) {
p_single <- make_single_plot(monthly_noon_data, pname, pal[[pname]])
print(p_single)
}
# I create a function in order to make a comparison between four different pollutants around noon on 1 st day of every month. I chose four pollutants and manipulated the information in a way that the values were visualised easily throughout the entire time frame. I started by creating a combined plot with all the four pollutants on the same plot using standardized values, which facilitates easy comparison of the trends of the four pollutants though in different units. Then i created four individual graphs, each representing every pollutant, thus allowing me to look at the pattern of each pollutant individually without being distracted by the rest. This graphs is easy to understand and demonstrate the change in the levels of pollutants with time.
# ==============================================================================
# Q.5. Create a function O3NOCorrelation() that calculates the correlation between ozone and nitrogen dioxide measurements for the entire data collection period. Plot the data used as a scatterplot.
# Ozone Nitrogen Dioxide Correlation
# ==============================================================================
library(tidyverse)
library(ggrepel)
O3NOCorrelation <- function(data_frame) {
cat("\n==============================================\n")
cat(" Ozone vs Nitrogen Dioxide\n")
cat(" Correlation and NO₂ Level Visualisation\n")
cat("==============================================\n\n")
# Extract and clean data
df <- data_frame %>%
select(Ozone, Nitrogen_Dioxide) %>%
drop_na()
cat("Cleaned dataset loaded.\n")
cat(paste("Number of valid data pairs:", nrow(df), "\n\n"))
# Calculate the correlation
r_value <- cor(df$Ozone, df$Nitrogen_Dioxide, method = "pearson")
cat("Pearson Correlation (r): ", round(r_value, 4), "\n")
# Create NO₂ Level Categories
df <- df %>%
mutate(
NO2_Level = case_when(
Nitrogen_Dioxide < quantile(Nitrogen_Dioxide, 0.33) ~ "Low NO₂",
Nitrogen_Dioxide < quantile(Nitrogen_Dioxide, 0.66) ~ "Medium NO₂",
TRUE ~ "High NO₂"
)
)
cat("O₂ Levels Categorised into: Low, Medium, High.\n\n")
# Scatterplot Showing NO2 Level
scatter_plot <- suppressWarnings(
ggplot(df, aes(x = Ozone, y = Nitrogen_Dioxide, colour = NO2_Level)) +
geom_point(size = 3, alpha = 0.8) +
geom_smooth(method = "lm", se = FALSE, colour = "black", linewidth = 1.3) +
# Correlation annotation
annotate("text",
x = Inf, y = Inf,
label = paste("r =", round(r_value, 3)),
hjust = 1.2, vjust = 1.5,
size = 6, fontface = "bold", colour = "#2C3E50") +
scale_colour_manual(values = c(
"Low NO₂" = "#1a9641", # green
"Medium NO₂" = "#fdae61", # orange
"High NO₂" = "#d7191c" # red
)) +
labs(
title = "Ozone vs Nitrogen Dioxide",
subtitle = "Points coloured by NO₂ Level (Low, Medium, High)",
x = "Ozone (µg/m³)",
y = "Nitrogen Dioxide (µg/m³)",
colour = "NO₂ Level"
) +
theme_minimal(base_size = 15) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5)
)
)
cat("scatterplot created with NO₂ level categories.\n")
print(scatter_plot)
return(r_value)
}
# Run the function
cor_value <- O3NOCorrelation(PortTalbotAQ)
##
## ==============================================
## Ozone vs Nitrogen Dioxide
## Correlation and NO₂ Level Visualisation
## ==============================================
##
## Cleaned dataset loaded.
## Number of valid data pairs: 83745
##
## Pearson Correlation (r): -0.4858
## O₂ Levels Categorised into: Low, Medium, High.
##
## scatterplot created with NO₂ level categories.
## `geom_smooth()` using formula = 'y ~ x'
cat("\n==============================================")
##
## ==============================================
cat("\n Final Correlation Value:", round(cor_value, 4))
##
## Final Correlation Value: -0.4858
cat("\n==============================================\n")
##
## ==============================================
# I observe how ozone and nitrogen dioxide are correlated with each other by establishing a function that determines the Pearson correlation between the two environmental pollutants.I started by choosing the two variables in the dataset and eliminating any non response cases. I followed it by finding the correlation coefficient to demonstrate the degree of relationship between the two pollutants. In order to visualise this relationship, I created a scatterplot in which every point would represent a matched combination of measurements of ozone and nitrogen dioxide. It also incorporates a fitted linear trend line and colour coded categories of the NO2 levels (low, medium, high) into the plot and, therefore, it is more convenient to interpret the pattern and the distribution of the data. This help me to give a numerical correlation as well as the visual representation of the relation between the two pollutants.
# ==============================================================================
# Q.6. Perform a t-test to identify whether there is a significant difference in ozone levels when wind speeds are below 5ms-1 and when they are equal to or greater than 5ms-1. You should state your experimental hypotheses, the result of the test, and whether you accept or reject the null hypothesis in your comments.
# T-Test for Does Ozone Differ Under Different Wind Speeds?
# ==============================================================================
# Prepare data means like select needed columns, remove missing values and define groups
t_test_data <- PortTalbotAQ %>%
dplyr::select(Ozone, Wind_Speed) %>%
tidyr::drop_na() %>%
dplyr::mutate(
Wind_Group = ifelse(Wind_Speed < 5, "<5 m/s", "≥5 m/s"),
Wind_Group = factor(Wind_Group, levels = c("<5 m/s", "≥5 m/s"))
)
# Descriptive statistics
group_summary <- t_test_data %>%
dplyr::group_by(Wind_Group) %>%
dplyr::summarise(
n = dplyr::n(),
mean_ozone = mean(Ozone),
sd_ozone = sd(Ozone),
.groups = "drop"
)
print(group_summary)
## # A tibble: 2 × 4
## Wind_Group n mean_ozone sd_ozone
## <fct> <int> <dbl> <dbl>
## 1 <5 m/s 42208 47.2 21.6
## 2 ≥5 m/s 41944 58.7 17.0
# Perform Welch two-sample t-test
t_test_result <- t.test(Ozone ~ Wind_Group, data = t_test_data,
alternative = "two.sided", var.equal = FALSE)
cat("\n================ T-TEST ================\n")
##
## ================ T-TEST ================
print(t_test_result)
##
## Welch Two Sample t-test
##
## data: Ozone by Wind_Group
## t = -86.531, df = 80001, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group <5 m/s and group ≥5 m/s is not equal to 0
## 95 percent confidence interval:
## -11.84095 -11.31642
## sample estimates:
## mean in group <5 m/s mean in group ≥5 m/s
## 47.15779 58.73648
# Interpretation at alpha = 0.05
alpha <- 0.05
pval <- t_test_result$p.value
cat("\n================ COMMENTARY ===================\n")
##
## ================ COMMENTARY ===================
cat("Hypotheses:\n")
## Hypotheses:
cat(" H0: Mean ozone is the same for <5 m/s and ≥5 m/s wind speeds.\n")
## H0: Mean ozone is the same for <5 m/s and ≥5 m/s wind speeds.
cat(" H1: Mean ozone differs between the two wind speed groups.\n\n")
## H1: Mean ozone differs between the two wind speed groups.
cat(sprintf("Group means:\n <5 m/s: %.3f\n ≥5 m/s: %.3f\n\n",
group_summary$mean_ozone[group_summary$Wind_Group == "<5 m/s"],
group_summary$mean_ozone[group_summary$Wind_Group == "≥5 m/s"]))
## Group means:
## <5 m/s: 47.158
## ≥5 m/s: 58.736
cat(sprintf("p-value = %.4f\n", pval))
## p-value = 0.0000
if (pval < alpha) {
cat("\nDecision: REJECT H0 (significant difference found).\n")
} else {
cat("\nDecision: FAIL TO REJECT H0 (no significant difference found).\n")
}
##
## Decision: REJECT H0 (significant difference found).
cat("================================================\n")
## ================================================
# I used to test the hypothesis that the level of ozone varies with the speed of wind. I split the dataset into two, the wind speeds with a wind speed that was less than 5 m/s and the wind speeds that were equal to or more than 5 m/s. I then took a mean of the ozone reading in each group after deleting the missing values, and used Welch two-sample t-test to compare the group mean. I explained the null hypothesis which is no difference in mean ozone between the two types of wind and the alternative hypothesis which puts that the means are dissimilar. The p-value of the t-test enabled me to reject or fail to reject the null hypothesis at the significance level of 0.05 and thus provided me with an opportunity to evaluate whether the speed of the wind had a statistically significant impact on the ozone concentrations.
# ==============================================================================
# Q.7. Perform an ANOVA test to identify if there are statistically significant differences in the levels of nitrogen dioxide during the three periods November-December 2016, November-December 2018, and November-December 2020. In each case the data should be taken for 17:00 on each day between November 1st and December 31st for the three years. If such a difference is found, perform the appropriate ad-hoc tests to identify where the differences lie. You should state your experimental hypotheses, the result of the test, and whether you accept or reject the null hypothesis in your comments. If you do post-hoc test then you should also explain what you are doing and why you are doing this in the comments.
# ANOVA on NO2 at 17:00 during Nov–Dec (2016 vs 2018 vs 2020)
# If ANOVA is significant then run Tukey HSD post-hoc
# ==============================================================================
# Packages
library(dplyr)
library(lubridate)
library(ggplot2)
runNO2_ANOVA_17h_NovDec <- function(data = PortTalbotAQ, alpha = 0.05) {
# Guard rails
req <- c("Nitrogen_Dioxide", "DateTime")
if (!all(req %in% names(data))) {
stop("Required columns missing: 'Nitrogen_Dioxide' and/or 'DateTime'.")
}
# Filter for Nov–Dec at 17:00 for years 2016, 2018, 2020
df <- data %>%
mutate(
Year = year(DateTime),
Month = month(DateTime),
Hour = hour(DateTime)
) %>%
filter(
Year %in% c(2016, 2018, 2020),
Month %in% c(11, 12),
Hour == 17
) %>%
transmute(
DateTime,
NO2 = as.numeric(Nitrogen_Dioxide),
Period = factor(
Year,
levels = c(2016, 2018, 2020),
labels = c("Nov–Dec 2016", "Nov–Dec 2018", "Nov–Dec 2020")
)
) %>%
filter(is.finite(NO2))
# Summary table per period
summary_table <- df %>%
group_by(Period) %>%
summarise(
N = n(),
Mean = round(mean(NO2), 3),
SD = round(sd(NO2), 3),
.groups = "drop"
)
cat("\n================ SUMMARY (NO2 at 17:00) ================\n")
print(summary_table)
# Check for enough data (≥2 per period)
n_ok <- all(summary_table$N[match(levels(df$Period), summary_table$Period)] >= 2)
if (!n_ok || nrow(summary_table) < 3) {
cat("\nNot enough data in at least one required period for ANOVA.\n",
"Need Nov–Dec of 2016, 2018, 2020 at 17:00 with ≥2 observations per period.\n",
"H0: mean(NO2) equal across periods; H1: at least one differs.\n",
"Decision: Cannot perform ANOVA due to insufficient data.\n")
return(invisible(list(summary = summary_table)))
}
# One way ANOVA
fit <- aov(NO2 ~ Period, data = df)
anova_tab <- summary(fit)
pval <- anova_tab[[1]][["Pr(>F)"]][1]
cat("\n===================== ANOVA =====================\n")
print(anova_tab)
# Visual with boxplot and jitter
p <- ggplot(df, aes(Period, NO2, fill = Period)) +
geom_boxplot(alpha = 0.8, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.35, size = 1) +
labs(
title = "NO₂ at 17:00 — Nov–Dec (2016 vs 2018 vs 2020)",
x = "Period",
y = "NO₂ (µg/m³)"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "none")
print(p)
# Commentary And Post-hoc
cat("\n================= COMMENTARY =================\n")
cat("Experimental hypotheses:\n")
cat(" H0: mean(NO2) Nov–Dec 2016 = mean(NO2) Nov–Dec 2018 = mean(NO2) Nov–Dec 2020\n")
cat(" H1: At least one period's mean is different.\n\n")
cat(sprintf("ANOVA p-value = %.4f\n", pval))
if (is.finite(pval) && pval < alpha) {
cat(sprintf("Decision at α = %.2f: REJECT H0. A significant difference exists.\n\n", alpha))
cat("Post-hoc rationale:\n")
cat(" Because the overall ANOVA is significant, apply Tukey's HSD to\n")
cat(" identify which specific period pairs differ while controlling the\n")
cat(" family-wise Type I error rate.\n\n")
tuk <- TukeyHSD(fit)
print(tuk)
tuk_df <- as.data.frame(tuk$Period)
tuk_df$Comparison <- rownames(tuk_df)
rownames(tuk_df) <- NULL
names(tuk_df) <- sub("p adj", "p.adj", names(tuk_df), fixed = TRUE)
tuk_df <- tuk_df[, c("Comparison", "diff", "lwr", "upr", if ("p.adj" %in% names(tuk_df)) "p.adj" else "p.value")]
names(tuk_df) <- c("Comparison", "Mean_Diff", "CI_Lower", "CI_Upper", "Adj_P")
tuk_df <- transform(
tuk_df,
Mean_Diff = round(Mean_Diff, 3),
CI_Lower = round(CI_Lower, 3),
CI_Upper = round(CI_Upper, 3),
Adj_P = signif(Adj_P, 4)
)
cat("\nTukey HSD pairwise differences with adjusted p-values:\n")
print(tuk_df)
return(invisible(list(summary = summary_table, anova = anova_tab, tukey = tuk_df)))
} else {
cat(sprintf("Decision at α = %.2f: FAIL TO REJECT H0. No significant difference detected.\n", alpha))
return(invisible(list(summary = summary_table, anova = anova_tab)))
}
}
# analysis
runNO2_ANOVA_17h_NovDec()
##
## ================ SUMMARY (NO2 at 17:00) ================
## # A tibble: 3 × 4
## Period N Mean SD
## <fct> <int> <dbl> <dbl>
## 1 Nov–Dec 2016 60 36.4 17.8
## 2 Nov–Dec 2018 59 23.7 14.7
## 3 Nov–Dec 2020 59 22.9 11.7
##
## ===================== ANOVA =====================
## Df Sum Sq Mean Sq F value Pr(>F)
## Period 2 6866 3433 15.35 7.2e-07 ***
## Residuals 175 39133 224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ================= COMMENTARY =================
## Experimental hypotheses:
## H0: mean(NO2) Nov–Dec 2016 = mean(NO2) Nov–Dec 2018 = mean(NO2) Nov–Dec 2020
## H1: At least one period's mean is different.
##
## ANOVA p-value = 0.0000
## Decision at α = 0.05: REJECT H0. A significant difference exists.
##
## Post-hoc rationale:
## Because the overall ANOVA is significant, apply Tukey's HSD to
## identify which specific period pairs differ while controlling the
## family-wise Type I error rate.
##
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = NO2 ~ Period, data = df)
##
## $Period
## diff lwr upr p adj
## Nov–Dec 2018-Nov–Dec 2016 -12.7239363 -19.204890 -6.242983 0.0000202
## Nov–Dec 2020-Nov–Dec 2016 -13.5170197 -19.997973 -7.036066 0.0000057
## Nov–Dec 2020-Nov–Dec 2018 -0.7930834 -7.301211 5.715044 0.9553023
##
##
## Tukey HSD pairwise differences with adjusted p-values:
## Comparison Mean_Diff CI_Lower CI_Upper Adj_P
## 1 Nov–Dec 2018-Nov–Dec 2016 -12.724 -19.205 -6.243 2.015e-05
## 2 Nov–Dec 2020-Nov–Dec 2016 -13.517 -19.998 -7.036 5.658e-06
## 3 Nov–Dec 2020-Nov–Dec 2018 -0.793 -7.301 5.715 9.553e-01
# I researched a question regarding the possibility of the presence of changes in the nitrogen dioxide (NO2) level during three selected periods of time(November to December 2016, 2018, and 2020. I filtered the dataset to take the measurements of NO2 at 17:00 in each day in the following three years. I conducted a one way ANOVA to test the null hypothesis that the data mean on NO2 at all three periods was the same. In case the ANOVA would result in a statistically significant difference, I used Tukey HSD post-hoc afterwards to understand which pairs of years would have significant differences, and why this post-hoc test was necessary. This step enabled me to identify the change in the concentration of NO2 at the three different time periods.