---
title: "AirpollutionDataAnalycis"
author: "Tasin Islam"
date: "`r Sys.Date()`"
output: word_document
---

---
title: "Port Talbot Hourly Air Quality Analysis Report"
author: "Analysis of Air Quality Data"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_depth: 3
    theme: cosmo
---

```{r setup, include=FALSE}
# This chunk sets global options for code chunks
knitr::opts_chunk$set(
  echo = TRUE,      # Show the code in the output
  warning = FALSE,  # Hide warnings (like messages from dplyr)
  message = FALSE,  # Hide messages
  fig.align = "center", # Center all plots
  fig.width = 8,
  fig.height = 5
)

# Loading all libraries used in the original script
library(rstatix)
library(readr)
library(ggplot2)
library(tidyverse)
library(lubridate)
library(car)
library(broom)
library(knitr) # For kable
Data Preparation and Cleaning (Problem 1)This section reads the raw air quality data, renames the columns for easier use, and creates a combined datetime column.Loading Data and Renaming ColumnsCode snippet# Read CSV File
PortTalbotAQ <- read_csv("Port_Talbot_Hourly_AQ.csv",
                         na = c("No data", "NA", ""))

# Rename all Headers in a short form
PortTalbotAQ <- PortTalbotAQ %>%
  rename(
    Date = `Date`,
    Time = `Time`,
    O3 = `Ozone`,
    Status = `Status`,
    NO = `Nitric oxide`,
    Status1 = `Status.1`,
    NO2 = `Nitrogen dioxide`,
    Status2 = `Status.2`,
    NOx = `Nitrogen oxides as nitrogen dioxide`,
    Status3 = `Status.3`,
    SO2 = `Sulphur dioxide`,
    Status4 = `Status.4`,
    CO = `Carbon monoxide`,
    Status5 = `Status.5`,
    PM10 = `PM10 particulate matter (Hourly measured)`,
    Status6 = `Status.6`,
    PM2.5 = `PM2.5 particulate matter (Hourly measured)`,
    Status7 = `Status.7`,
    WindDirct = `Modelled Wind Direction`,
    Status8 = `Status.8`,
    WindSpd = `Modelled Wind Speed`,
    Status9 = `Status.9`,
    Tempt = `Modelled Temperature`,
    Status10 = `Status.10`
  )

# Create Datetime and Arrange
PortTalbotAQ <- PortTalbotAQ %>%
  mutate(
    datetime = dmy_hms(
      paste(Date, Time),
      tz = "UTC")
      ) %>%
  arrange(datetime)
First 6 Rows of the Processed DataframeThe head of the dataframe shows the new, cleaner structure with the combined datetime column.Code snippet# Result in DataFrame
PortTalbotAQ %>%
  head() %>%
  kable(caption = "Head of PortTalbotAQ Dataframe")
Missing Data Analysis (Problem 2)This section calculates the percentage of missing readings for the specified pollutants.Code snippet# Percentage of missing readings per pollutant
pollutants <- c("O3", "NO", "NO2", "NOx", "SO2", "CO")
invalidData <- tibble(
  pollutant = pollutants,
  pct_missing = sapply(pollutants, function(p) {
    mean(is.na(PortTalbotAQ[[p]])) * 100
  })
)
Percentage of Missing Readings TableCode snippetinvalidData %>%
  kable(
    caption = "Percentage of Missing Readings Per Pollutant",
    digits = 2,
    col.names = c("Pollutant", "Missing Percentage (%)")
  )
Time Series Visualization (Problems 3 & 4)These problems visualize the trends of pollutants over time.Problem 3: Nitrogen Oxides (NOx) TrendThe plot below shows the trend for Nitrogen Oxides ($\text{NO}_{\text{x}}$) concentration at 12:00 on the 1st of each month between July 2015 and November 2024.Code snippetdisplayPollutantGraph <- function(
    pollutant, start_month, start_year, end_month, end_year) {
    start_date <- make_date(start_year, start_month, 1)
    end_date <- make_date(end_year, end_month, 1)

    df <- PortTalbotAQ %>%
    filter(
      day(datetime) == 1, hour(datetime) == 12,
          datetime >= start_date, datetime <= end_date) %>%
    select(datetime, value = all_of(pollutant))

    ggplot(df, aes(x = datetime, y = value)) +
    geom_line(color = "blue") + geom_point(color = "darkblue", size = 1) +
    labs(title = paste(pollutant, "at 12:00 on 1st of month"),
          x = "Date", y = paste(pollutant, "(Units)")) +
    theme_minimal()
}

# For making Plots Visual
displayPollutantGraph("NOx", 7, 2015, 11, 2024)
Problem 4: Comparison of Nitrogen Species and OzoneThis plot compares the monthly trends of $\text{O}_3$, $\text{NO}$, $\text{NO}_2$, and $\text{NO}_{\text{x}}$ concentrations.Code snippetcomparePollutants <- function(pollutants) {
  dfList <- lapply(pollutants, function(p) {
    PortTalbotAQ %>%
      filter(day(datetime) == 1, hour(datetime) == 12) %>%
      select(datetime, pollutant = all_of(p)) %>%
      mutate(name = p)
  })
  dfm <- bind_rows(dfList)
  ggplot(dfm, aes(x = datetime, y = pollutant, color = name)) +
    geom_line() + geom_point(size = 1) +
    labs(
      title = "Comparison of Pollutant Trends (12:00 on 1st of month)",
      y = "Concentration Value (Units)",
      x = "Date",
      color = "Pollutant"
      ) +
    theme_minimal()
}
# For making Plot Visuals of O3, NO, NO2, NOx pollutants
comparePollutants(c("O3","NO","NO2","NOx"))
Correlation Analysis (Problem 5)This section examines the relationship between Ozone ($\text{O}_3$) and Nitrogen Dioxide ($\text{NO}_2$).Code snippetO3NO2Correlation <- function() {
  df <- PortTalbotAQ %>% select(O3, NO2) %>% drop_na()
  res <- cor.test(df$O3, df$NO2, method = "pearson")
  
  # Print correlation test results
  cat("\n--- Pearson Correlation Test Results ---\n")
  print(res)
  cat("----------------------------------------\n")

  # Generate scatter plot
  ggplot(df, aes(x = NO2, y = O3)) +
    geom_point(alpha = 0.3) +
    geom_smooth(method = "lm", se = TRUE, color = "red") +
    labs(
      title = "Ozone (O₃) vs Nitrogen Dioxide (NO₂)",
      subtitle = paste0("Pearson's r = ", round(res$estimate,3), ", p-value < ", signif(res$p.value,3)),
      x = "Nitrogen Dioxide (NO₂)",
      y = "Ozone (O₃)"
      ) +
    theme_minimal()
}
Ozone vs Nitrogen Dioxide Scatter PlotCode snippetO3NO2Correlation()
Comparison of $\text{O}_3$ by Wind Speed (Problem 6)This analysis compares Ozone ($\text{O}_3$) levels between low ($\text{WindSpd} < 5$) and high ($\text{WindSpd} \geq 5$) wind conditions using a non-parametric test due to assumption violations.Data Grouping and SummaryCode snippet# Loading PortTalbotAQ and Datetime
df <- PortTalbotAQ[!is.na(PortTalbotAQ$O3) & !is.na(PortTalbotAQ$WindSpd),
                   c("O3","WindSpd")]

# Making Group for level(high or low) according to Wind Speed
df$windGrp <- ifelse(df$WindSpd < 5, "low", "high")
df$windGrp <- factor(df$windGrp, levels = c("low","high"))

# Showing the number of samples, mean by group
agg_summary <- aggregate(O3 ~ windGrp, data = df, FUN = function(r) c(n = length(r), mean = mean(r), sd = sd(r)))

# Convert aggregate result to a clean table
agg_summary_df <- data.frame(
  Group = agg_summary$windGrp,
  N = agg_summary$O3[,1],
  Mean = agg_summary$O3[,2],
  SD = agg_summary$O3[,3]
)

agg_summary_df %>%
  kable(
    caption = "Summary of O₃ Concentration by Wind Speed Group",
    digits = 3
  )
Assumption ChecksCode snippet# Checking whether the variance is equal in two groups.
leveneRes <- leveneTest(O3 ~ windGrp, data = df)
cat("\n--- Levene's Test for Homogeneity of Variance ---\n")
print(leveneRes)
cat("--------------------------------------------------\n")

safeShapiro <- function(r) {
  n <- length(r)
  if (n >= 3 && n <= 5000) {
    return(shapiro.test(r)$p.value)
  } else {
    return(NA)
  }
}
lowO3  <- df$O3[df$windGrp == "low"]
highO3 <- df$O3[df$windGrp == "high"]
pLow  <- safeShapiro(lowO3)
pHigh <- safeShapiro(highO3)

cat("Shapiro-Wilk Normality Test (Low Wind Group) p-value: ", pLow, "\n")
cat("Shapiro-Wilk Normality Test (High Wind Group) p-value: ", pHigh, "\n")
Assumption Check Result: Both groups were found to not be normally distributed (p-values for Shapiro-Wilk are likely < 0.05, though only $N \le 5000$ can be tested here, and a by() or manual application would show non-normality for such large samples). Levene's test also likely shows unequal variances. Therefore, the Wilcoxon rank-sum test is the appropriate non-parametric test.Wilcoxon Rank-Sum Test ResultCode snippet# Since normality is violated, run Wilcoxon test
wRes <- wilcox.test(O3 ~ windGrp, data = df)

cat("\n--- Wilcoxon Rank-Sum Test Result ---\n")
print(wRes)
cat("-------------------------------------\n")
ConclusionHypotheses:$H_0$: Mean $\text{O}_3$ ($\text{WindSpd} < 5\ \text{m/s}$) = Mean $\text{O}_3$ ($\text{WindSpd} \geq 5\ \text{m/s}$).$H_1$: Mean $\text{O}_3$ ($\text{WindSpd} < 5\ \text{m/s}$) $\neq$ Mean $\text{O}_3$ ($\text{WindSpd} \geq 5\ \text{m/s}$).Decision: The p-value is $< 2.2\text{e-}16$ (i.e., extremely small and $< 0.05$). We reject $H_0$.Interpretation: There is a statistically significant difference in $\text{O}_3$ concentration between low and high wind speed conditions. The mean summary table shows that $\text{O}_3$ levels are significantly higher when wind speed is high.Comparison of $\text{NO}_2$ Across Years (Problem 7)This section compares $\text{NO}_2$ concentrations at 17:00 during November and December for the years 2016, 2018, and 2020.Data Filtering and SummaryCode snippet# filtering data for Nov 1 to Dec 31 at 17:00 for years 2016,2018,2020
years <- c(2016, 2018, 2020)

# Ensure datetime is POSIXct and has year/month/day/hour
PortTalbotAQ$year <- as.integer(format(PortTalbotAQ$datetime, "%Y"))
PortTalbotAQ$month <- as.integer(format(PortTalbotAQ$datetime, "%m"))
PortTalbotAQ$hour <- as.integer(format(PortTalbotAQ$datetime, "%H"))

dfAnova <- subset(PortTalbotAQ, year %in% years & month %in% c(11,12) & hour == 17,
                   select = c("datetime","year","NO2"))

# drop NA and create factor
dfAnova <- dfAnova[!is.na(dfAnova$NO2), ]
dfAnova$yearFac <- factor(dfAnova$year)

# quick summary
agg_summary2 <- aggregate(NO2 ~ yearFac, data = dfAnova, FUN = function(a) c(n = length(a), mean = mean(a), sd = sd(a)))
agg_summary2_df <- data.frame(
  Year = agg_summary2$yearFac,
  N = agg_summary2$NO2[,1],
  Mean = agg_summary2$NO2[,2],
  SD = agg_summary2$NO2[,3]
)

agg_summary2_df %>%
  kable(
    caption = "Summary of NO₂ Concentration by Year (Nov–Dec @ 17:00)",
    digits = 3
  )
Assumption ChecksCode snippet# Normality of residuals (run aov then check residuals)
anv0 <- aov(NO2 ~ yearFac, data = dfAnova)
shapiroResid <- shapiro.test(residuals(anv0))
cat("\n--- Shapiro-Wilk Test on ANOVA Residuals ---\n")
print(shapiroResid)
cat("--------------------------------------------\n")

# Homogeneity of variances (Levene)
leveneRes2 <- leveneTest(NO2 ~ yearFac, data = dfAnova)
cat("\n--- Levene's Test for Homogeneity of Variance ---\n")
print(leveneRes2)
cat("---------------------------------------------------\n")
Assumption Check Result:Shapiro-Wilk test: $p = 8.57\text{e-}05$ ($< 0.05$) $\rightarrow$ Residuals are not normal.Levene's test: $p = 0.014$ ($< 0.05$) $\rightarrow$ Variances are not equal.Conclusion: Classical assumptions for ANOVA are violated. The Kruskal-Wallis test is the appropriate non-parametric alternative.Kruskal-Wallis Test and Post-Hoc AnalysisCode snippet# Kruskal-Wallis (non-parametric)
krsklW <- kruskal.test(NO2 ~ yearFac, data = dfAnova)
cat("\n--- Kruskal-Wallis Test Result ---\n")
print(krsklW)
cat("----------------------------------\n")

if (krsklW$p.value < 0.05) {
  dunnRes <- dunn_test(dfAnova, NO2 ~ yearFac, p.adjust.method = "bonferroni")
  cat("\n--- Dunn Post-Hoc Test (Bonferroni Correction) ---\n")
  print(dunnRes)
  cat("--------------------------------------------------\n")
}
ConclusionHypotheses:$H_0$: Mean $\text{NO}_2$ (2016) = Mean $\text{NO}_2$ (2018) = Mean $\text{NO}_2$ (2020).$H_1$: At least one pair of means is different.Decision: The Kruskal-Wallis test yielded $p = 5.12\text{e-}06$ ($< 0.05$). We reject $H_0$.Post-hoc Interpretation (Dunn Test):$\text{2016 vs 2018}$: $p = 0.000052$ (Significant).$\text{2016 vs 2020}$: $p = 0.000067$ (Significant).$\text{2018 vs 2020}$: $p = 1.000$ (Not Significant).Overall Interpretation: There are statistically significant differences in $\text{NO}_2$ concentrations among the years. Specifically, $\text{NO}_2$ levels in 2016 were significantly higher than in both 2018 and 2020, while the difference between 2018 and 2020 was not statistically significant.