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

# Loading necessary packahges
library(tidyverse)
library(lubridate)
library(janitor)
library(skimr)
library(corrplot)
library(rstatix)

# Step 1 - Import Dataset and Understand the Data

# Import the data
AirQuality <- read.csv("Port_Talbot_Hourly_AQ.csv",
                       header = TRUE,
                       stringsAsFactors = FALSE)

# Check original column names 
cat("\nOriginal column names:\n")
print(names(AirQuality))

# Check datatypes
cat("\nData types of column:\n")
str(AirQuality)

# Convert "No data" into missing values
AirQuality[AirQuality == "No data"] <- NA

# Convert date/time columns into proper date-time format.
AirQuality$DateTime <- dmy_hms(
  paste(AirQuality$Date, AirQuality$Time)
)

# Renaming the main analysis variables to concise, meaningful names. The original
# variables are retained only where needed; concise names make the analysis more readable 
AirQuality <- AirQuality %>%
  rename(
    ozone = 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.,
    pm25 = PM2.5.particulate.matter..Hourly.measured.,
    wind_direction =  Modelled.Wind.Direction,
    wind_speed = Modelled.Wind.Speed,
    temperature = Modelled.Temperature
  )

# Convert pollutant columns into numeric format
numeric_vars <- c(
  "ozone", "nitric_oxide", "nitrogen_dioxide", "nitrogen_oxides",
  "sulphur_dioxide", "carbon_monoxide", "pm10", "pm25",
  "wind_direction", "wind_speed", "temperature"
)

# Convert selected variables to numeric. This is a safeguard in case any column
# has been read as character due to non-numeric values in the raw file.
AirQuality <- AirQuality %>%
  mutate(across(all_of(numeric_vars), as.numeric))

# Check data types again after conversion
str(AirQuality)

cat("\nFirst few rows of the cleaned dataset:\n")
print(head(AirQuality))

cat("\nDataset dimensions: rows and columns\n")
print(dim(AirQuality))

cat("\nBasic summary statistics:\n")
print(summary(AirQuality[numeric_vars]))


# Step 2: Exploratory Visualisation

# Long-term time trend: PM2.5 and NO2 

# Hourly observations over ten years can make raw line plots overcrowded. Monthly
# averages are therefore used for time-series visualisation because they preserve
# the long-term trend while making the pattern easier to interpret.

monthly_pollution <- AirQuality %>%
  mutate(month_start = floor_date(DateTime, unit = "month")) %>%
  group_by(month_start) %>%
  summarise(
    mean_pm25 = mean(pm25, na.rm = TRUE),
    mean_no2 = mean(nitrogen_dioxide, na.rm = TRUE),
    mean_pm10 = mean(pm10, na.rm = TRUE),
    mean_wind_speed = mean(wind_speed, na.rm = TRUE),
    mean_temperature = mean(temperature, na.rm = TRUE),
    .groups = "drop"
  )

# A time-series line plot is used because the data are repeated observations over
# time. It helps identify long-term movement, seasonal cycles and unusual periods
# of high or low concentration for PM2.5.
p_pm25_time <- ggplot(monthly_pollution, aes(x = month_start, y = mean_pm25)) +
  geom_line(linewidth = 0.6, na.rm = TRUE) +
  geom_point(size = 0.8, na.rm = TRUE) +
  labs(
    title = "Monthly Mean PM2.5 Concentration in Port Talbot, 2015-2024",
    x = "Date",
    y = "Mean PM2.5 concentration"
  ) +
  theme_minimal()

print(p_pm25_time)

# Nitrogen dioxide is also plotted as a time-series because comparing two
# pollutants over time helps assess whether different pollutants show similar or
# contrasting temporal behaviour.
p_no2_time <- ggplot(monthly_pollution, aes(x = month_start, y = mean_no2)) +
  geom_line(linewidth = 0.6, na.rm = TRUE) +
  geom_point(size = 0.8, na.rm = TRUE) +
  labs(
    title = "Monthly Mean Nitrogen Dioxide Concentration in Port Talbot, 2015-2024",
    x = "Date",
    y = "Mean nitrogen dioxide concentration"
  ) +
  theme_minimal()

