setwd("D:/Users/Jervis Payne/Downloads/Acca7022 - file")
list.files(recursive = TRUE)

AirQuality <- read.csv("AirQualityUCI.csv", sep=";", dec=",", header=TRUE, stringsAsFactors=FALSE)
names(AirQuality)
str(AirQuality)
head(AirQuality, 5)
dim(AirQuality)
summary(AirQuality)

############################################################
# ACCA7022 Component 2 (March 2026)
# Script: AirQualityAnalysis.r
# Dataset: AirQualityUCI (UCI ML Repository)
# Student: 2127853
############################################################

# -------------------------------
# 1) Set working directory + check files
# -------------------------------
# Setting the working directory ensures the script can locate the dataset
# reliably using relative filenames, improving reproducibility when running the script.
setwd("D:/Users/Jervis Payne/Downloads/Acca7022 - file")

# Listing files provides quick evidence that the required dataset file is
# present in the working directory before import (helps troubleshooting).
# Optional: comment out once confirmed.
print(list.files(recursive = TRUE))

# -------------------------------
# 2) Import dataset into AirQuality
# -------------------------------
# Importing the dataset into a single dataframe called AirQuality is required
# and provides the base object used for all subsequent cleaning, visualisation, and testing.
# The separator (;) and decimal comma (,) match the formatting used by AirQualityUCI.csv.
AirQuality <- read.csv(
  "AirQualityUCI.csv",
  sep = ";",
  dec = ",",
  header = TRUE,
  stringsAsFactors = FALSE
)

# -------------------------------
# 3) Check column names + data types
# -------------------------------
# Checking column names and data types confirms the import was successful,
# identifies which variables are pollutants vs meteorological measures, and prevents
# incorrect statistical analysis due to wrong data types (e.g., characters vs numeric).
print(names(AirQuality))
str(AirQuality)

# -------------------------------
# 4) head, dimensions, descriptive stats
# -------------------------------
# Viewing the first rows verifies the dataset structure and typical values.
# Reporting dimensions confirms the number of observations and variables.
# Descriptive statistics provide baseline understanding and help detect missing/invalid codes.
print(head(AirQuality, 5))
print(dim(AirQuality))
print(summary(AirQuality))

# -------------------------------
# 5) Data cleaning + DateTime creation (FOUNDATIONAL FOR VALID ANALYSIS)
# -------------------------------
# The dataset contains non-informative columns (X, X.1) that are entirely NA and should be removed.
# It also contains invalid/missing sensor readings coded as -200. Converting -200 to NA prevents distorted summary
# statistics and ensures correlations/tests are performed only on valid measurements. Creating a DateTime variable
# enables correct time-series plots for hourly data.

# Remove columns that are entirely NA in this file
AirQuality <- AirQuality[, !names(AirQuality) %in% c("X", "X.1")]

# Convert -200 invalid codes to NA for measurement columns (do not apply to Date/Time)
numeric_cols <- setdiff(names(AirQuality), c("Date", "Time"))
for (col_name in numeric_cols) {
  AirQuality[[col_name]][AirQuality[[col_name]] == -200] <- NA
}

# Create a DateTime column (required for time trend visualisations)
# Justification: Combining Date and Time into a POSIXct DateTime supports correct chronological plotting.
AirQuality$DateTime <- as.POSIXct(
  paste(AirQuality$Date, AirQuality$Time),
  format = "%d/%m/%Y %H.%M.%S",
  tz = "UTC"
)

# -------------------------------
# PLAN 2: Data cleaning + DateTime creation (RUN BEFORE VALIDATION)
# -------------------------------
# X and X.1 are empty columns (all NA) and must be removed.
# Values of -200 represent invalid/missing readings in this file and must be treated as NA
# to avoid distorted summaries, plots, correlations, and hypothesis tests.

# Remove empty columns
AirQuality$X <- NULL
AirQuality$X.1 <- NULL

# Convert -200 invalid codes to NA for all measurement variables
# (Leave Date and Time as character for now)
measurement_cols <- setdiff(names(AirQuality), c("Date", "Time"))
for (col_name in measurement_cols) {
  AirQuality[[col_name]][AirQuality[[col_name]] == -200] <- NA
}

# Create DateTime (hourly time-series index)
# Justification: Required for time trend plots and seasonal grouping later.
AirQuality$DateTime <- as.POSIXct(
  paste(AirQuality$Date, AirQuality$Time),
  format = "%d/%m/%Y %H.%M.%S",
  tz = "UTC"
)
print(c("X" %in% names(AirQuality), "X.1" %in% names(AirQuality)))  # should be FALSE FALSE
print(min(AirQuality$CO.GT., na.rm = TRUE))                         # should NOT be -200



