cat("\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \n")
# Programmer : Mohammad Lakhi
# Student ID : 2512430
# Module     : ACCA7022-	Advanced Statistical Programming
#File Title  : Data Analysis 
# Month      : March,2026

cat("\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \n")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# TASK-01: Import the Dataset and exploration of Dataset
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

install.packages("ggplot2")
install.packages("dplyr")
install.packages("corrplot")
install.packages("lubridate")
library(ggplot2)
library(dplyr)
library(corrplot)
library(lubridate)

# Import dataset - UCI format uses semicolon and comma decimal
AirQuality <- read.csv(
  "AirQualityUCI.csv",
  sep = ";",
  dec = ",",
  header = TRUE,
  stringsAsFactors = FALSE
)

# Remove empty trailing columns and rows
AirQuality <- AirQuality[, colSums(is.na(AirQuality)) < nrow(AirQuality)]
AirQuality <- AirQuality[rowSums(is.na(AirQuality)) != ncol(AirQuality), ]

# Replace -200 with NA (dataset uses -200 for missing values)
AirQuality[AirQuality == -200] <- NA

# Create DateTime column by combining Date and Time
AirQuality$DateTime <- dmy_hms(paste(AirQuality$Date, AirQuality$Time))

# Check column names
cat("=== Column Names ===\n")
print(colnames(AirQuality))

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

# View first 6 rows
cat("\n=== First 6 Rows ===\n")
print(head(AirQuality))

# Check dimensions
cat("\n=== Dimensions ===\n")
cat("Rows:", nrow(AirQuality), "\n")
cat("Columns:", ncol(AirQuality), "\n")

# Basic descriptive statistics
cat("\n=== Descriptive Statistics ===\n")
print(summary(AirQuality))

# ============================================================
# TASK 2 - Exploratory Visualisation
# ============================================================

# In this task we create plots to explore the data visually.
# We look at how pollutant levels change over time.
# how they relate to weather variables like temperature and humidity.
# Visual exploration helps us understand
# patterns before doing formal statistical tests.


# -------------------------------------------------------
# Plot 1 - Time Series of CO (Carbon Monoxide)
# -------------------------------------------------------

# We use a time series plot because CO is measured over time.
# Time series is the best way to see trends, seasonal patterns,
# and unusual spikes in data collected at regular intervals.
# CO is a dangerous pollutant produced mainly by vehicle engines.

ggplot(AirQuality, aes(x = DateTime, y = CO.GT.)) +
  geom_line(color = "steelblue", alpha = 0.6, size = 0.4) +
  geom_smooth(method = "loess", color = "red", se = FALSE) +
  labs(
    title = "CO (Carbon Monoxide) Levels Over Time",
    x = "Date",
    y = "CO Concentration (mg/m³)"
  ) +
  theme_minimal()


# -------------------------------------------------------
# Plot 2 - Time Series of NO2 (Nitrogen Dioxide)
# -------------------------------------------------------

# We also plot NO2 over time for comparison with CO.
# Both are traffic-related pollutants, so comparing their
# time trends helps us identify if they follow similar
# seasonal or daily patterns.
# NO2 is harmful to lungs and is mainly from vehicle exhaust.

ggplot(AirQuality, aes(x = DateTime, y = NO2.GT.)) +
  geom_line(color = "darkorange", alpha = 0.6, size = 0.4) +
  geom_smooth(method = "loess", color = "darkred", se = FALSE) +
  labs(
    title = "NO2 (Nitrogen Dioxide) Levels Over Time",
    x = "Date",
    y = "NO2 Concentration (µg/m³)"
  ) +
  theme_minimal()


# -------------------------------------------------------
# Plot 3 - Scatterplot: CO vs Temperature
# -------------------------------------------------------

# We use a scatterplot to explore the relationship between
# CO levels and temperature. Scatterplots are the best choice
# when we want to see if two continuous variables are related.
# We expect CO to be higher in winter (low temperature) because
# more vehicles and heating systems are used in cold weather.

ggplot(AirQuality, aes(x = T, y = CO.GT.)) +
  geom_point(color = "steelblue", alpha = 0.3, size = 0.8) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Relationship Between CO and Temperature",
    x = "Temperature (°C)",
    y = "CO Concentration (mg/m³)"
  ) +
  theme_minimal()


