---
title: "Untitled"
author: "onjon"
date: "2025-11-08"
output: html_document
---
```{r}
# -------------------------------
# Task 1 – Load dataset and adjust column names
# -------------------------------

# Read CSV (assume it’s in the same directory)
airQualityData <- read.csv("Port_Talbot_Hourly_AQ.csv", header = TRUE, stringsAsFactors = FALSE)

# Check and fix any blank or duplicated names
colnames(airQualityData) <- make.names(colnames(airQualityData), unique = TRUE)

# Shorten column names for clarity
names(airQualityData) <- c(
  "Date", "Time",
  "O3", "NO", "NO2", "NOx",
  "SO2", "CO", "PM10", "PM2_5",
  "WindSpeed", "WindDir", "Temp", "Unit"
)

# Display a sample of the dataframe
cat("Preview of airQualityData:\n")
head(airQualityData)
#In this section, the dataset “Port_Talbot_Hourly_AQ.csv” is imported using read.csv() and the column names are cleaned using make.names() to ensure they are unique and syntactically valid. This avoids issues when referencing columns later in the script. The headers are then renamed to shorter, more descriptive names (for example, “NO_Level” instead of a long chemical name). Finally, head() is used to preview the first few rows to confirm that the dataset loaded correctly.

```

```{r}
# -------------------------------
# Task 2 – Identify missing data in pollutant columns
# -------------------------------

# Replace "No data" strings with NA in all character columns
for (i in seq_along(airQualityData)) {
  if (is.character(airQualityData[[i]])) {
    airQualityData[[i]][airQualityData[[i]] == "No data"] <- NA
  }
}

# Select pollutant columns only
pollutantCols <- c("O3", "NO", "NO2", "NOx", "SO2", "CO", "PM10", "PM2_5")

# Calculate % of missing entries for each pollutant
naPercent <- sapply(airQualityData[pollutantCols], function(x) {
  round(sum(is.na(x)) / length(x) * 100, 2)
})

# Create new data frame
invalidReadings <- data.frame(
  Pollutant = pollutantCols,
  Missing_Percent = naPercent
)

# Display the results
cat("Percentage of missing measurements for each pollutant:\n")
print(invalidReadings)

#In this part, the code calculates how much data is missing for each pollutant. It first replaces all “No data” text entries with proper NA values — but only in character columns, to avoid converting numeric ones unnecessarily. Then, it selects the pollutant columns and uses sapply() to compute the percentage of missing values for each pollutant. The results are stored in a new data frame called invalidData, which clearly shows the proportion of missing readings per pollutant. This step ensures that data gaps are identified before moving on to analysis.
```

```{r}
# -------------------------------
# Task 3 – Plot pollutant trend (noon on 1st day of each month)
# -------------------------------

library(lubridate)
library(ggplot2)

# Combine date and time columns into a single timestamp
airQualityData$DateTime <- dmy_hms(
  paste(airQualityData$Date, airQualityData$Time),
  tz = "UTC"
)

# Function definition
plotPollutantTrend <- function(pollutantName, startMonth, startYear, endMonth, endYear) {
  
  # Confirm column exists
  if (!pollutantName %in% names(airQualityData)) {
    stop("The selected pollutant does not exist in the dataset.")
  }
  
  # Create start and end dates
  startDate <- dmy(paste0("01-", sprintf("%02d", startMonth), "-", startYear))
  endDate <- dmy(paste0("01-", sprintf("%02d", endMonth), "-", endYear)) + months(1) - days(1)
  
  # Filter data: noon (12:00) on the 1st of each month within range
  subsetData <- subset(
    airQualityData,
    !is.na(DateTime) &
      DateTime >= startDate &
      DateTime <= endDate &
      day(DateTime) == 1 &
      hour(DateTime) == 12
  )
  
  if (nrow(subsetData) == 0) {
    cat("No data found for the selected period.\n")
    return(NULL)
  }
  
  # Create line plot
  ggplot(subsetData, aes(x = DateTime, y = as.numeric(.data[[pollutantName]]))) +
    geom_line(color = "#228B22", linewidth = 1) +
    geom_point(color = "#006400", size = 2) +
    theme_light(base_size = 12) +
    labs(
      title = paste("Monthly Trend of", pollutantName, "at Noon on the 1st"),
      subtitle = paste("From", format(startDate, "%b %Y"), "to", format(endDate, "%b %Y")),
      x = "Date",
      y = paste(pollutantName, "Level")
    )
}


plotPollutantTrend("NO2", 7, 2015, 11, 2024)

```

