# Student Name:chetanbharthi atit
# Student ID:2522450
# Module: ACCA7022 Advanced Statistical Programming
# Assignment: Data Analysis
# File: AirQualityAnalysis.r
# Purpose: Analyse Port Talbot hourly air-quality data using
#          descriptive statistics, visualisation, correlation
#          analysis and hypothesis testing.
# Date:10-07-2026

# Import Packages

library(readr)
library(dplyr)
library(ggplot2)
library(lubridate)
library(tidyr)
library(corrplot)
library(moments)

# 1st Step - Import and inspect data

# Providing datafile path
data_file <- "Port_Talbot_Hourly_AQ.csv"

# Importing the dataset in a dtaframe named AirQuality
AirQuality <- read_csv(data_file, show_col_types = FALSE)

# Check column names and data types.
cat("\n--- Column names ---\n")
print(names(AirQuality))

cat("\n--- Data structure ---\n")
str(AirQuality)

# Replaces "No data" with NA across the entire data frame
AirQuality[AirQuality == "No data"] <- NA

# Rename long pollutant and weather variable names into shorter, readable names.
# This improves code clarity and reduces errors caused by repeated long column names.
AirQuality <- AirQuality %>%
  rename(
    NO = `Nitric oxide`,
    NO2 = `Nitrogen dioxide`,
    NOx = `Nitrogen oxides as nitrogen dioxide`,
    SO2 = `Sulphur dioxide`,
    CO = `Carbon monoxide`,
    PM10 = `PM10 particulate matter (Hourly measured)`,
    PM25 = `PM2.5 particulate matter (Hourly measured)`,
    WindDirection = `Modelled Wind Direction`,
    WindSpeed = `Modelled Wind Speed`,
    Temperature = `Modelled Temperature`
  )

# Convert date/time columns into proper date-time format
# The dataset uses day/month/year format, so dmy() is used rather than mdy().
AirQuality$DateTime <- as.POSIXct(
  paste(AirQuality$Date, AirQuality$Time), 
  format = "%d/%m/%Y %H:%M:%S", 
  tz = "UTC"
)


# Convert pollutant and weather columns to numeric.
# This is necessary because columns containing "No data" may be imported as character.
numeric_vars <- c("Ozone", "NO", "NO2", "NOx", "SO2", "CO", "PM10", "PM25",
                  "WindDirection", "WindSpeed", "Temperature")

AirQuality <- AirQuality %>%
  mutate(across(all_of(numeric_vars), as.numeric))

cat("\n--- First few rows ---\n")
print(head(AirQuality))


cat("\n--- Data structure after cleaning ---\n")
str(AirQuality)

cat("\n--- Summary of selected numeric variables after cleaning ---\n")
print(summary(AirQuality[numeric_vars]))

# 2nd Step - Exploratory Visualisation

# Monthly aggregation is used for time-series plots because the original data are hourly
# over ten years. Plotting every hourly point would be overcrowded and very difficult to interpret.
monthly_trends <- AirQuality %>%
  mutate(MonthStart = floor_date(DateTime, unit = "month")) %>%
  group_by(MonthStart) %>%
  summarise(
    PM10_mean = mean(PM10, na.rm = TRUE),
    NO2_mean = mean(NO2, na.rm = TRUE),
    PM25_mean = mean(PM25, na.rm = TRUE),
    .groups = "drop"
  )

# Time-series plot for PM10.
# A line plot is suitable here because it shows how average pollutant levels change over time.
plot_pm10_time <- ggplot(monthly_trends, aes(x = MonthStart, y = PM10_mean)) +
  geom_line(linewidth = 0.6, na.rm = TRUE) +
  labs(
    title = "Monthly Mean PM10 Concentration in Margam, Port Talbot",
    x = "Year",
    y = "PM10 concentration"
  ) +
  theme_minimal()

print(plot_pm10_time)

