# ##############################################################################
# Student: Nate Fox 2311568
# Module: Statistical Programming (ACCB5021)
# File: airpollution.R
#
# This script will:
# 1. Load the air quality data set
# 2. Rename columns to make them easier to work with in R
# 3. Prepare the data for later analysis and visualization
# ##############################################################################

options(scipen=999) # Turns off scientific notation


install.packages ("readxl")  # Adding a package to read XLS file. 
library (readxl)

# ##############################################################################
# Task 1 Question 1 - Loading data and renaming column headers
# ##############################################################################

PortTalbotAQ<-read_excel("Port Talbot Hourly AQ.xlsx", sheet=1)

colnames(PortTalbotAQ) <- c(
  "date", "time", "ozone", "ozone_status", "nitric_oxide", "nitric_oxide_status",
  "nitrogen_dioxide", "nitrogen_dioxide_status", "nitrogen_oxides",
  "nitrogen_oxides_status", "sulphur_dioxide", "sulphur_dioxide_status",
  "carbon_monoxide", "carbon_monoxide_status", "PM10", "PM10_status", "PM2_5",
  "PM2_5_status", "wind_direction", "wind_direction_status", "wind_speed",
  "wind_speed_status", "temperature", "temperature_status"
)

# As the first 10 rows consist of non-usable data or now a duplication of heading will remove the first ten rows.

PortTalbotAQ <- PortTalbotAQ[-c(1:10), ] # Removing row 1-10 but keeping the columns.

# ##############################################################################
# Task 1 Question 2 - Creating new data frame showing the percentage of readings where no measurement has been recorded for each pollutant.
# ##############################################################################

# Identifying pollutant columns only, these are the measured air pollutants, excluding status and weather variables. 
pollutant_cols <- c(
  "ozone", "nitric_oxide", "nitrogen_dioxide", "nitrogen_oxides", "sulphur_dioxide", 
  "carbon_monoxide", "PM10", "PM2_5"
)

pollutant_cols # Check to see if worked and has an output 

invalid_percentages <- sapply(pollutant_cols, function(pollutant) { # Calculating percentage of "No data" readings for each pollutant
  no_data_count <- sum(PortTalbotAQ[[pollutant]] == "No data", na.rm = TRUE) # Counting the number of "No data" missing measurements
  (no_data_count / nrow(PortTalbotAQ)) * 100 # Calculate the percentage of invalid readings 
})

invalid_percentages # Check to see if worked and has an output

# Time to construct the data frame to store the results
invalidData <- data.frame(
  pollutant = names(invalid_percentages),
  percent_no_data = invalid_percentages
) # ran and data table appeared

invalidData # Check to see if worked and has although not very readable

invalidData$percent_no_data <- round(invalidData$percent_no_data,2) # Rounds the data to 2 decimal places (which needed to come before percentage)

invalidData$percent_no_data <- paste0(invalidData$percent_no_data, "%") # Adds the percentage symbol at the end of the data

invalidData # Checked to see if above worked. Rounded and percentage symbol included looks neater.
# The resulting data frame shows the proportion of missing measurements for each 
# pollutant across the full collection period. 

# ##############################################################################
# Task 1 Question 3 - Function to plot pollutant levels at noon on the 1st of each month
# for a given date range
# ##############################################################################