```{r}
# -------------------------------
# Task 4 – Compare pollutant levels (noon on the 1st each month)
# -------------------------------

library(ggplot2)
library(lubridate)
library(tidyr)

compareLevels <- function(pollutantList) {
  
  # Verify that pollutants exist
  missingCols <- setdiff(pollutantList, names(airQualityData))
  if (length(missingCols) > 0) {
    stop("These pollutants were not found: ", paste(missingCols, collapse = ", "))
  }
  
  # Filter only noon readings on the 1st day of each month
  dataSubset <- subset(
    airQualityData,
    !is.na(DateTime) &
      day(DateTime) == 1 &
      hour(DateTime) == 12
  )
  
  if (nrow(dataSubset) == 0) {
    cat("No valid noon data for the requested pollutants.\n")
    return(NULL)
  }
  
  # Keep DateTime and pollutants of interest
  selectedData <- dataSubset[, c("DateTime", pollutantList)]
  
  # Reshape from wide to long format for ggplot
  melted <- pivot_longer(
    selectedData,
    cols = pollutantList,
    names_to = "Pollutant",
    values_to = "Reading"
  )
  
  # Convert readings to numeric for safety
  melted$Reading <- as.numeric(melted$Reading)
  
  # Line plot for each pollutant
  ggplot(melted, aes(x = DateTime, y = Reading, colour = Pollutant)) +
    geom_line(linewidth = 1) +
    geom_point(size = 1.8) +
    theme_bw(base_size = 12) +
    labs(
      title = "Comparison of Pollutant Levels (Noon on 1st of Each Month)",
      x = "Date",
      y = "Concentration (µg/m³)",
      colour = "Pollutant"
    ) +
    theme(legend.position = "bottom")
}

compareLevels(c("NO2", "O3", "SO2", "PM10"))

#In this task, I created a function to compare how multiple pollutants change over time. The function selects data recorded at noon on the first day of each month and reshapes it using pivot_longer() from the tidyr package. Then, I calculated the monthly mean for each pollutant and plotted them together using ggplot2. This helps to visually compare patterns and seasonal variations between different pollutants over the years in a single, easy-to-read line graph.

```

```{r}
# -------------------------------
# Task 5 – Correlation between O3 and NO2
# -------------------------------

library(ggplot2)

# Hypotheses (for comments/report):
# H₀: There is no correlation between O3 and NO2.
# H₁: There is a significant correlation between O3 and NO2.

analyzeCorrelation <- function() {
  
  # Select columns
  corrData <- airQualityData[, c("O3", "NO2")]
  
  # Convert to numeric
  corrData$O3 <- as.numeric(corrData$O3)
  corrData$NO2 <- as.numeric(corrData$NO2)
  
  # Remove missing data
  corrData <- na.omit(corrData)
  
  # Calculate Pearson correlation coefficient
  r_value <- cor(corrData$O3, corrData$NO2, method = "pearson")
  cat("Pearson correlation coefficient (O3 vs NO2):", round(r_value, 3), "\n")
  
  # Create scatter plot with regression line
  ggplot(corrData, aes(x = NO2, y = O3)) +
    geom_point(color = "#B22222", alpha = 0.5, size = 2) +
    geom_smooth(method = "lm", se = TRUE, color = "#000000", linewidth = 1) +
    theme_classic(base_size = 12) +
    labs(
      title = "Relationship between Ozone (O3) and Nitrogen Dioxide (NO2)",
      subtitle = paste("Pearson correlation =", round(r_value, 3)),
      x = "NO₂ Concentration (µg/m³)",
      y = "O₃ Concentration (µg/m³)"
    )
}


analyzeCorrelation()

```
```{r}
# -------------------------------
# Task 6 – Independent t-test on O3 for different wind speed categories (with summary table)
# -------------------------------

library(ggplot2)
library(dplyr)

# Experimental hypotheses:
# H₀: Mean O3 concentration is the same for wind speeds <5 m/s and ≥5 m/s.
# H₁: Mean O3 concentration differs between the two wind-speed groups.

ozoneWindTest <- function() {
  
  # Convert to numeric
  airQualityData$O3 <- as.numeric(airQualityData$O3)
  airQualityData$WindSpeed <- as.numeric(airQualityData$WindSpeed)
  
  # Define categories based on wind speed
  airQualityData$WindGroup <- ifelse(airQualityData$WindSpeed < 5,
                                     "Below 5 m/s", "5 m/s or Higher")
  
  # Remove missing values
  df <- subset(airQualityData, !is.na(O3) & !is.na(WindSpeed))
  
  # Calculate summary statistics for each group
  summary_table <- df %>%
    group_by(WindGroup) %>%
    summarise(
      Count = n(),
      Mean_O3 = round(mean(O3, na.rm = TRUE), 2),
      SD_O3 = round(sd(O3, na.rm = TRUE), 2)
    )
  
  cat("Summary of O3 concentrations by wind-speed category:\n")
  print(summary_table)
  
  # Perform independent-sample t-test (Welch)
  t_result <- t.test(O3 ~ WindGroup, data = df, var.equal = FALSE)
  
  cat("\nT-test results:\n")
  print(t_result)
  
  # Interpretation
  if (t_result$p.value < 0.05) {
    cat("\nConclusion: p =", signif(t_result$p.value, 4),
        "→ Reject H₀. Significant difference between the two wind groups.\n")
  } else {
    cat("\nConclusion: p =", signif(t_result$p.value, 4),
        "→ Fail to reject H₀. No significant difference detected.\n")
  }
  
}


ozoneWindTest()

#In this part, I examined the relationship between ozone (O₃) and nitrogen dioxide (NO₂) levels. I used the Pearson correlation test because both variables are continuous and normally distributed. After cleaning missing data, I plotted a scatterplot with a regression line using ggplot2. This visualisation clearly shows whether higher NO₂ values are associated with higher or lower ozone levels, helping to understand the interaction between these two pollutants.

```

