--- title: "2434902" output: html_document --- ```{r} # 1.Setup and Data Loading # Install necessary library if (!requireNamespace("tidyverse", quietly = TRUE)) { install.packages("tidyverse") } if (!requireNamespace("lubridate", quietly = TRUE)) { install.packages("lubridate") } # Load the libraries for data manipulation and visualization library(tidyverse) library(lubridate) # Create a data frame PortTalbotAQ, ensuring that the headers are imported. # Also, combine Date and Time columns and rename headers immediately for clean workflow. PortTalbotAQ <- read.csv("Port_Talbot_Hourly_AQ.csv", header = TRUE, stringsAsFactors = FALSE, na.strings = "No data") %>% # Set "No data" to NA # Convert Date and Time to a single DateTime object mutate( DateTime = as.POSIXct(paste(Date, Time), format = "%d-%m-%Y %H:%M:%S", tz = "GMT") ) %>% select( DateTime, Ozone, Nitric_Oxide = `Nitric.oxide`, Nitrogen_Dioxide = `Nitrogen.dioxide`, Nitrogen_Oxides = `Nitrogen.oxides.as.nitrogen.dioxide`, # Renamed from 'Nitrogen oxides as nitrogen dioxide' Sulphur_Dioxide = `Sulphur.dioxide`, Carbon_Monoxide = `Carbon.monoxide`, PM10 = `PM10.particulate.matter..Hourly.measured.`, PM2.5 = `PM2.5.particulate.matter..Hourly.measured.`, Wind_Direction = `Modelled.Wind.Direction`, Wind_Speed = `Modelled.Wind.Speed`, Temperature = `Modelled.Temperature` ) # Coerce relevant column numeric_cols <- c("Ozone", "Nitric_Oxide", "Nitrogen_Dioxide", "Nitrogen_Oxides", "Sulphur_Dioxide", "Carbon_Monoxide", "PM10", "PM2.5", "Wind_Direction", "Wind_Speed", "Temperature") for (col in numeric_cols) { PortTalbotAQ[[col]] <- as.numeric(PortTalbotAQ[[col]]) } # Display data frame. print("1.3 Resulting PortTalbotAQ data frame (first 6 rows):") print(head(PortTalbotAQ)) ``` ```{r} # 2. Missing Data Analysis # Create a new data frame # readings for each pollution. pollutant_cols <- c("Ozone", "Nitric_Oxide", "Nitrogen_Dioxide", "Nitrogen_Oxides", "Sulphur_Dioxide", "Carbon_Monoxide", "PM10", "PM2.5") invalidData <- PortTalbotAQ %>% summarise(across(all_of(pollutant_cols), ~ mean(is.na(.)) * 100)) %>% pivot_longer(everything(), names_to = "Pollutant", values_to = "Percentage_Missing") # Display this data frame. print("\n2.2 InvalidData frame (Percentage of missing values per pollutant):") print(invalidData) ``` ```{r} # 3.Time-Series Plotting # Create a filtered dataset with only data monthly_noon_data <- PortTalbotAQ %>% mutate( Day = day(DateTime), Hour = hour(DateTime) ) %>% filter(Day == 1, Hour == 12) # Create a function displayPollutantGraph() displayPollutantGraph <- function(pollutant_name, start_month, start_year, end_month, end_year) { # Convert inputs to date objects for filtering start_date <- as.POSIXct(paste0("01-", start_month, "-", start_year, " 12:00:00"), format = "%d-%m-%Y %H:%M:%S", tz = "GMT") # Filter up to the end date of the last month/year specified. end_date <- as.POSIXct(paste0("01-", end_month, "-", end_year, " 12:00:00"), format = "%d-%m-%Y %H:%M:%S", tz = "GMT") # For inclusive end month, we use the start of the last month + 1 month. end_date_inclusive <- end_date + months(1) # Filter data for the specified pollutant and time range plot_data <- monthly_noon_data %>% filter(DateTime >= start_date, DateTime < end_date_inclusive) %>% # Select the date and the specified pollutant column select(DateTime, !!sym(pollutant_name)) %>% na.omit() # Plotting with ggplot2 p <- ggplot(plot_data, aes(x = DateTime, y = !!sym(pollutant_name))) + geom_line(color = "darkblue", size = 1) + geom_point(color = "red", size = 2) + labs( title = paste("Variation in", gsub("_", " ", pollutant_name), "Measured at Noon on the 1st of the Month"), subtitle = paste0("From ", start_month, "/", start_year, " to ", end_month, "/", end_year), x = "Date", y = paste(gsub("_", " ", pollutant_name), "Concentration") ) + theme_minimal() + theme( plot.title = element_text(hjust = 0.5, face = "bold"), plot.subtitle = element_text(hjust = 0.5), axis.text.x = element_text(angle = 45, hjust = 1) ) print(p) } # Display this using the data for nitrogen oxides as nitrogen dioxide between July 2015 and November 2024 inclusive. # The provided data file starts in Jan 2020, so the requested range of July 2015 to Nov 2024 is partially achievable. The graph will reflect the available subset. print("\n3.2 Graph for Nitrogen Oxides as Nitrogen Dioxide (July 2015 - Nov 2024):") # Pollutant name is 'Nitrogen_Oxides' (the renamed column) displayPollutantGraph("Nitrogen_Oxides", 7, 2015, 11, 2024) ``` ```{r} # 4.Time-Series Plotting - Multiple Pollutants # Create a function comparePollutants() comparePollutants <- function(pollutant_vector) { # Prepare data for plotting plot_data <- monthly_noon_data %>% # Select DateTime and the specified pollutant columns select(DateTime, all_of(pollutant_vector)) %>% # Reshape to long format for ggplot pivot_longer(cols = -DateTime, names_to = "Pollutant", values_to = "Concentration") %>% na.omit() # Plotting with ggplot p <- ggplot(plot_data, aes(x = DateTime, y = Concentration, color = Pollutant)) + geom_line(size = 1) + labs( title = "Comparison of Pollutant Levels Measured at Noon on the 1st of the Month", subtitle = "All Available Years", x = "Date", y = "Concentration (Mixed Units)", # Units are mixed (ugm-3, mgm-3), noting this limitation color = "Pollutant" ) + scale_color_brewer(palette = "Set1") + # Use a distinct color palette theme_minimal() + theme( plot.title = element_text(hjust = 0.5, face = "bold"), plot.subtitle = element_text(hjust = 0.5), axis.text.x = element_text(angle = 45, hjust = 1) ) print(p) } # Display this using an appropriate vector. selected_pollutants <- c("Ozone", "Nitrogen_Dioxide", "PM10", "PM2.5") print("\n4.2 Graph for Comparing 4 Pollutants:") comparePollutants(selected_pollutants) ``` ```{r} # 5.Correlation and Scatterplot # Create a function O3NOCorrelation() that calculates the correlation between ozone and nitrogen dioxide measurements for the entire data collection period. O3NOCorrelation <- function(data_frame) { # Select Ozone and Nitrogen Dioxide columns for the entire data collection period correlation_data <- data_frame %>% select(Ozone, Nitrogen_Dioxide) %>% # Remove rows with any NA values for correlation calculation and plotting na.omit() # Calculate the Pearson correlation coefficient correlation_value <- cor(correlation_data$Ozone, correlation_data$Nitrogen_Dioxide, method = "pearson") # Plot the data used as a scatterplot. p <- ggplot(correlation_data, aes(x = Ozone, y = Nitrogen_Dioxide)) + geom_point(alpha = 0.5, color = "darkgreen") + geom_smooth(method = "lm", col = "red", se = FALSE) + # Add linear regression line labs( title = paste("Scatterplot of Ozone vs. Nitrogen Dioxide"), subtitle = paste0("Pearson Correlation Coefficient (r): ", round(correlation_value, 4)), x = "Ozone Concentration", y = "Nitrogen Dioxide Concentration" ) + theme_minimal() + theme(plot.title = element_text(hjust = 0.5, face = "bold"), plot.subtitle = element_text(hjust = 0.5)) print("\n5.2 Scatterplot of Ozone vs. Nitrogen Dioxide:") print(p) # Return the correlation value return(correlation_value) } # Execute the function correlation_result <- O3NOCorrelation(PortTalbotAQ) print(paste0("\nCorrelation between Ozone and Nitrogen Dioxide: ", round(correlation_result, 4))) ``` ```{r} # 6.T-Test # Perform a t-test to identify whether there is a significant difference in ozone levels when wind speeds are below 5ms-1 and when they are equal to or greater than 5ms-1. # Filter data for the t-test and create the grouping factor t_test_data <- PortTalbotAQ %>% select(Ozone, Wind_Speed) %>% na.omit() %>% mutate( # Create a factor variable for the two groups Wind_Group = factor(ifelse(Wind_Speed < 5, "Below_5ms-1", "5ms-1_or_Greater")) ) # Perform the t-test # t.test is appropriate for comparing the means of two independent groups t_test_result <- t.test(Ozone ~ Wind_Group, data = t_test_data) print("\n6.2 T-Test Results:") print(t_test_result) # State your experimental hypotheses, the result of the test, and whether you accept or reject the null hypothesis in your comments. cat( " # T-Test # Experimental Hypotheses: # H0 (Null Hypothesis): There is no significant difference in mean ozone levels # when wind speeds are below 5ms-1 and when they are # equal to or greater than 5ms-1 (mu_low = mu_high). # Ha (Alternative Hypothesis): There is a significant difference in mean ozone # levels when wind speeds are below 5ms-1 and when # they are equal to or greater than 5ms-1 (mu_low != mu_high). # Result of the test: # The mean ozone concentration for Wind_Group '5ms-1_or_Greater' (", round(t_test_result$estimate[1], 4), ") is ", ifelse(t_test_result$estimate[1] > t_test_result$estimate[2], "higher", "lower"), " than for 'Below_5ms-1' (", round(t_test_result$estimate[2], 4), "). # The p-value from the t-test is ", round(t_test_result$p.value, 4), ". # Conclusion (Accept or Reject Null Hypothesis): # Assuming a significance level (alpha) of 0.05, since the p-value is ", ifelse(t_test_result$p.value < 0.05, "less than", "greater than or equal to"), " 0.05: # We ", ifelse(t_test_result$p.value < 0.05, "REJECT", "FAIL TO REJECT"), " the null hypothesis (H0). # Conclusion: There ", ifelse(t_test_result$p.value < 0.05, "IS", "IS NOT"), " a statistically significant difference in mean ozone levels based on the wind speed categories. " ) ``` ```{r} # TASK 7 (20 marks) # ANOVA: NO2 (Nitrogen_Dioxide) at 17:00 during Nov–Dec of 2016 vs 2018 vs 2020 # Tukey HSD only if ANOVA is significant runNO2_ANOVA_17h_NovDec <- function() { # Require columns req <- c("Nitrogen_Dioxide", "DateTime") if (!all(req %in% names(PortTalbotAQ))) { stop("Required columns missing: 'Nitrogen_Dioxide' and/or 'DateTime'.") } # Filter Nov–Dec, 17:00, years 2016/2018/2020 Year <- as.integer(format(PortTalbotAQ$DateTime, "%Y")) Month <- as.integer(format(PortTalbotAQ$DateTime, "%m")) Hour <- format(PortTalbotAQ$DateTime, "%H") keep <- Year %in% c(2016, 2018, 2020) & Month %in% c(11, 12) & Hour == "17" df <- data.frame( DateTime = PortTalbotAQ$DateTime[keep], NO2 = as.numeric(PortTalbotAQ$Nitrogen_Dioxide[keep]), Year = Year[keep] ) df <- df[is.finite(df$NO2), ] # Label periods df$Period <- factor( df$Year, levels = c(2016, 2018, 2020), labels = c("Nov–Dec 2016", "Nov–Dec 2018", "Nov–Dec 2020") ) # Summaries grp <- split(df$NO2, df$Period) summary_table <- data.frame( Period = names(grp), N = sapply(grp, length), Mean = round(sapply(grp, function(x) if(length(x)) mean(x) else NA_real_), 3), SD = round(sapply(grp, function(x) if(length(x)>1) sd(x) else NA_real_), 3), row.names = NULL ) print(summary_table) # Check groups: all three present and have at least 2 obs each ok_groups <- all(!is.na(match(levels(df$Period), summary_table$Period))) && all(summary_table$N[match(levels(df$Period), summary_table$Period)] >= 2) if (!ok_groups) { cat(" ANOVA Test Commentary Required periods: Nov–Dec 2016, Nov–Dec 2018, Nov–Dec 2020 at 17:00. At least one period is missing or has fewer than 2 observations. Therefore, a one-way ANOVA cannot be performed on the specified periods. H0: mean(NO2) is equal across the three periods. H1: At least one period's mean differs. Decision: Cannot test H0/H1 due to insufficient data in the required periods. ") return(invisible(list(summary = summary_table))) } # ANOVA fit <- aov(NO2 ~ Period, data = df) anova_tab <- summary(fit) print(anova_tab) pval <- anova_tab[[1]][["Pr(>F)"]][1] # Boxplot for a quick visual boxplot(NO2 ~ Period, data = df, xlab = "Period", ylab = "NO2 (as recorded)", main = "NO2 at 17:00 across Nov–Dec periods (2016 vs 2018 vs 2020)") grid(nx = NA, ny = NULL) cat("ANOVA Test Commentary Experimental hypotheses: H0: mean(NO2) Nov–Dec 2016 = mean(NO2) Nov–Dec 2018 = mean(NO2) Nov–Dec 2020 H1: At least one period's mean is different. Result: ANOVA p-value =", round(pval, 4), " ") if (is.finite(pval) && pval < 0.05) { cat("Decision (α = 0.05): REJECT H0. There IS a significant difference.\n") # Post-hoc: Tukey HSD to identify which pairs differ (controls family-wise error) tuk <- TukeyHSD(fit) print(tuk) tuk_df <- transform(as.data.frame(tuk[[1]]), Comparison = rownames(tuk[[1]])) rownames(tuk_df) <- NULL names(tuk_df) <- sub("^p adj$", "p.adj", names(tuk_df)) # normalize name tuk_df <- tuk_df[, c("Comparison", "diff", "lwr", "upr", if ("p.adj" %in% names(tuk_df)) "p.adj" else names(tuk_df)[grepl("^p", names(tuk_df))])] names(tuk_df) <- c("Comparison", "Mean_Diff", "CI_Lower", "CI_Upper", "Adj_P") tuk_df$Mean_Diff <- round(tuk_df$Mean_Diff, 3) tuk_df$CI_Lower <- round(tuk_df$CI_Lower, 3) tuk_df$CI_Upper <- round(tuk_df$CI_Upper, 3) tuk_df$Adj_P <- signif(tuk_df$Adj_P, 4) print(tuk_df) cat(" Post-hoc rationale: Because ANOVA found a significant overall effect, we applied Tukey's HSD to identify which specific period pairs differ while controlling family-wise Type I error. ") return(invisible(list(summary = summary_table, anova = anova_tab, tukey = tuk_df))) } else { cat("Decision (α = 0.05): FAIL TO REJECT H0. No significant difference detected. ") return(invisible(list(summary = summary_table, anova = anova_tab))) } } runNO2_ANOVA_17h_NovDec() ```