############################################################
# Student Name: Your Name
# Module: ACCA7022 Advanced Statistical Programming
# Assignment: Data Analysis
# File Name: AirQualityAnalysis.r
# Purpose: Air quality analysis of Port Talbot hourly dataset
############################################################

# -------------------------------
# 1. Load Required Packages
# -------------------------------

packages <- c("tidyverse", "lubridate", "corrplot", "GGally", "car")

install_if_missing <- function(pkg) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg)
    library(pkg, character.only = TRUE)
  }
}

lapply(packages, install_if_missing)

# -------------------------------
# 2. Import Dataset
# -------------------------------

AirQuality <- read.csv("Port_Talbot_Hourly_AQ.csv")

# -------------------------------
# 3. Initial Data Inspection
# -------------------------------

cat("Column Names:\n")
print(names(AirQuality))

cat("\nData Structure:\n")
str(AirQuality)

cat("\nFirst Few Rows:\n")
print(head(AirQuality))

cat("\nDataset Dimensions:\n")
print(dim(AirQuality))

cat("\nSummary Statistics:\n")
print(summary(AirQuality))

# -------------------------------
# 4. Rename Important Columns
# -------------------------------

library(dplyr)

AirQuality <- AirQuality %>%
  rename(
    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
  )


# -------------------------------
# 5. Convert Date and Time
# -------------------------------

# Use original date column if it exists, otherwise use Date
if ("Date_Original" %in% names(AirQuality)) {
  AirQuality$DateText <- as.character(AirQuality$Date_Original)
} else {
  AirQuality$DateText <- as.character(AirQuality$Date)
}

AirQuality$Time <- as.character(AirQuality$Time)

# Try both possible date formats
AirQuality$DateTime <- as.POSIXct(
  paste(AirQuality$DateText, AirQuality$Time),
  format = "%Y-%m-%d %H:%M:%S",
  tz = "UTC"
)

# If still NA, try dd/mm/yyyy format
if (all(is.na(AirQuality$DateTime))) {
  AirQuality$DateTime <- as.POSIXct(
    paste(AirQuality$DateText, AirQuality$Time),
    format = "%d/%m/%Y %H:%M:%S",
    tz = "UTC"
  )
}

AirQuality$Date <- as.Date(AirQuality$DateTime)
AirQuality$Month <- lubridate::month(AirQuality$DateTime, label = TRUE)
AirQuality$Year <- lubridate::year(AirQuality$DateTime)

cat("\nDateTime check:\n")
print(head(AirQuality$DateTime))
print(sum(is.na(AirQuality$DateTime)))

# -------------------------------
# 6. Convert Pollutants to Numeric
# -------------------------------

numeric_columns <- c(
  "Ozone",
  "Nitric_Oxide",
  "Nitrogen_Dioxide",
  "Nitrogen_Oxides",
  "Sulphur_Dioxide",
  "Carbon_Monoxide",
  "PM10",
  "PM25",
  "Wind_Direction",
  "Wind_Speed",
  "Temperature"
)

AirQuality[numeric_columns] <- lapply(
  AirQuality[numeric_columns],
  function(x) as.numeric(as.character(x))
)

cat("\nCleaned Data Structure:\n")
str(AirQuality)

cat("\nMissing Values:\n")
print(colSums(is.na(AirQuality)))

# -------------------------------
# 7. Descriptive Statistics
# -------------------------------

descriptive_stats <- AirQuality %>%
  select(all_of(numeric_columns)) %>%
  summarise_all(
    list(
      Mean = ~mean(., na.rm = TRUE),
      Median = ~median(., na.rm = TRUE),
      SD = ~sd(., na.rm = TRUE),
      Min = ~min(., na.rm = TRUE),
      Max = ~max(., na.rm = TRUE)
    )
  )

cat("\nDescriptive Statistics:\n")
print(descriptive_stats)

# -------------------------------
# 8. Time Series Visualisations
# -------------------------------

# Time series plots are used because the dataset contains hourly
# measurements across time. These plots help show long-term trends
# and changes in pollutant concentration.