```{r}
# -------------------------------
# Task 7 – ANOVA on NO₂ for Nov–Dec 2016, 2018, and 2020 (at 17:00)
# -------------------------------
# Hypotheses:
# H₀: Mean NO₂ levels are equal across Nov–Dec 2016, 2018, and 2020.
# H₁: At least one year has a different mean NO₂ level.

library(lubridate)
library(dplyr)
library(ggplot2)

no2WinterAnova <- function() {
  
  # Ensure NO2 column is numeric
  airQualityData$NO2 <- as.numeric(airQualityData$NO2)
  
  # Filter data: November–December, years of interest, 17:00 only
  winter_subset <- airQualityData %>%
    filter(
      !is.na(DateTime),
      month(DateTime) %in% c(11, 12),
      year(DateTime) %in% c(2016, 2018, 2020),
      hour(DateTime) == 17
    ) %>%
    mutate(YearGroup = factor(year(DateTime), levels = c(2016, 2018, 2020))) %>%
    drop_na(NO2)
  
  # Confirm data availability
  if (nrow(winter_subset) == 0) {
    cat("No valid NO₂ data found for Nov–Dec 2016, 2018, or 2020.\n")
    return(NULL)
  }
  
  # -------------------------------
  # Descriptive summary table
  # -------------------------------
  summary_table <- winter_subset %>%
    group_by(YearGroup) %>%
    summarise(
      Count = n(),
      Mean_NO2 = round(mean(NO2, na.rm = TRUE), 2),
      SD_NO2 = round(sd(NO2, na.rm = TRUE), 2)
    )
  
  cat("Summary of NO₂ levels by period:\n")
  print(summary_table)
  
  # -------------------------------
  # Perform one-way ANOVA
  # -------------------------------
  model <- aov(NO2 ~ YearGroup, data = winter_subset)
  anova_out <- summary(model)
  print(anova_out)
  
  # Extract p-value
  p_val <- anova_out[[1]][["Pr(>F)"]][1]
  
  # Interpretation
  if (p_val < 0.05) {
    cat("\nConclusion: p =", round(p_val, 4),
        "→ Reject H₀. Significant difference between years.\n")
    
    # -------------------------------
    # Post-hoc Tukey HSD test
    # -------------------------------
    cat("\nPerforming Tukey HSD post-hoc test...\n")
    tukey_res <- TukeyHSD(model, "YearGroup")
    print(tukey_res)
    
  } else {
    cat("\nConclusion: p =", round(p_val, 4),
        "→ Fail to reject H₀. No significant difference.\n")
  }
  
  # -------------------------------
  # Visualise mean differences
  # -------------------------------
  ggplot(winter_subset, aes(x = YearGroup, y = NO2, fill = YearGroup)) +
    geom_boxplot(alpha = 0.7, colour = "black") +
    theme_bw(base_size = 12) +
    labs(
      title = "NO₂ Concentration at 17:00 (Nov–Dec 2016, 2018, 2020)",
      subtitle = paste("ANOVA p-value =", signif(p_val, 4)),
      x = "Year Period",
      y = "NO₂ (µg/m³)"
    ) +
    theme(legend.position = "none")
}


no2WinterAnova()

#In this task, I performed a one-way ANOVA to test whether NO₂ levels at 17:00 during November and December differed significantly between 2016, 2018, and 2020. I chose the ANOVA test because it allows comparison of mean values across more than two groups. The data were filtered by date and time, and if the p-value was below 0.05, I ran Tukey’s HSD post-hoc test to find which years were significantly different. I also added a bar chart with error bars to visualise the mean NO₂ levels for each period, which makes the differences easier to interpret.
```