---
title: "component2: Data Analysis"
author: "Kandarp Sharda(2432187)"
output: word_document 
date: "2025-13-05"
---

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

```

## Important Points on Coding Conventions and Testing Structure

- Each task begins with a comment at the top of the code explaining its purpose.
- All variable and function names follow camelCase, except for `PortTalbotAQ`, which is written in full to show it contains the entire dataset and also it is metioned in task also that we have to create by using this variable name.
- A short explanation is provided at the end of each task to give an overview of how the code works.
- Visualisation explanations are included in the Reflection Summary, while normal outputs are explained within this documentation.
- For testing, normality and variance checks were carried out to decide whether to use parametric or non-parametric methods, and the appropriate test was then applied.
- For tasks that required logical filtering, I explained those parts of the code using small examples in comment, as shown in Task 3.
  


# Task1: Data Import, Cleaning, and Preprocessing 
```{r}
# ------------------------------------------------------------
# Task 1: Description
# ------------------------------------------------------------
# Create the PortTalbotAQ data frame by importing the dataset.
# Ensure headers are read correctly and rename them to shorter,
# clearer labels. Assume the data file is in the same folder
# as this script. Display the resulting cleaned data frame.
# ------------------------------------------------------------

#-----------------------------------------------------------
# Step 1: Import data and handle missing values
#-----------------------------------------------------------
# Import data and convert empty cells to NA
PortTalbotAQ <- read.csv("Port_Talbot_Hourly_AQ.csv",
                         stringsAsFactors = FALSE,
                         na.strings = c("No data", "", " "))

#-----------------------------------------------------------
# Step 2: Select only the necessary columns for analysis
#-----------------------------------------------------------
necessaryColumns  <- c("Date","Time","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")

# Filtering out necessary column
PortTalbotAqClean  <- PortTalbotAQ[, necessaryColumns]   

# Start from your trimmed table
PortTalbotAQ <- PortTalbotAqClean

#-----------------------------------------------------------
# Step 3: Rename columns and clean missing values
#-----------------------------------------------------------
# Rename to short, clear names
names(PortTalbotAQ) <- c(
  "date","time",
  "O3","NO","NO2","NOx",
  "SO2","CO",
  "PM10","PM2_5",
  "wind_dir","wind_speed","temp"
)


#-----------------------------------------------------------
# Step 4: Create datetime column (combine date and time)
#-----------------------------------------------------------
# Combine 'date' and 'time' into a single POSIXct datetime column
PortTalbotAQ$datetime <- as.POSIXct(
  paste(PortTalbotAQ$date, PortTalbotAQ$time),
  format = "%d-%m-%Y %H:%M",
  tz = "UTC"
)

#-----------------------------------------------------------
# Step 5: Reorder columns (move datetime to the first column)
#-----------------------------------------------------------
PortTalbotAQ <- PortTalbotAQ[, c("datetime", setdiff(names(PortTalbotAQ), "datetime"))]


#-----------------------------------------------------------
# Step 6: Convert all measurement columns to numeric
#-----------------------------------------------------------
numCols <- setdiff(names(PortTalbotAQ), c("datetime","date","time"))
for (nm in numCols) {
  # suppressWarnings() hides unnecessary conversion warnings in console output
  PortTalbotAQ[[nm]] <- suppressWarnings(as.numeric(PortTalbotAQ[[nm]]))
}

#-----------------------------------------------------------
# Step 7: Sort dataset chronologically by datetime
#-----------------------------------------------------------
PortTalbotAQ <- PortTalbotAQ[order(PortTalbotAQ$datetime), ]

#-----------------------------------------------------------
# Step 8: Quick data validation
#-----------------------------------------------------------
# Display structure and first few rows to confirm successful cleaning
PortTalbotAQ

## Code Explanation   
#In this task, I have prepared the “Port Talbot air quality” dataset for analysis. I used the “na.strings” option to handle missing values, selected only the necessary columns, and renamed them to make the data easier to read. For pollutants, I tried to use their chemical symbols for better clarity for column renaming. I also created a combined datetime column to make time-based filtering easier and converted all measurement fields to numeric values for consistency. Finally, I sorted the dataset in chronological order so that it can help in future visual representation.

```  


# Result explanation for task1
The resulting data frame shows the updated "PortTalbotA"` dataset with shorter header name. It confirms that the datetime column was created correctly and that all pollutant columns (such as O3, NO, NO2, NOx, SO2 and CO) were successfully converted to numeric values. The table also shows that the data is ordered chronologically, starting from 1 January 2015 at 01:00.


