---
title: "Airpollution"
author: "Dharmik"
date: "2025-11-08"
output: html_document
---

```{r}
# Clean PortTalbotAQ column names once and for all
names(PortTalbotAQ) <- make.names(names(PortTalbotAQ), unique = TRUE)

# Remove any unnamed or completely empty columns
PortTalbotAQ <- PortTalbotAQ[, !is.na(names(PortTalbotAQ)) & names(PortTalbotAQ) != ""]
PortTalbotAQ <- PortTalbotAQ[, colSums(!is.na(PortTalbotAQ)) > 0]

# Check all column names
print(names(PortTalbotAQ))

```

```{r}
# -------------------------------
# Question 1 – Import and rename headers
# -------------------------------

# Read the CSV file
PortTalbotAQ <- read.csv("Port_Talbot_Hourly_AQ.csv",
                         header = TRUE,
                         stringsAsFactors = FALSE)
# --- Clean and verify column names ---
names(PortTalbotAQ) <- make.names(names(PortTalbotAQ), unique = TRUE)

# Check if any names are empty or NA
print(names(PortTalbotAQ))


# Rename the headers to shorter but clear names
colnames(PortTalbotAQ) <- c(
  "DateRecorded", "TimeRecorded",
  "Ozone", "NO_Level", "NO2_Level", "NOx_Level",
  "SO2_Level", "CO_Level", "PM10", "PM2_5",
  "Wind_Speed", "Wind_Dir", "Temperature", "Unit"
)

# Display the first few rows
cat("Preview of PortTalbotAQ data frame:\n")
head(PortTalbotAQ)

#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}
# -------------------------------
# Question 2 – Calculate percentage of 'No data' readings
# -------------------------------

# Replace 'No data' with NA safely (only in character columns)
PortTalbotAQ[] <- lapply(PortTalbotAQ, function(col) {
  if (is.character(col)) {
    col[col == "No data"] <- NA
  }
  return(col)
})

# Select only the pollutant columns
pollutant_columns <- c(
  "Ozone", "NO_Level", "NO2_Level", "NOx_Level",
  "SO2_Level", "CO_Level", "PM10", "PM2_5"
)

# Compute % of missing values for each pollutant
missing_percentage <- sapply(PortTalbotAQ[pollutant_columns], function(x) {
  round(sum(is.na(x)) / length(x) * 100, 2)
})

# Create the invalidData data frame
invalidData <- data.frame(
  Pollutant = names(missing_percentage),
  MissingPercent = as.numeric(missing_percentage)
)

# Display the result
cat("Percentage of missing ('No data') readings per pollutant:\n")
print(invalidData)

#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}
# Remove all columns with no name
PortTalbotAQ <- PortTalbotAQ[, !is.na(names(PortTalbotAQ))]

# Standardise column names so they are all valid
names(PortTalbotAQ) <- make.names(names(PortTalbotAQ), unique = TRUE)

# Confirm clean names
print(names(PortTalbotAQ))

library(lubridate)

# Combine date and time columns correctly
PortTalbotAQ$DateTime <- dmy_hms(
  paste(PortTalbotAQ$DateRecorded, PortTalbotAQ$TimeRecorded),
  tz = "UTC"
)

# In this step, I cleaned and standardised the dataset before further analysis. I first removed any columns without names and used make.names() to ensure all column headers were valid and unique. Then, I combined the separate date and time columns into a single DateTime variable using dmy_hms() from the lubridate package. This was important because it allows accurate time-based filtering and plotting later in the analysis.
head(PortTalbotAQ$DateTime)

```

```{r}
library(ggplot2)

displayPollutantGraph <- function(pollutant, start_month, start_year, end_month, end_year) {
  if (!pollutant %in% names(PortTalbotAQ)) stop("Pollutant not found in dataset.")
  
  start_date <- dmy(paste0("01-", sprintf("%02d", start_month), "-", start_year))
  end_date   <- dmy(paste0("01-", sprintf("%02d", end_month), "-", end_year)) + months(1) - days(1)
  
  # Base R filtering (safe)
  valid_rows <- !is.na(PortTalbotAQ$DateTime) &
                PortTalbotAQ$DateTime >= start_date &
                PortTalbotAQ$DateTime <= end_date &
                day(PortTalbotAQ$DateTime) == 1 &
                hour(PortTalbotAQ$DateTime) == 12
  
  noon_data <- PortTalbotAQ[valid_rows, c("DateTime", pollutant)]
  
  if (nrow(noon_data) == 0) {
    cat("No noon data found between selected dates.\n")
    return(NULL)
  }
  
  ggplot(noon_data, aes(x = DateTime, y = as.numeric(.data[[pollutant]]))) +
    geom_line(color = "#0072B2", linewidth = 1) +
    geom_point(color = "#56B4E9", size = 2) +
    theme_minimal(base_size = 12) +
    labs(
      title = paste("Variation in", pollutant, "at Noon on the 1st of Each Month"),
      subtitle = paste("From", format(start_date, "%b %Y"), "to", format(end_date, "%b %Y")),
      x = "Date",
      y = paste(pollutant, "Level")
    )
}


displayPollutantGraph("NO2_Level", 7, 2015, 11, 2024)

```


