--- title: "Untitled" output: word_document: default html_document: default --- ```{r} # ER NO:- [2434357] # 1.Setup and Data Loading # Install necessary libraries if they are not already installed 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. # 'No data' entries are correctly interpreted as missing values (NA). # 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") ) %>% # Rename the headers to shorter but still understandable descriptors. # Select and rename columns in one step. select( Date, Time, 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 columns to numeric, converting any remaining non-numeric data to NA 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 the 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 invalidData that shows the percentage of missing # readings for each pollutant. 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 - Single Pollutant # Create a filtered dataset with only data from the 1st of each month at noon (12:00:00) 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 ggplot2 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 of 4 pollutants. 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 Commentary # 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 — Consolidated ANOVA script (All 3 versions together) # Requires: PortTalbotAQ with columns Nitrogen_Dioxide, DateTime (POSIXct) # Also uses: stats, graphics, utils; optionally dplyr, lubridate if available run_Task7_ALL <- function() { # ---- helpers ---- has_pkg <- function(p) { requireNamespace(p, quietly = TRUE) } `%||%` <- function(a, b) if (!is.null(a)) a else b stopifnot(all(c("Nitrogen_Dioxide", "DateTime") %in% names(PortTalbotAQ))) # Try to use dplyr/lubridate if available; else base R fallback use_dplyr <- has_pkg("dplyr") && has_pkg("lubridate") if (use_dplyr) { dplyr <- asNamespace("dplyr") lubri <- asNamespace("lubridate") } # Version 1: Your first function (robust; 2016 vs 2018 vs 2020; Tukey if sig) runNO2_ANOVA_17h_NovDec <- function() { req <- c("Nitrogen_Dioxide", "DateTime") if (!all(req %in% names(PortTalbotAQ))) { stop("Required columns missing: 'Nitrogen_Dioxide' and/or 'DateTime'.") } 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), ] df$Period <- factor( df$Year, levels = c(2016, 2018, 2020), labels = c("Nov–Dec 2016", "Nov–Dec 2018", "Nov–Dec 2020") ) 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 ) cat("\n================ V1: Summary ==================\n"); print(summary_table) 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(" V1 — 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))) } fit <- aov(NO2 ~ Period, data = df) anova_tab <- summary(fit) cat("\n================ V1: ANOVA ==================\n"); print(anova_tab) pval <- anova_tab[[1]][["Pr(>F)"]][1] boxplot(NO2 ~ Period, data = df, xlab = "Period", ylab = "NO2 (as recorded)", main = "V1: NO2 at 17:00 (Nov–Dec 2016 vs 2018 vs 2020)") grid(nx = NA, ny = NULL) cat(" V1 — 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), "\n") if (is.finite(pval) && pval < 0.05) { cat("Decision (α = 0.05): REJECT H0. There IS a significant difference.\n") tuk <- TukeyHSD(fit) cat("\n--------------- V1: Tukey HSD ---------------\n"); 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)) 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) cat("\nTukey HSD (tidy):\n"); 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))) } } # Version 2: Substitute years demo (2020, 2021, 2022) if earlier years missing run_Substitute_ANOVA <- function() { target_hour <- 17 target_months <- c(11, 12) subs_years <- c(2020, 2021, 2022) if (use_dplyr) { anova_data_sub <- PortTalbotAQ %>% dplyr$mutate( Month = lubri$month(DateTime), Year = lubri$year(DateTime), Hour = lubri$hour(DateTime) ) %>% dplyr$filter( Month %in% target_months, Year %in% subs_years, Hour == target_hour ) %>% dplyr$select(DateTime, Nitrogen_Dioxide, Year) %>% stats::na.omit() } else { dt <- PortTalbotAQ Month <- as.integer(format(dt$DateTime, "%m")) Year <- as.integer(format(dt$DateTime, "%Y")) Hour <- as.integer(format(dt$DateTime, "%H")) keep <- Month %in% target_months & Year %in% subs_years & Hour == target_hour anova_data_sub <- stats::na.omit(data.frame( DateTime = dt$DateTime[keep], Nitrogen_Dioxide = dt$Nitrogen_Dioxide[keep], Year = Year[keep] )) } if (length(unique(anova_data_sub$Year)) == 3) { anova_data_sub$Year_Factor <- factor(anova_data_sub$Year) aov_res <- aov(Nitrogen_Dioxide ~ Year_Factor, data = anova_data_sub) anova_sum <- summary(aov_res) pval <- anova_sum[[1]]$'Pr(>F)'[1] cat("\n================ V2: ANOVA (Substitute 2020/2021/2022) ================\n") print(anova_sum) cat(" # V2 — ANOVA Test Commentary (Substitute Data: 2020, 2021, 2022) # NOTE: If 2016 and 2018 are missing, this demonstrates the technique using # available Nov–Dec 17:00 data for 2020–2022. # H0: μ2020 = μ2021 = μ2022 # H1: At least one mean differs. # p-value:", round(pval, 4), "\n") boxplot(Nitrogen_Dioxide ~ Year_Factor, data = anova_data_sub, xlab = "Year", ylab = "NO2", main = "V2: NO2 at 17:00 (Nov–Dec 2020 vs 2021 vs 2022)") grid(nx = NA, ny = NULL) if (is.finite(pval) && pval < 0.05) { cat("Decision: REJECT H0; run Tukey HSD.\n") tuk <- TukeyHSD(aov_res) print(tuk) cat(" # Post-hoc: Tukey HSD controls family-wise error across year-pair comparisons. ") return(invisible(list(anova = anova_sum, tukey = tuk))) } else { cat("Decision: FAIL TO REJECT H0 (no significant difference).\n") return(invisible(list(anova = anova_sum))) } } else { cat("\n================ V2: ANOVA (Substitute) =================\n") cat("Cannot perform substitute ANOVA (2020–2022): insufficient groups/data.\n") return(invisible(NULL)) } } # ---------------------------------------------------------- # Version 3: Simple check + ANOVA (2016/2018/2020) mirroring your last block # ---------------------------------------------------------- run_Simple_ANOVA_Check <- function() { target_years <- c(2016, 2018, 2020) target_months <- c(11, 12) target_hour <- 17 if (use_dplyr) { anova_data <- PortTalbotAQ %>% dplyr$mutate( Year = lubri$year(DateTime), Month = lubri$month(DateTime), Hour = lubri$hour(DateTime) ) %>% dplyr$filter( Year %in% target_years, Month %in% target_months, Hour == target_hour ) %>% dplyr$select(Nitrogen_Dioxide, Year) %>% tidyr::drop_na() } else { dt <- PortTalbotAQ Year <- as.integer(format(dt$DateTime, "%Y")) Month <- as.integer(format(dt$DateTime, "%m")) Hour <- as.integer(format(dt$DateTime, "%H")) keep <- Year %in% target_years & Month %in% target_months & Hour == target_hour anova_data <- stats::na.omit(data.frame( Nitrogen_Dioxide = dt$Nitrogen_Dioxide[keep], Year = Year[keep] )) } group_counts <- table(anova_data$Year) if (length(group_counts) < 3) { cat(" ============================================================================== V3 — ANOVA Test Commentary ============================================================================== The dataset does not contain nitrogen dioxide measurements for all required years (2016, 2018, 2020) for Nov–Dec at 17:00. Because at least one year-group has no data, a one-way ANOVA cannot be performed. ============================================================================== ") return(invisible(NULL)) } else { anova_data$Year <- factor(anova_data$Year) aov_model <- aov(Nitrogen_Dioxide ~ Year, data = anova_data) anova_summary <- summary(aov_model) cat("\n================ V3: ANOVA (Simple) ==================\n"); print(anova_summary) p_value <- anova_summary[[1]]$"Pr(>F)"[1] cat(" ============================================================================== V3 — ANOVA Test Commentary ============================================================================== H0: μ2016 = μ2018 = μ2020 H1: At least one mean differs. Result: p-value =", round(p_value, 4), "\n") boxplot(Nitrogen_Dioxide ~ Year, data = anova_data, xlab = "Year", ylab = "NO2", main = "V3: NO2 at 17:00 (Nov–Dec 2016 vs 2018 vs 2020)") grid(nx = NA, ny = NULL) if (p_value < 0.05) { cat(" Since p < 0.05, we REJECT H0. Tukey HSD follows. ") tukey_result <- TukeyHSD(aov_model) print(tukey_result) cat(" ============================================================================== V3 — Post-hoc Test Commentary ============================================================================== Tukey HSD identifies which year-pairs differ with FWER control. ============================================================================== ") return(invisible(list(anova = anova_summary, tukey = tukey_result))) } else { cat(" Since p ≥ 0.05, we FAIL TO REJECT H0. No significant difference detected. ============================================================================== ") return(invisible(list(anova = anova_summary))) } } } # ---- RUN ALL THREE IN ORDER ---- cat("\n\n############################################################\n") cat("# TASK 7 — Consolidated Outputs: V1, V2 (Substitute), V3 #\n") cat("############################################################\n\n") res1 <- runNO2_ANOVA_17h_NovDec() res2 <- run_Substitute_ANOVA() res3 <- run_Simple_ANOVA_Check() invisible(list(V1 = res1, V2_substitute = res2, V3_simple = res3)) } # ---- Execute all three now ---- result <- run_Task7_ALL() ```