displayPollutantGraph <- function(pollutant, start_month, start_year, end_month, end_year){ # Defining a function with the pollutant
  # column name and the start and end of the date ranges
  
  PortTalbotAQ$date <- as.Date(PortTalbotAQ$date, format = "%Y-%m-%d") # Convert to date format for assist with plotting on a time axis later
  
  start_date <- as.Date(paste(start_year, sprintf("%02d", start_month), "01", sep = "-")) # Creating the start date of the requested range
  end_date <- as.Date(paste(end_year, sprintf("%02d", end_month), "01", sep ="-")) # Creating the end date of the requested range
  
  keep_rows <- PortTalbotAQ$date >= start_date & # Creating a logical vector identifying start date
    PortTalbotAQ$date <= end_date & # Identifying end date
    format(PortTalbotAQ$date, "%d") == "01" & # Identifying the range will start on the first of the month
    PortTalbotAQ$time == "12:00:00" # Identifying what was recorded precisely at noon (12:00:00)
  
  df <- PortTalbotAQ[keep_rows, c("date", "time", pollutant)] # Subset the original data frame to keep only rows that meet above conditions
  
  df[[pollutant]][df[[pollutant]] == "No data"] <- NA # Replacing any "No data" entries in the pollution column with NA
  df[[pollutant]] <- as.numeric(df[[pollutant]]) # Converts pollutant values from character to numeric
  df <- df[!is.na(df[[pollutant]]), ] # Removal of NA's to ensure any data pulled are valid measurements 
  
  if (nrow(df) == 0) { # Adding so if no valid data remains after filtering it will stop execution. 
    stop(paste("No valid", pollutant, "data found for the specified range")) # Display informative error message
  }
  
  plot(
    df$date,
    df[[pollutant]],
    type = "o", # Controls how the points are drawn (o = points and lines)
    pch = 16, # Controls how the points are visualized (16 = filled circle)
    xlab = "date",
    ylab = paste(pollutant, "level"),
    main = paste("variation in", pollutant, "at noon on the 1st of each month")
    )
}

displayPollutantGraph("nitrogen_oxides", 7, 2015, 11, 2024) # After June 2022 there are no valid nitrogen_oxides readings at  
# noon on the 1st of the month in this dataset.

# ##############################################################################
# Task 1 Question 4 - Creating a function to compare multiple pollutants on the 
# same graph (noon on the 1st of each month across 2015-2024)
# ##############################################################################

comparePollutants <- function(pollutants){
  PortTalbotAQ$date <- as.Date(PortTalbotAQ$date, format ="%Y-%m-%d") # Ensuring the column is date format.
  start_date <- as.Date("2015-01-01") # Specifying start date
  end_date <- as.Date("2024-12-31") # specifying end date
  
  keep_rows <- PortTalbotAQ$date >= start_date & # This section is setting the parameters of the requested conditions
    PortTalbotAQ$date <= end_date &
    format(PortTalbotAQ$date, "%d") == "01" &
    PortTalbotAQ$time == "12:00:00"
  
  df <- PortTalbotAQ[keep_rows, c("date", pollutants)] # Ensure the function only keeps date & requested pollutant columns
  
  for (p in pollutants) {
    df[[p]][df[[p]] == "No data"] <- NA
    df[[p]] <- suppressWarnings(as.numeric(df[[p]]))
  } # This section is converting "No data" to NA and each pollutant column to numeric
  
  df <- df[rowSums(is.na(df[pollutants])) < length(pollutants), ] # Removing rows where all selected pollutants are missing
  
  if (nrow(df) == 0) {
    stop("No valid data found for the selected pollutants at noon on the 1st of the month (2015-2024).")
  }
  
  y_min <- min(sapply(pollutants, function(p) df[[p]]), na.rm =TRUE)
  y_max <- max(sapply(pollutants, function(p) df[[p]]), na.rm =TRUE)
  # The above determines the parameters for the selected pollutants 
  
  cols <- rep(c("blue", "red", "darkgreen", "purple"), length.out = length(pollutants)) # Adding colours to the plot lines to distinguish pollutants apart
  
  plot( # Arranging plot similar to Task 1 Question 3.
    df$date,
    df[[pollutants[1]]],
    type = "o",
    pch = 16,
    col = cols[1], # Adding colours in the plot
    xlab ="Date",
    ylab = "Pollutant level",
    main = "Comparison of pollutant levels (Noon on 1st of month, 2015-2024)",
    ylim = c(y_min, y_max) # This is being used as plotting many unlike Task 1 question 3 where we only plotted one.
  )
  
  if (length(pollutants) > 1) { # Adding the remaining pollutants as lines/points on the same graph
    for (i in 2:length(pollutants)) {
      lines(df$date, 
            df[[pollutants[i]]], 
            type = "o", 
            pch =16,
            col = cols[i] # Adding colours into the legend
            )
    }
  }
  
  legend("topright", 
         legend = pollutants, 
         col = cols,
         lty = 1, 
         pch = 16, 
         bty = "n"
         ) # lty = line type (1 = solid line), bty = box type (n = no box)
}

