# 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:

# Step 1: Import Dataset and Understand the dataset

# Import Libraries
library(tidyverse)
library(lubridate)
library(corrplot)

# Import Dataset
AirQuality <- read.csv("Port_Talbot_Hourly_AQ.csv",
                       stringsAsFactors = FALSE)
# Check column names
names(AirQuality)

# Check data types
str(AirQuality)

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

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

# Convert pollutant columns into numeric format
pollutants <- c(
  "Ozone",
  "Nitric.oxide",
  "Nitrogen.dioxide",
  "Nitrogen.oxides.as.nitrogen.dioxide",
  "Sulphur.dioxide",
  "Carbon.monoxide",
  "PM10.particulate.matter..Hourly.measured.",
  "PM2.5.particulate.matter..Hourly.measured."
)

weather_vars <- c(
  "Modelled.Wind.Direction",
  "Modelled.Wind.Speed",
  "Modelled.Temperature"
)

# Convert pollutant and weather columns to numeric
AirQuality[pollutants] <- lapply(AirQuality[pollutants], as.numeric)
AirQuality[weather_vars] <- lapply(AirQuality[weather_vars], as.numeric)

# statistics
head(AirQuality)
dim(AirQuality)
summary(AirQuality[pollutants])
summary(AirQuality[weather_vars])

# Step 2: Exploratory Visualisation

# Monthly averages are used because hourly data across several years  can be visually crowded. Aggregation makes long-term patterns clearer.

monthly_data <- AirQuality %>%
  mutate(Month = floor_date(DateTime, "month")) %>%
  group_by(Month) %>%
  summarise(
    PM25 = mean(`PM2.5.particulate.matter..Hourly.measured.`, na.rm = TRUE),
    NO2 = mean(`Nitrogen.dioxide`, na.rm = TRUE),
    PM10 = mean(`PM10.particulate.matter..Hourly.measured.`, na.rm = TRUE),
    WindSpeed = mean(`Modelled.Wind.Speed`, na.rm = TRUE),
    Temperature = mean(`Modelled.Temperature`, na.rm = TRUE)
  )

# Visualizing Monthly Average PM2.5 Levels
ggplot(monthly_data, aes(x = Month, y = PM25)) +
  geom_line() +
  labs(
    title = "Monthly Average PM2.5 Levels in Port Talbot",
    x = "Date",
    y = "PM2.5 concentration"
  ) +
  theme_minimal()

# Visualizing Monthly Average Nitrogen Dioxide Levels
ggplot(monthly_data, aes(x = Month, y = NO2)) +
  geom_line() +
  labs(
    title = "Monthly Average Nitrogen Dioxide Levels in Port Talbot",
    x = "Date",
    y = "Nitrogen dioxide concentration"
  ) +
  theme_minimal()


# Pollutant vs Wind Speed and Pollutant vs Temperature
# Scatterplots are suitable because both variables are continuous.
# They allow the relationship between pollutant levels and weather conditions to be visually assessed.