ggplot(AirQuality, aes(x = DateTime, y = PM25)) +
  geom_line(na.rm = TRUE) +
  labs(
    title = "Time Trend of PM2.5 in Port Talbot",
    x = "Date",
    y = "PM2.5 Concentration"
  ) +
  theme_minimal()

ggplot(AirQuality, aes(x = DateTime, y = Nitrogen_Dioxide)) +
  geom_line(na.rm = TRUE) +
  labs(
    title = "Time Trend of Nitrogen Dioxide in Port Talbot",
    x = "Date",
    y = "Nitrogen Dioxide Concentration"
  ) +
  theme_minimal()

# -------------------------------
# 9. Pollutants and Weather Relationships
# -------------------------------

# Scatterplots are suitable because they show whether two continuous
# variables have a visible relationship.

ggplot(AirQuality, aes(x = Wind_Speed, y = PM25)) +
  geom_point(alpha = 0.3, na.rm = TRUE) +
  geom_smooth(method = "lm", se = FALSE, na.rm = TRUE) +
  labs(
    title = "Relationship Between Wind Speed and PM2.5",
    x = "Wind Speed",
    y = "PM2.5 Concentration"
  ) +
  theme_minimal()

ggplot(AirQuality, aes(x = Temperature, y = Nitrogen_Dioxide)) +
  geom_point(alpha = 0.3, na.rm = TRUE) +
  geom_smooth(method = "lm", se = FALSE, na.rm = TRUE) +
  labs(
    title = "Relationship Between Temperature and Nitrogen Dioxide",
    x = "Temperature",
    y = "Nitrogen Dioxide Concentration"
  ) +
  theme_minimal()

ggplot(AirQuality, aes(x = Temperature, y = Ozone)) +
  geom_point(alpha = 0.3, na.rm = TRUE) +
  geom_smooth(method = "lm", se = FALSE, na.rm = TRUE) +
  labs(
    title = "Relationship Between Temperature and Ozone",
    x = "Temperature",
    y = "Ozone Concentration"
  ) +
  theme_minimal()

# -------------------------------
# 10. Correlation Analysis
# -------------------------------

pollutants <- AirQuality %>%
  select(
    Ozone,
    Nitric_Oxide,
    Nitrogen_Dioxide,
    Nitrogen_Oxides,
    Sulphur_Dioxide,
    Carbon_Monoxide,
    PM10,
    PM25
  )

correlation_matrix <- cor(
  pollutants,
  use = "complete.obs",
  method = "pearson"
)

cat("\nCorrelation Matrix:\n")
print(round(correlation_matrix, 2))

corrplot(
  correlation_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  tl.cex = 0.8,
  number.cex = 0.7,
  title = "Correlation Heatmap of Air Pollutants",
  mar = c(0, 0, 2, 0)
)

# Pairwise scatterplot matrix
GGally::ggpairs(
  pollutants %>% select(Ozone, Nitrogen_Dioxide, PM10, PM25),
  title = "Pairwise Relationships Between Selected Pollutants"
)

# -------------------------------
# 11. Hypothesis Testing
# -------------------------------

# Research question:
# Does PM2.5 concentration differ between low wind speed and high wind speed?

# Null hypothesis H0:
# There is no significant difference in PM2.5 levels between low and high wind speed.

# Alternative hypothesis H1:
# There is a significant difference in PM2.5 levels between low and high wind speed.

median_wind <- median(AirQuality$Wind_Speed, na.rm = TRUE)

AirQuality <- AirQuality %>%
  mutate(
    Wind_Group = ifelse(Wind_Speed >= median_wind, "High Wind Speed", "Low Wind Speed")
  )

hypothesis_data <- AirQuality %>%
  select(PM25, Wind_Group) %>%
  drop_na()

cat("\nPM2.5 by Wind Group Summary:\n")
print(
  hypothesis_data %>%
    group_by(Wind_Group) %>%
    summarise(
      Mean = mean(PM25),
      Median = median(PM25),
      SD = sd(PM25),
      Count = n()
    )
)

