# Air Quality Analysis for Port Talbot Hourly AQ Dataset
# This script imports, cleans, analyses and visualises air quality data collected between 2015 and 2024.

# Load required libraries
library(readr)
library(dplyr)
library(ggplot2)
library(corrplot)
library(car)

# Set the working directory
setwd(getwd())

# Read the dataset
AirQuality <- read_csv("Port_Talbot_Hourly_AQ.csv", show_col_types = FALSE)

# Replace "No data" values with NA
AirQuality[AirQuality == "No data"] <- NA

# Remove all status columns as they are not required for analysis
AirQuality <- AirQuality[, !grepl("^Status", names(AirQuality))]

# Convert the Date column into Date format
AirQuality$Date <- as.Date(
  AirQuality$Date,
  tryFormats = c("%d/%m/%Y", "%d-%m-%Y")
)

# Identify all numeric columns
num_cols <- setdiff(names(AirQuality), c("Date", "Time"))

# Convert character variables into numeric values
AirQuality[num_cols] <- lapply(
  AirQuality[num_cols],
  function(x) as.numeric(trimws(as.character(x)))
)

# Combine Date and Time into one DateTime variable
AirQuality$DateTime <- as.POSIXct(
  paste(AirQuality$Date, AirQuality$Time),
  format = "%Y-%m-%d %H:%M:%S"
)

# Extract month from each observation
m <- as.numeric(format(AirQuality$Date, "%m"))

# Create a Season variable for seasonal analysis
AirQuality$Season <- factor(
  ifelse(m %in% c(12,1,2),"Winter",
         ifelse(m %in% c(3,4,5),"Spring",
                ifelse(m %in% c(6,7,8),"Summer","Autumn"))),
  levels = c("Winter","Spring","Summer","Autumn")
)

# Write output to both the console and a text file
sink("Analysis_Output.txt", split = TRUE)

cat("Air Quality Analysis\n\n")

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

# Display data structure
cat("\nData types\n")
str(AirQuality)

# Display first six observations
cat("\nFirst six rows\n")
print(head(AirQuality))

# Display dataset dimensions
cat("\nDataset dimensions\n")
print(dim(AirQuality))

# Produce descriptive statistics
cat("\nDescriptive statistics\n")
print(summary(AirQuality))


# Time series plots show long-term trends in pollutant levels
# Plot ozone concentration over time
p1 <- ggplot(AirQuality, aes(DateTime, Ozone)) +
  geom_line(colour = "blue") +
  labs(title = "Ozone over Time",
       x = "Date",
       y = "Ozone") +
  theme_minimal()

ggsave("TimeSeries_Ozone.png", p1, width = 10, height = 5)

# Plot PM2.5 concentration over time
p2 <- ggplot(
  AirQuality,
  aes(DateTime, `PM2.5 particulate matter (Hourly measured)`)
) +
  geom_line(colour = "red") +
  labs(title = "PM2.5 over Time",
       x = "Date",
       y = "PM2.5") +
  theme_minimal()

ggsave("TimeSeries_PM25.png", p2, width = 10, height = 5)


# Scatter plots examine relationships between variables
# Relationship between ozone and temperature
p3 <- ggplot(
  AirQuality,
  aes(`Modelled Temperature`, Ozone)
) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", se = FALSE, colour = "red") +
  theme_minimal() +
  labs(title = "Ozone vs Temperature",
       x = "Modelled Temperature",
       y = "Ozone")

ggsave("Scatter_Ozone_Temperature.png", p3, width = 7, height = 6)

# Relationship between PM2.5 and wind speed
p4 <- ggplot(
  AirQuality,
  aes(`Modelled Wind Speed`,
      `PM2.5 particulate matter (Hourly measured)`)
) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", se = FALSE, colour = "blue") +
  theme_minimal() +
  labs(title = "PM2.5 vs Wind Speed",
       x = "Modelled Wind Speed",
       y = "PM2.5")

ggsave("Scatter_PM25_WindSpeed.png", p4, width = 7, height = 6)


# Correlation analysis
# Select pollutant variables for correlation analysis
corr_data <- AirQuality %>%
  select(
    Ozone,
    `Nitrogen dioxide`,
    `PM2.5 particulate matter (Hourly measured)`,
    `PM10 particulate matter (Hourly measured)`
  )

# Calculate correlation matrix
corr_matrix <- cor(
  corr_data,
  use = "pairwise.complete.obs"
)

cat("\nCorrelation matrix\n")
print(round(corr_matrix, 3))

# Produce correlation heatmap
png("Correlation_Heatmap.png", 900, 700)

corrplot(
  corr_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  tl.col = "black"
)

dev.off()

# Produce scatterplot matrix
png("Scatterplot_Matrix.png", 900, 900)

pairs(corr_data, main = "Scatterplot Matrix")

dev.off()

# Provide interpretation of correlation values
cat("\nCorrelation interpretation\n")
cat("PM2.5 and PM10 showed a moderate positive correlation (r = 0.605).\n")
cat("Ozone and nitrogen dioxide showed a moderate negative correlation (r = -0.486).\n")
cat("Nitrogen dioxide and PM2.5 showed a weak positive correlation (r = 0.315).\n")
cat("Ozone and PM10 showed a very weak positive correlation (r = 0.100).\n")


# Independent samples t-test
# Calculate the median wind speed
median_ws <- median(
  AirQuality$`Modelled Wind Speed`,
  na.rm = TRUE
)

# Divide observations into high and low wind speed groups
AirQuality$WindGroup <- ifelse(
  AirQuality$`Modelled Wind Speed` >= median_ws,
  "High",
  "Low"
)

cat("\nHypothesis test\n")
cat("Null hypothesis: Mean ozone concentration is equal under high and low wind speed.\n")
cat("Alternative hypothesis: Mean ozone concentration differs under high and low wind speed.\n")

# Perform independent samples t-test
tt <- t.test(Ozone ~ WindGroup, data = AirQuality)

print(tt)

# Interpret the hypothesis test
if (tt$p.value < 0.05) {
  cat("Interpretation: Reject the null hypothesis. Ozone concentrations differ significantly between wind speed groups.\n")
} else {
  cat("Interpretation: Fail to reject the null hypothesis.\n")
}


# One-way ANOVA
cat("\nSeasonal comparison\n")
cat("Null hypothesis: Mean ozone concentration is equal across all seasons.\n")
cat("Alternative hypothesis: At least one season has a different mean ozone concentration.\n")

# Perform one-way ANOVA
fit <- aov(Ozone ~ Season, data = AirQuality)

print(summary(fit))

# Perform Tukey post hoc comparisons
cat("\nTukey HSD post hoc test\n")
print(TukeyHSD(fit))

# Produce seasonal boxplot
p5 <- ggplot(
  AirQuality,
  aes(Season, Ozone, fill = Season)
) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Ozone by Season",
       x = "Season",
       y = "Ozone")

ggsave("Boxplot_Ozone_Season.png", p5, width = 8, height = 6)

# Interpret ANOVA results
if (summary(fit)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("Interpretation: Reject the null hypothesis. Ozone concentrations vary significantly across seasons.\n")
} else {
  cat("Interpretation: Fail to reject the null hypothesis.\n")
}
sink()