# Time-series plot for NO2.
# NO2 is selected as a second pollutant to compare whether gaseous pollution shows
# a different long-term pattern from particulate pollution.
plot_no2_time <- ggplot(monthly_trends, aes(x = MonthStart, y = NO2_mean)) +
  geom_line(linewidth = 0.6, na.rm = TRUE) +
  labs(
    title = "Monthly Mean NO2 Concentration in Margam, Port Talbot",
    x = "Year",
    y = "NO2 concentration"
  ) +
  theme_minimal()

print(plot_no2_time)


# Scatterplot of PM10 and wind speed.
# A scatterplot is appropriate because both variables are continuous, and the smoother
# helps assess whether pollutant concentration changes as wind speed increases.
plot_pm10_wind <- ggplot(AirQuality, aes(x = WindSpeed, y = PM10)) +
  geom_point(alpha = 0.15, na.rm = TRUE) +
  geom_smooth(method = "loess", se = TRUE, na.rm = TRUE) +
  labs(
    title = "Relationship between PM10 and Wind Speed",
    x = "Modelled wind speed",
    y = "PM10 concentration"
  ) +
  theme_minimal()

print(plot_pm10_wind)

# Scatterplot of NO2 and temperature.
# This plot explores whether a weather-related variable is associated with a gaseous pollutant.
plot_no2_temp <- ggplot(AirQuality, aes(x = Temperature, y = NO2)) +
  geom_point(alpha = 0.15, na.rm = TRUE) +
  geom_smooth(method = "loess", se = TRUE, na.rm = TRUE) +
  labs(
    title = "Relationship between NO2 and Temperature",
    x = "Modelled temperature",
    y = "NO2 concentration"
  ) +
  theme_minimal()

print(plot_no2_temp)


# 3rd Step - Correlation analysis

# Five pollutants have been selected here for the correlation analysis. 
pollutants <- c("PM10", "PM25", "NO2", "NOx", "Ozone")

correlation_data <- AirQuality %>%
  select(all_of(pollutants)) %>%
  drop_na()

# Check skewness before selecting the correlation method.
# Air-quality data often have high positive skew because most observations are moderate,
# but occasional pollution episodes create extreme high values.
skewness_results <- correlation_data %>%
  summarise(across(everything(), ~ moments::skewness(.x, na.rm = TRUE)))

cat("\n--- Skewness values for selected pollutants ---\n")
print(round(skewness_results, 3))

# Since four of the five pollutants show skewness above 1, Spearman correlation is selected.
# Spearman is rank-based and is more robust for skewed data and extreme pollution episodes
# than Pearson correlation.
cor_matrix <- cor(correlation_data, method = "spearman")

cat("\n--- Spearman correlation matrix for selected pollutants ---\n")
print(round(cor_matrix, 3))

# Display a correlogram of the Spearman correlation matrix.
# A heatmap provides a clear visual summary of the strength and direction
# of pairwise relationships between pollutants.

corrplot(
  cor_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  number.cex = 0.75,
  tl.cex = 0.9,
  tl.col = "black",
  diag = FALSE,
  title = "Spearman Correlations between Selected Pollutants",
  mar = c(0, 0, 2, 0)
)

# Scatterplot matrix to visually examine pairwise relationships
# between the selected pollutants.

pairs(
  correlation_data,
  main = "Scatterplot Matrix of Selected Pollutants",
  pch = 16,
  cex = 0.25
)

# Interpretation:
# - PM10 and PM2.5 show a strong positive correlation, suggesting similar emission sources.
# - NO2 and NOx have a very strong positive correlation, indicating they often increase together.
# - Ozone is moderately negatively correlated with NO2 and NOx, suggesting an inverse relationship.
# - PM10 has weak positive correlations with NO2 and NOx, indicating limited association.
# - Correlation shows association only and does not imply causation.

# 4th Step - Hypothesis Test — Two Conditions

#    PM10 under low vs high wind speed