print(p_no2_time)

# Relationship between pollutants and weather variables

# Scatterplots are appropriate here because both PM2.5 and wind speed are
# continuous variables. The smoothing line helps reveal the overall relationship
# without assuming that the relationship is perfectly linear.
p_pm25_wind <- AirQuality %>%
  filter(!is.na(pm25), !is.na(wind_speed)) %>%
  ggplot(aes(x = wind_speed, y = pm25)) +
  geom_point(alpha = 0.15, size = 0.6) +
  geom_smooth(method = "loess", se = TRUE, linewidth = 0.8) +
  labs(
    title = "Relationship between Wind Speed and PM2.5",
    x = "Modelled wind speed (m/s)",
    y = "PM2.5 concentration"
  ) +
  theme_minimal()

print(p_pm25_wind)

# This weather-related scatterplot explores whether nitrogen dioxide varies
# with local temperature. This is suitable because both variables are continuous,
# and the plot can show whether the relationship is positive, negative or weak.
p_no2_temp <- AirQuality %>%
  filter(!is.na(nitrogen_dioxide), !is.na(temperature)) %>%
  ggplot(aes(x = temperature, y = nitrogen_dioxide)) +
  geom_point(alpha = 0.15, size = 0.6) +
  geom_smooth(method = "loess", se = TRUE, linewidth = 0.8) +
  labs(
    title = "Relationship between Temperature and Nitrogen Dioxide",
    x = "Modelled temperature (degrees C)",
    y = "Nitrogen dioxide concentration"
  ) +
  theme_minimal()

print(p_no2_temp)

# Step 3: Correlation Analysis
# Six pollutants are compared: PM2.5, PM10, NO2,NO, Ozone, and SO2
pollutant_vars <- c("pm25", "pm10", "nitrogen_dioxide", "nitric_oxide", "ozone", "sulphur_dioxide")

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

# Check distribution shape to choose Pearson vs Spearman
# Histograms are used to visually check for skewness. Air-quality
# concentrations are typically right-skewed (many low readings, occasional
# high pollution episodes), which violates the approximate-normality
# assumption behind Pearson correlation.
par(mfrow = c(3, 3))
for (v in pollutant_vars) {
  hist(correlation_data[[v]], main = paste("Distribution of", v),
       xlab = v, col = "lightgrey", breaks = 40)
}
par(mfrow = c(3, 3))

# Spearman correlation is selected because air-quality data is majorly right-skewed
# Spearman correlation is rank-based and is therefore less sensitive to outliers than Pearson correlation.
spearman_cor <- cor(correlation_data, method = "spearman")

cat("\nSpearman correlation matrix for selected pollutants:\n")
print(round(spearman_cor, 3))

# Visualise the correlation matrix as a heatmap
# A heatmap is used because it displays the full correlation matrix at once
# and uses colour intensity to make the strength of each relationship
# immediately visible, which is clearer than reading a table of numbers.
heatmap(spearman_cor, symm = TRUE, Rowv = NA, Colv = NA,
        col = colorRampPalette(c("blue", "white", "red"))(50),
        main = "Spearman Correlation Heatmap: Pollutants")

# A scatterplot matrix (pairs plot) is also produced to visually confirm the
# direction and monotonic shape of each pairwise relationship, supporting
# the correlation coefficients reported above.
pairs(correlation_data[sample(seq_len(nrow(correlation_data)), 3000), ],
      pch = 16, col = rgb(0, 0, 0, 0.1),
      main = "Scatterplot Matrix of Pollutant Concentrations")

# PM2.5 and PM10 show a strong positive correlation (r = 0.688),
# indicating that concentrations of fine and coarse particulate matter
# tend to increase and decrease together.

