# 1. SETUP & DATA PRE-PROCESSING

# Load required libraries
library(readxl)      # Data ingestion
library(tidyverse)   # Data manipulation and visualization (ggplot2, dplyr)
library(lubridate)   # Date/Time engineering
library(imputeTS)    # Advanced time-series imputation
library(GGally)      # Advanced multivariate visualizations

cat("\n 1. DATA INGESTION & CLEANING \n")

# Import Data (Ensure this file path matches your computer)
file_path <- "D:/Data Sets/air+quality/AirQualityUCI.xlsx"
AirQuality <- read_excel(file_path)

# Handle Missing Data
# Convert sensor error codes (-200) to NA, then drop empty Excel rows
AirQuality[AirQuality == -200] <- NA
AirQuality <- AirQuality %>% filter(!is.na(Date))

# Advanced Imputation
# Use linear interpolation to estimate missing hours, preserving the chronological timeline
AirQuality$`CO(GT)` <- na_interpolation(AirQuality$`CO(GT)`)
AirQuality$`NO2(GT)` <- na_interpolation(AirQuality$`NO2(GT)`)
AirQuality$`NOx(GT)` <- na_interpolation(AirQuality$`NOx(GT)`)
AirQuality$`C6H6(GT)` <- na_interpolation(AirQuality$`C6H6(GT)`)
AirQuality$`PT08.S5(O3)` <- na_interpolation(AirQuality$`PT08.S5(O3)`)
AirQuality$T <- na_interpolation(AirQuality$T)
AirQuality$RH <- na_interpolation(AirQuality$RH)

# Feature Engineering
# Create Timestamp, Hour, Month, Season, and Weekend variables for grouping
AirQuality$Timestamp <- ymd_hms(paste(format(AirQuality$Date, "%Y-%m-%d"),
                                      format(AirQuality$Time, "%H:%M:%S")))
AirQuality$Hour <- as.numeric(format(AirQuality$Timestamp, "%H"))
AirQuality$Month <- as.numeric(format(AirQuality$Timestamp, "%m"))
AirQuality$DayOfWeek <- factor(weekdays(AirQuality$Timestamp),
                               levels = c("Monday", "Tuesday", "Wednesday",
                                          "Thursday", "Friday", "Saturday", "Sunday"))
AirQuality$IsWeekend <- ifelse(AirQuality$DayOfWeek %in% c("Saturday", "Sunday"),
                               "Weekend", "Weekday")

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

# Order seasons chronologically for logical plotting
AirQuality$Season <- factor(AirQuality$Season, levels = c("Spring", "Summer", "Autumn", "Winter"))

cat("Data Pre-processing Complete. Dimensions:\n")
print(dim(AirQuality))

# 2. EXPLORATORY VISUALISATION
cat("\n 2. GENERATING EXPLORATORY VISUALISATIONS \n")

# 2.1 The Diurnal Cycle
# RATIONALE: Line graphs optimally display chronological trends like the 24-hour cycle.
Hourly_Trend <- AirQuality %>%
  group_by(Hour) %>%
  summarise(Mean_NO2 = mean(`NO2(GT)`, na.rm = TRUE))

plot_time <- ggplot(Hourly_Trend, aes(x = Hour, y = Mean_NO2)) +
  geom_line(color = "firebrick", linewidth = 1.2) +
  geom_point(size = 3) +
  scale_x_continuous(breaks = 0:23) +
  labs(title = "Diurnal Cycle: Average Hourly NO2 Concentration",
       x = "Hour (0-23)", y = "Mean NO2") +
  theme_minimal()
print(plot_time)

# 2.2 Relational Scatterplot (Temp vs O3)
# RATIONALE: Scatterplots with regression lines illustrate the correlation between two continuous variables.
plot_scatter <- ggplot(AirQuality, aes(x = T, y = `PT08.S5(O3)`)) +
  geom_point(color = "steelblue", alpha = 0.2) +
  geom_smooth(method = "lm", color = "darkorange", linewidth = 1.2) +
  labs(title = "Relational Analysis: Temperature vs O3 Sensor",
       subtitle = "Weak relationship between temperature and O3 sensor response (r = -0.089)",
       x = "Temperature (C)", y = "PT08.S5 (O3 Indicator)") +
  theme_minimal()
print(plot_scatter)

# 2.3 Seasonal Comparison of Multiple Pollutants
# RATIONALE: Facet plots allow side-by-side comparison of different gases with varying measurement units.
pollutants_seasonal <- AirQuality %>%
  filter(!is.na(Season)) %>%
  select(Season, `CO(GT)`, `NO2(GT)`, `NOx(GT)`, `C6H6(GT)`) %>%
  pivot_longer(cols = -Season, names_to = "Pollutant", values_to = "Concentration")

plot_facet <- ggplot(pollutants_seasonal, aes(x = Season, y = Concentration, fill = Season)) +
  geom_boxplot(outlier.alpha = 0.2) +
  facet_wrap(~ Pollutant, scales = "free_y") +
  scale_fill_brewer(palette = "Set3") +
  labs(title = "Seasonal Distribution of Major Pollutants",
       x = "Season", y = "Concentration") +
  theme_minimal() + theme(legend.position = "none")