# Research focus:
# Does PM10 differ between low and high wind-speed conditions?
#
# Null hypothesis H0:
# There is no statistically significant difference in PM10 levels between
# low and high wind-speed conditions.
#
# Alternative hypothesis H1:
# There is a statistically significant difference in PM10 levels between
# low and high wind-speed conditions.

# Select the variables required for the analysis and remove missing values.
pm10_wind_test <- AirQuality %>%
  select(PM10, WindSpeed) %>%
  drop_na()

# Use the median wind speed as a data-driven threshold to create two groups.
# This avoids selecting an arbitrary cut-off value.
median_wind <- median(pm10_wind_test$WindSpeed, na.rm = TRUE)

pm10_wind_test <- pm10_wind_test %>%
  mutate(
    WindGroup = if_else(
      WindSpeed <= median_wind,
      "Low wind speed",
      "High wind speed"
    )
  )

# Display descriptive statistics for each wind-speed group.
cat("\n--- PM10 Summary by Wind-Speed Group ---\n")

print(
  pm10_wind_test %>%
    group_by(WindGroup) %>%
    summarise(
      n = n(),
      Mean = mean(PM10),
      Median = median(PM10),
      SD = sd(PM10),
      .groups = "drop"
    )
)


# Check normality before selecting the statistical test


# The Shapiro-Wilk test assesses whether PM10 follows a normal distribution.
# Because the dataset is very large, a random sample (maximum 5000 observations)
# is used since the Shapiro-Wilk test becomes overly sensitive for large samples.

set.seed(123)

sample_pm10 <- pm10_wind_test %>%
  group_by(WindGroup) %>%
  group_modify(~ slice_sample(.x, n = min(nrow(.x), 5000))) %>%
  ungroup()

cat("\n--- Shapiro-Wilk Normality Test ---\n")

shapiro_results <- by(
  sample_pm10$PM10,
  sample_pm10$WindGroup,
  shapiro.test
)

print(shapiro_results)

# The Shapiro-Wilk tests were significant for both wind-speed groups
# (p < 0.05), so the normality assumption was violated. Therefore,
# the Wilcoxon rank-sum test was appropriate.



# Wilcoxon Rank-Sum Test

pm10_wilcox <- wilcox.test(
  PM10 ~ WindGroup,
  data = pm10_wind_test,
  exact = FALSE
)

cat("\n--- Wilcoxon Rank-Sum Test ---\n")
print(pm10_wilcox)

# Interpretation:
# The Wilcoxon test was statistically significant (W = 909754218,
# p < 2.2e-16). Therefore, H0 is rejected.
# PM10 concentrations differed significantly between low and high
# wind-speed conditions. The median PM10 value was higher during
# high wind speed (18.4) than low wind speed (15.4).



# Boxplot

# A boxplot is used because it clearly compares the median,
# spread and potential outliers of PM10 between the two
# independent wind-speed groups.

plot_pm10_wind_group <- ggplot(
  pm10_wind_test,
  aes(x = WindGroup, y = PM10)
) +
  geom_boxplot(outlier.alpha = 0.2) +
  labs(
    title = "PM10 Levels under Low and High Wind-Speed Conditions",
    subtitle = paste(
      "Wind-speed groups split at the median (",
      round(median_wind, 2),
      ")",
      sep = ""
    ),
    x = "Wind-Speed Group",
    y = "PM10 Concentration"
  ) +
  theme_minimal()

print(plot_pm10_wind_group)

# This suggests that, in this dataset, higher wind speed is associated
# with higher PM10 levels.

# 5th Step - Seasonal Group Comparison


# Seasons are defined using standard meteorological groupings:
# Winter = Dec, Jan, Feb; Spring = Mar, Apr, May;
# Summer = Jun, Jul, Aug; Autumn = Sep, Oct, Nov.
AirQuality <- AirQuality %>%
  mutate(
    Season = case_when(
      month(DateTime) %in% c(12, 1, 2) ~ "Winter",
      month(DateTime) %in% c(3, 4, 5) ~ "Spring",
      month(DateTime) %in% c(6, 7, 8) ~ "Summer",
      month(DateTime) %in% c(9, 10, 11) ~ "Autumn",
      TRUE ~ NA_character_
    ),
    Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Autumn"))
  )