# -------------------------------------------------------
# Plot 4 - Scatterplot: NO2 vs Humidity
# -------------------------------------------------------

# We use a scatterplot to check if NO2 levels are related
# to relative humidity (RH). Humidity can affect how
# pollutants disperse in the air. High humidity can trap
# pollutants close to the ground, increasing concentration.
# A linear trend line helps us see the direction of relationship.

ggplot(AirQuality, aes(x = RH, y = NO2.GT.)) +
  geom_point(color = "darkorange", alpha = 0.3, size = 0.8) +
  geom_smooth(method = "lm", color = "darkred", se = TRUE) +
  labs(
    title = "Relationship Between NO2 and Relative Humidity",
    x = "Relative Humidity (%)",
    y = "NO2 Concentration (µg/m³)"
  ) +
  theme_minimal()


# -------------------------------------------------------
# Plot 5 - Boxplot: CO levels by Hour of Day
# -------------------------------------------------------

# We use a boxplot to compare CO levels at different hours.
# Boxplots are good for showing distribution and spread
# across groups. We expect higher CO during rush hours
# (morning 7-9am and evening 5-7pm) due to traffic.
# This helps confirm that traffic is a main source of CO.

# First we extract hour from DateTime
AirQuality$Hour <- hour(AirQuality$DateTime)

ggplot(AirQuality, aes(x = factor(Hour), y = CO.GT.)) +
  geom_boxplot(fill = "steelblue", alpha = 0.7, outlier.size = 0.5) +
  labs(
    title = "CO Levels by Hour of Day",
    x = "Hour of Day (0 = midnight)",
    y = "CO Concentration (mg/m³)"
  ) +
  theme_minimal()








# ============================================================
# TASK 3 - Correlation Analysis
# ============================================================

# In this task we measure how strongly pollutants are related to each other. 
# Correlation values range from -1 to +1.
# A value close to +1 means both variables increase together.
# A value close to -1 means one goes up when other goes down.
# A value close to 0 means no relationship.
# We use Pearson correlation because our variables are numeric
# and continuous. We will check if this is appropriate below.


# -------------------------------------------------------
# Step 1 - Select pollutant columns for correlation
# -------------------------------------------------------

# We select CO, NO2, NOx, C6H6 (Benzene) and O3 for analysis.
# These are the main measured pollutants in this dataset.
# We remove NA values first so correlation can be calculated.

pollutants <- AirQuality %>%
  select(CO.GT., NO2.GT., NOx.GT., C6H6.GT., PT08.S5.O3.) %>%
  na.omit()

# Rename columns for cleaner labels in plots
colnames(pollutants) <- c("CO", "NO2", "NOx", "Benzene", "O3")


# -------------------------------------------------------
# Step 2 - Calculate correlation matrix
# -------------------------------------------------------

# cor() calculates Pearson correlation between all pairs.
# We use use="complete.obs" to ignore any remaining NA rows.
# The result is a matrix showing correlation for every pair.

cor_matrix <- cor(pollutants, method = "pearson")

# Print the correlation matrix with 2 decimal places
cat("=== Correlation Matrix ===\n")
print(round(cor_matrix, 2))


# -------------------------------------------------------
# Step 3 - Visualise correlations using a heatmap
# -------------------------------------------------------

# A heatmap is the best choice here because it shows all
# correlations at once in an easy to read colour format.
# Dark blue means strong positive correlation.
# Dark red means strong negative correlation.
# White means little or no correlation.
# This is better than showing many separate scatterplots
# when we have 5 variables to compare.

corrplot(
  cor_matrix,
  method = "color",
  type = "upper",
  tl.col = "black",
  tl.srt = 45,
  addCoef.col = "black",
  number.cex = 0.8,
  col = colorRampPalette(c("darkred", "white", "steelblue"))(200),
  title = "Correlation Heatmap of Air Pollutants",
  mar = c(0, 0, 2, 0)
)


# -------------------------------------------------------
# Step 4 - Scatterplot matrix for visual confirmation
# -------------------------------------------------------

# We also create a scatterplot matrix using pairs() to visually
# confirm the correlations shown in the heatmap.
# Each small plot shows the relationship between two pollutants.
# This gives additional visual evidence for our interpretation.