# Task 2: Calculate percentage of missing values for each pollutant
```{r}

# ------------------------------------------------------------
# Task 2: Description
# ------------------------------------------------------------
# Create a new data frame called invalidData that calculates the
# percentage of missing readings for each pollutant. Display the
# resulting table to show how much data is unrecorded for each measure.
# ------------------------------------------------------------

#-----------------------------------------------------------
# Step 1: Identify pollutant columns
#-----------------------------------------------------------
# Exclude non-numeric fields ('datetime', 'date', and 'time')
# to focus only on pollutant measurement columns.
pollutantCols <- setdiff(names(PortTalbotAQ), c("datetime", "date", "time"))

#-----------------------------------------------------------
# Step 2: Get the total number of records
#-----------------------------------------------------------
totalRows <- nrow(PortTalbotAQ)

#-----------------------------------------------------------
# Step 3: Calculate missing data percentage for each pollutant
#-----------------------------------------------------------
# Create a new dataframe that stores each pollutant name
# along with the percentage of missing values in that column
invalidData <- data.frame(
  Pollutant = pollutantCols,
  Missing_Percentage = sapply(PortTalbotAQ[, pollutantCols], function(x) {
    sum(is.na(x)) / totalRows * 100
  })
)

#-----------------------------------------------------------
# Step 4: Round and display results
#-----------------------------------------------------------
# Display nicely rounded percentages
invalidData$Missing_Percentage <- round(invalidData$Missing_Percentage, 2)

# Display the final missing value summary table
invalidData

## Code Explanation  
#In this task, I have calculated the percentage of missing data for each pollutant in the dataset. I first excluded non-numeric fields like #“datetime”, “date” and “time” to focus only on pollutant measurements. I have then used the sapply() function to efficiently apply the #is.na() check across all pollutant columns and calculate the percentage of missing values, storing the results in a new dataframe.

```
# Result explanation for task2
The invalidData table shows the percentage of missing readings for each pollutant. Pollutants like O3, NO, NO2 and wind variables have very low missing values, while PM10 and PM2.5 have slightly higher gaps. 


# Task 3 – Create a function to display monthly noon pollutant trends (displayPollutantGraph)   
```{r}   
# ------------------------------------------------------------
# Task 3: Description
# ------------------------------------------------------------
# Create a function displayPollutantGraph() that takes a pollutant name,
# a start month/year and an end month/year. It should extract readings
# taken at noon on the 1st of each month and plot how the pollutant
# changes over time. Use this function to plot NOx levels from
# July 2015 to November 2024.
# ------------------------------------------------------------

displayPollutantGraph <- function(data, pollutantName, startMonth, startYear, endMonth, endYear){
  
  #-----------------------------------------------------------
  # Step 1: Create time boundaries for the selected period
  #-----------------------------------------------------------
  # It Represents the beginning of the start month (1st day at 00:00).
  start <- as.POSIXct(sprintf("01-%02d-%04d 00:00", startMonth, startYear),
                      format = "%d-%m-%Y %H:%M", tz = "UTC")
  
  # It marks the beginning of the end month
  endMonthStart <- as.POSIXct(sprintf("01-%02d-%04d 00:00", endMonth, endYear),
                                format = "%d-%m-%Y %H:%M", tz = "UTC")
  
  # It represents the last moment of the end month (23:59 on the final day)
  end <- seq(endMonthStart, length = 2, by = "1 month")[2] - 1
  
  # Example:
  # Suppose the function is called with:
  #   startMonth = 7,  start_year = 2015
  #   endMonth   = 11, end_year   = 2024
  #
  # Then:
  #   start will be "2015-07-01 00:00:00 UTC" -> beginning of July 2015
  #   endMonthStart will be "2024-11-01 00:00:00 UTC" -> beginning of November 2024
  #   end will be "2024-11-30 23:59:59 UTC" -> gets the last moment of the end month by taking 
  # the start of the next month and subtracting one second

  
  #-----------------------------------------------------------
  # Step 2: Filter data for 1st day of each month at 12:00 (noon)
  #-----------------------------------------------------------
  idx <- format(data$datetime, "%d") == "01" &
         format(data$datetime, "%H") == "12" &
         format(data$datetime, "%M") == "00" &
         data$datetime >= start & data$datetime <= end

  #-----------------------------------------------------------
  # Step 3: Subset data and rename columns for clarity
  #-----------------------------------------------------------
  df <- data[idx, c("datetime", pollutantName)]
  names(df) <- c("datetime", "value")
  
  # Check if any data is available within the specified range
  if (nrow(df) == 0) {
    message("No data found for the chosen range/time filter.")
    return(invisible(NULL))
  }
  
  #-----------------------------------------------------------
  # Step 4: Plot variation in pollutant levels
  #-----------------------------------------------------------
  plot(df$datetime, df$value, type = "l",
         xlab = "Date (1st of month at 12:00)",
         ylab = pollutantName,
         main = sprintf("%s – Monthly Noon Levels (%02d/%04d to %02d/%04d)",
                        pollutantName, startMonth, startYear, endMonth, endYear))
  points(df$datetime, df$value, pch = 16, cex = 0.6)
  grid()
}

displayPollutantGraph(
  data = PortTalbotAQ,
  pollutantName = "NOx",
  startMonth = 7, startYear = 2015,
  endMonth = 11, endYear = 2024
  
)

## Code Explanation 
#In this task, I created a function to display pollutant variations over a specific time period. I built start and end time boundaries using #as.POSIXct() and filtered the dataset to include readings taken at noon on the first day of each month. While doing this task I also learned #to create formatted strings using the sprint(), similar to C syntax.

```

