---
title: "Untitled"
author: "Tasin Islam"
date: "2025-11-12"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#Analysis of Port Talbot Hourly AQ Dataset
#==========================================


#loading all libraries
#-----------------------
  library(rstatix)
  library(readr)
  library(ggplot2)
  library(tidyverse)
  library(lubridate)
  library(car)
  library(broom)

#Problem-1
#==========

#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`
  )
#Result in DataFrame
PortTalbotAQ <- PortTalbotAQ %>%
  mutate(
    datetime = dmy_hms(
      paste(Date, Time),
      tz = "UTC")
      ) %>%
  arrange(datetime)

print(head(PortTalbotAQ))

#Problem-2
#==========

#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
  })
)
print(invalidData)

#Problem-3
#==========

#Displaying nitrogen oxides as nitrogen dioxide(NO) between July 2015 and November 2024
displayPollutantGraph <- 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() + geom_point() +
    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
#==========

#Comparing vector of 4 pollutants
comparePollutants <- 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() + labs(y = "Value", x = "Date") + theme_minimal()
}
# For making Plot Visuals of O3, NO, NO2, NOx pollutants
comparePollutants(c("O3","NO","NO2","NOx"))

#Problem-5
#==========

#Correlation between ozone and nitrogen dioxide
O3NO2Correlation <- function() {
  df <- PortTalbotAQ %>% select(O3, NO2) %>% drop_na()
  res <- cor.test(df$O3, df$NO2, method = "pearson")
  print(res)
  ggplot(df, aes(x = NO2, y = O3)) +
    geom_point() + geom_smooth(method = "lm", se = TRUE) +
    labs(title = "O3 vs NO2", subtitle = paste0("r = ", round(res$estimate,3),
                                                ", p = ", signif(res$p.value,3)))
}
O3NO2Correlation()

#Problem-6
#==========

# ==== 1. ডাটা প্রস্তুতি ====
# নিশ্চিত করুন PortTalbotAQ লোডেড এবং datetime টিকা আছে।
df <- PortTalbotAQ[!is.na(PortTalbotAQ$O3) & !is.na(PortTalbotAQ$WindSpd),
                   c("O3","WindSpd")]

# গ্রুপ বানানো
df$windGrp <- ifelse(df$WindSpd < 5, "low", "high")
df$windGrp <- factor(df$windGrp, levels = c("low","high"))

# সামারি (নমুনা সংখ্যা, গড়, স্ট্যান্ডার্ড ডেভিয়েশন)
aggregate(O3 ~ windGrp, data = df, FUN = function(r) c(n = length(r), mean = mean(r), sd = sd(r)))

# ==== 2. আনুমানিক চেক ====# ==== 2FUN = . আনুমানিক চেক ====
# Shapiro (normality) -- প্রতিটি গ্রুপে
by(df$O3, df$windGrp, function(r) {
  n <- length(r)
  if (n >= 3 && n <= 5000) {
    shapiro.test(r)
  } else {
    cat("Group size", n, "is not valid for Shapiro test. Returning NA.\n")
    NA
  }
})

# Levene test (variance equality)
leveneRes <- leveneTest(O3 ~ windGrp, data = df)
print(leveneRes)

# ==== 3. উপযুক্ত টেস্ট নির্বাচন ও চালানো ====
# যদি দুই গ্রুপই প্রায় normal: use Welch t-test (safe)
# যদি normality খারাপ হলে: use wilcox.test (nonparametric)
# Define helper function for safe Shapiro test
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("Low group Shapiro p =", pLow, "\n")
cat("High group Shapiro p =", pHigh, "\n")

# Then safely create your normality flag
nrmlOk <- (!is.na(pLow) && pLow > 0.05) && (!is.na(pHigh) && pHigh > 0.05)
if (nrmlOk) {
  tRes <- t.test(O3 ~ windGrp, data = df, var.equal = FALSE) # Welch
  print(tRes)
  # effect size
  d <- cohen.d(O3 ~ windGrp, data = df)
  print(d)
} else {
  wRes <- wilcox.test(O3 ~ windGrp, data = df)
  print(wRes)
}


#Problem-7
#==========

# filter for Nov 1 - Dec 31 at 17:00 for each year

# 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
dfAnova <- dfAnova[!is.na(dfAnova$NO2), ]

# create factor
dfAnova$yearFac <- factor(dfAnova$year)

# quick summary
table(dfAnova$yearFac)
aggregate(NO2 ~ yearFac, data = dfAnova, FUN = function(a) c(n = length(a), mean = mean(a), sd = sd(a)))

# ==== 2. Check ANOVA assumptions ====
# 2a. Normality of residuals (run ANOVA then check residuals)
anv0 <- aov(NO2 ~ yearFac, data = dfAnova)
shapiroResid <- shapiro.test(residuals(anv0))
print(shapiroResid)

# 2b. Homogeneity of variances (Levene)
leveneRes2 <- leveneTest(NO2 ~ yearFac, data = dfAnova)
print(leveneRes2)

# ==== 3. Run ANOVA or alternative ====
if (shapiroResid$p.value > 0.05 && leveneRes2$`Pr(>F)`[1] > 0.05) {
  # Use classical ANOVA
  summary(anv0)
  # effect size (eta squared)
  anvRes <- anovaTest(data = dfAnova, dv = NO2, between = yearFac)
  print(anvRes)
  eta <- dfAnova %>% anovaTest(NO2 ~ yearFac) %>% getAnvTbl() %>% .[["ges"]]
  # post-hoc Tukey
  tukeyRes <- TukeyHSD(anv0)
  print(tukeyRes)
} else {
  # Kruskal-Wallis (non-parametric)
  krsklW <- kruskal.test(NO2 ~ yearFac, data = dfAnova)
  print(krsklW)
  if (krsklW$p.value < 0.05) {
    dunnRes <- dunn_test(dfAnova, NO2 ~ yearFac, p.adjust.method = "bonferroni")
    print(dunnRes)
  }
}