# Nitrogen dioxide (NO2) and nitric oxide (NO) have a strong positive
# correlation (r = 0.626), suggesting they are likely emitted from similar
# combustion-related sources.

# Sulphur dioxide (SO2) has moderate positive correlations with PM10,
# PM2.5, NO2 and NO (r = 0.359–0.489), indicating these pollutants may
# occur together under similar environmental conditions.

# Ozone (O3) is moderately negatively correlated with NO2 (r = -0.418)
# and weakly negatively correlated with NO (r = -0.281). This suggests
# that higher nitrogen oxide concentrations are generally associated with
# lower ozone levels.

# Step 4: Hypothesis Test — Two Conditions

# nitrogen_dioxide is compared between low and high wind-speed conditions, because
# wind speed influences the atmospheric dispersion of gaseous pollutants.
#
# H0: There is no statistically significant difference in nitrogen_dioxide levels
#     between low and high wind-speed conditions.
# H1: There is a statistically significant difference in nitrogen_dioxide levels
#     between low and high wind-speed conditions.

# The median wind speed is used as a data-driven threshold to split the
# data into two groups, rather than an arbitrary cut-off value.
wind_median <- median(AirQuality$wind_speed, na.rm = TRUE)

wind_test_data <- AirQuality %>%
  filter(!is.na(nitrogen_dioxide), !is.na(wind_speed)) %>%
  mutate(
    wind_group = if_else(wind_speed <= wind_median, "Low wind speed", "High wind speed"),
    wind_group = factor(wind_group, levels = c("Low wind speed", "High wind speed"))
  )

cat("\nMedian wind speed used as threshold (m/s):\n")
print(wind_median)

cat("\nnitrogen_dioxide summary by wind-speed group:\n")
wind_group_summary <- wind_test_data %>%
  group_by(wind_group) %>%
  summarise(
    n = n(),
    mean_nitrogen_dioxide = mean(nitrogen_dioxide, na.rm = TRUE),
    median_nitrogen_dioxide = median(nitrogen_dioxide, na.rm = TRUE),
    sd_nitrogen_dioxide = sd(nitrogen_dioxide, na.rm = TRUE),
    iqr_nitrogen_dioxide = IQR(nitrogen_dioxide, na.rm = TRUE),
    .groups = "drop"
  )
print(wind_group_summary)

# Check distribution shape using Shapiro-Wilk tests on a sample from each group.
# A fixed seed makes the sampling reproducible.
set.seed(7022)
normality_wind <- wind_test_data %>%
  group_by(wind_group) %>%
  slice_head(n = 5000) %>%
  shapiro_test(nitrogen_dioxide)

cat("\nShapiro-Wilk normality checks for nitrogen_dioxide by wind-speed group:\n")
print(normality_wind)

# Both p-values are far below 0.05, indicating that the NO₂ data are not normally distributed. 
# Therefore, choosing the Wilcoxon rank-sum test instead of an independent samples t-test is appropriate.

wind_wilcox <- wilcox.test(nitrogen_dioxide ~ wind_group, data = wind_test_data, exact = FALSE)

cat("\nWilcoxon rank-sum test comparing nitrogen_dioxide by wind-speed group:\n")
print(wind_wilcox)

# A boxplot is added to support interpretation of the test. It shows differences
# in median, spread and outliers between the two wind-speed conditions.
p_nitrogen_dioxide_wind_box <- ggplot(wind_test_data, aes(x = wind_group, y = nitrogen_dioxide)) +
  geom_boxplot(outlier.alpha = 0.15) +
  labs(
    title = "Nitrogen Dioxide Concentration under Low and High Wind-Speed Conditions",
    subtitle = paste0("Groups defined using median wind speed = ", round(wind_median, 2), " m/s"),
    x = "Wind-speed condition",
    y = "Nitrogen dioxide concentration (µg/m³)"
  ) +
  theme_minimal()

print(p_nitrogen_dioxide_wind_box)