# Task 4 – Create a function to display monthly noon pollutant trends (comparePollutants)   
```{r}
# ------------------------------------------------------------
# Task 4: Description
# ------------------------------------------------------------
# Create a function comparePollutants() that takes a vector of pollutant
# names and plots how each pollutant changes over time on the same graph.
# Use noon readings on the 1st of each month across all ten years.
# Demonstrate this using a suitable vector of four pollutants.
# -----------------------------------------------------------

comparePollutants <- function(data, pollutants){

  #-----------------------------------------------------------
  # Step 1: Filter records for 1st day of each month at 12:00 (noon)
  #-----------------------------------------------------------
  # This ensures the comparison is made using consistent time points
  # across all years (same day and hour each month).
  idx <- format(data$datetime, "%d") == "01" &
         format(data$datetime, "%H") == "12" &
         format(data$datetime, "%M") == "00" 

  #-----------------------------------------------------------
  # Step 2: Subset data to keep datetime and selected pollutants
  #-----------------------------------------------------------
  df <- data[idx, c("datetime", pollutants)]  
  
  
  #-----------------------------------------------------------
  # Step 3: Auto-generate distinct colours for each pollutant line
  #-----------------------------------------------------------   
  colors <- rainbow(length(pollutants))
  
  #-----------------------------------------------------------
  # Step 4: Plot the first pollutant as the base line chart
  #-----------------------------------------------------------
  # Calculate the range of all pollutant values to determine appropriate y-axis limits
  # This prevents data points from being cut off when the plot area is adjusted
  yMax <- max(df[, pollutants], na.rm = TRUE)
  yMin <- min(df[, pollutants], na.rm = TRUE)

  # Adjust plot margins to make room for legend on the right
  # mar = c(bottom, left, top, right); increased right margin from 2 to 8
  par(mar = c(5, 4, 4, 8), xpd = TRUE)

  # Create base plot with adjusted y-axis limits (5% padding at top)
  plot(df$datetime, df[[pollutants[1]]], type = "l",
     col = colors[1], lwd = 2,
     xlab = "Date", ylab = "Concentration level",
     main = "Pollutant Comparison at Noon", 
     ylim = c(yMin, yMax * 1.05))

  
  #-----------------------------------------------------------
  # Step 5: Add lines for the remaining pollutants
  #-----------------------------------------------------------
  # Each pollutant is plotted in a different colour
  # to allow visual comparison over the same timeline.
  for (i in 2:length(pollutants)) {
      lines(df$datetime, df[[pollutants[i]]], col = colors[i], lwd = 2)
  }
    
  #-----------------------------------------------------------
  # Step 6: Add legend to identify each pollutant
  #-----------------------------------------------------------
  # The legend is positioned at the top-left corner for visibility. 
  legend(x = "topright", inset = c(-0.25, 0),
         legend = pollutants, col = colors, lwd = 2, bty = "n")
  
  # Reset margins
  par(mar = c(5, 4, 4, 2) + 0.1, xpd = FALSE)

}

pollutants <- c("O3", "NO2", "PM10", "PM2_5")
comparePollutants(PortTalbotAQ, pollutants)

##Code Explanation  
#In this task, I created a function O3NOCorrelation() to compare how multiple pollutants changed over a set time period. I filtered the data #to include noon readings from the first day of each month and used a multi-line plot to show all trends together. This made seasonal patterns #and differences between pollutants easier to understand(explained in the visualisation section). I also learned to adjust the legend using #the par() function so it did not overlap with the graph.
```