pm10_season <- AirQuality %>%
  select(PM10, Season) %>%
  drop_na()

cat("\n--- PM10 summary by season ---\n")
print(
  pm10_season %>%
    group_by(Season) %>%
    summarise(
      n = n(),
      mean_PM10 = mean(PM10, na.rm = TRUE),
      median_PM10 = median(PM10, na.rm = TRUE),
      sd_PM10 = sd(PM10, na.rm = TRUE),
      IQR_PM10 = IQR(PM10, na.rm = TRUE),
      .groups = "drop"
    )
)

# Boxplots are used because they compare distributions across multiple groups,
# showing medians, spread and outliers more clearly than a table alone.
plot_pm10_season <- ggplot(pm10_season, aes(x = Season, y = PM10)) +
  geom_boxplot(outlier.alpha = 0.2, na.rm = TRUE) +
  labs(
    title = "Seasonal Distribution of PM10 Concentrations",
    x = "Season",
    y = "PM10 concentration"
  ) +
  theme_minimal()

print(plot_pm10_season)


# Seasonal hypotheses:
#
# Null hypothesis H0:
# There is no statistically significant difference in PM10 levels across
# Winter, Spring, Summer and Autumn.
#
# Alternative hypothesis H1:
# At least one season has a statistically different PM10 level.

# The Shapiro-Wilk test is used to check whether PM10 is normally
# distributed within each season. A random sample of up to 5000
# observations per season is used because the dataset is large.

set.seed(123)

sample_pm10_season <- pm10_season %>%
  group_by(Season) %>%
  group_modify(~ slice_sample(.x, n = min(nrow(.x), 5000))) %>%
  ungroup()

cat("\n--- Shapiro-Wilk Normality Test for PM10 by Season ---\n")

shapiro_season_results <- by(
  sample_pm10_season$PM10,
  sample_pm10_season$Season,
  shapiro.test
)

print(shapiro_season_results)

# Kruskal-Wallis is selected because there are four independent seasonal groups
# and the Shapiro-Wilk tests indicate that PM10 is not normally distributed
# within the seasonal groups. This is the non-parametric alternative to
# one-way ANOVA.

pm10_kruskal <- kruskal.test(PM10 ~ Season, data = pm10_season)

cat("\n--- Kruskal-Wallis test: PM10 by season ---\n")
print(pm10_kruskal)

# The Kruskal-Wallis test showed a statistically significant
# difference in PM10 concentrations between seasons
# (χ² = 1439.6, df = 3, p < 2.2e-16). Therefore, H0 is rejected.

# The seasonal summaries indicate that Spring had the highest
# median PM10 concentration (18.5), followed by Summer (18.0),
# Winter (16.9) and Autumn (14.5).

# Since the Kruskal-Wallis test was statistically significant,
# post-hoc pairwise Wilcoxon tests were performed to identify
# which seasons differed in PM10 concentrations. Holm adjustment
# was applied to control the family-wise error rate.

cat("\n--- Post-hoc Pairwise Wilcoxon Tests (Holm Adjustment) ---\n")

posthoc_season <- pairwise.wilcox.test(
  pm10_season$PM10,
  pm10_season$Season,
  p.adjust.method = "holm",
  exact = FALSE
)

print(posthoc_season)

# Interpretation:
# The overall Kruskal-Wallis test was statistically significant
# (χ² = 1439.6, df = 3, p < 2.2e-16), so H0 was rejected.

# The post-hoc results showed that all pairwise seasonal comparisons were
# statistically significant after Holm adjustment. The difference
# between Spring and Summer was also significant (adjusted p = 0.035),
# although it was smaller than the differences involving Autumn or Winter.

# Overall, PM10 concentrations varied significantly across seasons,
# with the highest median concentration observed in Spring and the
# lowest in Autumn.