# Clear environment
rm(list = ls())

# Load required libraries
library(tidyverse)
library(lubridate)
library(ggplot2)
library(corrplot)

# ------------------------------------------------------------------------------
# SECTION 1: DATA IMPORT
# ------------------------------------------------------------------------------

AirQuality <- read.csv("Port_Talbot_Hourly_AQ.csv", stringsAsFactors = FALSE)

cat("--- Section 1: Data Import ---\n\n")

# Check dimensions
cat("Dataset dimensions:",
    nrow(AirQuality),
    "rows x",
    ncol(AirQuality),
    "columns\n\n")

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

# Check data types
cat("Data types:\n")
str(AirQuality)
cat("\n")

# First few rows
cat("First 5 rows:\n")
print(head(AirQuality, 5))
cat("\n")

# Basic descriptive statistics
cat("Basic descriptive statistics:\n")
print(summary(AirQuality))

# ------------------------------------------------------------------------------
# SECTION 2: MISSING DATA ANALYSIS
# ------------------------------------------------------------------------------

cat("\n--- Section 2: Missing Data Analysis ---\n\n")

# Identify frequency of "No data" strings to assess data quality
count_no_data <- function(x) {
  sum(x == "No data", na.rm = TRUE)
}

missing_counts <- sapply(AirQuality, count_no_data)

# Create a summary table to quantify gaps in the dataset
missing_summary <- data.frame(
  Column = names(AirQuality),
  Missing_Count = missing_counts,
  Percent_Missing = round((missing_counts / nrow(AirQuality)) * 100, 2)
)

# Filter to show only columns requiring attention and sort by impact
missing_summary <- missing_summary[missing_summary$Missing_Count > 0, ]
missing_summary <- missing_summary[order(-missing_summary$Missing_Count), ]

if (nrow(missing_summary) > 0) {
  cat("Columns with missing data:\n")
  print(missing_summary, row.names = FALSE)
} else {
  cat("No missing data identified.\n")
}

cat("\nSummary of missing data findings:\n")
cat("Total observations: ", nrow(AirQuality), "\n")
if (nrow(missing_summary) > 0) {
  cat("Max missing count:  ",
      max(missing_summary$Missing_Count),
      "\n")
  cat("Min missing count:  ",
      min(missing_summary$Missing_Count),
      "\n\n")
}

# ------------------------------------------------------------------------------
# SECTION 3: DATA CLEANING & CONVERSION
# ------------------------------------------------------------------------------

cat("--- Section 3: Data Cleaning & Conversion ---\n\n")

# Helper function to convert columns to numeric, treating "No data" as NA
extract_numeric <- function(x) {
  x[x == "No data"] <- NA
  as.numeric(sub(" .*", "", x))
}

# Adjust timestamps: Convert "24:00:00" to "00:00:00" of the next day
# to allow correct parsing into a POSIXct format.
AirQuality$Date_parsed <- as.Date(AirQuality$Date, format = "%d/%m/%Y")
is_hour24 <- AirQuality$Time == "24:00:00"

AirQuality$Date_parsed[is_hour24] <- AirQuality$Date_parsed[is_hour24] + 1
AirQuality$Time[is_hour24] <- "00:00:00"
AirQuality$Date <- format(AirQuality$Date_parsed, "%d/%m/%Y")
AirQuality$Date_parsed <- NULL

cat("Corrected",
    sum(is_hour24),
    "entries with 24:00:00 timestamp.\n\n")