# ===============================
# PLAN 3: Exploratory Visualisation
# ===============================
# Examine time trends in at least two pollutants and explore relationships
# between pollutants and weather variables (T, RH, AH). Plot choices must be justified in code comments.

# -------------------------------
# 3A) Time trends (at least two pollutants)
# -------------------------------

# A line plot is appropriate for hourly time-series data because it shows
# temporal patterns (trends, seasonality, peaks/spikes) in pollutant concentrations over time.
plot(AirQuality$DateTime, AirQuality$CO.GT., type = "l",
     xlab = "DateTime", ylab = "CO (mg/m^3)", main = "Time Trend: CO (CO.GT.)")

# Using a second pollutant time-series supports comparison of how different pollutants
# behave over time (e.g., co-movement, different seasonal profiles), strengthening EDA conclusions.
plot(AirQuality$DateTime, AirQuality$NO2.GT., type = "l",
     xlab = "DateTime", ylab = "NO2 (µg/m^3)", main = "Time Trend: NO2 (NO2.GT.)")


# -------------------------------
# 3B) Pollutants vs weather variables (relationships)
# -------------------------------

#  A scatter plot is appropriate for exploring relationships between two continuous variables
# (pollutant vs weather) and for revealing non-linear patterns and outliers that may affect later test choice.
plot(AirQuality$T, AirQuality$CO.GT.,
     xlab = "Temperature (°C)", ylab = "CO (mg/m^3)", main = "CO vs Temperature")

# LOWESS summarises the overall relationship, but requires complete (non-missing) x/y pairs.
cc1 <- complete.cases(AirQuality$T, AirQuality$CO.GT.)
lines(lowess(AirQuality$T[cc1], AirQuality$CO.GT.[cc1]), lwd = 2)

# Relative humidity can influence pollutant dispersion/chemistry; a scatter plot tests whether
# higher humidity is associated with different pollutant levels.
plot(AirQuality$RH, AirQuality$NO2.GT.,
     xlab = "Relative Humidity (%)", ylab = "NO2 (µg/m^3)", main = "NO2 vs Relative Humidity")

# LOWESS provides a robust trend summary while excluding missing observations.
cc2 <- complete.cases(AirQuality$RH, AirQuality$NO2.GT.)
lines(lowess(AirQuality$RH[cc2], AirQuality$NO2.GT.[cc2]), lwd = 2)

# Absolute humidity is another weather-related variable; using it demonstrates broader exploration
# of meteorological relationships rather than relying on a single factor.
plot(AirQuality$AH, AirQuality$CO.GT.,
     xlab = "Absolute Humidity (AH)", ylab = "CO (mg/m^3)", main = "CO vs Absolute Humidity")

# LOWESS trend line summarises the overall association using complete cases only.
cc3 <- complete.cases(AirQuality$AH, AirQuality$CO.GT.)
lines(lowess(AirQuality$AH[cc3], AirQuality$CO.GT.[cc3]), lwd = 2)


# ===============================
# PLAN 4: Correlations between pollutants
# ===============================
# Calculate correlations between at least three pollutants and present a correlation visual
# (e.g., scatterplots, correlogram, heatmap) with an interpretation in context.

# Selected at least three pollutant variables (using the ground-truth pollutant columns in this dataset)
pollutant_cols <- c("CO.GT.", "NOx.GT.", "NO2.GT.", "C6H6.GT.")  # 4 pollutants (≥3 required)
pollutants <- AirQuality[, pollutant_cols]

# Pairwise complete observations avoid discarding whole rows when some variables are missing (NA).
corr_mat <- cor(pollutants, use = "pairwise.complete.obs", method = "pearson")
print(corr_mat)

# -------------------------------
# Correlation visualisation (heatmap)
# -------------------------------
# A heatmap summarises correlation strength and direction across multiple pollutant pairs
# in one figure, making patterns easy to identify and compare.
heatmap(corr_mat, symm = TRUE, main = "Correlation Heatmap (Pollutants)")

# -------------------------------
# Interpretation
# -------------------------------
# Interpretation: Replace the example text below with your actual observations from corr_mat.
# Example interpretation template:
# - Identify the strongest positive correlation (which pollutant pair and approximate r value).
# - Identify any weak/near-zero correlation.
# - Provide a short air-quality explanation (e.g., shared traffic emissions, different chemical behaviour).
#
# Example: "CO and NOx show a positive correlation, consistent with shared combustion/traffic sources.
# NO2 may correlate differently due to atmospheric conversion processes and local dispersion conditions."