ggplot(hypothesis_data, aes(x = Wind_Group, y = PM25)) +
  geom_boxplot(na.rm = TRUE) +
  labs(
    title = "PM2.5 Levels Under Low and High Wind Speed",
    x = "Wind Speed Group",
    y = "PM2.5 Concentration"
  ) +
  theme_minimal()

# Normality test using sample because large datasets can make
# Shapiro-Wilk overly sensitive.

set.seed(123)

sample_low <- hypothesis_data %>%
  filter(Wind_Group == "Low Wind Speed") %>%
  pull(PM25) %>%
  sample(min(5000, length(.)))

sample_high <- hypothesis_data %>%
  filter(Wind_Group == "High Wind Speed") %>%
  pull(PM25) %>%
  sample(min(5000, length(.)))

cat("\nShapiro-Wilk Test for Low Wind Speed PM2.5:\n")
print(shapiro.test(sample_low))

cat("\nShapiro-Wilk Test for High Wind Speed PM2.5:\n")
print(shapiro.test(sample_high))

# Wilcoxon test is used because air pollution data are often skewed
# and may not satisfy normality assumptions.

wilcox_result <- wilcox.test(PM25 ~ Wind_Group, data = hypothesis_data)

cat("\nWilcoxon Test Result for PM2.5 by Wind Speed Group:\n")
print(wilcox_result)

# Optional t-test for comparison
t_test_result <- t.test(PM25 ~ Wind_Group, data = hypothesis_data)

cat("\nIndependent Samples t-test Result:\n")
print(t_test_result)

# -------------------------------
# 12. Seasonal Group Comparison
# -------------------------------

AirQuality <- AirQuality %>%
  mutate(
    Season = case_when(
      month(Date) %in% c(12, 1, 2) ~ "Winter",
      month(Date) %in% c(3, 4, 5) ~ "Spring",
      month(Date) %in% c(6, 7, 8) ~ "Summer",
      month(Date) %in% c(9, 10, 11) ~ "Autumn",
      TRUE ~ NA_character_
    ),
    Season = factor(
      Season,
      levels = c("Winter", "Spring", "Summer", "Autumn")
    )
  )

season_data <- AirQuality %>%
  select(PM25, Season) %>%
  drop_na()

cat("\nPM2.5 by Season Summary:\n")
print(
  season_data %>%
    group_by(Season) %>%
    summarise(
      Mean = mean(PM25),
      Median = median(PM25),
      SD = sd(PM25),
      Count = n()
    )
)

ggplot(season_data, aes(x = Season, y = PM25)) +
  geom_boxplot(na.rm = TRUE) +
  labs(
    title = "Seasonal Comparison of PM2.5 Levels",
    x = "Season",
    y = "PM2.5 Concentration"
  ) +
  theme_minimal()

# Hypotheses:
# H0: PM2.5 levels are the same across all seasons.
# H1: At least one season has significantly different PM2.5 levels.

# ANOVA test
anova_result <- aov(PM25 ~ Season, data = season_data)

cat("\nANOVA Result for PM2.5 Across Seasons:\n")
print(summary(anova_result))

cat("\nTukey Post-Hoc Test:\n")
print(TukeyHSD(anova_result))

# Kruskal-Wallis test as a non-parametric alternative
kruskal_result <- kruskal.test(PM25 ~ Season, data = season_data)

cat("\nKruskal-Wallis Test Result:\n")
print(kruskal_result)

# -------------------------------
# 13. Monthly Boxplot
# -------------------------------

ggplot(AirQuality, aes(x = Month, y = PM25)) +
  geom_boxplot(na.rm = TRUE) +
  labs(
    title = "Monthly Distribution of PM2.5 Levels",
    x = "Month",
    y = "PM2.5 Concentration"
  ) +
  theme_minimal()

# -------------------------------
# 14. Save Cleaned Data and Outputs
# -------------------------------

write.csv(AirQuality, "Cleaned_Port_Talbot_AirQuality.csv", row.names = FALSE)

write.csv(
  round(correlation_matrix, 3),
  "Pollutant_Correlation_Matrix.csv",
  row.names = TRUE
)

cat("\nAnalysis completed successfully.\n")