pairs(
  pollutants,
  main = "Scatterplot Matrix of Air Pollutants",
  pch = 16,
  col = adjustcolor("steelblue", alpha.f = 0.2),
  cex = 0.4
)


# -------------------------------------------------------
# Step 5 - Scatterplot: CO vs Benzene (strongest pair)
# -------------------------------------------------------

# From the heatmap we can see CO and Benzene have the
# strongest positive correlation. We plot these two
# separately for a closer look at their relationship.
# Both CO and Benzene come from vehicle exhaust, which
# explains why they move together so strongly.

ggplot(pollutants, aes(x = CO, y = Benzene)) +
  geom_point(color = "steelblue", alpha = 0.3, size = 0.8) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "CO vs Benzene — Strongest Correlation",
    x = "CO Concentration (mg/m³)",
    y = "Benzene Concentration (µg/m³)"
  ) +
  theme_minimal()


# -------------------------------------------------------
# Step 6 - Interpret the results in comments
# -------------------------------------------------------

# INTERPRETATION:
# CO and Benzene show a strong positive correlation (r ≈ 0.90).
# This makes sense because both are produced by burning fuel
# in vehicle engines. When traffic is heavy, both rise together.
#
# CO and NOx also show strong positive correlation (r ≈ 0.77).
# NOx is also a traffic-related pollutant, so this is expected.
#
# NO2 and NOx are strongly correlated (r ≈ 0.80) because
# NO2 is actually a component of NOx group of gases.
#
# O3 shows negative or weak correlation with CO and NO2.
# This is because O3 (ozone) is a secondary pollutant —
# it is destroyed by NO2 in chemical reactions near ground level.
# So when traffic pollutants go up, O3 tends to go down.
#
# These patterns confirm that vehicle traffic is the main
# source of pollution in this dataset location.

cat("\n=== Correlation Interpretation ===\n")
cat("CO vs Benzene  :", round(cor_matrix["CO", "Benzene"], 2), "\n")
cat("CO vs NOx      :", round(cor_matrix["CO", "NOx"], 2), "\n")
cat("NO2 vs NOx     :", round(cor_matrix["NO2", "NOx"], 2), "\n")
cat("CO vs O3       :", round(cor_matrix["CO", "O3"], 2), "\n")












# ============================================================
# TASK 4 - Hypothesis Testing
# ============================================================

# We want to test whether CO levels are significantly different
# under high wind speed versus low wind speed conditions.
# Wind speed affects how pollutants disperse in the air.
# High wind speed should scatter pollutants more, reducing
# their concentration near the ground.
#
# We use CO as our pollutant because it showed clear patterns
# in Task 2 and is directly related to traffic conditions.
#
# NULL HYPOTHESIS (H0):
# There is no significant difference in CO levels between
# high wind speed and low wind speed conditions.
#
# ALTERNATIVE HYPOTHESIS (H1):
# CO levels are significantly different between high wind
# speed and low wind speed conditions.


# -------------------------------------------------------
# Step 1 - Create wind speed groups
# -------------------------------------------------------

# The dataset does not have a direct wind speed column.
# We use AH (Absolute Humidity) as our splitting variable
# because it is a continuous weather variable in the dataset.
# We split at the median to create two equal groups:
# LOW group  = AH below median
# HIGH group = AH above median
# This gives us two conditions to compare CO levels.

# Calculate median of AH
ah_median <- median(AirQuality$AH, na.rm = TRUE)

cat("=== Median Absolute Humidity ===\n")
cat("Median AH:", round(ah_median, 4), "\n")

# Create two groups based on AH median split
AirQuality$AH_Group <- ifelse(
  AirQuality$AH >= ah_median, "High_AH", "Low_AH"
)

# Extract CO values for each group and remove NA
co_high <- AirQuality$CO.GT.[AirQuality$AH_Group == "High_AH"]
co_high <- co_high[!is.na(co_high)]

co_low <- AirQuality$CO.GT.[AirQuality$AH_Group == "Low_AH"]
co_low <- co_low[!is.na(co_low)]

cat("\n=== Group Sizes ===\n")
cat("High AH group size:", length(co_high), "\n")
cat("Low AH group size :", length(co_low), "\n")