```{r}
# -------------------------------
# Question 4 – Compare multiple pollutants across 10 years
# -------------------------------

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

comparePollutants <- function(pollutant_vector) {
  
  # Validate that all pollutants exist
  missing <- setdiff(pollutant_vector, names(PortTalbotAQ))
  if (length(missing) > 0) stop("Missing pollutants: ", paste(missing, collapse = ", "))
  
  #   Filter: only noon on the 1st day
  df <- PortTalbotAQ[
    !is.na(PortTalbotAQ$DateTime) &
      day(PortTalbotAQ$DateTime) == 1 &
      hour(PortTalbotAQ$DateTime) == 12,
  ]
  
  if (nrow(df) == 0) {
    cat("No matching rows found for 12:00 on 1st of each month.\n")
    return(NULL)
  }
  
  #   Extract year–month for grouping
  df$YearMonth <- format(df$DateTime, "%Y-%m")
  
  # Keep only relevant columns and reshape to long form
  tidy_df <- df[, c("YearMonth", pollutant_vector)]
  long_df <- pivot_longer(tidy_df, cols = pollutant_vector,
                          names_to = "Pollutant", values_to = "Value")
  
  #   Calculate mean value per month per pollutant
  summary_df <- long_df %>%
    group_by(YearMonth, Pollutant) %>%
    summarise(MeanValue = mean(as.numeric(Value), na.rm = TRUE), .groups = "drop")
  
  #   Convert YearMonth to a proper date (first of each month)
  summary_df$YearMonth <- as.Date(paste0(summary_df$YearMonth, "-01"))
  
  #   Plot – one smooth line per pollutant across the years
  ggplot(summary_df, aes(x = YearMonth, y = MeanValue, color = Pollutant)) +
    geom_line(linewidth = 1) +
    geom_point(size = 2) +
    theme_minimal(base_size = 12) +
    labs(
      title = "Variation in Pollutant Levels at Noon on the 1st (2015–2024)",
      x = "Year",
      y = "Average Concentration",
      color = "Pollutant"
    ) +
    theme(legend.position = "bottom")
}

comparePollutants(c("NO2_Level", "Ozone", "PM10", "SO2_Level"))


#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}
# -------------------------------
# Question 5 – Correlation between Ozone and NO2
# -------------------------------

library(ggplot2)

O3NOCorrelation <- function() {
  
  # Select only required columns
  df <- PortTalbotAQ[, c("Ozone", "NO2_Level")]
  
  # Convert to numeric (some may still be character)
  df$Ozone <- as.numeric(df$Ozone)
  df$NO2_Level <- as.numeric(df$NO2_Level)
  
  # Remove missing values
  df <- na.omit(df)
  
  # Calculate correlation coefficient (Pearson)
  correlation <- cor(df$Ozone, df$NO2_Level, method = "pearson")
  
  cat("Correlation between Ozone and NO2:", round(correlation, 3), "\n")
  
  # Scatterplot with regression line
  p <- ggplot(df, aes(x = NO2_Level, y = Ozone)) +
    geom_point(color = "#0072B2", alpha = 0.6, size = 2) +
    geom_smooth(method = "lm", se = TRUE, color = "#D55E00", linewidth = 1) +
    theme_minimal(base_size = 12) +
    labs(
      title = "Relationship between Ozone and Nitrogen Dioxide Levels",
      subtitle = paste("Pearson correlation =", round(correlation, 3)),
      x = "NO₂ Level",
      y = "O₃ Level"
    )
  
  print(p)
  
  # Return the numeric correlation value
  return(correlation)
}


O3NOCorrelation()

```