# Interpretation:
# he p-value is below 0.05; therefore, the null hypothesis is rejected. 
# Nitrogen dioxide concentrations differ significantly between low and high wind-speed conditions. 
# The higher median concentration observed under low wind-speed conditions suggests 
# that reduced atmospheric dispersion may contribute to pollutant accumulation.
# However, this statistical association does not establish a causal relationship.\n")

# Step 5: Seasonal Group Comparison

# Research question: Do nitrogen dioxide (NO2) concentrations differ across
# Winter, Spring, Summer and Autumn?
#
# Null hypothesis (H0): There is no statistically significant difference in NO2
# levels across seasons.
#
# Alternative hypothesis (H1): At least one season has a statistically different
# NO2 level.
#
# A season variable is created using standard meteorological seasons:
# Winter = December, January, February
# Spring = March, April, May
# Summer = June, July, August
# Autumn = September, October, November

season_data <- AirQuality %>%
  filter(!is.na(nitrogen_dioxide), !is.na(DateTime)) %>%
  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"))
  ) %>%
  filter(!is.na(season))

cat("\nNitrogen dioxide summary by season:\n")
season_summary <- season_data %>%
  group_by(season) %>%
  summarise(
    n = n(),
    mean_NO2 = mean(nitrogen_dioxide, na.rm = TRUE),
    median_NO2 = median(nitrogen_dioxide, na.rm = TRUE),
    sd_NO2 = sd(nitrogen_dioxide, na.rm = TRUE),
    iqr_NO2 = IQR(nitrogen_dioxide, na.rm = TRUE),
    .groups = "drop"
  )
print(season_summary)

# Boxplots are used because the task requires a group comparison and boxplots
# clearly show median, spread, skewness and outliers for each season.
p_NO2_season_box <- ggplot(season_data, aes(x = season, y = nitrogen_dioxide)) +
  geom_boxplot(outlier.alpha = 0.12) +
  labs(
    title = "Seasonal Distribution of Nitrogen Dioxide Concentration in Port Talbot",
    x = "Season",
    y = expression(NO[2]~concentration~(mu*g/m^3))
  ) +
  theme_minimal()

print(p_NO2_season_box)

# Normality and variance checks are included before choosing the seasonal test.
# With large hourly air-quality data and likely skewness, a non-parametric test is
# more defensible than assuming ANOVA conditions automatically hold.
set.seed(7022)
normality_season <- season_data %>%
  group_by(season) %>%
  slice_head(n = 5000) %>%
  shapiro_test(nitrogen_dioxide)

cat("\nShapiro-Wilk normality checks for nitrogen dioxide by season:\n")
print(normality_season)

# Shapiro tests are significant, so normality is not met and Kruskal-Wallis test is ideal for this task. 

variance_check <- levene_test(season_data, nitrogen_dioxide ~ season)

cat("\nLevene test for equality of variance across seasons:\n")
print(variance_check)

# Levene test is significant, so equal variance is not met.

# Kruskal-Wallis is selected because it compares more than two independent groups
# without requiring normality. It tests whether at least one seasonal distribution
# differs from the others.
season_kruskal <- kruskal.test(nitrogen_dioxide ~ season, data = season_data)

cat("\nKruskal-Wallis test comparing nitrogen dioxide across seasons:\n")
print(season_kruskal)

# here, P-value is less than 0.05 so, the Kruskal-Wallis test is significant
# Therefore, pairwise Wilcoxon tests identify
# which seasonal pairs differ. Holm adjustment controls the family-wise error rate
# across multiple comparisons.
season_posthoc <- pairwise.wilcox.test(
x = season_data$nitrogen_dioxide,
g = season_data$season,
p.adjust.method = "holm",
exact = FALSE
)
cat("\nPost-hoc pairwise Wilcoxon tests with Holm p-value adjustment:\n")
print(season_posthoc)

# interpretation:
# The p-value is below 0.05, so the null hypothesis is rejected.
# Nitrogen dioxide concentrations differ significantly across seasons. 
# Post-hoc Wilcoxon tests show all seasonal pairs differ significantly.