print(plot_facet)

# 2.4 The Multivariate Matrix
# RATIONALE: GGally summarises multicollinearity, densities, and scatter relationships in one matrix.
matrix_data <- AirQuality %>% select(`CO(GT)`, `NO2(GT)`, T, RH) %>% drop_na()
plot_matrix <- ggpairs(matrix_data, title = "Multivariate Matrix: Pollutants vs Weather",
  lower = list(continuous = wrap("points", alpha = 0.1, color = "darkblue")),
  diag  = list(continuous = wrap("densityDiag", fill = "steelblue"))) +
  theme_minimal()
print(plot_matrix)

# 3. ASSUMPTION CHECKING & HYPOTHESIS TESTING
cat("\n 3. ASSUMPTION CHECKING: NORMALITY \n")

# 3.1 Visual Proof of Normality (Q-Q Plot)
# RATIONALE: A Q-Q plot visually assesses whether data follows a theoretical normal distribution.
plot_qq <- ggplot(AirQuality, aes(sample = `NO2(GT)`)) +
  stat_qq(color = "steelblue", alpha = 0.4) +
  stat_qq_line(color = "red", linewidth = 1.2) +
  labs(title = "Q-Q Plot: Checking Normality of NO2",
       subtitle = "Deviations from the red line indicate non-normal data",
       x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()
print(plot_qq)

# 3.2 Mathematical Proof of Normality (Shapiro-Wilk)
set.seed(123)  # Ensure reproducibility
sample_no2 <- sample(AirQuality$`NO2(GT)`[!is.na(AirQuality$`NO2(GT)`)], 5000)
print(shapiro.test(sample_no2))
cat("Interpretation: p < 0.05 confirms non-normal distribution. Non-parametric tests required.\n")

cat("\n INFERENCE 1: PREDICTIVE MODELING (Linear Regression) \n")
# H0: Temp and Humidity do not predict NO2.
model_no2 <- lm(`NO2(GT)` ~ T + RH, data = AirQuality)
print(summary(model_no2))

cat("\n INFERENCE 2: TWO-SAMPLE COMPARISON (Welch T-Test) \n")
# H0: Mean CO is equal on Weekdays and Weekends.
print(t.test(`CO(GT)` ~ IsWeekend, data = AirQuality))

# 4. SEASONAL GROUP COMPARISON
cat("\n 4. SEASONAL GROUP COMPARISON (NO2) \n")
# RATIONALE FOR TEST: Because the Shapiro-Wilk test confirmed non-normal distribution,
# we use the Kruskal-Wallis Rank Sum Test instead of a standard ANOVA.
# H0: The median NO2 concentration is equal across all four seasons.
# H1: At least one season has a significantly different median NO2 concentration.
kw_test_result <- kruskal.test(`NO2(GT)` ~ Season, data = AirQuality)
print(kw_test_result)
cat("Interpretation: Reject H0 if p < 0.05, confirming significant seasonal variance in NO2.\n")

plot_seasonal_no2 <- ggplot(AirQuality %>% filter(!is.na(Season)),
                            aes(x = Season, y = `NO2(GT)`, fill = Season)) +
  geom_boxplot(outlier.colour = "red", outlier.alpha = 0.5) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Distribution of NO2 Levels Across Seasons",
       subtitle = "Kruskal-Wallis confirms statistically significant variance",
       x = "Meteorological Season", y = "NO2 Concentration (microg/m^3)") +
  theme_minimal() + theme(legend.position = "none")
print(plot_seasonal_no2)

# 5. ADVANCED PROGRAMMING (Automated Function)
cat("\n 5. AUTOMATED ANALYSIS \n")

# Custom function to calculate stats and plot distributions for any chosen gas dynamically
analyze_pollutant <- function(df, column_name) {
  if (!column_name %in% names(df)) stop("Error: Pollutant column not found.")
  clean_data <- na.omit(df[[column_name]])
  stats <- list(
    Pollutant      = column_name,
    Valid_Readings = length(clean_data),
    Mean           = round(mean(clean_data), 2),
    Median         = round(median(clean_data), 2),
    SD             = round(sd(clean_data), 2)
  )
  p <- ggplot(df, aes(x = .data[[column_name]])) +
    geom_histogram(fill = "seagreen", color = "white", bins = 30) +
    geom_vline(aes(xintercept = mean(.data[[column_name]], na.rm = TRUE)),
               color = "red", linetype = "dashed", linewidth = 1.2) +
    labs(title = paste("Automated Distribution Analysis:", column_name),
         subtitle = "Red dashed line indicates the mean concentration",
         x = column_name, y = "Frequency") +
    theme_minimal()
  print(p)
  return(stats)
}

# Execute automation
cat("\nResults for Benzene:\n")
print(analyze_pollutant(AirQuality, "C6H6(GT)"))
cat("\nResults for Nitrogen Oxides:\n")
print(analyze_pollutant(AirQuality, "NOx(GT)"))