# Interpretation (based on corr_mat):
# CO and C6H6 show the strongest positive correlation (r ≈ 0.93), suggesting these pollutants often increase
# and decrease together, consistent with shared combustion/traffic-related emission behaviour in urban air.
# CO also correlates strongly with NOx (r ≈ 0.80) and moderately with NO2 (r ≈ 0.68), indicating that periods
# of elevated combustion emissions tend to coincide with higher nitrogen oxides.
# NOx and NO2 are strongly correlated (r ≈ 0.76), which is expected because NO2 is chemically related to NOx,
# but the correlation is not perfect, likely reflecting atmospheric conversion (NO ↔ NO2), dispersion, and
# time-of-day effects. The weakest relationship here is NO2 with C6H6 (r ≈ 0.61), implying that despite shared
# emission influences, NO2 variability is also driven by chemistry and meteorological conditions.


# A scatterplot matrix complements the heatmap by showing the raw bivariate relationships,
# potential non-linearity, and outliers that may not be visible from correlation coefficients alone.
pairs(pollutants, main = "Scatterplot Matrix (Pollutants)")


print(c("RH" %in% names(AirQuality), "NO2.GT." %in% names(AirQuality)))

rh_median <- median(AirQuality$RH, na.rm = TRUE)
AirQuality$RH_Group <- ifelse(AirQuality$RH >= rh_median, "High_RH", "Low_RH")
table(AirQuality$RH_Group, useNA = "ifany")



# ===============================
# PLAN 5: Two-condition hypothesis test
# ===============================
# Choose one pollutant and compare it under two conditions (e.g., high vs low wind speed).
# This dataset does not contain wind variables, so i use a defensible meteorological split using RH (humidity):
# High_RH vs Low_RH based on the median. This creates two conditions for hypothesis testing.

# -------------------------------
# 5A) Define hypotheses (explicit)
# -------------------------------
# H0 (null): The central tendency (median) of NO2.GT. is the same for High_RH and Low_RH.
# H1 (alternative): The central tendency (median) of NO2.GT. differs between High_RH and Low_RH.

# -------------------------------
# 5B) Choose test 
# -------------------------------
# Pollutant data are often skewed and include outliers. The Wilcoxon rank-sum test
# (Mann–Whitney U) is a robust non-parametric alternative to the two-sample t-test because it does not
# assume normality. It tests whether the distributions (and typically medians) differ between groups.

# Run Wilcoxon test (automatically excludes NA values)
wilcox_NO2_RH <- wilcox.test(NO2.GT. ~ RH_Group, data = AirQuality, exact = FALSE)
print(wilcox_NO2_RH)

# -------------------------------
# 5C) Report group summaries 
# -------------------------------
# Reporting medians/IQRs supports interpretation by showing the direction and practical size
# of the difference, not only the p-value.
NO2_high <- AirQuality$NO2.GT.[AirQuality$RH_Group == "High_RH"]
NO2_low  <- AirQuality$NO2.GT.[AirQuality$RH_Group == "Low_RH"]

print(c(
  Median_High_RH = median(NO2_high, na.rm = TRUE),
  Median_Low_RH  = median(NO2_low,  na.rm = TRUE)
))

print(c(
  IQR_High_RH = IQR(NO2_high, na.rm = TRUE),
  IQR_Low_RH  = IQR(NO2_low,  na.rm = TRUE)
))

print(c(Median_Difference_LowMinusHigh =
          median(NO2_low, na.rm=TRUE) - median(NO2_high, na.rm=TRUE)))


# Interpretation (write in comments after seeing output):
# - If p-value < 0.05: reject H0; conclude NO2 differs between High_RH and Low_RH.
# - Use the medians to state which condition has higher NO2 and give a plausible explanation.


# Interpretation:
# The Wilcoxon rank-sum test shows a statistically significant difference in NO2 between humidity groups
# (p-value < 2.2e-16), so H0 is rejected at the 5% significance level.
# The median NO2 is lower under High_RH (median = 102) than under Low_RH (median = 114), indicating that
# typical NO2 concentrations are higher during lower-humidity conditions in this dataset.
# The IQRs are similar (High_RH = 64, Low_RH = 61), suggesting comparable variability, but with a shifted
# central tendency. A plausible explanation is that humidity is linked to atmospheric mixing and chemistry,
# so changes in RH can coincide with conditions that alter measured NO2 levels (e.g., dispersion and conversion).