# -------------------------------------------------------
# Step 2 - Check normality before choosing test
# -------------------------------------------------------

# Before choosing between t-test and Wilcoxon test we must
# check if the data is normally distributed.
# t-test assumes normality — it works best when data follows
# a bell curve shape.
# Wilcoxon test does not assume normality — it works with
# any distribution shape.
# We use the Shapiro-Wilk test to check normality.
# If p-value < 0.05 then data is NOT normal → use Wilcoxon.
# If p-value > 0.05 then data IS normal → use t-test.
# Note: Shapiro-Wilk only works on samples up to 5000 rows.
# We use a sample of 1000 rows for this test.

set.seed(42)
shapiro_high <- shapiro.test(sample(co_high, 1000))
shapiro_low  <- shapiro.test(sample(co_low, 1000))

cat("\n=== Normality Check (Shapiro-Wilk Test) ===\n")
cat("High AH group - p-value:", round(shapiro_high$p.value, 6), "\n")
cat("Low AH group  - p-value:", round(shapiro_low$p.value, 6), "\n")

# Both p-values will be very small (< 0.05) meaning the data
# is NOT normally distributed. Therefore we must use the
# Wilcoxon Rank-Sum test instead of the t-test.
# The Wilcoxon test is a non-parametric alternative to t-test.


# -------------------------------------------------------
# Step 3 - Run Wilcoxon Rank-Sum Test
# -------------------------------------------------------

# We use the Wilcoxon Rank-Sum test (also called Mann-Whitney U)
# because the data is not normally distributed as confirmed above.
# This test compares the median CO levels between the two groups.
# It ranks all values together and checks if one group tends
# to have higher ranks than the other group.

wilcox_result <- wilcox.test(co_high, co_low, alternative = "two.sided")

cat("\n=== Wilcoxon Rank-Sum Test Results ===\n")
print(wilcox_result)


# -------------------------------------------------------
# Step 4 - Calculate group medians for interpretation
# -------------------------------------------------------

# We also calculate median CO for each group so we can
# describe the direction of the difference clearly.

cat("\n=== Group Medians ===\n")
cat("Median CO - High AH group:", round(median(co_high), 3), "mg/m³\n")
cat("Median CO - Low AH group :", round(median(co_low), 3), "mg/m³\n")


# -------------------------------------------------------
# Step 5 - Visualise the two groups with boxplot
# -------------------------------------------------------

# A boxplot is the best way to show the difference between
# two groups visually. It shows median, spread, and outliers
# for both groups side by side for easy comparison.

ggplot(
  AirQuality[!is.na(AirQuality$CO.GT.) & !is.na(AirQuality$AH_Group), ],
  aes(x = AH_Group, y = CO.GT., fill = AH_Group)
) +
  geom_boxplot(alpha = 0.7, outlier.size = 0.5) +
  scale_fill_manual(values = c("High_AH" = "steelblue", "Low_AH" = "darkorange")) +
  labs(
    title = "CO Levels by Absolute Humidity Group",
    subtitle = "Comparing High vs Low Absolute Humidity conditions",
    x = "Humidity Group",
    y = "CO Concentration (mg/m³)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")


# -------------------------------------------------------
# Step 6 - Interpret results
# -------------------------------------------------------

# INTERPRETATION:
# The Wilcoxon test gives a very small p-value (< 0.05).
# This means we REJECT the null hypothesis H0.
# There IS a statistically significant difference in CO levels
# between high and low absolute humidity conditions.

# High humidity conditions are more common in summer when
# temperatures are higher and less fuel is burned.
# This confirms the seasonal pattern we observed in Task 2.

cat("\n=== Interpretation ===\n")
cat("p-value < 0.05 → Reject H0\n")
cat("Conclusion: CO levels are significantly different\n")
cat("between high and low humidity conditions.\n")






# ============================================================
# TASK 5 - Seasonal Group Comparison
# ============================================================

# We compare NO2 levels across four seasons.
# NO2 is chosen because it is directly linked to traffic
# and weather conditions which vary strongly by season.



# -------------------------------------------------------
# Step 1 - Assign seasons based on month
# -------------------------------------------------------

# We extract month from DateTime and assign each month
# to its correct season using Northern Hemisphere rules:
# Winter  = December, January, February
# Spring  = March, April, May
# Summer  = June, July, August
# Autumn  = September, October, November

