# ==========================================================
# Air Quality Analysis Using R
# ==========================================================


# ==========================================================
# Import Packages
# ==========================================================

library(tidyverse)
library(lubridate)

library(corrplot)


# ==========================================================
# Import Dataset
# ==========================================================

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


# ==========================================================
# Initial Dataset Inspection
# ==========================================================

colnames(AirQuality)

str(AirQuality)


# ==========================================================
# Data Cleaning
# ==========================================================

# Remove status columns because they contain metadata
# rather than numerical measurements.

AirQuality <- AirQuality %>%
  select(-starts_with("Status"))


# Replace "No data" values with NA

AirQuality[AirQuality == "No data"] <- NA


# Convert Date format

AirQuality$Date <- dmy(AirQuality$Date)


# Convert Time format

AirQuality$Time <- format(
  strptime(
    AirQuality$Time,
    format="%H:%M:%S"
  ),
  "%H:%M:%S"
)


# Convert numerical variables

numeric_columns <- c(
  "Ozone",
  "Nitric oxide",
  "Nitrogen dioxide",
  "Nitrogen oxides as nitrogen dioxide",
  "Sulphur dioxide",
  "Carbon monoxide",
  "PM10 particulate matter (Hourly measured)",
  "PM2.5 particulate matter (Hourly measured)",
  "Modelled Wind Direction",
  "Modelled Wind Speed",
  "Modelled Temperature"
)


AirQuality[numeric_columns] <-
  lapply(
    AirQuality[numeric_columns],
    as.numeric
  )


# Create complete date-time variable
# Required because measurements are hourly.

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


AirQuality <- AirQuality %>%
  arrange(DateTime)



# ==========================================================
# TASK 1: Initial Data Exploration
# ==========================================================


# First six observations
# This confirms successful data import and cleaning.

head(AirQuality,6)



# Dataset dimensions

cat(
  "Cleaned Dataset Dimensions:",
  nrow(AirQuality),
  "rows and",
  ncol(AirQuality),
  "columns"
)



# Descriptive statistics

analysis_variables <- AirQuality %>%
  select(
    Ozone,
    `Nitrogen dioxide`,
    `PM10 particulate matter (Hourly measured)`,
    `PM2.5 particulate matter (Hourly measured)`,
    `Modelled Wind Speed`,
    `Modelled Temperature`
  )



descriptive_statistics <- data.frame(
  
  Variable = names(analysis_variables),
  
  Valid_Observations =
    sapply(
      analysis_variables,
      function(x) sum(!is.na(x))
    ),
  
  Missing_Values =
    sapply(
      analysis_variables,
      function(x) sum(is.na(x))
    ),
  
  Mean =
    sapply(
      analysis_variables,
      function(x) mean(x,na.rm=TRUE)
    ),
  
  Standard_Deviation =
    sapply(
      analysis_variables,
      function(x) sd(x,na.rm=TRUE)
    ),
  
  Median =
    sapply(
      analysis_variables,
      function(x) median(x,na.rm=TRUE)
    ),
  
  Minimum =
    sapply(
      analysis_variables,
      function(x) min(x,na.rm=TRUE)
    ),
  
  Maximum =
    sapply(
      analysis_variables,
      function(x) max(x,na.rm=TRUE)
    )
  
)


descriptive_statistics[,4:8] <-
  round(
    descriptive_statistics[,4:8],
    2
  )


print(descriptive_statistics)



# Missing values

colSums(is.na(AirQuality))



# ==========================================================
# TASK 2: Exploratory Visualisation
# ==========================================================


# A time-series line plot was selected because PM2.5 measurements
# were collected hourly across a ten-year period. A line plot is
# suitable for ordered temporal data because it connects observations
# chronologically and allows identification of long-term trends,
# seasonal changes, pollution peaks and unusual events.
# A bar chart was not selected because the dataset contains thousands
# of continuous observations where individual bars would reduce clarity.

ggplot(
  AirQuality,
  aes(
    x=DateTime,
    y=`PM2.5 particulate matter (Hourly measured)`
  )
)+
  geom_line(
    colour="steelblue"
  )+
  labs(
    title="PM2.5 Concentration Over Time",
    x="Date",
    y="PM2.5 (µg/m³)"
  )+
  theme_minimal()



# A second time-series plot was created for nitrogen dioxide to compare
# another pollutant's behaviour over the same time period.
# Using the same visual format allows comparison between pollutants
# and helps identify whether changes occur simultaneously or follow
# different temporal patterns.


ggplot(
  AirQuality,
  aes(
    x=DateTime,
    y=`Nitrogen dioxide`
  )
)+
  geom_line(
    colour="darkgreen"
  )+
  labs(
    title="Nitrogen Dioxide Concentration Over Time",
    x="Date",
    y="Nitrogen Dioxide (µg/m³)"
  )+
  theme_minimal()



# A scatterplot was selected because PM2.5 concentration and wind speed
# are both continuous numerical variables.
# Scatterplots are appropriate for investigating possible relationships
# between two continuous variables because they show direction,
# strength, variability and potential outliers.
# The fitted linear regression line provides an indication of the
# overall relationship pattern.


ggplot(
  AirQuality,
  aes(
    x=`Modelled Wind Speed`,
    y=`PM2.5 particulate matter (Hourly measured)`
  )
)+
  geom_point(
    alpha=0.4
  )+
  geom_smooth(
    method="lm",
    se=FALSE
  )+
  labs(
    title="Relationship Between PM2.5 and Wind Speed",
    x="Wind Speed (m/s)",
    y="PM2.5 (µg/m³)"
  )+
  theme_minimal()



# A scatterplot was selected to examine the relationship between
# nitrogen dioxide concentration and temperature because both variables
# are continuous measurements.
# This plot type allows assessment of whether temperature changes are
# associated with increases or decreases in nitrogen dioxide levels,
# while also showing variability and extreme observations.