# Add 4 pollutants together
pollutants4 <- c("ozone", "nitrogen_dioxide", "PM10", "PM2_5")

comparePollutants(pollutants4) # Checking to see if worked

# ##############################################################################
# Task 1 Question 5 - Creating a function to calculate correlation between ozone
# and nitrogen dioxide and plot the data as a scatterplot.
# ##############################################################################

O3NOCorrelation <- function(p1,p2) {    # Creating the function
  
  df <- PortTalbotAQ[, c(p1, p2)] # Creating a smaller data frame with only the two columns needed
  
  df$ozone[df$ozone == "No data"] <- NA
  df$nitrogen_dioxide[df$nitrogen_dioxide == "No data"] <- NA # Replace "No data" with NA for both variables
  
  # Time to convert both columns to numeric (suppress warning caused by NA coercion)
  df$ozone <- suppressWarnings(as.numeric(df$ozone))
  df$nitrogen_dioxide <- suppressWarnings(as.numeric(df$nitrogen_dioxide))
  
  df <- df[!is.na(df$ozone) & !is.na(df$nitrogen_dioxide), ] # Keeps only rows where both variables required are valid numbers
  
  if (nrow(df) == 0) {
    stop("No valid paired ozone and nitorgen_dioxide measurements available to calculate correlation.") # Safety check in cvase nothing remains
  }
  
  corr_value <- cor(df$ozone, df$nitrogen_dioxide) # Calculate Pearson correlation
  
  plot(                    # Plot the paired data as a scatterplot
    df$ozone,
    df$nitrogen_dioxide,
    pch = 16,
    col = rgb(0, 0, 0, 0.2), # Helps show density patterns
    xlab = "Ozone",
    ylab = "Nitrogen dioxide",
    main = paste("Ozone vs Nitrogen dioxide (Correlation =", round(corr_value, 3), ")") # rounded to 3 decimal places and rounded
  )
  
  return(round(corr_value, 3))
}

O3NOCorrelation("ozone", "nitrogen_dioxide") # Returns output as required

# rm(O3NOCorrelation)

# ##############################################################################
# Task 1 Question 6 - T-test: Ozone levels when wind speed < 5 m/s vs >= 5 m/s
# ##############################################################################

# Experimental hypotheses:
# H0 (null): The mean ozone level is the same when wind speed < 5 m/s and when wind speed >= 5 m/s.
# H1 (alternative): The mean ozone level is different between the two wind speed groups.

df <- PortTalbotAQ[, c("ozone", "wind_speed")] # Creating a smaller data frame with only required columns needed

df$ozone[df$ozone == "No data"] <- NA
df$wind_speed[df$wind_speed == "No data"] <- NA # Similar to previous code replacing No Data with NA for both data types

df$ozone <- suppressWarnings(as.numeric(df$ozone))
df$wind_speed <- suppressWarnings(as.numeric(df$wind_speed)) #Converting both columns to numeric (suppress warning from coercion)

df <- df[!is.na(df$ozone) & !is.na(df$wind_speed), ] # Keeps only rows where both ozone and wind_speed are valid

# Splitting ozone into two groups to assist further calculations (based on wind speed)
ozone_low_wind <- df$ozone[df$wind_speed < 5] # Wind speed below 5 m/s
ozone_high_wind <- df$ozone[df$wind_speed >= 5] # wind speed 5 m/s or higher

if (length(ozone_low_wind) == 0 | length(ozone_high_wind) == 0) {
  stop("Not enough valid data in one or both wind speed groups to perform the t-test.") # This ensures both groups hold data and if no data is present this comment will be outputted. 
}

ttest_result <- t.test(ozone_low_wind, ozone_high_wind, alternative = "two.sided") # Two-sample t-test (Welch)

ttest_result # Provides output of above Welch T-test