# Process and restructure the dataset for analysis
AirQuality <- AirQuality %>%
  mutate(
    DateTime = as.POSIXct(paste(Date, Time), format = "%d/%m/%Y %H:%M:%S", tz = "UTC"),
    Date = as.Date(DateTime),
    Hour = hour(DateTime),
    Month = month(DateTime),
    Year = year(DateTime),
    DayOfWeek = wday(DateTime, label = TRUE),
    Season = case_when(
      Month %in% c(12, 1, 2) ~ "Winter",
      Month %in% c(3, 4, 5) ~ "Spring",
      Month %in% c(6, 7, 8) ~ "Summer",
      Month %in% c(9, 10, 11) ~ "Autumn"
    ),
    # Convert pollutant and weather variables to numeric format
    Ozone = extract_numeric(Ozone),
    Nitric_Oxide = extract_numeric(Nitric.oxide),
    NO2 = extract_numeric(Nitrogen.dioxide),
    NOx = extract_numeric(Nitrogen.oxides.as.nitrogen.dioxide),
    SO2 = extract_numeric(Sulphur.dioxide),
    CO = extract_numeric(Carbon.monoxide),
    PM10 = extract_numeric(PM10.particulate.matter..Hourly.measured.),
    PM25 = extract_numeric(PM2.5.particulate.matter..Hourly.measured.),
    Wind_Direction = extract_numeric(Modelled.Wind.Direction),
    Wind_Speed = extract_numeric(Modelled.Wind.Speed),
    Temperature = extract_numeric(Modelled.Temperature)
  ) %>%
  # Select only the cleaned variables for the final analysis data frame
  select(
    DateTime,
    Date,
    Hour,
    Month,
    Year,
    DayOfWeek,
    Season,
    Ozone,
    Nitric_Oxide,
    NO2,
    NOx,
    SO2,
    CO,
    PM10,
    PM25,
    Wind_Direction,
    Wind_Speed,
    Temperature
  )

cat("Data cleaning complete.\n")
cat("Dimensions:",
    nrow(AirQuality),
    "rows x",
    ncol(AirQuality),
    "columns.\n\n")

# Verify processing results
cat("Missing DateTime values:", sum(is.na(AirQuality$DateTime)), "\n")
cat("First 5 rows of cleaned dataset:\n")
print(head(AirQuality, 5))

# ------------------------------------------------------------------------------
# SECTION 4: DESCRIPTIVE STATISTICS
# ------------------------------------------------------------------------------

cat("--- Section 4: Descriptive Statistics ---\n\n")

# Calculate summary statistics for all pollutants at once
pollutant_summary <- AirQuality %>%
  select(PM25, PM10, NO2, SO2, Ozone, CO) %>%
  summarise(across(
    everything(),
    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)
    )
  )) %>%
  pivot_longer(everything(),
               names_to = c("Pollutant", "Stat"),
               names_sep = "_") %>%
  pivot_wider(names_from = Stat, values_from = value)

cat("Pollutant Summary Statistics (µg/m³):\n")
print(as.data.frame(pollutant_summary), row.names = FALSE)

# Weather variables summary
cat("\nWeather Variables:\n")
weather_summary <- AirQuality %>%
  summarise(
    Temp_Mean = mean(Temperature, na.rm = TRUE),
    Temp_SD   = sd(Temperature, na.rm = TRUE),
    Wind_Mean = mean(Wind_Speed, na.rm = TRUE),
    Wind_SD   = sd(Wind_Speed, na.rm = TRUE)
  )
print(as.data.frame(weather_summary))