ggplot(
  AirQuality,
  aes(
    x=`Modelled Temperature`,
    y=`Nitrogen dioxide`
  )
)+
  geom_point(
    alpha=0.4
  )+
  geom_smooth(
    method="lm",
    se=FALSE
  )+
  labs(
    title="Relationship Between Nitrogen Dioxide and Temperature",
    x="Temperature (°C)",
    y="Nitrogen Dioxide (µg/m³)"
  )+
  theme_minimal()



# ==========================================================
# TASK 3: Correlation Analysis
# ==========================================================


pollutants <- AirQuality %>%
  select(
    Ozone,
    `Nitrogen dioxide`,
    `PM10 particulate matter (Hourly measured)`,
    `PM2.5 particulate matter (Hourly measured)`
  )


# Calculate Pearson correlation matrix.
#
# Pearson correlation was selected because the variables are
# continuous numerical pollutant measurements.
# The correlation coefficient ranges from -1 to +1:
# positive values indicate that variables increase together,
# negative values indicate an inverse relationship,
# and values close to zero indicate weak linear association.


correlation_matrix <- cor(
  pollutants,
  use="complete.obs",
  method="pearson"
)



print(correlation_matrix)



corrplot(
  correlation_matrix,
  method="color",
  type="upper",
  addCoef.col="black",
  tl.col="black",
  tl.srt=45
)


# Interpretation of correlation results:
#
# PM10 and PM2.5 show the strongest positive relationship.
# This was expected because both represent particulate matter
# and may increase during similar pollution episodes.
#
# Ozone and nitrogen dioxide show a negative relationship.
# This may indicate different atmospheric behaviour because
# these pollutants are influenced by different chemical processes.
#
# Nitrogen dioxide and PM2.5 show a weaker positive relationship,
# suggesting that some emission sources may affect both pollutants,
# but they are also influenced by other factors such as weather
# conditions and pollutant dispersion.
#
# Correlation values represent association only and cannot prove
# that one pollutant causes changes in another.
# Additional analysis considering season, wind direction,
# temperature and emission sources would provide deeper insight.



# ==========================================================
# TASK 4: Hypothesis Testing
# ==========================================================


# Hypothesis testing:
# The research question is whether PM2.5 concentrations differ
# between low and high wind speed conditions.
#
# Welch Two-Sample t-test was selected because:
# 1. PM2.5 is a continuous numerical variable.
# 2. Wind speed conditions create two independent groups.
# 3. Welch's test compares the means of two groups without assuming
#    equal variance between groups.
#
# This method is appropriate because environmental datasets often
# contain unequal variability between different weather conditions.


median_wind <- median(
  AirQuality$`Modelled Wind Speed`,
  na.rm=TRUE
)



AirQuality$Wind_Group <- ifelse(
  AirQuality$`Modelled Wind Speed` <= median_wind,
  "Low Wind Speed",
  "High Wind Speed"
)



# Welch two sample t-test selected because:
# - PM2.5 is continuous
# - Wind group contains two independent categories
# - Welch test does not assume equal variance between groups.


# Null hypothesis (H0):
# There is no significant difference in PM2.5 concentrations
# between low and high wind speed conditions.

# Alternative hypothesis (H1):
# There is a significant difference in PM2.5 concentrations
# between low and high wind speed conditions.



t_test_result <- t.test(
  `PM2.5 particulate matter (Hourly measured)` ~ Wind_Group,
  data=AirQuality
)


t_test_result



# ==========================================================
# Effect Size Calculation
# ==========================================================


high_pm <- AirQuality %>%
  filter(Wind_Group=="High Wind Speed") %>%
  pull(`PM2.5 particulate matter (Hourly measured)`)


low_pm <- AirQuality %>%
  filter(Wind_Group=="Low Wind Speed") %>%
  pull(`PM2.5 particulate matter (Hourly measured)`)


high_pm <- high_pm[!is.na(high_pm)]
low_pm <- low_pm[!is.na(low_pm)]


pooled_sd <- sqrt(
  (
    (length(high_pm)-1)*var(high_pm)+
      (length(low_pm)-1)*var(low_pm)
  )/
    (length(high_pm)+length(low_pm)-2)
)


cohens_d <-
  (mean(high_pm)-mean(low_pm))/pooled_sd


cat(
  "Cohen's d:",
  round(cohens_d,3)
)



# ==========================================================
# TASK 5: Seasonal Analysis
# ==========================================================


AirQuality$Month <-
  month(AirQuality$Date)



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",
  
  TRUE
  ~ "Autumn"
  
)



AirQuality$Season <-
  factor(
    AirQuality$Season,
    levels=c(
      "Winter",
      "Spring",
      "Summer",
      "Autumn"
    )
  )


# A boxplot was selected because season is a categorical variable
# containing four groups, while PM2.5 concentration is a continuous
# numerical variable.
# Boxplots are suitable for comparing distributions because they show
# median values, interquartile ranges, variability and potential
# outliers between seasons.
# This visualisation also supports the later ANOVA analysis by showing
# whether seasonal differences appear visually meaningful.


ggplot(
  AirQuality,
  aes(
    x=Season,
    y=`PM2.5 particulate matter (Hourly measured)`,
    fill=Season
  )
)+
  geom_boxplot()+
  labs(
    title="PM2.5 Concentrations Across Seasons",
    x="Season",
    y="PM2.5 (µg/m³)"
  )+
  theme_minimal()



# ANOVA test

anova_result <- aov(
  `PM2.5 particulate matter (Hourly measured)` ~ Season,
  data=AirQuality
)


summary(anova_result)



# Post hoc comparison

TukeyHSD(anova_result)