# The above worked and provided an output. Time to add the decision rule and conclusion
# If p-value < 0.05, reject H0 (significant difference in mean ozone between groups)
# If p-value >= 0.05, fail to reject H0 (no significant difference detected)

alpha <- 0.05

if (ttest_result$p.value < alpha) {
  message ("Conclusion: p-value = ", round(ttest_result$p.value, 4),
           " < ", alpha, ". Reject H0. There is a statistically significant difference in mean ozone levels between wind speed groups.") # Reject H0
} else {
  message ("Conclusion: p-value = ", round(ttest_result$p.value, 4),
           " >= ", alpha, ". Fail to reject H0. There is no statistically significant difference in mean ozone levels between wind speed groups.") # Fail to reject H0
}

# When the above has been run it has outputted the Reject H0 option which shows the code is working.

# ##############################################################################
# Task 1 Question 7 - ANOVA: Nitrogen dioxide levels across selected periods
# (Nov-Dec 2016, 2018, 2020 at 17:00)
# ##############################################################################

# Experimental Hypotheses again similar to Task 1 Question 6:
# H0 (null): Mean nitrogen dioxide levels are equal across all three periods required
# H1 (alternative): At least one period has a different mean nitrogen dioxide level

# As the question is requiring specific dates will intially ensure the date column is date format.
PortTalbotAQ$date <- as.Date(PortTalbotAQ$date, format = "%Y-%m-%d")

df <- PortTalbotAQ[, c("date", "time", "nitrogen_dioxide")] # Picking only relevant columns

df$nitrogen_dioxide[df$nitrogen_dioxide == "No data"] <- NA # Replacing "No data with NA like before

df$nitrogen_dioxide <- suppressWarnings(as.numeric(df$nitrogen_dioxide)) # Converting to numeric

df <- df[!is.na(df$nitrogen_dioxide), ] # Removing rows with missing values

df <- df[df$time == "17:00:00", ] # Setting the exact time in which the data need to be pulled

df$year <- as.numeric(format(df$date, "%Y")) # Creating year & month variables
df$month <- as.numeric(format(df$date, "%m"))

df <- df[df$month %in% c(11, 12), ] # Keeping only Nov & December data

df <- df[df$year %in% c(2016, 2018, 2020), ] # Keeping only the three years of interest

df$period <- factor(df$year,
                    levels = c(2016, 2018, 2020),
                    labels = c("Nov-Dec 2016", "Nov-Dec 2018", "Nov-Dec 2020")) # Creating a factor representing the three comparison periods

if (length(unique(df$period)) < 3) {
  stop("Not enough data available for all three periods to perform ANOVA.") # Automated response if not enough available data to carry out request.
}

# Main code has now been added along with a safety check time to add the ANOVA code

anova_result <- aov(nitrogen_dioxide ~ period, data = df) # Run one-way ANOVA

summary(anova_result) # Displaying result

# Post-Hoc tests if significant

Alpha <- 0.05
anova_p <- summary(anova_result)[[1]]["period", "Pr(>F)"]

if (anova_p < alpha) {
  message("ANOVA result: p-value = ", round(anova_p, 4),
          " < ", alpha, ". Reject H0. At least one period differs significantly.")
  # Post-Hoc explanation: ANOVA tells us a difference exists overall, but not which groups differ
  # while Tukey HSD compares all pairs while controlling for multiple comparisons.
  TukeyHSD(anova_result)
} else {
  message("ANOVA result: p-value = ", round(anova_p, 4),
          " >= ", alpha, ". Fail to reject H0. No signficant differences detected.")
}

# ANOVA results indicated a statistically significant difference in mean nitrogen dioxide
# levels across the three periods (p < 0.05), therefore the null hypothesis was rejected.
#
# A Tukey HSD post-hoc test was conducted to identify where the differences lay.
# The results showed that nitrogen dioxide levels in Nov-Dec 2016 were significantly
# higher than in both Nov-Dec 2018 and Nov-Dec 2020.
#
# No statistically significant difference was found between the Nov-Dec 2018 and
# Nov-Dec 2020 periods.
# ##############################################################################