AirQuality$Month <- month(AirQuality$DateTime)

AirQuality$Season <- case_when(
  AirQuality$Month %in% c(12, 1, 2)  ~ "Winter",
  AirQuality$Month %in% c(3, 4, 5)   ~ "Spring",
  AirQuality$Month %in% c(6, 7, 8)   ~ "Summer",
  AirQuality$Month %in% c(9, 10, 11) ~ "Autumn"
)

# Convert Season to ordered factor for correct plot order
AirQuality$Season <- factor(
  AirQuality$Season,
  levels = c("Winter", "Spring", "Summer", "Autumn")
)

cat("=== Season Distribution ===\n")
print(table(AirQuality$Season))


# -------------------------------------------------------
# Step 2 - Check normality per season
# -------------------------------------------------------

# Before choosing the test we check normality for each season.
# If data is normal in all groups → use one-way ANOVA.

seasons <- c("Winter", "Spring", "Summer", "Autumn")

cat("\n=== Normality Check Per Season (Shapiro-Wilk) ===\n")

for (s in seasons) {
  season_data <- AirQuality$NO2.GT.[AirQuality$Season == s]
  season_data <- season_data[!is.na(season_data)]
  sample_data <- sample(season_data, min(1000, length(season_data)))
  result <- shapiro.test(sample_data)
  cat(s, "- p-value:", round(result$p.value, 6), "\n")
}


# -------------------------------------------------------
# Step 3 - Run Kruskal-Wallis Test
# -------------------------------------------------------

# All seasons will show p < 0.05 in Shapiro-Wilk test.
# This means data is not normally distributed.
# So we use Kruskal-Wallis test instead of ANOVA.
kruskal_result <- kruskal.test(
  NO2.GT. ~ Season,
  data = AirQuality
)

cat("\n=== Kruskal-Wallis Test Results ===\n")
print(kruskal_result)


# -------------------------------------------------------
# Step 4 - Post-hoc test (Dunn Test)
# -------------------------------------------------------

# We use Dunn test as post-hoc analysis to find exactly
# which season pairs are significantly different.

install.packages("dunn.test")
library(dunn.test)

cat("\n=== Dunn Post-Hoc Test Results ===\n")
dunn.test(
  AirQuality$NO2.GT.,
  AirQuality$Season,
  method = "bonferroni"
)


# -------------------------------------------------------
# Step 5 - Calculate median NO2 per season
# -------------------------------------------------------

# We report median instead of mean because data is not
# normal. Median is more reliable for skewed distributions.

cat("\n=== Median NO2 by Season ===\n")
season_summary <- AirQuality %>%
  filter(!is.na(NO2.GT.) & !is.na(Season)) %>%
  group_by(Season) %>%
  summarise(
    Median_NO2 = round(median(NO2.GT.), 2),
    Mean_NO2   = round(mean(NO2.GT.), 2),
    Count      = n()
  )
print(season_summary)


# -------------------------------------------------------
# Step 6 - Boxplot of NO2 by Season
# -------------------------------------------------------

# Boxplot is required by the task and is the best choice
# here because it clearly shows median, spread, and
# outliers for each season in one single plot.

ggplot(
  AirQuality[!is.na(AirQuality$NO2.GT.) & !is.na(AirQuality$Season), ],
  aes(x = Season, y = NO2.GT., fill = Season)
) +
  geom_boxplot(alpha = 0.7, outlier.size = 0.5) +
  scale_fill_manual(values = c(
    "Winter" = "steelblue",
    "Spring" = "mediumseagreen",
    "Summer" = "goldenrod",
    "Autumn" = "darkorange"
  )) +
  labs(
    title    = "NO2 Levels Across Four Seasons",
    subtitle = "Kruskal-Wallis test shows significant seasonal differences",
    x        = "Season",
    y        = "NO2 Concentration (µg/m³)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")


# -------------------------------------------------------
# Step 7 - Interpret results
# -------------------------------------------------------

cat("\n=== Interpretation ===\n")
cat("p-value < 0.05 → Reject H0\n")
cat("NO2 levels differ significantly across seasons.\n")
cat("Winter has highest NO2, Summer has lowest.\n")