# PM2.5 statistics segmented by season
cat("\nPM2.5 Summary by Season:\n")
seasonal_stats <- AirQuality %>%
  group_by(Season) %>%
  summarise(
    N      = n(),
    Mean   = mean(PM25, na.rm = TRUE),
    Median = median(PM25, na.rm = TRUE),
    SD     = sd(PM25, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  arrange(factor(Season, levels = c("Winter", "Spring", "Summer", "Autumn")))

print(as.data.frame(seasonal_stats), row.names = FALSE)

# ------------------------------------------------------------------------------
# SECTION 5: EXPLORATORY VISUALIZATIONS
# ------------------------------------------------------------------------------
# This section explores pollutant patterns, identifies trends over time,
# and examines relationships between pollution levels and weather conditions.

cat("--- Section 5: Exploratory Visualizations ---\n\n")

# PM2.5 variation across seasons.
# The boxplot compares pollution levels between seasons and highlights
# differences in distribution, spread, and potential outliers.

p1 <- ggplot(AirQuality, aes(x = PM25)) +
  geom_histogram(
    bins = 50,
    fill = "#2E86AB",
    color = "black",
    alpha = 0.7
  ) +
  labs(title = "Distribution of PM2.5 Concentrations", x = "PM2.5 (µg/m³)", y = "Frequency (hours)") +
  theme_minimal()
ggsave(
  "01_PM25_Distribution.png",
  p1,
  width = 10,
  height = 6,
  dpi = 150
)
cat("✓ Plot 1: Histogram generated\n")

# Comparing pollution levels across the seasons.
# This helps us pick up on any seasonal spikes and keeps an eye
# on potential outliers for different times of the year.
p2 <- ggplot(AirQuality, aes(x = Season, y = PM25, fill = Season)) +
  geom_boxplot(alpha = 0.7, outlier.size = 1) +
  scale_fill_manual(
    values = c(
      "Winter" = "#1F77B4",
      "Spring" = "#FF7F0E",
      "Summer" = "#2CA02C",
      "Autumn" = "#D62728"
    )
  ) +
  labs(title = "PM2.5 Levels by Season", x = "Season", y = "PM2.5 (µg/m³)") +
  theme_minimal() + theme(legend.position = "none")
ggsave(
  "02_PM25_Seasonal.png",
  p2,
  width = 10,
  height = 6,
  dpi = 150
)
cat("✓ Plot 2: Seasonal boxplot generated\n")

# Time trends for PM2.5 and NO2.
# By looking at these together on a monthly basis, we can see if they share
# the same long-term patterns, which might point to similar sources.
monthly_avg <- AirQuality %>%
  mutate(YearMonth = floor_date(Date, "month")) %>%
  group_by(YearMonth) %>%
  summarise(
    PM25_avg = mean(PM25, na.rm = TRUE),
    NO2_avg = mean(NO2, na.rm = TRUE),
    .groups = 'drop'
  )

p3 <- ggplot(monthly_avg, aes(x = YearMonth)) +
  geom_line(aes(y = PM25_avg, color = "PM2.5"), linewidth = 1) +
  geom_line(aes(y = NO2_avg, color = "NO2"), linewidth = 1) +
  scale_color_manual(values = c("PM2.5" = "#2E86AB", "NO2" = "#A23B72")) +
  labs(title = "Monthly Average Trends: PM2.5 and NO2", x = "Date", y = "Concentration (µg/m³)") +
  theme_minimal()

ggsave(
  "03_Timeseries_Trends.png",
  p3,
  width = 10,
  height = 6,
  dpi = 150
)
cat("✓ Plot 3: Time series generated\n")

# PM2.5 relationship with wind speed.
# The scatterplot examines whether wind speed influences pollutant levels.
# A negative relationship would suggest that stronger winds help disperse PM2.5.

p4 <- ggplot(AirQuality, aes(x = Wind_Speed, y = PM25)) +
  geom_point(alpha = 0.2,
             size = 1,
             color = "#2E86AB") +
  geom_smooth(
    method = "lm",
    color = "#A23B72",
    fill = "#A23B72",
    alpha = 0.2
  ) +
  labs(title = "Weather Relationship: PM2.5 vs Wind Speed", x = "Wind Speed (m/s)", y = "PM2.5 (µg/m³)") +
  theme_minimal()

ggsave(
  "04_PM25_WindSpeed.png",
  p4,
  width = 10,
  height = 6,
  dpi = 150
)
cat("✓ Plot 4: Wind dispersion scatterplot generated\n")

# Hourly PM2.5 patterns.
# Breaking this down by the hour of the day is a great way to spot
# spikes that might line up with rush hour traffic.
hourly_avg <- AirQuality %>%
  group_by(Hour) %>%
  summarise(Mean_PM25 = mean(PM25, na.rm = TRUE), .groups = 'drop')

p5 <- ggplot(hourly_avg, aes(x = Hour, y = Mean_PM25)) +
  geom_col(fill = "#2E86AB", alpha = 0.7) +
  scale_x_continuous(breaks = seq(0, 23, 2)) +
  labs(title = "Average PM2.5 by Hour of Day", x = "Hour (0-23)", y = "Mean PM2.5 (µg/m³)") +
  theme_minimal()

ggsave(
  "05_PM25_HourlyPattern.png",
  p5,
  width = 10,
  height = 6,
  dpi = 150
)
cat("✓ Plot 5: Diurnal cycle bar chart generated\n")

# Temperature vs PM2.5.
# Finally, let’s see if temperature plays a role in these pollution
# levels. The seasonal coloring adds a bit more context to the relationship.
p6 <- ggplot(AirQuality, aes(x = Temperature, y = PM25, color = Season)) +
  geom_point(alpha = 0.1, size = 1) +
  geom_smooth(
    method = "lm",
    fill = NA,
    linewidth = 1,
    se = FALSE
  ) +
  scale_color_manual(
    values = c(
      "Winter" = "#1F77B4",
      "Spring" = "#FF7F0E",
      "Summer" = "#2CA02C",
      "Autumn" = "#D62728"
    )
  ) +
  labs(title = "Weather Relationship: Temperature vs PM2.5", x = "Temperature (°C)", y = "PM2.5 (µg/m³)") +
  theme_minimal()

ggsave(
  "06_PM25_Temperature.png",
  p6,
  width = 10,
  height = 6,
  dpi = 150
)
cat("✓ Plot 6: Seasonal temperature impact scatterplot generated\n\n")

# ------------------------------------------------------------------------------
# SECTION 6: CORRELATION ANALYSIS
# ------------------------------------------------------------------------------
# This section examines relationships between the main pollutants.
# Spearman correlation is used because air quality variables are often skewed
# and the method is less affected by non-normal distributions and outliers.

cat("--- Section 6: Correlation Analysis ---\n\n")

# Filtering for complete data to ensure our correlation matrix is accurate.
pollutant_data <- AirQuality %>%
  select(PM25, PM10, NO2, SO2, Ozone, CO) %>%
  filter(complete.cases(.))

cat("Rows used for correlation:",
    nrow(pollutant_data),
    "/",
    nrow(AirQuality),
    "\n\n")

# Calculating the Spearman correlation matrix.
spearman_matrix <- cor(pollutant_data, method = "spearman", use = "complete.obs")

cat("Spearman Correlation Matrix:\n")
print(round(spearman_matrix, 3))

# Breaking down the key relationships so they're easy to interpret.
cat("\nINTERPRETATION OF KEY RELATIONSHIPS:\n")
cat(
  "- PM2.5 & PM10: r =",
  round(spearman_matrix["PM25", "PM10"], 3),
  "(Strong positive, likely sharing the same source)\n"
)
cat(
  "- PM2.5 & NO2:  r =",
  round(spearman_matrix["PM25", "NO2"], 3),
  "(Moderate positive, common in traffic-heavy areas)\n"
)
cat(
  "- NO2 & Ozone:  r =",
  round(spearman_matrix["NO2", "Ozone"], 3),
  "(Moderate negative, a sign of typical urban photochemistry)\n\n"
)

# ------------------------------------------------------------------------------
# VISUALIZATION: CORRELATION HEATMAP
# ------------------------------------------------------------------------------
# A heatmap is great for getting a "big picture" view of how everything correlates.

png(
  "07_Pollutant_Correlations_Heatmap.png",
  width = 800,
  height = 700,
  res = 100
)
corrplot(
  spearman_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  diag = FALSE,
  tl.col = "black",
  tl.srt = 45,
  title = "Pollutant Correlation Matrix",
  mar = c(2, 2, 3, 2)
)
dev.off()
cat("✓ Saved: 07_Pollutant_Correlations_Heatmap.png\n\n")

# ------------------------------------------------------------------------------
# VISUALIZATION: CORRELATION SCATTERPLOTS
# ------------------------------------------------------------------------------
# Scatterplots are used to visualise selected pollutant relationships.
# Trend lines help identify the direction of association between variables.

# Helper function to describe the strength of the relationship
classify_r <- function(r) {
  if (abs(r) >= 0.65)
    return(ifelse(r > 0, "Strong Positive", "Strong Negative"))
  if (abs(r) >= 0.4)
    return(ifelse(r > 0, "Moderate Positive", "Moderate Negative"))
  if (abs(r) < 0.05)
    return("No Correlation")
  return(ifelse(r > 0, "Weak Positive", "Weak Negative"))
}

plot_pairs <- list(c("PM25", "PM10"), c("PM25", "NO2"), c("NO2", "Ozone"))
filenames <- c(
  "08_Correlation_PM25_PM10.png",
  "09_Correlation_PM25_NO2.png",
  "10_Correlation_NO2_Ozone.png"
)

for (i in 1:length(plot_pairs)) {
  x_v <- plot_pairs[[i]][1]
  y_v <- plot_pairs[[i]][2]
  r_val <- spearman_matrix[x_v, y_v]
  
  p <- ggplot(AirQuality, aes_string(x = x_v, y = y_v)) +
    geom_point(alpha = 0.2,
               size = 1,
               color = "#2E86AB") +
    geom_smooth(
      method = "lm",
      color = "#A23B72",
      fill = "#A23B72",
      alpha = 0.2,
      linewidth = 1.2
    ) +
    labs(
      title = sprintf("%s vs %s (r=%.3f)", x_v, y_v, r_val),
      subtitle = paste("Relationship:", classify_r(r_val))
    ) +
    theme_minimal()
  
  ggsave(filenames[i],
         p,
         width = 10,
         height = 6,
         dpi = 150)
  cat("✓ Saved:", filenames[i], "\n")
}

# ------------------------------------------------------------------------------
# SECTION 7: HYPOTHESIS TESTING
# ------------------------------------------------------------------------------

# This section examines whether PM2.5 concentrations differ between low and
# high wind speed conditions, and across seasons, using non-parametric tests.

cat("--- Section 7: Hypothesis Testing ---\n\n")

# Pre-filter for analysis
AQ_clean <- AirQuality %>%
  filter(PM25 >= 0, Wind_Speed >= 0, !is.na(PM25), !is.na(Wind_Speed))

# ------------------------------------------------------------------------------
# SECTION 7.1: WIND SPEED EFFECT ON PM2.5
# ------------------------------------------------------------------------------

cat("Hypothesis 1: Effect of Wind Speed on PM2.5\n")
cat("Null Hypothesis (H0): PM2.5 is the same in high vs low wind conditions\n")
cat("Alternative Hypothesis (H1): PM2.5 differs between high and low wind conditions\n\n")

wind_q25 <- quantile(AQ_clean$Wind_Speed, 0.25, na.rm = TRUE)
wind_q75 <- quantile(AQ_clean$Wind_Speed, 0.75, na.rm = TRUE)

AQ_wind <- AQ_clean %>%
  mutate(
    Wind_Category = case_when(
      Wind_Speed <= wind_q25 ~ "Low",
      Wind_Speed >= wind_q75 ~ "High",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(Wind_Category))

low_wind_pm25  <- AQ_wind %>% filter(Wind_Category == "Low")  %>% pull(PM25)
high_wind_pm25 <- AQ_wind %>% filter(Wind_Category == "High") %>% pull(PM25)

# Check normality before choosing test method (sample if >5000, Shapiro-Wilk limit)
if (length(low_wind_pm25) > 5000) {
  sample_low  <- sample(low_wind_pm25, 5000)
  sample_high <- sample(high_wind_pm25, 5000)
} else {
  sample_low  <- low_wind_pm25
  sample_high <- high_wind_pm25
}

shapiro_low  <- shapiro.test(sample_low)
shapiro_high <- shapiro.test(sample_high)

cat("Normality Check (Shapiro-Wilk):\n")
cat("Low Wind:  p-value =",
    format.pval(shapiro_low$p.value, digits = 3),
    "\n")
cat("High Wind: p-value =",
    format.pval(shapiro_high$p.value, digits = 3),
    "\n")
cat("Conclusion: Data is not normally distributed → Mann-Whitney U test used.\n\n")

mw_test <- wilcox.test(PM25 ~ Wind_Category, data = AQ_wind)

cat("Mann-Whitney U Test Results:\n")
cat("P-value: ", format.pval(mw_test$p.value, digits = 3), "\n")
cat(
  "Conclusion: ",
  ifelse(
    mw_test$p.value < 0.05,
    "Statistically significant.",
    "Not significant."
  ),
  "\n\n"
)

# ------------------------------------------------------------------------------
# VISUALIZATION: Wind Speed Effect Boxplot
# ------------------------------------------------------------------------------
# Visualizes the distribution of PM2.5 concentrations across seasons.
# Matches the seasonal pattern shown in Figure 2 (Section 5).
# Supports the statistical results from the Kruskal–Wallis test.

p_wind <- ggplot(AQ_wind, aes(x = Wind_Category, y = PM25, fill = Wind_Category)) +
  geom_boxplot(alpha = 0.7) +
  scale_fill_manual(values = c("Low" = "#2E86AB", "High" = "#A23B72")) +
  labs(title = "PM2.5 Distribution by Wind Speed Category", x = "Wind Speed Category", y = "PM2.5 (µg/m³)") +
  theme_minimal()

ggsave(
  "11_Wind_PM25_Boxplot.png",
  p_wind,
  width = 8,
  height = 6,
  dpi = 150
)
cat("✓ Plot 11: Wind speed effect boxplot generated\n\n")

# ------------------------------------------------------------------------------
# SECTION 7.2: SEASONAL DIFFERENCES IN PM2.5
# ------------------------------------------------------------------------------

cat("Hypothesis 2: Seasonal Differences in PM2.5\n")
cat("Null Hypothesis (H0): PM2.5 distribution is identical across seasons\n")
cat("Alternative Hypothesis (H1): At least one season differs significantly\n\n")

# Kruskal-Wallis is the appropriate non-parametric alternative to ANOVA
# for comparing >2 groups with non-normal data.

kw_test <- kruskal.test(PM25 ~ Season, data = AQ_clean)

cat("Kruskal-Wallis Test Results:\n")
cat(
  "H-statistic: ",
  round(kw_test$statistic, 2),
  ", p-value: ",
  format.pval(kw_test$p.value, digits = 3),
  "\n\n"
)

# Post-hoc pairwise comparisons with Bonferroni correction
# We multiply p-values by the number of pairs to control Type I error inflation
cat("[Post-Hoc] Seasonal Pairwise Wilcoxon Tests (Bonferroni corrected):\n")

seasons <- unique(as.character(AQ_clean$Season))
combos <- combn(seasons, 2, simplify = FALSE)
num_tests <- length(combos)

pairwise_results <- lapply(combos, function(pair) {
  sub_data <- AQ_clean %>% filter(Season %in% pair)
  test <- wilcox.test(PM25 ~ Season, data = sub_data)
  data.frame(
    Comparison = paste(pair, collapse = " vs "),
    P_Value = min(1, round(test$p.value * num_tests, 4)) # Bonferroni adjustment
  )
}) %>% bind_rows()

print(pairwise_results, row.names = FALSE)
cat("\n")

# ------------------------------------------------------------------------------
# VISUALIZATION: Seasonal Comparison Boxplot
# ------------------------------------------------------------------------------

# WHY: Boxplot visualizes distribution of PM2.5 across seasons.
# References Figure 2 (Section 5) which shows same seasonal pattern.
# Confirms statistical findings from Kruskal-Wallis test.

p_season <- ggplot(AQ_clean, aes(x = Season, y = PM25, fill = Season)) +
  geom_boxplot(alpha = 0.7) +
  scale_fill_manual(
    values = c(
      "Winter" = "#1F77B4",
      "Spring" = "#FF7F0E",
      "Summer" = "#2CA02C",
      "Autumn" = "#D62728"
    )
  ) +
  labs(title = "PM2.5 Distribution by Season", x = "Season", y = "PM2.5 (µg/m³)") +
  theme_minimal()

ggsave(
  "12_Seasonal_PM25_Boxplot.png",
  p_season,
  width = 10,
  height = 6,
  dpi = 150
)
cat("✓ Plot 12: Seasonal comparison boxplot generated\n\n")

# ------------------------------------------------------------------------------
# Summary of Methodological Choices
# ------------------------------------------------------------------------------

test_summary <- data.frame(
  Analysis = c("Correlation", "Wind Effect", "Seasonal Effect"),
  Method   = c("Spearman", "Mann-Whitney U", "Kruskal-Wallis"),
  Reason   = c(
    "Non-linear/Skewed",
    "2 Groups, Non-normal",
    ">2 Groups, Non-normal"
  )
)

cat("Summary of Statistical Methods:\n")
print(test_summary, row.names = FALSE)
cat("\n")
cat("--- ANALYSIS COMPLETE ---\n")