# ===============================
# PLAN 6: Seasonal comparison
# ===============================
# Compared one pollutant across seasons (Winter/Spring/Summer/Autumn),
# produce boxplots, run an appropriate test (ANOVA or Kruskal–Wallis), and interpret results.

# -------------------------------
# 6A) Create Season variable from Date
# -------------------------------
# Seasons are derived from month to allow comparison of pollutant behaviour across
# Winter/Spring/Summer/Autumn (expected to differ due to weather, emission patterns, and chemistry).
AirQuality$DateParsed <- as.Date(AirQuality$Date, format = "%d/%m/%Y")
AirQuality$Month <- as.integer(format(AirQuality$DateParsed, "%m"))

AirQuality$Season <- NA
AirQuality$Season[AirQuality$Month %in% c(12, 1, 2)]  <- "Winter"
AirQuality$Season[AirQuality$Month %in% c(3, 4, 5)]   <- "Spring"
AirQuality$Season[AirQuality$Month %in% c(6, 7, 8)]   <- "Summer"
AirQuality$Season[AirQuality$Month %in% c(9, 10, 11)] <- "Autumn"
AirQuality$Season <- factor(AirQuality$Season, levels = c("Winter", "Spring", "Summer", "Autumn"))

# Evidence check: confirm season counts
print(table(AirQuality$Season, useNA = "ifany"))

# -------------------------------
# 6B) Choose pollutant for seasonal comparison
# -------------------------------
# NO2 is a key pollutant with strong links to traffic/combustion and can vary seasonally
# due to atmospheric mixing and photochemical effects.
season_pollutant <- AirQuality$NO2.GT.

# -------------------------------
# 6C) Boxplot by season (visual comparison)
# -------------------------------
# Boxplots are appropriate for comparing distributions across multiple groups (four seasons),
# showing medians, spread (IQR), and outliers in a compact form.
boxplot(season_pollutant ~ AirQuality$Season,
        xlab = "Season", ylab = "NO2 (µg/m^3)",
        main = "Seasonal Comparison of NO2 (NO2.GT.)")

# -------------------------------
# 6D) Hypotheses
# -------------------------------
# H0 (null): The central tendency of NO2 is the same across Winter, Spring, Summer, and Autumn.
# H1 (alternative): At least one season differs in NO2 central tendency.

# -------------------------------
# 6E) Choosen test
# -------------------------------
# Pollutant measurements commonly show skewness/outliers and may violate ANOVA assumptions.
# The Kruskal–Wallis test is a robust non-parametric method for comparing more than two groups without
# assuming normality. It tests whether at least one group distribution differs.
kruskal_season_NO2 <- kruskal.test(season_pollutant ~ AirQuality$Season)
print(kruskal_season_NO2)

# -------------------------------
# 6F) Seasonal summaries
# -------------------------------
# Reporting medians/IQRs per season supports interpretation by showing direction and practical size.
season_levels <- levels(AirQuality$Season)

for (s in season_levels) {
  vals <- season_pollutant[AirQuality$Season == s]
  cat("\nSeason:", s,
      "| n =", sum(!is.na(vals)),
      "| median =", median(vals, na.rm = TRUE),
      "| IQR =", IQR(vals, na.rm = TRUE))
}
cat("\n")

# -------------------------------
# 6G) Interpretation
# -------------------------------
# Interpretation:
# - If p-value < 0.05: reject H0; conclude NO2 differs by season (at least one season differs).
# - Use the boxplot + seasonal medians to state which season(s) have higher/lower typical NO2.
# - Provide a plausible explanation (e.g., winter inversion/heating, summer photochemistry/dispersion).


# Interpretation:
# The Kruskal–Wallis test indicates a statistically significant seasonal difference in NO2
# (chi-squared = 955.64, df = 3, p-value < 2.2e-16), so H0 is rejected at the 5% level.
# Seasonal medians show the highest typical NO2 in Winter (median = 135), followed by Spring (106.5)
# and Autumn (105), with the lowest median in Summer (88.5). The IQRs indicate Winter and Autumn also
# show wider variability (IQR = 66 and 65) than Spring (52), suggesting more spread in those seasons.
# A plausible explanation is that winter conditions can reduce atmospheric mixing (stable air/inversions)
# and increase combustion-related emissions (e.g., heating and traffic), raising NO2, while summer has
# stronger dispersion and photochemical processes that can reduce observed NO2 levels.


# Practical differences (median shifts):
print(c(
  Winter_minus_Summer = 135 - 88.5,
  Winter_minus_Spring = 135 - 106.5,
  Winter_minus_Autumn = 135 - 105
))