# Task-5 - Analyse correlation between Ozone (O3) and Nitrogen Dioxide (NO2)
```{r}

# ------------------------------------------------------------
# Task 5: Description
# ------------------------------------------------------------
# Create a function O3NOCorrelation() that calculates the correlation
# between ozone (O3) and nitrogen dioxide (NO2) for the full dataset.
# Use the cleaned data, handle missing values, and choose an appropriate
# correlation method. Plot the data as a scatterplot to show the
# relationship between O3 and NO2.
# ------------------------------------------------------------

O3NOCorrelation <- function(data){
  
  #-----------------------------------------------------------
  # Step 1: Clean the data (keep only complete O3 & NO2 rows)
  #-----------------------------------------------------------
  df <- data[, c("O3", "NO2")]
  df <- df[complete.cases(df), ]
  

  #-----------------------------------------------------------
  # Step 2: Check normality (Shapiro-Wilk test)
  # The Shapiro–Wilk test only supports sample sizes between 3 and 5000.
  # Therefore, if there are more than 5000 observations, we take a random
  # subset of 5000 values to make the test valid and still representative.
  #-----------------------------------------------------------
  
  # Reason for sampling it because Shapiro-Wilk  works perfectly till 5000 observations
  sampleO3 <- if (length(df$O3) > 5000)  sample(df$O3, 5000) else df$O3  
  sampleNO2 <- if (length(df$NO2) > 5000)  sample(df$NO2, 5000) else df$NO2
  
  pO3 <- shapiro.test(sampleO3)$p.value 
  pNO2 <- shapiro.test(sampleNO2)$p.value 
  
  # Decide correlation method
  usePearson  <- (pO3 >= 0.05 && pNO2 >= 0.05)
  method <- if (usePearson) "pearson" else "spearman"
  
  # Explain the decision
  if (usePearson) {
    cat("Both O3 and NO2 are normally distributed (p >= 0.05),",
      "so using Pearson correlation (parametric).\n")
  } else {
    cat("At least one variable is not normally distributed (p < 0.05),",
      "so using Spearman correlation (non-parametric).\n")
  }
  
  #-----------------------------------------------------------
  # Step 3: Plot sprintter + trend lines
  #-----------------------------------------------------------
  plot(df$NO2, df$O3,
       main = paste("O3 vs NO2 -", toupper(method)),
       xlab = "Nitrogen Dioxide (NO2)",
       ylab = "Ozone (O3)",
       pch = 19, col = rgb(0.1, 0.4, 0.8, 0.5))
  
  # Add LOWESS (non-linear) trend line
  lines(lowess(df$NO2, df$O3), col = "red", lwd = 2)
  
  # Add linear regression line
  abline(lm(O3 ~ NO2, data = df), col = "darkgreen", lty = 2, lwd = 2)
  
  
  #-----------------------------------------------------------
  # Step 4: Show correlations
  #-----------------------------------------------------------
  corTest <- cor.test(df$O3, df$NO2, method = method, exact = FALSE)
  correlationValue <- unname(corTest$estimate)
  correlationP     <- corTest$p.value
  
  cat(sprintf("Correlation coefficient (method = %s): %.3f\n", method, correlationValue))
  
  if (correlationValue > 0) {
      cat("Positive correlation: As NO2 increases, O3 also increases.\n")
  } else if (correlationValue < 0) {
      cat("Negative correlation: As NO2 increases, O3 decreases.\n")
  } 
  
  legend("topright",
       legend = c(
         # Show which correlation method and value were used
         sprintf("%s correlation = %.3f", toupper(method), correlationValue),
         
         # Explain what the red and green lines mean on the plot
         "Red line = LOWESS (non-linear trend)",
         "Green dashed line = Linear regression"
       ),
       
       # Match colours and line styles to each description
       col = c( "red", "darkgreen"),
       lwd = c(NA, 2, 2),
       lty = c(NA, 1, 2),
       
       # Remove legend box and adjust text size
       bty = "n",
       cex = 0.9)
  
  #-----------------------------------------------------------
  # Step 5: Results table
  #-----------------------------------------------------------
  resultsTable <- data.frame(
    Metric = c(
      "Shapiro-Wilk(O3)",
      "Shapiro-Wilk(NO2)",
      "Correlation method",
      "Correlation coefficient",
      "Correlation p-value",
      "Hypotheses Decision"
    ),
    Value = c(
      pO3,
      pNO2,
      toupper(method),
      correlationValue,
      correlationP,
      ifelse(correlationP > 0.05,
             "Fail to reject H0 - No significant association between O3 and NO2",
             "Reject H0 - Significant association between O3 and NO2")
    ),
    stringsAsFactors = FALSE
  )
  
  resultsTable
  
  }

O3NOCorrelation(PortTalbotAQ)

## Code Explanation   
#For this task, I kept only complete pairs of NO2 and O3, to avoid bias from missing values. I tested normality using the Shapiro-Wilk test, #using random sampling when the sample exceeded 5000. Based on the result, I used Pearson for normal data and Spearman otherwise. I then #plotted a scatter graph with linear and LOWESS fits and created a small results table showing the p-values, chosen test, and outcome. Here I #learned that the Shapiro-Wilk test cannot handle samples beyond 5,000, which led me to use random sampling to work within its limitations.
```
# Result explanation for task5
The Shapiro-Wilk test for both O3 and NO2 produced extremely small p-values, meaning both variables strongly deviate from a normal distribution. Because normality was not met, the Spearman method was used for correlation. The Spearman coefficient is -0.412, indicating a moderate negative relationship between O3 and NO2. The correlation p-value is shown as 0, meaning it is so small that R rounds it down to zero. This makes the result highly significant. Therefore, we reject the null hypothesis and conclude that there is a significant association between O3 and NO2.