# PM2.5 vs Modelled Wind Speed
ggplot(AirQuality, aes(x = `Modelled.Wind.Speed`,
                       y = `PM2.5.particulate.matter..Hourly.measured.`)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(
    title = "Relationship Between PM2.5 and Wind Speed",
    x = "Modelled wind speed",
    y = "PM2.5 concentration"
  ) +
  theme_minimal()

# Nitrogen Dioxide vs Modelled Temperature
ggplot(AirQuality, aes(x = `Modelled.Temperature`,
                       y = `Nitrogen.dioxide`)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(
    title = "Relationship Between Nitrogen Dioxide and Temperature",
    x = "Modelled temperature",
    y = "Nitrogen dioxide concentration"
  ) +
  theme_minimal()

# Step 3: Correlation Analysis

# Selecting Pollutants
correlation_pollutants <- c(
  "PM2.5.particulate.matter..Hourly.measured.",
  "PM10.particulate.matter..Hourly.measured.",
  "Nitrogen.dioxide",
  "Sulphur.dioxide"
)

# Removing missing values
cor_data <- AirQuality %>%
  select(all_of(correlation_pollutants)) %>%
  drop_na()

# Check normality before selecting the correlation method

# Q-Q plots are used to assess whether pollutant concentrations
# are approximately normally distributed.

par(mfrow = c(2,2))

qqnorm(sample(cor_data$`PM2.5.particulate.matter..Hourly.measured.`,5000),
       main = "PM2.5")
qqline(sample(cor_data$`PM2.5.particulate.matter..Hourly.measured.`,5000))

qqnorm(sample(cor_data$`PM10.particulate.matter..Hourly.measured.`,5000),
       main = "PM10")
qqline(sample(cor_data$`PM10.particulate.matter..Hourly.measured.`,5000))

qqnorm(sample(cor_data$Nitrogen.dioxide,5000),
       main = "NO2")
qqline(sample(cor_data$Nitrogen.dioxide,5000))

qqnorm(sample(cor_data$Sulphur.dioxide,5000),
       main = "SO2")
qqline(sample(cor_data$Sulphur.dioxide,5000))

par(mfrow=c(1,1))

# Spearman correlation is used because air-quality contain skewness and extreme pollution peaks. 
# Spearman correlation is less sensitive to outliers than Pearson correlation.

# creating correlation matrix
cor_matrix <- cor(cor_data, method = "spearman")

print(cor_matrix)

# creating correlogram
# A correlogram provides a clear visual summary of the strength and direction
# of relationships between pollutants.

corrplot(
  cor_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  tl.cex = 0.8,
  number.cex = 0.7,
  title = "Spearman Correlations Between Pollutants",
  mar = c(0, 0, 2, 0)
)

# interpretation:
# The strongest correlation was observed between PM2.5 and PM10 (ρ = 0.69),
# suggesting that concentrations of fine and coarse particulate matter tend to increase together.

# PM10 and Sulphur dioxide showed a moderate positive correlation (ρ = 0.49),
# while the remaining pollutant pairs displayed weak positive correlations
# (ρ = 0.26–0.39). These weaker relationships suggest that, although the pollutants may share some common sources, they are also influenced by different atmospheric processes and meteorological conditions.

# Step 4: Hypothesis Test — Two Conditions

# Research question:
# Are PM2.5 levels different under low and high wind-speed conditions?

# H0: There is no statistically significant difference in PM2.5 levels between low and high wind-speed conditions.
# H1: There is a statistically significant difference in PM2.5 levels between low and high wind-speed conditions.

# The median wind speed is used as a data-driven threshold 

median_wind <- median(AirQuality$`Modelled.Wind.Speed`, na.rm = TRUE)

wind_test_data <- AirQuality %>%
  mutate(
    WindGroup = ifelse(
      `Modelled.Wind.Speed` <= median_wind,
      "Low wind speed",
      "High wind speed"
    )
  ) %>%
  select(`PM2.5.particulate.matter..Hourly.measured.`, WindGroup) %>%
  drop_na()

# Summarise PM2.5 levels by wind-speed group
wind_group_summary <- wind_test_data %>%
  group_by(WindGroup) %>%
  summarise(
    n = n(),
    mean_PM25 = mean(`PM2.5.particulate.matter..Hourly.measured.`),
    median_PM25 = median(`PM2.5.particulate.matter..Hourly.measured.`),
    sd_PM25 = sd(`PM2.5.particulate.matter..Hourly.measured.`)
  )

print(wind_group_summary)

qqnorm(sample(wind_test_data$`PM2.5.particulate.matter..Hourly.measured.`, 5000))
qqline(sample(wind_test_data$`PM2.5.particulate.matter..Hourly.measured.`, 5000))


# A Wilcoxon rank-sum test is used because PM2.5 concentrations are
# not normally distributed. This non-parametric test compares the
# distributions of two independent groups without assuming normality.
# This test compares distributions between two independent groups without requiring normality.

wind_test <- wilcox.test(
  `PM2.5.particulate.matter..Hourly.measured.` ~ WindGroup,
  data = wind_test_data
)

print(wind_test)

# Interpretation:
# The mean and median PM2.5 concentrations are very similar under both high and low wind-speed conditions, suggesting only small differences in pollutant levels between the two groups.

# The Wilcoxon rank-sum test produced a p-value of 0.1533, which is greater than the 0.05 significance level. 
# Therefore, the null hypothesis is not rejected. 
# This indicates that there is no statistically significant difference in PM2.5 concentrations between low and high wind-speed

# Step 5: Seasonal Group Comparison

# Research question:
# Do PM2.5 levels vary across Winter, Spring, Summer and Autumn?

# H0: There is no statistically significant difference in PM2.5 levels
# across seasons.
# H1: At least one season has a statistically different PM2.5 level.

season_data <- AirQuality %>%
  mutate(
    MonthNumber = month(DateTime),
    Season = case_when(
      MonthNumber %in% c(12, 1, 2) ~ "Winter",
      MonthNumber %in% c(3, 4, 5) ~ "Spring",
      MonthNumber %in% c(6, 7, 8) ~ "Summer",
      MonthNumber %in% c(9, 10, 11) ~ "Autumn"
    ),
    Season = factor(Season,
                    levels = c("Winter", "Spring", "Summer", "Autumn"))
  ) %>%
  select(`PM2.5.particulate.matter..Hourly.measured.`, Season) %>%
  drop_na()

# Boxplots are used to compare pollutant distributions across seasons because they show medians, spread and outliers clearly for each group.

ggplot(season_data,
       aes(x = Season,
           y = `PM2.5.particulate.matter..Hourly.measured.`)) +
  geom_boxplot() +
  labs(
    title = "PM2.5 Levels Across Seasons in Port Talbot",
    x = "Season",
    y = "PM2.5 concentration"
  ) +
  theme_minimal()

# Summarise PM2.5 levels by season
season_summary <- season_data %>%
  group_by(Season) %>%
  summarise(
    n = n(),
    mean_PM25 = mean(`PM2.5.particulate.matter..Hourly.measured.`),
    median_PM25 = median(`PM2.5.particulate.matter..Hourly.measured.`),
    sd_PM25 = sd(`PM2.5.particulate.matter..Hourly.measured.`)
  )

print(season_summary)

qqnorm(sample(season_data$`PM2.5.particulate.matter..Hourly.measured.`,5000))
qqline(sample(season_data$`PM2.5.particulate.matter..Hourly.measured.`,5000))

# The Q-Q plot shows  that PM2.5 concentrations are not normally distributed. The data also appear positively skewed with several extreme observations.
# A Kruskal-Wallis test is used because there are four independent seasonal groups and pollutant data are not normally distributed.

season_test <- kruskal.test(
  `PM2.5.particulate.matter..Hourly.measured.` ~ Season,
  data = season_data
)

print(season_test)

# As the Kruskal-Wallis test is significant, pairwise Wilcoxon tests are used to identify which seasons differ. 
# Bonferroni adjustment controls the risk of false positives from multiple comparisons.

pairwise_season_test <- pairwise.wilcox.test(
  season_data$`PM2.5.particulate.matter..Hourly.measured.`,
  season_data$Season,
  p.adjust.method = "bonferroni"
)

print(pairwise_season_test)

# Interpretation:
# The Kruskal-Wallis test indicated a statistically significant difference
# in PM2.5 concentrations across the four seasons
# (χ² = 1343.8, df = 3, p < 0.001).
# Therefore, the null hypothesis was rejected, indicating that at least one
# season has different PM2.5 concentrations.

# Summary statistics suggest that PM2.5 concentrations were highest during
# Spring (mean = 11.10 µg/m³) and lowest during Autumn (mean = 7.96 µg/m³),
# while Winter and Summer showed intermediate concentrations.

# Pairwise Wilcoxon tests with Bonferroni adjustment showed statistically
# significant differences between every pair of seasons (all adjusted
# p-values < 0.05). This suggests that seasonal variation has a measurable
# influence on PM2.5 concentrations in Port Talbot.