```{r}
# -------------------------------
# Question 6 – t-test on Ozone levels vs Wind Speed
# -------------------------------

library(ggplot2)

# Experimental hypotheses (for your report / comments):
# H0 (Null): Mean Ozone level is the same for wind speeds <5 and ≥5.
# H1 (Alternative): Mean Ozone level differs between the two wind-speed groups.

tTestOzoneWind <- function() {
  
  # Ensure numeric data
  PortTalbotAQ$Ozone <- as.numeric(PortTalbotAQ$Ozone)
  PortTalbotAQ$Wind_Speed <- as.numeric(PortTalbotAQ$Wind_Speed)
  
  # Create two groups
  low_wind  <- PortTalbotAQ$Ozone[PortTalbotAQ$Wind_Speed < 5 & !is.na(PortTalbotAQ$Ozone)]
  high_wind <- PortTalbotAQ$Ozone[PortTalbotAQ$Wind_Speed >= 5 & !is.na(PortTalbotAQ$Ozone)]
  
  # Perform independent-samples t-test (Welch’s version, unequal variances)
  result <- t.test(low_wind, high_wind, var.equal = FALSE)
  
  # Display results
  print(result)
  
  # Interpretation guidance
  if (result$p.value < 0.05) {
    cat("\nConclusion: p-value <", result$p.value,
        "→ Reject H0. There is a significant difference in mean Ozone levels.\n")
  } else {
    cat("\nConclusion: p-value =", result$p.value,
        "→ Fail to reject H0. No significant difference detected.\n")
  }
  
  # Optional simple boxplot
  df_plot <- data.frame(
    Group = c(rep("Wind < 5 m/s", length(low_wind)),
              rep("Wind ≥ 5 m/s", length(high_wind))),
    Ozone = c(low_wind, high_wind)
  )
  
  ggplot(df_plot, aes(x = Group, y = Ozone, fill = Group)) +
    geom_boxplot(alpha = 0.7) +
    theme_minimal(base_size = 12) +
    labs(
      title = "Ozone Levels under Different Wind Speed Conditions",
      subtitle = paste("p-value =", signif(result$p.value, 3)),
      x = "",
      y = "Ozone concentration"
    ) +
    theme(legend.position = "none")
}


tTestOzoneWind()

#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}
# -------------------------------
# Question 7 – ANOVA on NO₂ levels (Nov–Dec 2016, 2018, 2020 at 17:00)
# -------------------------------
# Experimental Hypotheses:
# H₀ (Null): The mean NO₂ levels are the same during Nov–Dec 2016, 2018, and 2020.
# H₁ (Alternative): At least one of these periods has a different mean NO₂ level.

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

anova_NO2 <- function() {
  
  # Ensure numeric data
  PortTalbotAQ$NO2_Level <- as.numeric(PortTalbotAQ$NO2_Level)
  
  # Filter the dataset:
  # Include only records from November and December of 2016, 2018, and 2020 at 17:00.
  winter_data <- PortTalbotAQ %>%
    filter(!is.na(DateTime),
           month(DateTime) %in% c(11, 12),
           year(DateTime) %in% c(2016, 2018, 2020),
           hour(DateTime) == 17) %>%
    mutate(Period = factor(year(DateTime), 
                           levels = c(2016, 2018, 2020))) %>%
    select(DateTime, NO2_Level, Period) %>%
    drop_na(NO2_Level)
  
  # Check there is data
  if (nrow(winter_data) == 0) {
    cat("No NO₂ data available for the selected dates.\n")
    return(NULL)
  }
  
  # -------------------------------
  # Perform One-Way ANOVA
  # -------------------------------
  model <- aov(NO2_Level ~ Period, data = winter_data)
  summary_result <- summary(model)
  print(summary_result)
  
  # Extract the p-value from the ANOVA summary
  p_value <- summary_result[[1]][["Pr(>F)"]][1]
  
  # -------------------------------
  # Interpret the result
  # -------------------------------
  if (p_value < 0.05) {
    cat("\nConclusion: p =", round(p_value, 4),
        "→ Reject H₀. There is a statistically significant difference in mean NO₂ levels between years.\n")
  } else {
    cat("\nConclusion: p =", round(p_value, 4),
        "→ Fail to reject H₀. No significant difference between years.\n")
  }
  
  # -------------------------------
  # Post-hoc test (if significant)
  # -------------------------------
  # Explanation:
  # If ANOVA shows a difference, we perform Tukey’s HSD test to find which pairs of years differ.
  if (p_value < 0.05) {
    cat("\nRunning Tukey HSD post-hoc test to locate the differences...\n")
    tukey <- TukeyHSD(model, "Period")
    print(tukey)
  }
  
  # -------------------------------
  # Plot mean NO₂ levels for visual comparison
  # -------------------------------
  ggplot(winter_data, aes(x = Period, y = NO2_Level, fill = Period)) +
    stat_summary(fun = mean, geom = "bar", color = "black", alpha = 0.7) +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2) +
    theme_minimal(base_size = 12) +
    labs(
      title = "NO₂ Levels at 17:00 (Nov–Dec 2016, 2018, 2020)",
      subtitle = paste("p-value:", signif(p_value, 4)),
      x = "Year (Period)",
      y = "NO₂ concentration (µg/m³)"
    ) +
    theme(legend.position = "none")
}

anova_NO2()

#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.

```