# Task-6  Test whether Ozone (O3) levels differ across wind-speed groups
```{r}
# ------------------------------------------------------------
# Task 6: Description
# ------------------------------------------------------------
# Perform a t-test to check whether ozone (O3) levels differ
# between two wind-speed groups:
#   - Group 1: wind speed < 5 m/s
#   - Group 2: wind speed >= 5 m/s
#
# State the hypotheses clearly, run the test, and report whether
# the result leads to accepting or rejecting the null hypothesis.
# ------------------------------------------------------------

library(car)

# ------------------------------------------------------------
# Hypotheses
# ------------------------------------------------------------
# H0: The mean O3 (ozone) concentration is equal between the two wind-speed groups
# that mean wind speed < 5 m/s and wind speed >= 5 m/s.
#
# H1: The mean O3 (ozone) concentration differs significantly between
#     the two wind-speed groups.
# ------------------------------------------------------------


# ------------------------------------------------------------
# Step 1: Split the data into two groups based on wind speed
# ------------------------------------------------------------
o3LowWind  <- PortTalbotAQ[PortTalbotAQ$wind_speed < 5, "O3"]
o3HighWind <- PortTalbotAQ[PortTalbotAQ$wind_speed >= 5, "O3"]


# Remove missing values
o3LowWind  <- na.omit(o3LowWind)
o3HighWind <- na.omit(o3HighWind)


# ------------------------------------------------------------
# Step 2: Visualise the distributions for both groups
# ------------------------------------------------------------
plot(density(o3LowWind), 
     col = "blue", 
     lwd = 2,
     main = "Normal Distribution Curves of O3 Levels by Wind Speed",
     xlab = "O3 Concentration",
     ylab = "Density")

# Add the second group curve
lines(density(o3HighWind), col = "red", lwd = 2, lty = 2)

legend("topright",
       legend = c(
        "o3LowWind", "o3HighWind"
       ),
       col = c("blue", "red"),
       lwd = c(2, 3),
       lty = c(2, 1),
       bty = "n",
       cex = 0.9
    )


# ------------------------------------------------------------
# Step 3: Normality Test (Shapiro-Wilk)
# ------------------------------------------------------------
# The Shapiro-Wilk test checks if data follows a normal distribution.
# It is valid for sample sizes between 3 and 5000 observations.

# If sample size > 5000, randomly sample 5000 values for validity.
sampleLow  <- if (length(o3LowWind)  > 5000) sample(o3LowWind,  5000) else o3LowWind
sampleHigh <- if (length(o3HighWind) > 5000) sample(o3HighWind, 5000) else o3HighWind

pLow  <- shapiro.test(sampleLow)$p.value
pHigh <- shapiro.test(sampleHigh)$p.value

cat(sprintf("Shapiro-Wilk p-values:\n  Low Wind: %.10f\n  High Wind: %.10f\n\n", pLow, pHigh))      

# ------------------------------------------------------------
# Step4 : Test equal variance    
# ------------------------------------------------------------
# #Perform Levene’s Test for equal variances
tmp <- data.frame(
    value = c(o3LowWind, o3HighWind),
    group = factor(c(rep("Low", length(o3LowWind)),
                     rep("High", length(o3HighWind))))
  )
  
levRes   <- leveneTest(value ~ group, data = tmp)
pLevene  <- levRes$`Pr(>F)`[1]
cat("Levene’s Test p-value:", round(pLevene, 6), "\n")

# ------------------------------------------------------------
# Step 5: Decide on test type based on normality
# ------------------------------------------------------------
# If both groups are normally distributed (p >= 0.05), use a parametric test (t-test).
# Otherwise, use a non-parametric test (Wilcoxon).
if (pLow >= 0.05 && pHigh >= 0.05 && pLevene >= 0.05) {
  
  cat("Both groups are normally distributed — using independent t-test (parametric).\n\n")
  testResult <- t.test(o3LowWind, o3HighWind, var.equal = FALSE)
  
} else {
  
  cat("At least one group is not normally distributed — using Wilcoxon rank-sum test (non-parametric).\n\n")
  testResult <- wilcox.test(o3LowWind, o3HighWind)
}

# ------------------------------------------------------------
# Step 5: Display test result and interpretation
# ------------------------------------------------------------
resultsTable <- data.frame(
  Metric = c(
    "Shapiro-Wilk (Low Wind)",
    "Shapiro-Wilk (High Wind)",
    "Levene’s Test (Equal Variance)",
    "Main Test p-value",
    "Test Used",
    "Hypotheses Decision"
  ),
  
  Value = c(
    pLow,
    pHigh,
    pLevene,
    testResult$p.value,
    testResult$method,
    ifelse(testResult$p.value > 0.05,
           "Fail to reject H0: O3 concentrations are similar across wind-speed groups (no significant difference)",
           "Reject H0: O3 concentrations differ significantly between the two wind-speed groups")
  ),
  
  stringsAsFactors = FALSE
)

resultsTable

## Code explanationn    
#In this task, I split the O3 data into low-wind and high-wind groups and removed missing values. I used density plots to compare their #distributions, then checked normality with the Shapiro-Wilk test and variances with Levene’s test. Depending on the results, I used either a #parametric t-test or a non-parametric Wilcoxon test. I then created a clear results table showing the p-values, chosen test, and conclusion.
```
# Result explanation for task6
The Shapiro-Wilk tests for both low-wind and high-wind groups gave extremely small p-values, so neither group is normally distributed. Levene’s test also gave a p-value reported as 0, which does not mean the value is exactly zero, but that it is so small R rounds it down. This shows the variances are clearly unequal. Because normality and equal-variance assumptions were both violated, a non-parametric test was needed. The Wilcoxon test was used, and its p-value is also shown as 0 (again meaning “extremely small”). This indicates a highly significant difference between the two groups. Therefore, we reject the null hypothesis and conclude that O3 concentrations differ significantly between low-wind and high-wind conditions.



# Task 7 – Comparing Nitrogen Dioxide (NO2) Levels Across Years(2016, 2018, 2020) Using ANOVA / Kruskal–Wallis
```{r}   
# ------------------------------------------------------------
# Task 7: Description
# ------------------------------------------------------------
# Perform an ANOVA test to check whether NO2 levels differ between:
#   - Nov–Dec 2016
#   - Nov–Dec 2018
#   - Nov–Dec 2020
#
# Use data recorded at 17:00 each day between 1 Nov and 31 Dec
# for all three years.
#
# If a significant difference is found, run the appropriate post-hoc
# test (e.g., Tukey or pairwise Wilcoxon) to identify which years
# differ and explain why this post-hoc test is chosen.
#
# Clearly state the hypotheses, report the ANOVA result, and indicate
# whether the null hypothesis is accepted or rejected.

library(car)

# ------------------------------------------------------------
# Hypotheses
# ------------------------------------------------------------
# H0: The mean NO2 concentration at 17:00 (Nov-Dec) is equal across 2016, 2018, and 2020.
# H1: At least one year’s mean NO2 concentration differs.
# ------------------------------------------------------------

pollutantName <- "NO2"

# ------------------------------------------------------------
# Step 1: Define helper function to fetch pollutant data
# ------------------------------------------------------------
fetchPollutantData <- function(data, pollutantName,
                               startMonth, startYear,
                               endMonth, endYear) {
  
  # Build start and end time bounds
  start <- as.POSIXct(sprintf("01-%02d-%04d 00:00",
                              startMonth, startYear),
                      format = "%d-%m-%Y %H:%M", tz = "UTC")
  
  endMonthStart <- as.POSIXct(sprintf("01-%02d-%04d 00:00",
                                      endMonth, endYear),
                              format = "%d-%m-%Y %H:%M", tz = "UTC")
  
  # Last moment of the end month
  end <- seq(endMonthStart, length = 2, by = "1 month")[2] - 1
  
  # Keep all records at 17:00 between the start and end dates
  idx <- format(data$datetime, "%H") == "17" &
         format(data$datetime, "%M") == "00" &
         data$datetime >= start & data$datetime <= end
  
  # Subset and rename
  df <- data[idx, c("datetime", pollutantName)]
  names(df) <- c("datetime", "value")
  
  if (nrow(df) == 0) {
    message("No data found for the chosen range/time filter.")
    return(invisible(NULL))
  }
  
  return(df)
}

# ------------------------------------------------------------
# Step 2: Fetch NO2 data for the required years
# ------------------------------------------------------------
no2_2016 <- fetchPollutantData(PortTalbotAQ, pollutantName, 11, 2016, 12, 2016)
no2_2018 <- fetchPollutantData(PortTalbotAQ, pollutantName, 11, 2018, 12, 2018)
no2_2020 <- fetchPollutantData(PortTalbotAQ, pollutantName, 11, 2020, 12, 2020)

# ------------------------------------------------------------
# Step 3: Normality test (Shapiro-Wilk)
# ------------------------------------------------------------
# Shapiro-Wilk works for 3 <= n <= 5000, and NA values are not allowed.
# Each group has <5000 points, so no sampling required.

pNo2_2016 <- shapiro.test(na.omit(no2_2016$value))
pNo2_2018 <- shapiro.test(na.omit(no2_2018$value))
pNo2_2020 <- shapiro.test(na.omit(no2_2020$value))

# Display p-values
cat("Shapiro–Wilk Normality p-values:\n")
cat("  2016:", pNo2_2016$p.value, "\n")
cat("  2018:", pNo2_2018$p.value, "\n")
cat("  2020:", pNo2_2020$p.value, "\n\n")


# ------------------------------------------------------------
# Step 4: Combine all years into one data frame
# ------------------------------------------------------------
# ANOVA / Kruskal require a single column of values
# and a grouping column (value, year).

no2_2016$year <- "2016"
no2_2018$year <- "2018"
no2_2020$year <- "2020"

combinedDf <- rbind(no2_2016, no2_2018, no2_2020)
combinedDf$year <- factor(combinedDf$year, levels = c("2016", "2018", "2020"))

# ------------------------------------------------------------
# Step :5 Test equal variance    
# ------------------------------------------------------------
# #Perform Levene’s Test for equal variances
leveneRes <- leveneTest(value ~ year, data = combinedDf)
pLevene <- leveneRes$`Pr(>F)`[1]

if (pLevene > 0.05) {
  cat("Variances are equal across groups (p >", round(pLevene, 3), ") \n")
} else {
  cat("Variances are significantly different (p =", round(pLevene, 3), ")\n")
}

# ------------------------------------------------------------
# Step 6: Choose test type based on normality and variance equality   
# ------------------------------------------------------------
# If all groups are normally distributed (p > 0.05) use ANOVA.
# Otherwise use Kruskal-Wallis (non-parametric equivalent).
# ------------------------------------------------------------

if (pNo2_2016$p.value >= 0.05 &
    pNo2_2018$p.value >= 0.05 &
    pNo2_2020$p.value >= 0.05 & pLevene >= 0.05) {
  
  cat("All groups approximately normal using One-Way ANOVA.\n\n")
  testUsed <- "ANOVA"
  
  # -------------------------------
  # Perform One-Way ANOVA
  # -------------------------------
  result <- aov(value ~ year, data = combinedDf)
  print(summary(result))
  
  # Extract ANOVA p-value
  mainP <- summary(result)[[1]][["Pr(>F)"]][1]
  
  # ----------------------------------------------------------
  # Step 7: Post-hoc analysis (Tukey HSD) if ANOVA significant
  # ----------------------------------------------------------
  # ANOVA only shows that at least one mean differs.
  # Tukey HSD identifies which specific pairs of years differ.
  # ----------------------------------------------------------
  
  if (!is.na(pAnova) && pAnova < 0.05) {
    cat("\nConclusion: Reject H0 - At least one year's mean NO2 differs.\n")
    cat("Post-hoc Test: Tukey HSD (pairwise mean comparisons)\n\n")
    print(TukeyHSD(result))
  } else {
    cat("\nConclusion: Fail to reject H0 - Means appear similar across years.\n")
  }
  
} else {
  
  cat("At least one group not normally distributed so  using Kruskal-Wallis test.\n\n")
  testUsed <- "Kruskal-Wallis"

  
  # -------------------------------
  # Perform Kruskal-Wallis Test
  # -------------------------------
  result <- kruskal.test(value ~ year, data = combinedDf)
  mainP <- result$p.value
  
  # ----------------------------------------------------------
  # Step :7 Post-hoc analysis (Pairwise Wilcoxon) if significant
  # ----------------------------------------------------------
  # Kruskal-Wallis compares medians across groups.
  # If significant, we perform pairwise Wilcoxon rank-sum tests
  # to find which year pairs differ.
  # The Benjamini-Hochberg (BH) adjustment corrects
  # for multiple comparisons to limit false positives.
  # ----------------------------------------------------------
  
  if (!is.na(result$p.value) && result$p.value < 0.05) {
    ("\nConclusion: Reject H0 - At least one year's median NO2 differs.\n")
    cat("Post-hoc Test: Pairwise Wilcoxon (BH-adjusted p-values)\n\n")
    print(pairwise.wilcox.test(combinedDf$value,
                               combinedDf$year,
                               p.adjust.method = "BH"))
  } else {
    print("\nConclusion: Fail to reject H0 - Medians appear similar across years.\n")
  }
}

# ------------------------------------------------------------
# Step 7: Results summary table (p-values + chosen test)
# ------------------------------------------------------------
resultsTable <- data.frame(
  Metric = c("Shapiro 2016", 
             "Shapiro 2018", 
             "Shapiro 2020", 
             "Levene’s Test (Equal Variance)", 
             "Main Test p-value", 
             "Test Used", 
             "Hypotheses Decision"),
  
  Value = c(
    round(pNo2_2016$p.value, 6),
    round(pNo2_2018$p.value, 6),
    round(pNo2_2020$p.value, 6),
    round(pLevene, 6),
    round(mainP, 6),
    testUsed,
    ifelse(mainP > 0.05,
           "Fail to reject H0: Mean NO2 concentrations are similar (no significant difference)",
           "Reject H0: At least one year's mean NO2 concentration differs significantly")
  ),
  
  stringsAsFactors = FALSE
)


# ------------------------------------------------------------
# Step 8: Boxplot visualisation for clarity
# ------------------------------------------------------------
boxplot(value ~ year, data = combinedDf,
        main = "NO2 Concentration at 17:00 (Nov–Dec, 2016–2020)",
        xlab = "Year",
        ylab = "NO2 Concentration",
        col = c("skyblue", "salmon", "lightgreen"),
        border = "gray40")
grid()
resultsTable

## Code explatnation  
#In this task, I compared NO2 levels across 2016, 2018 and 2020 (November–December) using ANOVA and Kruskal-Wallis tests. I first checked #normality with the Shapiro-Wilk test and variances with Levene’s test to decide the correct method. When assumptions were met, I used one-way #ANOVA with Tukey’s HSD; otherwise, I used the Kruskal-Wallis test with a pairwise Wilcoxon post-hoc test. 

```
# Result explanation for task7
The Shapiro-Wilk tests showed that NO2 values for 2018 and 2020 were not normally distributed (p < 0.05), while 2016 was borderline. Levene’s test also produced a significant result (p = 0.014), indicating unequal variances across the three years. Because the assumptions for ANOVA were not met, the Kruskal-Wallis test was used. The test returned a very small p-value, meaning at least one year’s NO2 concentration is significantly different. Therefore, we reject the null hypothesis.

A post-hoc pairwise Wilcoxon test was then carried out to identify where these differences occur. The results show significant differences between 2016 vs 2018 and 2016 vs 2020. However, 2018 vs 2020 shows no significant difference. This means that 2016 is the year that differs most, while 2018 and 2020 have similar NO2 levels.